Supporting both AD and non-AD variables through templating

Materials

Consider the GenericConstantMaterial.md material. The header file defines a templated class GenericConstantMaterialTempl<is_ad> with a material property of GenericMaterialProperty<Real, is_ad which evaluates to the correct type depending on the value of is_ad.

#pragma once

#include "Material.h"

/**
 * This material automatically declares as material properties whatever is passed to it
 * through the parameters 'prop_names' and uses the values from 'prop_values' as the values
 * for those properties.
 *
 * This is not meant to be used in a production capacity... and instead is meant to be used
 * during development phases for ultimate flexibility.
 */
template <bool is_ad>
class GenericConstantMaterialTempl : public Material
{
public:
  static InputParameters validParams();

  GenericConstantMaterialTempl(const InputParameters & parameters);

protected:
  virtual void initQpStatefulProperties() override;
  virtual void computeQpProperties() override;

  const std::vector<std::string> & _prop_names;
  const std::vector<Real> & _prop_values;

  unsigned int _num_props;

  std::vector<GenericMaterialProperty<Real, is_ad> *> _properties;
};

typedef GenericConstantMaterialTempl<false> GenericConstantMaterial;
typedef GenericConstantMaterialTempl<true> ADGenericConstantMaterial;

Note the typedef's at the end of the header: when GenericConstantMaterial is used in an input file, this class is instantiated with is_ad = false, while when ADGenericConstantMaterial is used, is_ad = true.

The corresponding source file with templated methods is

#include "GenericConstantMaterial.h"

registerMooseObject("MooseApp", GenericConstantMaterial);
registerMooseObject("MooseApp", ADGenericConstantMaterial);

template <bool is_ad>
InputParameters
GenericConstantMaterialTempl<is_ad>::validParams()
{
  InputParameters params = Material::validParams();
  params.addClassDescription(
      "Declares material properties based on names and values prescribed by input parameters.");
  params.addRequiredParam<std::vector<std::string>>(
      "prop_names", "The names of the properties this material will have");
  params.addRequiredParam<std::vector<Real>>("prop_values",
                                             "The values associated with the named properties");
  params.set<MooseEnum>("constant_on") = "SUBDOMAIN";
  params.declareControllable("prop_values");
  return params;
}

template <bool is_ad>
GenericConstantMaterialTempl<is_ad>::GenericConstantMaterialTempl(
    const InputParameters & parameters)
  : Material(parameters),
    _prop_names(getParam<std::vector<std::string>>("prop_names")),
    _prop_values(getParam<std::vector<Real>>("prop_values"))
{
  unsigned int num_names = _prop_names.size();
  unsigned int num_values = _prop_values.size();

  if (num_names != num_values)
    mooseError(
        "Number of prop_names must match the number of prop_values for a GenericConstantMaterial!");

  _num_props = num_names;

  _properties.resize(num_names);

  for (unsigned int i = 0; i < _num_props; i++)
    _properties[i] = &declareGenericProperty<Real, is_ad>(_prop_names[i]);
}

template <bool is_ad>
void
GenericConstantMaterialTempl<is_ad>::initQpStatefulProperties()
{
  computeQpProperties();
}

template <bool is_ad>
void
GenericConstantMaterialTempl<is_ad>::computeQpProperties()
{
  for (unsigned int i = 0; i < _num_props; i++)
    (*_properties[i])[_qp] = _prop_values[i];
}

template class GenericConstantMaterialTempl<false>;
template class GenericConstantMaterialTempl<true>;

Note that both GenericConstantMaterial and ADGenericConstantMaterial are registered to the app, and the material property is declared using the templated declareGenericMaterialProperty<T, is_ad>(mat_prop).

Kernels

Other classes can be templated in a similar fashion. Consider a Kernel example. In this case, the class is derived from the GenericKernel base class, to ensure that the computeQpResidual() method returns the correct type for both AD and non-AD variables (note that it is of type GenericReal<is_ad>). Otherwise, the basic concept is the same as before.

#pragma once

#include "GenericKernel.h"

class Function;

/**
 * This kernel implements a generic functional
 * body force term:
 * $ - c \cdof f \cdot \phi_i $
 *
 * The coefficient and function both have defaults
 * equal to 1.0.
 */
template <bool is_ad>
class BodyForceTempl : public GenericKernel<is_ad>
{
public:
  static InputParameters validParams();

  BodyForceTempl(const InputParameters & parameters);

protected:
  virtual GenericReal<is_ad> computeQpResidual() override;

