Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "ArrayKernel.h"
      11             : 
      12             : #include "Assembly.h"
      13             : #include "MooseVariableFE.h"
      14             : #include "MooseVariableScalar.h"
      15             : #include "SubProblem.h"
      16             : #include "NonlinearSystem.h"
      17             : #include "FEProblemBase.h"
      18             : 
      19             : #include "libmesh/threads.h"
      20             : #include "libmesh/quadrature.h"
      21             : 
      22             : InputParameters
      23       24993 : ArrayKernel::validParams()
      24             : {
      25       24993 :   InputParameters params = KernelBase::validParams();
      26       24993 :   params.registerBase("ArrayKernel");
      27       24993 :   return params;
      28           0 : }
      29             : 
      30        1841 : ArrayKernel::ArrayKernel(const InputParameters & parameters)
      31             :   : KernelBase(parameters),
      32             :     MooseVariableInterface<RealEigenVector>(this,
      33             :                                             false,
      34             :                                             "variable",
      35             :                                             Moose::VarKindType::VAR_SOLVER,
      36             :                                             Moose::VarFieldType::VAR_FIELD_ARRAY),
      37        3682 :     _var(*mooseVariable()),
      38        1841 :     _test(_var.phi()),
      39        1841 :     _grad_test(_var.gradPhi()),
      40        1841 :     _array_grad_test(_var.arrayGradPhi()),
      41        1841 :     _phi(_assembly.phi(_var)),
      42        1841 :     _grad_phi(_assembly.gradPhi(_var)),
      43        1841 :     _u(_is_implicit ? _var.sln() : _var.slnOld()),
      44        1841 :     _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld()),
      45        1841 :     _count(_var.count()),
      46        7364 :     _work_vector(_count)
      47             : {
      48        1841 :   addMooseVariableDependency(mooseVariable());
      49             : 
      50        1841 :   _save_in.resize(_save_in_strings.size());
      51        1841 :   _diag_save_in.resize(_diag_save_in_strings.size());
      52             : 
      53        1854 :   for (unsigned int i = 0; i < _save_in_strings.size(); i++)
      54             :   {
      55          13 :     ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _save_in_strings[i]);
      56             : 
      57          13 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_save_in_strings[i]))
      58           0 :       paramError("save_in", "cannot use solution variable as save-in variable");
      59             : 
      60          13 :     if (var->feType() != _var.feType())
      61           0 :       paramError(
      62             :           "save_in",
      63             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
      64           0 :           moose::internal::incompatVarMsg(*var, _var));
      65             : 
      66          13 :     _save_in[i] = var;
      67          13 :     var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
      68          13 :     addMooseVariableDependency(var);
      69             :   }
      70             : 
      71        1841 :   _has_save_in = _save_in.size() > 0;
      72             : 
      73        1854 :   for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
      74             :   {
      75          13 :     ArrayMooseVariable * var = &_subproblem.getArrayVariable(_tid, _diag_save_in_strings[i]);
      76             : 
      77          13 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_diag_save_in_strings[i]))
      78           0 :       paramError("diag_save_in", "cannot use solution variable as diag save-in variable");
      79             : 
      80          13 :     if (var->feType() != _var.feType())
      81           0 :       paramError(
      82             :           "diag_save_in",
      83             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
      84           0 :           moose::internal::incompatVarMsg(*var, _var));
      85             : 
      86          13 :     _diag_save_in[i] = var;
      87          13 :     var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
      88          13 :     addMooseVariableDependency(var);
      89             :   }
      90             : 
      91        1841 :   _has_diag_save_in = _diag_save_in.size() > 0;
      92        1841 : }
      93             : 
      94             : void
      95     4945023 : ArrayKernel::computeResidual()
      96             : {
      97     4945023 :   prepareVectorTag(_assembly, _var.number());
      98             : 
      99     4945023 :   precalculateResidual();
     100    27372747 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     101             :   {
     102    22427724 :     initQpResidual();
     103   133059756 :     for (_i = 0; _i < _test.size(); _i++)
     104             :     {
     105   110632032 :       _work_vector.setZero();
     106   110632032 :       computeQpResidual(_work_vector);
     107             :       mooseAssert(_work_vector.size() == _count,
     108             :                   "Size of local residual is not equal to the number of array variable compoments");
     109   110632032 :       _work_vector *= _JxW[_qp] * _coord[_qp];
     110   110632032 :       _assembly.saveLocalArrayResidual(_local_re, _i, _test.size(), _work_vector);
     111             :     }
     112             :   }
     113             : 
     114     4945023 :   accumulateTaggedLocalResidual();
     115             : 
     116     4945023 :   if (_has_save_in)
     117             :   {
     118         768 :     for (const auto & var : _save_in)
     119             :     {
     120         384 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     121         384 :       if (avar)
     122         384 :         avar->addSolution(_local_re);
     123             :       else
     124           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     125             :     }
     126             :   }
     127     4945023 : }
     128             : 
     129             : void
     130      471476 : ArrayKernel::computeJacobian()
     131             : {
     132      471476 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     133             : 
     134      471476 :   precalculateJacobian();
     135     2362564 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     136             :   {
     137     1891088 :     initQpJacobian();
     138     9456112 :     for (_i = 0; _i < _test.size(); _i++)
     139    38199040 :       for (_j = 0; _j < _phi.size(); _j++)
     140             :       {
     141    30634016 :         _work_vector = computeQpJacobian() * _JxW[_qp] * _coord[_qp];
     142   122536064 :         _assembly.saveDiagLocalArrayJacobian(
     143    30634016 :             _local_ke, _i, _test.size(), _j, _phi.size(), _var.number(), _work_vector);
     144             :       }
     145             :   }
     146             : 
     147      471476 :   accumulateTaggedLocalMatrix();
     148             : 
     149      471476 :   if (_has_diag_save_in)
     150             :   {
     151           0 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     152           0 :     for (const auto & var : _diag_save_in)
     153             :     {
     154           0 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     155           0 :       if (avar)
     156           0 :         avar->addSolution(diag);
     157             :       else
     158           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     159             :     }
     160           0 :   }
     161      471476 : }
     162             : 
     163             : RealEigenVector
     164    22556848 : ArrayKernel::computeQpJacobian()
     165             : {
     166    22556848 :   return RealEigenVector::Zero(_var.count());
     167             : }
     168             : 
     169             : void
     170      703890 : ArrayKernel::computeOffDiagJacobian(const unsigned int jvar_num)
     171             : {
     172      703890 :   const auto & jvar = getVariable(jvar_num);
     173             : 
     174      703890 :   bool same_var = (jvar_num == _var.number());
     175             : 
     176      703890 :   prepareMatrixTag(_assembly, _var.number(), jvar_num);
     177             : 
     178      703890 :   precalculateOffDiagJacobian(jvar_num);
     179     3555750 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     180             :   {
     181     2851860 :     initQpOffDiagJacobian(jvar);
     182    14575116 :     for (_i = 0; _i < _test.size(); _i++)
     183    61750200 :       for (_j = 0; _j < _phi.size(); _j++)
     184             :       {
     185    50026944 :         _work_matrix = computeQpOffDiagJacobian(jvar) * _JxW[_qp] * _coord[_qp];
     186   200107776 :         _assembly.saveFullLocalArrayJacobian(
     187    50026944 :             _local_ke, _i, _test.size(), _j, _phi.size(), _var.number(), jvar_num, _work_matrix);
     188             :       }
     189             :   }
     190             : 
     191      703890 :   accumulateTaggedLocalMatrix();
     192             : 
     193      703890 :   if (_has_diag_save_in && same_var)
     194             :   {
     195         256 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     196         512 :     for (const auto & var : _diag_save_in)
     197             :     {
     198         256 :       auto * avar = dynamic_cast<ArrayMooseVariable *>(var);
     199         256 :       if (avar)
     200         256 :         avar->addSolution(diag);
     201             :       else
     202           0 :         mooseError("Save-in variable for an array kernel must be an array variable");
     203             :     }
     204         256 :   }
     205      703890 : }
     206             : 
     207             : RealEigenMatrix
     208    46777768 : ArrayKernel::computeQpOffDiagJacobian(const MooseVariableFEBase & jvar)
     209             : {
     210    46777768 :   if (jvar.number() == _var.number())
     211    44982792 :     return computeQpJacobian().asDiagonal();
     212             :   else
     213     1794976 :     return RealEigenMatrix::Zero(_var.count(), jvar.count());
     214             : }
     215             : 
     216             : void
     217           0 : ArrayKernel::computeOffDiagJacobianScalar(unsigned int jvar)
     218             : {
     219           0 :   MooseVariableScalar & jv = _sys.getScalarVariable(_tid, jvar);
     220           0 :   prepareMatrixTag(_assembly, _var.number(), jvar);
     221             : 
     222           0 :   for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     223           0 :     for (_i = 0; _i < _test.size(); _i++)
     224             :     {
     225           0 :       _work_matrix = computeQpOffDiagJacobianScalar(jv) * _JxW[_qp] * _coord[_qp];
     226           0 :       _assembly.saveFullLocalArrayJacobian(
     227           0 :           _local_ke, _i, _test.size(), 0, 1, _var.number(), jvar, _work_matrix);
     228             :     }
     229             : 
     230           0 :   accumulateTaggedLocalMatrix();
     231           0 : }
     232             : 
     233             : RealEigenMatrix
     234           0 : ArrayKernel::computeQpOffDiagJacobianScalar(const MooseVariableScalar & jvar)
     235             : {
     236           0 :   return RealEigenMatrix::Zero(_var.count(), (unsigned int)jvar.order() + 1);
     237             : }