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 "ParsedODEKernel.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseVariableScalar.h"
      14             : #include "SystemBase.h"
      15             : #include "MooseApp.h"
      16             : 
      17             : #include "libmesh/fparser_ad.hh"
      18             : 
      19             : registerMooseObject("MooseApp", ParsedODEKernel);
      20             : 
      21             : InputParameters
      22       15423 : ParsedODEKernel::validParams()
      23             : {
      24       15423 :   InputParameters params = ODEKernel::validParams();
      25       15423 :   params += FunctionParserUtils<false>::validParams();
      26       15423 :   params.addClassDescription("Parsed expression ODE kernel.");
      27             : 
      28       15423 :   params.addRequiredCustomTypeParam<std::string>(
      29             :       "function", "FunctionExpression", "function expression");
      30       15423 :   params.deprecateParam("function", "expression", "02/07/2024");
      31       15423 :   params.addCoupledVar("args", "Scalar variables coupled in the parsed expression.");
      32       15423 :   params.deprecateCoupledVar("args", "coupled_variables", "02/07/2024");
      33       15423 :   params.addParam<std::vector<std::string>>(
      34             :       "constant_names", {}, "Vector of constants used in the parsed expression");
      35       15423 :   params.addParam<std::vector<std::string>>(
      36             :       "constant_expressions",
      37             :       {},
      38             :       "Vector of values for the constants in constant_names (can be an FParser expression)");
      39       15423 :   params.addParam<std::vector<PostprocessorName>>(
      40             :       "postprocessors", {}, "Vector of postprocessor names used in the function expression");
      41             : 
      42       15423 :   return params;
      43           0 : }
      44             : 
      45         575 : ParsedODEKernel::ParsedODEKernel(const InputParameters & parameters)
      46             :   : ODEKernel(parameters),
      47             :     FunctionParserUtils(parameters),
      48         575 :     _function(getParam<std::string>("expression")),
      49         575 :     _nargs(isCoupledScalar("coupled_variables") ? coupledScalarComponents("coupled_variables") : 0),
      50         575 :     _args(_nargs),
      51         575 :     _arg_names(_nargs),
      52         575 :     _func_dFdarg(_nargs),
      53         575 :     _number_of_nl_variables(_sys.nVariables()),
      54        1725 :     _arg_index(_number_of_nl_variables, -1)
      55             : {
      56             :   // build variables argument (start with variable the kernel is operating on)
      57         575 :   std::string variables = _var.name();
      58             : 
      59             :   // add additional coupled variables
      60         823 :   for (unsigned int i = 0; i < _nargs; ++i)
      61             :   {
      62         248 :     _arg_names[i] = getScalarVar("coupled_variables", i)->name();
      63         248 :     variables += "," + _arg_names[i];
      64         248 :     _args[i] = &coupledScalarValue("coupled_variables", i);
      65             : 
      66             :     // populate number -> arg index lookup table skipping aux variables
      67         248 :     unsigned int number = coupledScalar("coupled_variables", i);
      68         248 :     if (number < _number_of_nl_variables)
      69         248 :       _arg_index[number] = i;
      70             :   }
      71             : 
      72             :   // add postprocessors
      73         575 :   auto pp_names = getParam<std::vector<PostprocessorName>>("postprocessors");
      74         575 :   _pp.resize(pp_names.size());
      75         588 :   for (unsigned int i = 0; i < pp_names.size(); ++i)
      76             :   {
      77          13 :     variables += "," + pp_names[i];
      78          13 :     _pp[i] = &getPostprocessorValueByName(pp_names[i]);
      79             :   }
      80             : 
      81             :   // Note: we do not use the FunctionParsedUtils::parsedFunctionSetup because we are building
      82             :   // multiple expressions, and we can share the MPI barriers by keeping the code here.
      83             :   // base function object
      84         575 :   _func_F = std::make_shared<SymFunction>();
      85             : 
      86             :   // set FParser interneal feature flags
      87         575 :   setParserFeatureFlags(_func_F);
      88             : 
      89             :   // add the constant expressions
      90         575 :   addFParserConstants(_func_F,
      91             :                       getParam<std::vector<std::string>>("constant_names"),
      92             :                       getParam<std::vector<std::string>>("constant_expressions"));
      93             : 
      94             :   // parse function
      95         575 :   if (_func_F->Parse(_function, variables) >= 0)
      96           0 :     mooseError("Invalid function\n",
      97           0 :                _function,
      98             :                "\nin ParsedODEKernel ",
      99           0 :                name(),
     100             :                ".\n",
     101           0 :                _func_F->ErrorMsg());
     102             : 
     103             :   // on-diagonal derivative
     104         575 :   _func_dFdu = std::make_shared<SymFunction>(*_func_F);
     105             : 
     106             :   // let rank 0 do the work first (diff and JIT) to populate caches
     107         575 :   if (_communicator.rank() != 0)
     108         162 :     _communicator.barrier();
     109             : 
     110         575 :   if (_func_dFdu->AutoDiff(_var.name()) != -1)
     111           0 :     mooseError("Failed to take first derivative w.r.t. ", _var.name());
     112             : 
     113             :   // off-diagonal derivatives
     114         823 :   for (unsigned int i = 0; i < _nargs; ++i)
     115             :   {
     116         248 :     _func_dFdarg[i] = std::make_shared<SymFunction>(*_func_F);
     117             : 
     118         248 :     if (_func_dFdarg[i]->AutoDiff(_arg_names[i]) != -1)
     119           0 :       mooseError("Failed to take first derivative w.r.t. ", _arg_names[i]);
     120             :   }
     121             : 
     122             :   // optimize
     123         575 :   if (!_disable_fpoptimizer)
     124             :   {
     125         575 :     _func_F->Optimize();
     126         575 :     _func_dFdu->Optimize();
     127         823 :     for (unsigned int i = 0; i < _nargs; ++i)
     128         248 :       _func_dFdarg[i]->Optimize();
     129             :   }
     130             : 
     131             :   // just-in-time compile
     132         575 :   if (_enable_jit)
     133             :   {
     134         575 :     _func_F->JITCompile();
     135         575 :     _func_dFdu->JITCompile();
     136         823 :     for (unsigned int i = 0; i < _nargs; ++i)
     137         248 :       _func_dFdarg[i]->JITCompile();
     138             :   }
     139             : 
     140             :   // wait for ranks > 0 to catch up
     141         575 :   if (_communicator.rank() == 0)
     142         413 :     _communicator.barrier();
     143             : 
     144             :   // reserve storage for parameter passing buffer
     145         575 :   _func_params.resize(1 + _nargs + _pp.size());
     146         575 : }
     147             : 
     148             : void
     149       83922 : ParsedODEKernel::updateParams()
     150             : {
     151       83922 :   _func_params[0] = _u[_i];
     152             : 
     153      111738 :   for (unsigned int j = 0; j < _nargs; ++j)
     154       27816 :     _func_params[j + 1] = (*_args[j])[_i];
     155       84195 :   for (unsigned int j = 0; j < _pp.size(); ++j)
     156         273 :     _func_params[j + 1 + _nargs] = *_pp[j];
     157       83922 : }
     158             : 
     159             : Real
     160       63957 : ParsedODEKernel::computeQpResidual()
     161             : {
     162       63957 :   updateParams();
     163       63957 :   return evaluate(_func_F);
     164             : }
     165             : 
     166             : Real
     167       15790 : ParsedODEKernel::computeQpJacobian()
     168             : {
     169       15790 :   updateParams();
     170       15790 :   return evaluate(_func_dFdu);
     171             : }
     172             : 
     173             : Real
     174        5281 : ParsedODEKernel::computeQpOffDiagJacobianScalar(unsigned int jvar)
     175             : {
     176        5281 :   int i = _arg_index[jvar];
     177        5281 :   if (i < 0)
     178        1106 :     return 0.0;
     179             : 
     180        4175 :   updateParams();
     181        4175 :   return evaluate(_func_dFdarg[i]);
     182             : }