ParsedODEKernel.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://mooseframework.inl.gov
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "ParsedODEKernel.h"
 
// MOOSE includes
#include "MooseVariableScalar.h"
#include "SystemBase.h"
#include "MooseApp.h"
 
#include "libmesh/fparser_ad.hh"
 
registerMooseObject("MooseApp", ParsedODEKernel);
 
InputParameters
ParsedODEKernel::validParams()
{
  InputParameters params = ODEKernel::validParams();
  params += FunctionParserUtils<false>::validParams();
  params.addClassDescription("Parsed expression ODE kernel.");
 
  params.addRequiredCustomTypeParam<std::string>(
      "expression", "FunctionExpression", "function expression");
  params.addCoupledVar("coupled_variables", "Scalar variables coupled in the parsed expression.");
  params.addParam<std::vector<std::string>>(
      "constant_names", {}, "Vector of constants used in the parsed expression");
  params.addParam<std::vector<std::string>>(
      "constant_expressions",
      {},
      "Vector of values for the constants in constant_names (can be an FParser expression)");
  params.addParam<std::vector<PostprocessorName>>(
      "postprocessors", {}, "Vector of postprocessor names used in the function expression");
 
  return params;
}
 
ParsedODEKernel::ParsedODEKernel(const InputParameters & parameters)
  : ODEKernel(parameters),
    FunctionParserUtils(parameters),
    _function(getParam<std::string>("expression")),
    _nargs(isCoupledScalar("coupled_variables") ? coupledScalarComponents("coupled_variables") : 0),
    _args(_nargs),
    _arg_names(_nargs),
    _func_dFdarg(_nargs),
    _number_of_nl_variables(_sys.nVariables()),
    _arg_index(_number_of_nl_variables, -1)
{
  // build variables argument (start with variable the kernel is operating on)
  std::string variables = _var.name();
 
  // add additional coupled variables
  for (unsigned int i = 0; i < _nargs; ++i)
  {
    _arg_names[i] = getScalarVar("coupled_variables", i)->name();
    variables += "," + _arg_names[i];
    _args[i] = &coupledScalarValue("coupled_variables", i);
 
    // populate number -> arg index lookup table skipping aux variables
    unsigned int number = coupledScalar("coupled_variables", i);
    if (number < _number_of_nl_variables)
      _arg_index[number] = i;
  }
 
  // add postprocessors
  auto pp_names = getParam<std::vector<PostprocessorName>>("postprocessors");
  _pp.resize(pp_names.size());
  for (unsigned int i = 0; i < pp_names.size(); ++i)
  {
    variables += "," + pp_names[i];
    _pp[i] = &getPostprocessorValueByName(pp_names[i]);
  }
 
  // Note: we do not use the FunctionParsedUtils::parsedFunctionSetup because we are building
  // multiple expressions, and we can share the MPI barriers by keeping the code here.
  // base function object
  _func_F = std::make_shared<SymFunction>();
 
  // set FParser interneal feature flags
  setParserFeatureFlags(_func_F);
 
  // add the constant expressions
  addFParserConstants(_func_F,
                      getParam<std::vector<std::string>>("constant_names"),
                      getParam<std::vector<std::string>>("constant_expressions"));
 
  // parse function
  if (_func_F->Parse(_function, variables) >= 0)
    mooseError("Invalid function\n",
               _function,
               "\nin ParsedODEKernel ",
               name(),
               ".\n",
               _func_F->ErrorMsg());
 
  // on-diagonal derivative
  _func_dFdu = std::make_shared<SymFunction>(*_func_F);
 
  // let rank 0 do the work first (diff and JIT) to populate caches
  if (_communicator.rank() != 0)
    _communicator.barrier();
 
  if (_func_dFdu->AutoDiff(_var.name()) != -1)
    mooseError("Failed to take first derivative w.r.t. ", _var.name());
 
  // off-diagonal derivatives
  for (unsigned int i = 0; i < _nargs; ++i)
  {
    _func_dFdarg[i] = std::make_shared<SymFunction>(*_func_F);
 
    if (_func_dFdarg[i]->AutoDiff(_arg_names[i]) != -1)
      mooseError("Failed to take first derivative w.r.t. ", _arg_names[i]);
  }
 
  // optimize
  if (!_disable_fpoptimizer)
  {
    _func_F->Optimize();
    _func_dFdu->Optimize();
    for (unsigned int i = 0; i < _nargs; ++i)
      _func_dFdarg[i]->Optimize();
  }
 
  // just-in-time compile
  if (_enable_jit)
  {
    _func_F->JITCompile();
    _func_dFdu->JITCompile();
    for (unsigned int i = 0; i < _nargs; ++i)
      _func_dFdarg[i]->JITCompile();
  }
 
  // wait for ranks > 0 to catch up
  if (_communicator.rank() == 0)
    _communicator.barrier();
 
  // reserve storage for parameter passing buffer
  _func_params.resize(1 + _nargs + _pp.size());
}
 
void
ParsedODEKernel::updateParams()
{
  _func_params[0] = _u[_i];
 
  for (unsigned int j = 0; j < _nargs; ++j)
    _func_params[j + 1] = (*_args[j])[_i];
  for (unsigned int j = 0; j < _pp.size(); ++j)
    _func_params[j + 1 + _nargs] = *_pp[j];
}
 
Real
ParsedODEKernel::computeQpResidual()
{
  updateParams();
  return evaluate(_func_F);
}
 
Real
ParsedODEKernel::computeQpJacobian()
{
  updateParams();
  return evaluate(_func_dFdu);
}
 
Real
ParsedODEKernel::computeQpOffDiagJacobianScalar(unsigned int jvar)
{
  int i = _arg_index[jvar];
  if (i < 0)
    return 0.0;
 
  updateParams();
  return evaluate(_func_dFdarg[i]);
}