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       57624 : EigenKernel::validParams()
      25             : {
      26       57624 :   InputParameters params = Kernel::validParams();
      27      172872 :   params.addParam<bool>(
      28      115248 :       "eigen", true, "Use for eigenvalue problem (true) or source problem (false)");
      29      172872 :   params.addParam<PostprocessorName>(
      30      115248 :       "eigen_postprocessor", 1.0, "The name of the postprocessor that provides the eigenvalue.");
      31       57624 :   params.registerBase("EigenKernel");
      32       57624 :   return params;
      33           0 : }
      34             : 
      35         378 : EigenKernel::EigenKernel(const InputParameters & parameters)
      36             :   : Kernel(parameters),
      37         378 :     _eigen(getParam<bool>("eigen")),
      38         378 :     _eigen_sys(
      39         378 :         dynamic_cast<MooseEigenSystem *>(&_fe_problem.getNonlinearSystemBase(_sys.number()))),
      40         378 :     _eigenvalue(NULL)
      41             : {
      42             :   // The name to the postprocessor storing the eigenvalue
      43         378 :   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         378 :   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         378 :     EigenExecutionerBase * exec = dynamic_cast<EigenExecutionerBase *>(_app.getExecutioner());
      55         378 :     if (exec)
      56         338 :       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         378 :   if (eigen_pp_name.empty())
      62          40 :     _eigenvalue = &getPostprocessorValue("eigen_postprocessor");
      63             : 
      64             :   // If the name does exist, then use the postprocessor value
      65             :   else
      66             :   {
      67         338 :     if (_is_implicit)
      68         169 :       _eigenvalue = &getPostprocessorValueByName(eigen_pp_name);
      69             :     else
      70             :     {
      71         169 :       EigenExecutionerBase * exec = dynamic_cast<EigenExecutionerBase *>(_app.getExecutioner());
      72         169 :       if (exec)
      73         169 :         _eigenvalue = &exec->eigenvalueOld();
      74             :       else
      75           0 :         _eigenvalue = &getPostprocessorValueOldByName(eigen_pp_name);
      76             :     }
      77             :   }
      78         378 : }
      79             : 
      80             : void
      81      789588 : EigenKernel::computeResidual()
      82             : {
      83      789588 :   prepareVectorTag(_assembly, _var.number());
      84             : 
      85      789588 :   precalculateResidual();
      86             : 
      87             :   mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
      88      789588 :   Real one_over_eigen = 1.0 / *_eigenvalue;
      89     3947076 :   for (_i = 0; _i < _test.size(); _i++)
      90    15785712 :     for (_qp = 0; _qp < _qrule->n_points(); _qp++)
      91    12628224 :       _local_re(_i) += _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpResidual();
      92             : 
      93      789588 :   accumulateTaggedLocalResidual();
      94      789588 :   if (_has_save_in)
      95             :   {
      96           0 :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
      97           0 :     for (const auto & var : _save_in)
      98           0 :       var->sys().solution().add_vector(_local_re, var->dofIndices());
      99           0 :   }
     100      789588 : }
     101             : 
     102             : void
     103       22408 : EigenKernel::computeJacobian()
     104             : {
     105       22408 :   if (!_is_implicit)
     106           0 :     return;
     107             : 
     108       22408 :   prepareMatrixTag(_assembly, _var.number(), _var.number());
     109             : 
     110       22408 :   precalculateJacobian();
     111             :   mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
     112       22408 :   Real one_over_eigen = 1.0 / *_eigenvalue;
     113      111752 :   for (_i = 0; _i < _test.size(); _i++)
     114      446144 :     for (_j = 0; _j < _phi.size(); _j++)
     115     1782848 :       for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     116     1426048 :         _local_ke(_i, _j) += _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpJacobian();
     117             : 
     118       22408 :   accumulateTaggedLocalMatrix();
     119             : 
     120       22408 :   if (_has_diag_save_in && !_sys.computingScalingJacobian())
     121             :   {
     122         750 :     DenseVector<Number> diag = _assembly.getJacobianDiagonal(_local_ke);
     123         750 :     Threads::spin_mutex::scoped_lock lock(Threads::spin_mtx);
     124        1500 :     for (const auto & var : _diag_save_in)
     125         750 :       var->sys().solution().add_vector(diag, var->dofIndices());
     126         750 :   }
     127             : }
     128             : 
     129             : void
     130       11264 : EigenKernel::computeOffDiagJacobian(const unsigned int jvar_num)
     131             : {
     132       11264 :   if (!_is_implicit)
     133           0 :     return;
     134             : 
     135       11264 :   if (jvar_num == _var.number())
     136        8264 :     computeJacobian();
     137             :   else
     138             :   {
     139        3000 :     const auto & jvar = getVariable(jvar_num);
     140             : 
     141        3000 :     prepareMatrixTag(_assembly, _var.number(), jvar_num);
     142             : 
     143             :     // This (undisplaced) jvar could potentially yield the wrong phi size if this object is acting
     144             :     // on the displaced mesh
     145        3000 :     auto phi_size = jvar.dofIndices().size();
     146             :     mooseAssert(
     147             :         phi_size * jvar.count() == _local_ke.n(),
     148             :         "The size of the phi container does not match the number of local Jacobian columns");
     149             : 
     150        3000 :     if (_local_ke.m() != _test.size())
     151           0 :       return;
     152             : 
     153        3000 :     precalculateOffDiagJacobian(jvar_num);
     154             :     mooseAssert(*_eigenvalue != 0.0, "Can't divide by zero eigenvalue in EigenKernel!");
     155        3000 :     Real one_over_eigen = 1.0 / *_eigenvalue;
     156        3000 :     if (jvar.count() == 1)
     157             :     {
     158       15000 :       for (_i = 0; _i < _test.size(); _i++)
     159       60000 :         for (_j = 0; _j < phi_size; _j++)
     160      240000 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     161      192000 :             _local_ke(_i, _j) +=
     162      192000 :                 _JxW[_qp] * _coord[_qp] * one_over_eigen * computeQpOffDiagJacobian(jvar.number());
     163             :     }
     164             :     else
     165             :     {
     166           0 :       unsigned int n = phi_size;
     167           0 :       for (_i = 0; _i < _test.size(); _i++)
     168           0 :         for (_j = 0; _j < n; _j++)
     169           0 :           for (_qp = 0; _qp < _qrule->n_points(); _qp++)
     170             :           {
     171           0 :             RealEigenVector v = _JxW[_qp] * _coord[_qp] * one_over_eigen *
     172           0 :                                 computeQpOffDiagJacobianArray(static_cast<ArrayMooseVariable &>(
     173           0 :                                     const_cast<MooseVariableFieldBase &>(jvar)));
     174           0 :             for (unsigned int k = 0; k < v.size(); ++k)
     175           0 :               _local_ke(_i, _j + k * n) += v(k);
     176           0 :           }
     177             :     }
     178             : 
     179        3000 :     accumulateTaggedLocalMatrix();
     180             :   }
     181             : }
     182             : 
     183             : bool
     184        3452 : EigenKernel::enabled() const
     185             : {
     186        3452 :   bool flag = MooseObject::enabled();
     187        3452 :   if (_eigen)
     188             :   {
     189        3220 :     if (!_eigen_sys)
     190           0 :       mooseError("Eigen kernel ",
     191           0 :                  name(),
     192             :                  " requires a MooseEigenSystem and was designed to work with old eigenvalue",
     193             :                  " executioners such as 'NonlinearEigen'.  It is suggested to use the new",
     194             :                  " eigenvalue executioner 'Eigenvalue' along with kernel tagging");
     195             : 
     196        3220 :     if (_is_implicit)
     197        1610 :       return flag && (!_eigen_sys->activeOnOld());
     198             :     else
     199        1610 :       return flag && _eigen_sys->activeOnOld();
     200             :   }
     201             :   else
     202         232 :     return flag;
     203             : }