Add ADReaction (idaholab#21009)

framework/doc/content/source/kernels/Reaction.md

@@ -1,10 +1,10 @@
-# Reaction
+# Reaction / ADReaction
 
 ## Description
 
-`Reaction` implements a simple first-order reaction term where the rate of
-reaction is directly proportional to the governing variable $u$. Its weak form
-is given by
+`Reaction` (and its automatic differentiation version, `ADReaction`) implements
+a simple first-order reaction term where the rate of reaction is directly proportional
+to the governing variable $u$. Its weak form is given by
 \begin{equation}
 (\psi_i, \lambda u_h)
 \end{equation}
@@ -21,9 +21,8 @@ parameters. An example block is shown below for a diffusion-reaction equation:
 !listing test/tests/fvkernels/fv_coupled_var/coupled.i block=Kernels
 
 !alert note
-There is no AD (automatic differentiation) or FV (finite volume) version of `Reaction`.
-If you wish to use AD / FV, use [ADCoupledForce](/ADCoupledForce.md) /
-[FVCoupledForce](/FVCoupledForce.md) respectively.
+There is no FV (finite volume) version of `Reaction`. If you wish to use FV, use
+[/FVCoupledForce.md].
 
 !syntax parameters /Kernels/Reaction
 

framework/include/kernels/Reaction.h

@@ -9,18 +9,29 @@
 
 #pragma once
 
-#include "Kernel.h"
+#include "GenericKernel.h"
 
-class Reaction : public Kernel
+/**
+ *  Implements a simple consuming reaction term with weak form $(\\psi_i, \\lambda u_h)$.
+ */
+
+template <bool is_ad>
+class ReactionTempl : public GenericKernel<is_ad>
 {
 public:
   static InputParameters validParams();
 
-  Reaction(const InputParameters & parameters);
+  ReactionTempl(const InputParameters & parameters);
 
 protected:
-  virtual Real computeQpResidual() override;
+  virtual GenericReal<is_ad> computeQpResidual() override;
   virtual Real computeQpJacobian() override;
 
+  /// Scalar coefficient representing the relative amount consumed per unit time
   const Real & _rate;
+
+  usingGenericKernelMembers;
 };
+
+typedef ReactionTempl<false> Reaction;
+typedef ReactionTempl<true> ADReaction;

framework/src/kernels/Reaction.C

@@ -10,31 +10,42 @@
 #include "Reaction.h"
 
 registerMooseObject("MooseApp", Reaction);
+registerMooseObject("MooseApp", ADReaction);
 
+template <bool is_ad>
 InputParameters
-Reaction::validParams()
+ReactionTempl<is_ad>::validParams()
 {
-  InputParameters params = Kernel::validParams();
+  InputParameters params = GenericKernel<is_ad>::validParams();
   params.addClassDescription(
       "Implements a simple consuming reaction term with weak form $(\\psi_i, \\lambda u_h)$.");
   params.addParam<Real>(
       "rate", 1.0, "The $(\\lambda)$ multiplier, the relative amount consumed per unit time.");
   return params;
 }
 
-Reaction::Reaction(const InputParameters & parameters)
-  : Kernel(parameters), _rate(getParam<Real>("rate"))
+template <bool is_ad>
+ReactionTempl<is_ad>::ReactionTempl(const InputParameters & parameters)
+  : GenericKernel<is_ad>(parameters), _rate(this->template getParam<Real>("rate"))
 {
 }
 
-Real
-Reaction::computeQpResidual()
+template <bool is_ad>
+GenericReal<is_ad>
+ReactionTempl<is_ad>::computeQpResidual()
 {
   return _test[_i][_qp] * _rate * _u[_qp];
 }
 
+template <bool is_ad>
 Real
-Reaction::computeQpJacobian()
+ReactionTempl<is_ad>::computeQpJacobian()
 {
+  // This function will never be called for the AD version. But because C++ does
+  // not support an optional function declaration based on a template parameter,
+  // we must keep this this template for all cases.
   return _test[_i][_qp] * _rate * _phi[_j][_qp];
 }
+
+template class ReactionTempl<false>;
+template class ReactionTempl<true>;

test/tests/kernels/ad_reaction/ad_reaction.i

@@ -0,0 +1,46 @@
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 5
+  ny = 5
+[]
+
+[Variables]
+  [u]
+  []
+[]
+
+[Kernels]
+  [diffusion]
+    type = ADDiffusion
+    variable = u
+  []
+
+  [reaction]
+    type = ADReaction
+    variable = u
+  []
+
+  [force]
+    type = ADBodyForce
+    variable = u
+  []
+[]
+
+[BCs]
+  [left]
+    type = ADDirichletBC
+    boundary = left
+    variable = u
+    value = 0
+  []
+[]
+
+[Executioner]
+  type = Steady
+  solve_type = NEWTON
+[]
+
+[Outputs]
+  exodus = true
+[]

test/tests/kernels/ad_reaction/tests

@@ -0,0 +1,17 @@
+[Tests]
+  [ad_reaction]
+    type = 'Exodiff'
+    input = 'ad_reaction.i'
+    exodiff = 'ad_reaction_out.e'
+    requirement = 'The system shall contain a consuming reaction term object whose Jacobian is calculated via forward automatic differentiation.'
+    design = 'Reaction.md'
+    issues = '#21009'
+  []
+  [ad_reaction_jac]
+    type = PetscJacobianTester
+    input = 'ad_reaction.i'
+    requirement = 'The system shall produce a perfect Jacobian for a consuming reaction term object using forward automatic differentiation.'
+    design = 'Reaction.md'
+    issues = '#21009'
+  []
+[]