  /// Scale factor
  const Real & _scale;

  /// Optional function value
  const Function & _function;

  /// Optional Postprocessor value
  const PostprocessorValue & _postprocessor;

  // AD/non-AD version of the quadrature point coordinates
  const MooseArray<Moose::GenericType<Point, is_ad>> * _generic_q_point;

  usingGenericKernelMembers;
};

typedef BodyForceTempl<false> BodyForce;
typedef BodyForceTempl<true> ADBodyForce;

The corresponding source file is

#include "BodyForce.h"

// MOOSE
#include "Function.h"
#include "Assembly.h"

registerMooseObject("MooseApp", BodyForce);
registerMooseObject("MooseApp", ADBodyForce);

template <bool is_ad>
InputParameters
BodyForceTempl<is_ad>::validParams()
{
  InputParameters params = GenericKernel<is_ad>::validParams();
  params.addClassDescription("Demonstrates the multiple ways that scalar values can be introduced "
                             "into kernels, e.g. (controllable) constants, functions, and "
                             "postprocessors. Implements the weak form $(\\psi_i, -f)$.");
  params.addParam<Real>("value", 1.0, "Coefficient to multiply by the body force term");
  params.addParam<FunctionName>("function", "1", "A function that describes the body force");
  params.addParam<PostprocessorName>(
      "postprocessor", 1, "A postprocessor whose value is multiplied by the body force");
  params.declareControllable("value");
  return params;
}

template <bool is_ad>
BodyForceTempl<is_ad>::BodyForceTempl(const InputParameters & parameters)
  : GenericKernel<is_ad>(parameters),
    _scale(this->template getParam<Real>("value")),
    _function(getFunction("function")),
    _postprocessor(getPostprocessorValue("postprocessor")),
    _generic_q_point(this->_use_displaced_mesh ? &this->_assembly.template genericQPoints<is_ad>()
                                               : nullptr)
{
}

template <bool is_ad>
GenericReal<is_ad>
BodyForceTempl<is_ad>::computeQpResidual()
{
  if (_generic_q_point)
    return -_test[_i][_qp] * _scale * _postprocessor *
           _function.value(_t, (*_generic_q_point)[_qp]);
  else
    return -_test[_i][_qp] * _scale * _postprocessor * _function.value(_t, _q_point[_qp]);
}

template class BodyForceTempl<false>;
template class BodyForceTempl<true>;

One important observation in this case is that this class (BodyForceTempl) derives from a templated base class (GenericKernel<is_ad). This is slightly more complicated than the material example above. Members of the base class are not available in this derived class without including them with the using declaration in the header file (usingGenericKernelMembers), which is defined in the GenericKernel header

#pragma once

#include "Kernel.h"
#include "ADKernel.h"

template <bool is_ad>
class GenericKernel : public Kernel
{
public:
  static InputParameters validParams() { return Kernel::validParams(); };
  GenericKernel(const InputParameters & parameters) : Kernel(parameters) {}
};

template <>
class GenericKernel<true> : public ADKernel
{
public:
  static InputParameters validParams() { return ADKernel::validParams(); };
  GenericKernel(const InputParameters & parameters) : ADKernel(parameters) {}
};

#define usingGenericKernelMembers                                                                  \
  usingFunctionInterfaceMembers;                                                                   \
  usingPostprocessorInterfaceMembers;                                                              \
  usingMooseObjectMembers;                                                                         \
  usingTransientInterfaceMembers;                                                                  \
  usingTaggingInterfaceMembers;                                                                    \
  usingBlockRestrictableMembers;                                                                   \
  using GenericKernel<is_ad>::_qp;                                                                 \
  using GenericKernel<is_ad>::_i;                                                                  \
  using GenericKernel<is_ad>::_j;                                                                  \
  using GenericKernel<is_ad>::_u;                                                                  \
  using GenericKernel<is_ad>::_phi;                                                                \
  using GenericKernel<is_ad>::_test;                                                               \
  using GenericKernel<is_ad>::_q_point;                                                            \
  using GenericKernel<is_ad>::_var;                                                                \
  using GenericKernel<is_ad>::getVar;                                                              \
  using Coupleable::coupled;                                                                       \
  using Coupleable::coupledComponents

In addition, templated methods used in a class derived from a templated base class (like in example above) must be prefixed with this->template to avoid compiler ambiguity, for example, the use of getParam<T> in the BodyForce kernel above

_scale(this->template getParam<Real>("value")),