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       16051 : ParsedODEKernel::validParams()
      23             : {
      24       16051 :   InputParameters params = ODEKernel::validParams();
      25       16051 :   params += FunctionParserUtils<false>::validParams();
      26       32102 :   params.addClassDescription("Parsed expression ODE kernel.");
      27             : 
      28       96306 :   params.addRequiredCustomTypeParam<std::string>(
      29             :       "expression", "FunctionExpression", "function expression");
      30       64204 :   params.addCoupledVar("coupled_variables", "Scalar variables coupled in the parsed expression.");
      31       64204 :   params.addParam<std::vector<std::string>>(
      32             :       "constant_names", {}, "Vector of constants used in the parsed expression");
      33       64204 :   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       48153 :   params.addParam<std::vector<PostprocessorName>>(
      38             :       "postprocessors", {}, "Vector of postprocessor names used in the function expression");
      39             : 
      40       16051 :   return params;
      41           0 : }
      42             : 
      43         588 : ParsedODEKernel::ParsedODEKernel(const InputParameters & parameters)
      44             :   : ODEKernel(parameters),
      45             :     FunctionParserUtils(parameters),
      46         588 :     _function(getParam<std::string>("expression")),
      47        1530 :     _nargs(isCoupledScalar("coupled_variables") ? coupledScalarComponents("coupled_variables") : 0),
      48        1176 :     _args(_nargs),
      49        1176 :     _arg_names(_nargs),
      50        1176 :     _func_dFdarg(_nargs),
      51         588 :     _number_of_nl_variables(_sys.nVariables()),
      52        2352 :     _arg_index(_number_of_nl_variables, -1)
      53             : {
      54             :   // build variables argument (start with variable the kernel is operating on)
      55         588 :   std::string variables = _var.name();
      56             : 
      57             :   // add additional coupled variables
      58         836 :   for (unsigned int i = 0; i < _nargs; ++i)
      59             :   {
      60         496 :     _arg_names[i] = getScalarVar("coupled_variables", i)->name();
      61         248 :     variables += "," + _arg_names[i];
      62         496 :     _args[i] = &coupledScalarValue("coupled_variables", i);
      63             : 
      64             :     // populate number -> arg index lookup table skipping aux variables
      65         496 :     unsigned int number = coupledScalar("coupled_variables", i);
      66         248 :     if (number < _number_of_nl_variables)
      67         248 :       _arg_index[number] = i;
      68             :   }
      69             : 
      70             :   // add postprocessors
      71        1176 :   auto pp_names = getParam<std::vector<PostprocessorName>>("postprocessors");
      72         588 :   _pp.resize(pp_names.size());
      73         601 :   for (unsigned int i = 0; i < pp_names.size(); ++i)
      74             :   {
      75          13 :     variables += "," + pp_names[i];
      76          13 :     _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         588 :   _func_F = std::make_shared<SymFunction>();
      83             : 
      84             :   // set FParser interneal feature flags
      85         588 :   setParserFeatureFlags(_func_F);
      86             : 
      87             :   // add the constant expressions
      88        2940 :   addFParserConstants(_func_F,
      89             :                       getParam<std::vector<std::string>>("constant_names"),
      90             :                       getParam<std::vector<std::string>>("constant_expressions"));
      91             : 
      92             :   // parse function
      93         588 :   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         588 :   _func_dFdu = std::make_shared<SymFunction>(*_func_F);
     103             : 
     104             :   // let rank 0 do the work first (diff and JIT) to populate caches
     105         588 :   if (_communicator.rank() != 0)
     106         165 :     _communicator.barrier();
     107             : 
     108         588 :   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         836 :   for (unsigned int i = 0; i < _nargs; ++i)
     113             :   {
     114         248 :     _func_dFdarg[i] = std::make_shared<SymFunction>(*_func_F);
     115             : 
     116         248 :     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         588 :   if (!_disable_fpoptimizer)
     122             :   {
     123         588 :     _func_F->Optimize();
     124         588 :     _func_dFdu->Optimize();
     125         836 :     for (unsigned int i = 0; i < _nargs; ++i)
     126         248 :       _func_dFdarg[i]->Optimize();
     127             :   }
     128             : 
     129             :   // just-in-time compile
     130         588 :   if (_enable_jit)
     131             :   {
     132         588 :     _func_F->JITCompile();
     133         588 :     _func_dFdu->JITCompile();
     134         836 :     for (unsigned int i = 0; i < _nargs; ++i)
     135         248 :       _func_dFdarg[i]->JITCompile();
     136             :   }
     137             : 
     138             :   // wait for ranks > 0 to catch up
     139         588 :   if (_communicator.rank() == 0)
     140         423 :     _communicator.barrier();
     141             : 
     142             :   // reserve storage for parameter passing buffer
     143         588 :   _func_params.resize(1 + _nargs + _pp.size());
     144         588 : }
     145             : 
     146             : void
     147       84006 : ParsedODEKernel::updateParams()
     148             : {
     149       84006 :   _func_params[0] = _u[_i];
     150             : 
     151      111822 :   for (unsigned int j = 0; j < _nargs; ++j)
     152       27816 :     _func_params[j + 1] = (*_args[j])[_i];
     153       84279 :   for (unsigned int j = 0; j < _pp.size(); ++j)
     154         273 :     _func_params[j + 1 + _nargs] = *_pp[j];
     155       84006 : }
     156             : 
     157             : Real
     158       64013 : ParsedODEKernel::computeQpResidual()
     159             : {
     160       64013 :   updateParams();
     161      192039 :   return evaluate(_func_F);
     162             : }
     163             : 
     164             : Real
     165       15818 : ParsedODEKernel::computeQpJacobian()
     166             : {
     167       15818 :   updateParams();
     168       47454 :   return evaluate(_func_dFdu);
     169             : }
     170             : 
     171             : Real
     172        5281 : ParsedODEKernel::computeQpOffDiagJacobianScalar(unsigned int jvar)
     173             : {
     174        5281 :   int i = _arg_index[jvar];
     175        5281 :   if (i < 0)
     176        1106 :     return 0.0;
     177             : 
     178        4175 :   updateParams();
     179       12525 :   return evaluate(_func_dFdarg[i]);
     180             : }