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