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 "EigenKernel.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "Assembly.h"
      14             : #include "EigenExecutionerBase.h"
      15             : #include "Executioner.h"
      16             : #include "MooseApp.h"
      17             : #include "MooseEigenSystem.h"
      18             : #include "MooseVariableFE.h"
      19             : #include "FEProblemBase.h"
      20             : 
      21             : #include "libmesh/quadrature.h"
      22             : 
      23             : InputParameters
      24       12742 : EigenKernel::validParams()
      25             : {
      26       12742 :   InputParameters params = Kernel::validParams();
      27       38226 :   params.addParam<bool>(
      28       25484 :       "eigen", true, "Use for eigenvalue problem (true) or source problem (false)");
      29       38226 :   params.addParam<PostprocessorName>(
      30       25484 :       "eigen_postprocessor", 1.0, "The name of the postprocessor that provides the eigenvalue.");
      31       12742 :   params.registerBase("EigenKernel");
      32       12742 :   return params;
      33           0 : }
      34             : 
      35         339 : EigenKernel::EigenKernel(const InputParameters & parameters)
      36             :   : Kernel(parameters),
      37         339 :     _eigen(getParam<bool>("eigen")),
      38         339 :     _eigen_sys(
      39         339 :         dynamic_cast<MooseEigenSystem *>(&_fe_problem.getNonlinearSystemBase(_sys.number()))),
      40         339 :     _eigenvalue(NULL)
      41             : {
      42             :   // The name to the postprocessor storing the eigenvalue
      43         339 :   std::string eigen_pp_name;
      44             : 
      45             :   // If the "eigen_postprocessor" is given, use it. The isParamValid does not work here because of
      46             :   // the default value, which
      47             :   // you don't want to use if an EigenExecutioner exists.
      48        1017 :   if (!isDefaultPostprocessorValue("eigen_postprocessor"))
      49           0 :     eigen_pp_name = getPostprocessorName("eigen_postprocessor");
      50             : 
      51             :   // Attempt to extract the eigenvalue postprocessor from the Executioner
      52             :   else
      53             :   {
      54         339 :     EigenExecutionerBase * exec = dynamic_cast<EigenExecutionerBase *>(_app.getExecutioner());
      55         339 :     if (exec)
      56         906 :       eigen_pp_name = exec->getParam<PostprocessorName>("bx_norm");
      57             :   }
      58             : 
      59             :   // If the postprocessor name was not provided and an EigenExecutionerBase is not being used,
      60             :   // use the default value from the "eigen_postprocessor" parameter
      61         339 :   if (eigen_pp_name.empty())
      62         111 :     _eigenvalue = &getPostprocessorValue("eigen_postprocessor");
      63             : 
      64             :   // If the name does exist, then use the postprocessor value
      65             :   else
      66             :   {
      67         302 :     if (_is_implicit)
      68         151 :       _eigenvalue = &getPostprocessorValueByName(eigen_pp_name);
      69             :     else
      70             :     {
      71         151 :       EigenExecutionerBase * exec = dynamic_cast<EigenExecutionerBase *>(_app.getExecutioner());
      72         151 :       if (exec)
      73         151 :         _eigenvalue = &exec->eigenvalueOld();
      74             :       else
      75           0 :         _eigenvalue = &getPostprocessorValueOldByName(eigen_pp_name);
      76             :     }
      77             :   }
      78         339 : }
      79             : 
      80             : void
      81      695180 : EigenKernel::computeResidual()
      82             : {
      83      695180 :   prepareVectorTag(_assembly, _var.number());
      84             : 
      85      695180 :   precalculateResidual();
      86             : 
      87             :   mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
      88      695180 :   Real one_over_eigen = 1.0 / *_eigenvalue;
      89     3475132 :   for (_i = 0; _i < _test.size(); _i++)
      90    13898224 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      91    11118272 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpResidual();
      92             : 
      93      695180 :   accumulateTaggedLocalResidual();
      94      695180 :   if (_has_save_in)
      95           0 :     for (const auto & var : _save_in)
      96           0 :       var->sys().solution().add_vector(_local_re, var->dofIndices());
      97      695180 : }
      98             : 
      99             : void
     100       18372 : EigenKernel::computeJacobian()
     101             : {
     102       18372 :   if (!_is_implicit)
     103           0 :     return;
     104             : 
     105       18372 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     106             : 
     107       18372 :   precalculateJacobian();
     108             :   mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
     109       18372 :   Real one_over_eigen = 1.0 / *_eigenvalue;
     110       91604 :   for (_i = 0; _i < _test.size(); _i++)
     111      365648 :     for (_j = 0; _j < _phi.size(); _j++)
     112     1461056 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     113     1168640 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpJacobian();
     114             : 
     115       18372 :   accumulateTaggedLocalMatrix();
     116             : 
     117       18372 :   if (_has_diag_save_in && !_sys.computingScalingJacobian())
     118             :   {
     119         675 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     120        1350 :     for (const auto & var : _diag_save_in)
     121         675 :       var->sys().solution().add_vector(diag, var->dofIndices());
     122         675 :   }
     123             : }
     124             : 
     125             : void
     126        9368 : EigenKernel::computeOffDiagJacobian(const unsigned int jvar_num)
     127             : {
     128        9368 :   if (!_is_implicit)
     129           0 :     return;
     130             : 
     131        9368 :   if (jvar_num == _var.number())
     132        6668 :     computeJacobian();
     133             :   else
     134             :   {
     135        2700 :     const auto & jvar = getVariable(jvar_num);
     136             : 
     137        2700 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     138             : 
     139             :     // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting
     140             :     // on the displaced mesh
     141        2700 :     auto phi_size = jvar.dofIndices().size();
     142             :     mooseAssert(
     143             :         phi_size * jvar.count() == _local_ke.n(),
     144             :         "The size of the phi container does not match the number of local Jacobian columns");
     145             : 
     146        2700 :     if (_local_ke.m() != _test.size())
     147           0 :       return;
     148             : 
     149        2700 :     precalculateOffDiagJacobian(jvar_num);
     150             :     mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
     151        2700 :     Real one_over_eigen = 1.0 / *_eigenvalue;
     152        2700 :     if (jvar.count() == 1)
     153             :     {
     154       13500 :       for (_i = 0; _i < _test.size(); _i++)
     155       54000 :         for (_j = 0; _j < phi_size; _j++)
     156      216000 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     157      172800 :             _local_ke(_i, _j) +=
     158      172800 :                 _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpOffDiagJacobian(jvar.number());
     159             :     }
     160             :     else
     161             :     {
     162           0 :       unsigned int n = phi_size;
     163           0 :       for (_i = 0; _i < _test.size(); _i++)
     164           0 :         for (_j = 0; _j < n; _j++)
     165           0 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     166             :           {
     167           0 :             RealEigenVector v = _JxW[_qp] * _coord[_qp] * one_over_eigen *
     168           0 :                                 computeQpOffDiagJacobianArray(static_cast<ArrayMooseVariable &>(
     169           0 :                                     const_cast<MooseVariableFieldBase &>(jvar)));
     170           0 :             for (unsigned int k = 0; k < v.size(); ++k)
     171           0 :               _local_ke(_i, _j + k * n) += v(k);
     172           0 :           }
     173             :     }
     174             : 
     175        2700 :     accumulateTaggedLocalMatrix();
     176             :   }
     177             : }
     178             : 
     179             : bool
     180        3082 : EigenKernel::enabled() const
     181             : {
     182        3082 :   bool flag = MooseObject::enabled();
     183        3082 :   if (_eigen)
     184             :   {
     185        2868 :     if (!_eigen_sys)
     186           0 :       mooseError("Eigen kernel ",
     187           0 :                  name(),
     188             :                  " requires a MooseEigenSystem and was designed to work with old eigenvalue",
     189             :                  " executioners such as 'NonlinearEigen'.  It is suggested to use the new",
     190             :                  " eigenvalue executioner 'Eigenvalue' along with kernel tagging");
     191             : 
     192        2868 :     if (_is_implicit)
     193        1434 :       return flag && (!_eigen_sys->activeOnOld());
     194             :     else
     195        1434 :       return flag && _eigen_sys->activeOnOld();
     196             :   }
     197             :   else
     198         214 :     return flag;
     199             : }