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     1686219 : Kernel::validParams()
      25             : {
      26     1686219 :   InputParameters params = KernelBase::validParams();
      27     1686219 :   params.registerBase("Kernel");
      28     1686219 :   return params;
      29           0 : }
      30             : 
      31       82658 : 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      165212 :     _var(*mooseVariable()),
      39       82606 :     _test(_var.phi()),
      40       82606 :     _grad_test(_var.gradPhi()),
      41       82606 :     _phi(_assembly.phi(_var)),
      42       82606 :     _grad_phi(_assembly.gradPhi(_var)),
      43       82606 :     _u(_is_implicit ? _var.sln() : _var.slnOld()),
      44      165264 :     _grad_u(_is_implicit ? _var.gradSln() : _var.gradSlnOld())
      45             : {
      46       82606 :   addMooseVariableDependency(mooseVariable());
      47       82606 :   _save_in.resize(_save_in_strings.size());
      48       82606 :   _diag_save_in.resize(_diag_save_in_strings.size());
      49             : 
      50       82772 :   for (unsigned int i = 0; i < _save_in_strings.size(); i++)
      51             :   {
      52         170 :     MooseVariable * var = &_subproblem.getStandardVariable(_tid, _save_in_strings[i]);
      53             : 
      54         170 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_save_in_strings[i]))
      55           4 :       paramError("save_in", "cannot use solution variable as save-in variable");
      56             : 
      57         166 :     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         166 :     _save_in[i] = var;
      64         166 :     var->sys().addVariableToZeroOnResidual(_save_in_strings[i]);
      65         166 :     addMooseVariableDependency(var);
      66             :   }
      67             : 
      68       82602 :   _has_save_in = _save_in.size() > 0;
      69             : 
      70       82714 :   for (unsigned int i = 0; i < _diag_save_in_strings.size(); i++)
      71             :   {
      72         116 :     MooseVariable * var = &_subproblem.getStandardVariable(_tid, _diag_save_in_strings[i]);
      73             : 
      74         116 :     if (_fe_problem.getNonlinearSystemBase(_sys.number()).hasVariable(_diag_save_in_strings[i]))
      75           4 :       paramError("diag_save_in", "cannot use solution variable as diag save-in variable");
      76             : 
      77         112 :     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         112 :     _diag_save_in[i] = var;
      84         112 :     var->sys().addVariableToZeroOnJacobian(_diag_save_in_strings[i]);
      85         112 :     addMooseVariableDependency(var);
      86             :   }
      87             : 
      88       82598 :   _has_diag_save_in = _diag_save_in.size() > 0;
      89       82598 : }
      90             : 
      91             : void
      92   488748355 : Kernel::computeResidual()
      93             : {
      94   488748355 :   if (!hasVectorTags())
      95          72 :     return;
      96             : 
      97   488748283 :   prepareVectorTag(_assembly, _var.number());
      98             : 
      99   488748283 :   precalculateResidual();
     100  2709212744 :   for (_i = 0; _i < _test.size(); _i++)
     101 13903854642 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     102 11683390181 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * computeQpResidual();
     103             : 
     104   488748239 :   accumulateTaggedLocalResidual();
     105             : 
     106   488748239 :   if (_has_save_in)
     107             :   {
     108        5670 :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     109       12060 :     for (const auto & var : _save_in)
     110        6390 :       var->sys().solution().add_vector(_local_re, var->dofIndices());
     111        5670 :   }
     112             : }
     113             : 
     114             : void
     115   120120151 : Kernel::computeJacobian()
     116             : {
     117   120120151 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     118             : 
     119   120120151 :   precalculateJacobian();
     120   638239217 :   for (_i = 0; _i < _test.size(); _i++)
     121  2995772069 :     for (_j = 0; _j < _phi.size(); _j++)
     122 16743841951 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     123 14266188969 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpJacobian();
     124             : 
     125   120120130 :   accumulateTaggedLocalMatrix();
     126             : 
     127   120120130 :   if (_has_diag_save_in && !_sys.computingScalingJacobian())
     128             :   {
     129         716 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     130         716 :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     131        1500 :     for (const auto & var : _diag_save_in)
     132         784 :       var->sys().solution().add_vector(diag, var->dofIndices());
     133         716 :   }
     134   120120130 : }
     135             : 
     136             : void
     137    52723324 : Kernel::computeOffDiagJacobian(const unsigned int jvar_num)
     138             : {
     139    52723324 :   if (jvar_num == _var.number())
     140    33707969 :     computeJacobian();
     141             :   else
     142             :   {
     143    19015355 :     const auto & jvar = getVariable(jvar_num);
     144    19015355 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     145             : 
     146             :     // phi_size may not be equal to _phi.size(), e.g. when jvar is a vector variable
     147    19015355 :     const auto phi_size = jvar.dofIndices().size();
     148             :     mooseAssert(
     149             :         phi_size * jvar.count() == _local_ke.n(),
     150             :         "The size of the phi container does not match the number of local Jacobian columns");
     151             : 
     152             :     mooseAssert(_local_ke.m() == _test.size(),
     153             :                 "If these are not the same, then we shouldn't even be calling this method");
     154             : 
     155    19015355 :     precalculateOffDiagJacobian(jvar_num);
     156    19015355 :     if (jvar.count() == 1)
     157             :     {
     158    93256525 :       for (_i = 0; _i < _test.size(); _i++)
     159   374242384 :         for (_j = 0; _j < phi_size; _j++)
     160  1523418544 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     161  1223418482 :             _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobian(jvar.number());
     162             :     }
     163             :     else
     164             :     {
     165        1152 :       const auto n = cast_int<unsigned int>(phi_size);
     166        5760 :       for (_i = 0; _i < _test.size(); _i++)
     167       23040 :         for (_j = 0; _j < phi_size; _j++)
     168       92160 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     169             :           {
     170             :             const RealEigenVector v =
     171       73728 :                 _JxW[_qp] * _coord[_qp] *
     172      147456 :                 computeQpOffDiagJacobianArray(
     173       73728 :                     static_cast<ArrayMooseVariable &>(const_cast<MooseVariableFieldBase &>(jvar)));
     174      221184 :             for (unsigned int k = 0; k < v.size(); ++k)
     175      147456 :               _local_ke(_i, _j + k * n) += v(k);
     176       73728 :           }
     177             :     }
     178             : 
     179    19015355 :     accumulateTaggedLocalMatrix();
     180             :   }
     181    52723324 : }
     182             : 
     183             : void
     184           0 : Kernel::computeOffDiagJacobianScalar(const unsigned int jvar)
     185             : {
     186           0 :   MooseVariableScalar & jv = _sys.getScalarVariable(_tid, jvar);
     187           0 :   prepareMatrixTag(_assembly, _var.number(), jvar);
     188             : 
     189           0 :   for (_i = 0; _i < _test.size(); _i++)
     190           0 :     for (_j = 0; _j < jv.order(); _j++)
     191           0 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     192           0 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * computeQpOffDiagJacobianScalar(jvar);
     193             : 
     194           0 :   accumulateTaggedLocalMatrix();
     195           0 : }
     196             : 
     197             : void
     198         546 : Kernel::computeResidualAndJacobian()
     199             : {
     200         546 :   computeResidual();
     201             : 
     202        1533 :   for (const auto & [ivariable, jvariable] : _fe_problem.couplingEntries(_tid, _sys.number()))
     203             :   {
     204         987 :     const unsigned int ivar = ivariable->number();
     205         987 :     const unsigned int jvar = jvariable->number();
     206             : 
     207         987 :     if (ivar != _var.number())
     208         294 :       continue;
     209             : 
     210         693 :     if (_is_implicit)
     211             :     {
     212         693 :       prepareShapes(jvar);
     213         693 :       computeOffDiagJacobian(jvar);
     214             :     }
     215             :   }
     216             : 
     217             :   /// TODO: add nonlocal Jacobians and scalar Jacobians
     218         546 : }