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 "Kernel.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "MooseVariableFE.h"
      15             : #include "MooseVariableScalar.h"
      16             : #include "SubProblem.h"
      17             : #include "NonlinearSystem.h"
      18             : #include "FEProblemBase.h"
      19             : 
      20             : #include "libmesh/threads.h"
      21             : #include "libmesh/quadrature.h"
      22             : 
      23             : InputParameters
      24      474212 : Kernel::validParams()
      25             : {
      26      474212 :   InputParameters params = KernelBase::validParams();
      27      474212 :   params.registerBase("Kernel");
      28      474212 :   return params;
      29           0 : }
      30             : 
      31       76403 : Kernel::Kernel(const InputParameters & parameters)
      32             :   : KernelBase(parameters),
      33             :     MooseVariableInterface<Real>(this,
      34             :                                  false,
      35             :                                  "variable",
      36             :                                  Moose::VarKindType::VAR_SOLVER,
      37             :                                  Moose::VarFieldType::VAR_FIELD_STANDARD),
      38      152728 :     _var(*mooseVariable()),
      39       76364 :     _test(_var.phi()),
      40       76364 :     _grad_test(_var.gradPhi()),
      41       76364 :     _phi(_assembly.phi(_var)),
      42       76364 :     _grad_phi(_assembly.gradPhi(_var)),
      43       76364 :     _u(_is_implicit ? _var.sln() : _var.slnOld()),
      44      305495 :     _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld())
      45             : {
      46       76364 :   addMooseVariableDependency(mooseVariable());
      47       76364 :   _save_in.resize(_save_in_strings.size());
      48       76364 :   _diag_save_in.resize(_diag_save_in_strings.size());
      49             : 
      50       76516 :   for (unsigned int i = 0; i < _save_in_strings.size(); i++)
      51             :   {
      52         155 :     MooseVariable * var = &_subproblem.getStandardVariable(_tid, _save_in_strings[i]);
      53             : 
      54         155 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_save_in_strings[i]))
      55           6 :       paramError("save_in", "cannot use solution variable as save-in variable");
      56             : 
      57         152 :     if (var->feType() != _var.feType())
      58           0 :       paramError(
      59             :           "save_in",
      60             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
      61           0 :           moose::internal::incompatVarMsg(*var, _var));
      62             : 
      63         152 :     _save_in[i] = var;
      64         152 :     var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
      65         152 :     addMooseVariableDependency(var);
      66             :   }
      67             : 
      68       76361 :   _has_save_in = _save_in.size() > 0;
      69             : 
      70       76465 :   for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
      71             :   {
      72         107 :     MooseVariable * var = &_subproblem.getStandardVariable(_tid, _diag_save_in_strings[i]);
      73             : 
      74         107 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_diag_save_in_strings[i]))
      75           6 :       paramError("diag_save_in", "cannot use solution variable as diag save-in variable");
      76             : 
      77         104 :     if (var->feType() != _var.feType())
      78           0 :       paramError(
      79             :           "diag_save_in",
      80             :           "saved-in auxiliary variable is incompatible with the object's nonlinear variable: ",
      81           0 :           moose::internal::incompatVarMsg(*var, _var));
      82             : 
      83         104 :     _diag_save_in[i] = var;
      84         104 :     var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
      85         104 :     addMooseVariableDependency(var);
      86             :   }
      87             : 
      88       76358 :   _has_diag_save_in = _diag_save_in.size() > 0;
      89       76358 : }
      90             : 
      91             : void
      92   433189067 : Kernel::computeResidual()
      93             : {
      94   433189067 :   if (!hasVectorTags())
      95       19464 :     return;
      96             : 
      97   433169603 :   prepareVectorTag(_assembly, _var.number());
      98             : 
      99   433169603 :   precalculateResidual();
     100  2411737307 :   for (_i = 0; _i < _test.size(); _i++)
     101 12488582297 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     102 10510014593 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
     103             : 
     104   433169564 :   accumulateTaggedLocalResidual();
     105             : 
     106   433169564 :   if (_has_save_in)
     107       10720 :     for (const auto & var : _save_in)
     108        5680 :       var->sys().solution().add_vector(_local_re, var->dofIndices());
     109             : }
     110             : 
     111             : void
     112   103795575 : Kernel::computeJacobian()
     113             : {
     114   103795575 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     115             : 
     116   103795575 :   precalculateJacobian();
     117   553449810 :   for (_i = 0; _i < _test.size(); _i++)
     118  2618667135 :     for (_j = 0; _j < _phi.size(); _j++)
     119 14883476250 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     120 12714463369 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
     121             : 
     122   103795556 :   accumulateTaggedLocalMatrix();
     123             : 
     124   103795556 :   if (_has_diag_save_in && !_sys.computingScalingJacobian())
     125             :   {
     126         636 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     127        1332 :     for (const auto & var : _diag_save_in)
     128         696 :       var->sys().solution().add_vector(diag, var->dofIndices());
     129         636 :   }
     130   103795556 : }
     131             : 
     132             : void
     133    45368856 : Kernel::computeOffDiagJacobian(const unsigned int jvar_num)
     134             : {
     135    45368856 :   if (jvar_num == _var.number())
     136    29023259 :     computeJacobian();
     137             :   else
     138             :   {
     139    16345597 :     const auto & jvar = getVariable(jvar_num);
     140    16345597 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     141             : 
     142             :     // It's possible that this variable has not been requested for coupling anywhere for our
     143             :     // associated system/assembly. E.g. in one ALE simulation, the displacement diffusion kernel
     144             :     // runs on the reference mesh and all the velocity variable couplings are on the displaced mesh
     145    16345597 :     if (jvar.dofIndices().empty())
     146      389751 :       return;
     147             : 
     148             :     // phi_size may not be equal to _phi.size(), e.g. when jvar is a vector variable
     149    15955846 :     const auto phi_size = jvar.phiSize();
     150             :     mooseAssert(
     151             :         phi_size * jvar.count() == _local_ke.n(),
     152             :         "The size of the phi container does not match the number of local Jacobian columns");
     153             : 
     154             :     mooseAssert(_local_ke.m() == _test.size(),
     155             :                 "If these are not the same, then we shouldn't even be calling this method");
     156             : 
     157    15955846 :     precalculateOffDiagJacobian(jvar_num);
     158    15955846 :     if (jvar.count() == 1)
     159             :     {
     160    79352835 :       for (_i = 0; _i < _test.size(); _i++)
     161   325867761 :         for (_j = 0; _j < phi_size; _j++)
     162  1393476434 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     163  1131006686 :             _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar.number());
     164             :     }
     165             :     else
     166             :     {
     167        1024 :       const auto n = cast_int<unsigned int>(phi_size);
     168        5120 :       for (_i = 0; _i < _test.size(); _i++)
     169       20480 :         for (_j = 0; _j < phi_size; _j++)
     170       81920 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     171             :           {
     172             :             const RealEigenVector v =
     173       65536 :                 _JxW[_qp] * _coord[_qp] *
     174      131072 :                 computeQpOffDiagJacobianArray(
     175       65536 :                     static_cast<ArrayMooseVariable &>(const_cast<MooseVariableFieldBase &>(jvar)));
     176      196608 :             for (unsigned int k = 0; k < v.size(); ++k)
     177      131072 :               _local_ke(_i, _j + k * n) += v(k);
     178       65536 :           }
     179             :     }
     180             : 
     181    15955846 :     accumulateTaggedLocalMatrix();
     182             :   }
     183             : }
     184             : 
     185             : void
     186           0 : Kernel::computeOffDiagJacobianScalar(const unsigned int jvar)
     187             : {
     188           0 :   MooseVariableScalar & jv = _sys.getScalarVariable(_tid, jvar);
     189           0 :   prepareMatrixTag(_assembly, _var.number(), jvar);
     190             : 
     191           0 :   for (_i = 0; _i < _test.size(); _i++)
     192           0 :     for (_j = 0; _j < jv.order(); _j++)
     193           0 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     194           0 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobianScalar(jvar);
     195             : 
     196           0 :   accumulateTaggedLocalMatrix();
     197           0 : }
     198             : 
     199             : void
     200         485 : Kernel::computeResidualAndJacobian()
     201             : {
     202         485 :   computeResidual();
     203             : 
     204        1369 :   for (const auto & [ivariable, jvariable] : _fe_problem.couplingEntries(_tid, _sys.number()))
     205             :   {
     206         884 :     const unsigned int ivar = ivariable->number();
     207         884 :     const unsigned int jvar = jvariable->number();
     208             : 
     209         884 :     if (ivar != _var.number())
     210         266 :       continue;
     211             : 
     212         618 :     if (_is_implicit)
     213             :     {
     214         618 :       prepareShapes(jvar);
     215         618 :       computeOffDiagJacobian(jvar);
     216             :     }
     217             :   }
     218             : 
     219             :   /// TODO: add nonlocal Jacobians and scalar Jacobians
     220         485 : }