Add a Physics for radioactive decay

doc/content/source/kernels/ReleasingKernel.md

@@ -0,0 +1,12 @@
+# ReleasingKernel
+
+!syntax description /Kernels/ReleasingKernel
+
+This kernel is the volumetric equivalent of the [ReleasingNodalKernel.md]. See this nodal
+kernel documentation for information on the releasing mechanism.
+
+!syntax parameters /Kernels/ReleasingKernel
+
+!syntax inputs /Kernels/ReleasingKernel
+
+!syntax children /Kernels/ReleasingKernel

doc/content/source/kernels/TrappingKernel.md

@@ -0,0 +1,12 @@
+# TrappingKernel
+
+!syntax description /Kernels/TrappingKernel
+
+This kernel is the volumetric equivalent of the [TrappingNodalKernel.md]. See this nodal
+kernel documentation for information on the trapping mechanism.
+
+!syntax parameters /Kernels/TrappingKernel
+
+!syntax inputs /Kernels/TrappingKernel
+
+!syntax children /Kernels/TrappingKernel

include/kernels/ReleasingKernel.h

@@ -0,0 +1,41 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2024 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#pragma once
+
+#include "Kernel.h"
+
+/**
+ * Implements the contribution to the residual and Jacobian of a species release using a volumetric
+ * integration
+ */
+class ReleasingKernel : public Kernel
+{
+public:
+  ReleasingKernel(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  Real computeQpResidual() override;
+  Real computeQpJacobian() override;
+  Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+  /// Release coefficient
+  const Real _alpha_r;
+  /// Energy from detrapping
+  const Real _detrapping_energy;
+  /// Local temperature
+  const VariableValue & _temperature;
+  /// Species concentration
+  const VariableValue & _v;
+  /// Index of the species variable
+  const unsigned int _v_index;
+  /// Whether the v variable is the kernel's variable parameter variable
+  const bool _v_is_u;
+};

include/kernels/TrappingKernel.h

@@ -0,0 +1,44 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2024 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#pragma once
+
+#include "ADKernel.h"
+
+class Function;
+
+/**
+ * Implements the contribution to the residual and Jacobian using a volumetric integration and
+ * automatic differentiation
+ */
+class TrappingKernel : public ADKernel
+{
+public:
+  TrappingKernel(const InputParameters & parameters);
+
+  static InputParameters validParams();
+
+protected:
+  ADReal computeQpResidual() override;
+
+  /// Trapping coefficient
+  const Real _alpha_t;
+  /// Trapping energy
+  const Real _trapping_energy;
+  const Real _N;
+  const Function & _Ct0;
+  /// Concentration of the mobile species
+  const ADVariableValue & _mobile_concentration;
+  /// Number of other species competing for the same traps
+  unsigned int _n_other_concs;
+  /// Concentration of the other species
+  std::vector<const ADVariableValue *> _trapped_concentrations;
+  const Real _trap_per_free;
+  /// Local temperature
+  const ADVariableValue & _temperature;
+};

include/nodal_kernels/ReleasingNodalKernel.h

@@ -20,8 +20,24 @@ class ReleasingNodalKernel : public NodalKernel
 protected:
   Real computeQpResidual() override;
   Real computeQpJacobian() override;
+  Real computeQpOffDiagJacobian(unsigned int jvar) override;
 
+  /// Release coefficient
   const Real _alpha_r;
+  /// Energy from detrapping
   const Real _detrapping_energy;
+  /// Local temperature
   const VariableValue & _temperature;
+  /// Species concentration
+  const VariableValue & _v;
+  /// Index of the species variable
+  const unsigned int _v_index;
+  /// Whether the v variable is the kernel's variable parameter variable
+  const bool _v_is_u;
+  /// Whether the kernels are mass lumped to make it compatible
+  const bool _mass_lumped;
+  /// Local node mass
+  const VariableValue & _nodal_mass;
+  /// An array with ones for convenience
+  const VariableValue _one;
 };

include/nodal_kernels/TrappingNodalKernel.h

@@ -44,6 +44,12 @@ class TrappingNodalKernel : public NodalKernel
   const Real _trap_per_free;
   const VariableValue & _temperature;
   LocalDN _jacobian;
+  /// Whether the kernels are mass lumped to make it compatible
+  const bool _mass_lumped;
+  /// Local node mass
+  const VariableValue & _nodal_mass;
+  /// An array with ones for convenience
+  const VariableValue _one;
 
 private:
   void ADHelper();

include/physics/SpeciesRadioactiveDecay.h

@@ -0,0 +1,51 @@
+/********************************************************/
+/*             DO NOT MODIFY THIS HEADER                */
+/* TMAP8: Tritium Migration Analysis Program, Version 8 */
+/*                                                      */
+/*    Copyright 2021 Battelle Energy Alliance, LLC      */
+/*               ALL RIGHTS RESERVED                    */
+/********************************************************/
+
+#pragma once
+
+#include "SpeciesPhysicsBase.h"
+
+class ActionComponent;
+
+/**
+ * Creates all the objects needed to solve for the radioactive decay of the species
+ */
+class SpeciesRadioactiveDecay : public SpeciesPhysicsBase
+{
+public:
+  static InputParameters validParams();
+
+  SpeciesRadioactiveDecay(const InputParameters & parameters);
+
+  void addComponent(const ActionComponent & component) override;
+
+protected:
+  /// Return the name of the species variable
+  /// @param c_i index of the component
+  /// @param s_j index of  the species
+  VariableName getSpeciesVariableName(unsigned int c_i, unsigned int s_j) const;
+
+  /// Decay products on each component, for each species, for each decay reaction
+  std::vector<std::vector<std::vector<std::vector<VariableName>>>> _decay_products;
+  /// Decay constants on each component, for each species, for each decay reaction
+  std::vector<std::vector<std::vector<Real>>> _decay_constants;
+  /// Branching rations on each component, for each species, for each decay reaction, for each product species
+  std::vector<std::vector<std::vector<Real>>> _branching_ratios;
+
+  /// Whether to define a single variable for each species for all components, or a different one for each component
+  const bool _single_variable_set;
+  /// Whether to define initial conditions for the decaying species
+  const bool _add_initial_conditions;
+  /// Whether to use a nodal or a volumetric formulation for the definition of kernels, indexed per component
+  std::vector<bool> _use_nodal_formulations;
+
+private:
+  virtual void addSolverVariables() override;
+  virtual void addInitialConditions() override;
+  virtual void addFEKernels() override;
+};

include/physics/SpeciesTrappingPhysics.h

@@ -31,6 +31,8 @@ class SpeciesTrappingPhysics : public SpeciesPhysicsBase
   /// @param s_j index of  the species
   VariableName getSpeciesVariableName(unsigned int c_i, unsigned int s_j) const;
 
+  /// The discretization method for the trapping equations
+  const MooseEnum _discretization;
   /// The mobile species of interest
   std::vector<std::vector<VariableName>> _mobile_species_names;
 
@@ -54,6 +56,7 @@ class SpeciesTrappingPhysics : public SpeciesPhysicsBase
   const bool _single_variable_set;
 
 private:
+  bool isNodal() const { return _discretization == "nodal"; }
   virtual void addSolverVariables() override;
   virtual void addInitialConditions() override;
   virtual void addFEKernels() override;

src/base/TMAP8App.C

@@ -42,6 +42,8 @@ TMAP8App::registerAll(Factory & f, ActionFactory & af, Syntax & syntax)
   // TMAP8 specific Physics
   registerSyntax("SorptionExchangePhysics", "Physics/SorptionExchange/*");
   registerSyntax("SpeciesTrappingPhysics", "Physics/SpeciesTrapping/*");
+  registerSyntax("SpeciesDiffusionReactionCG", "Physics/SpeciesDiffusionReaction/*");
+  registerSyntax("SpeciesRadioactiveDecay", "Physics/SpeciesRadioactiveDecay/*");
 
   // Shorter syntax for MOOSE Physics used by TMAP8
   registerSyntax("DiffusionCG", "Physics/Diffusion/*");

src/kernels/ReleasingKernel.C

@@ -0,0 +1,69 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2024 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#include "ReleasingKernel.h"
+
+registerMooseObject("TMAP8App", ReleasingKernel);
+
+InputParameters
+ReleasingKernel::validParams()
+{
+  InputParameters params = Kernel::validParams();
+  params.addClassDescription(
+      "Implements a residual describing the release of trapped species in a material.");
+  params.addRequiredParam<Real>("alpha_r", "The release rate coefficient (1/s)");
+  params.addCoupledVar("trapped_concentration",
+                       "If specified, variable to compute the release rate with. "
+                       "Else, uses the 'variable' argument");
+  params.addParam<Real>("detrapping_energy", 0, "The detrapping energy (K)");
+  params.addRequiredCoupledVar("temperature", "The temperature (K)");
+  return params;
+}
+
+ReleasingKernel::ReleasingKernel(const InputParameters & parameters)
+  : Kernel(parameters),
+    _alpha_r(getParam<Real>("alpha_r")),
+    _detrapping_energy(getParam<Real>("detrapping_energy")),
+    _temperature(coupledValue("temperature")),
+    _v(isParamValid("trapped_concentration") ? coupledValue("trapped_concentration") : _u),
+    _v_index(coupled("trapped_concentration")),
+    _v_is_u(!isCoupled("trapped_concentration") ||
+            (coupled("trapped_concentration") == variable().number()))
+{
+}
+
+Real
+ReleasingKernel::computeQpResidual()
+{
+  // std::cout << "Volumetric "
+  //           << _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] <<
+  //           std::endl;
+  return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _v[_qp] * _test[_i][_qp];
+}
+
+Real
+ReleasingKernel::computeQpJacobian()
+{
+  if (_v_is_u)
+    return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _phi[_j][_qp] *
+           _test[_i][_qp];
+  else
+    return 0;
+}
+
+Real
+ReleasingKernel::computeQpOffDiagJacobian(unsigned int jvar)
+{
+  if (!_v_is_u && jvar == _v_index)
+    return _alpha_r * std::exp(-_detrapping_energy / _temperature[_qp]) * _phi[_j][_qp] *
+           _test[_i][_qp];
+  else
+    return 0;
+
+  // TODO: add temperature off-diagonal term
+}

src/kernels/TrappingKernel.C

@@ -0,0 +1,77 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2024 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#include "TrappingKernel.h"
+#include "Function.h"
+
+registerMooseObject("TMAP8App", TrappingKernel);
+
+InputParameters
+TrappingKernel::validParams()
+{
+  InputParameters params = ADKernel::validParams();
+  params.addClassDescription(
+      "Implements a residual describing the trapping of a species in a material.");
+  params.addRequiredParam<Real>("alpha_t",
+                                "The trapping rate coefficient. This has units of 1/s (e.g. no "
+                                "number densities are involved)");
+  params.addParam<Real>("trapping_energy", 0, "The trapping energy (K)");
+  params.addRequiredParam<Real>("N", "The atomic number density of the host material (1/m^3)");
+  params.addRequiredParam<FunctionName>(
+      "Ct0", "The fraction of host sites that can contribute to trapping as a function (-)");
+  params.addParam<Real>(
+      "trap_per_free",
+      1.,
+      "An estimate for the ratio of the concentration magnitude of trapped species to free "
+      "species. Setting a value for this can be helpful in producing a well-scaled matrix (-)");
+  params.addRequiredCoupledVar(
+      "mobile_concentration",
+      "The variable representing the mobile concentration of solute particles (1/m^3)");
+  params.addCoupledVar("trapped_concentration",
+                       "The variable representing the trapped concentration of solute particles "
+                       "(1/m^3). If unspecified, defaults to the 'variable' parameter");
+  params.addCoupledVar("other_trapped_concentration_variables",
+                       "Other variables representing trapped particle concentrations.");
+  params.addRequiredCoupledVar("temperature", "The temperature (K)");
+  return params;
+}
+
+TrappingKernel::TrappingKernel(const InputParameters & parameters)
+  : ADKernel(parameters),
+    _alpha_t(getParam<Real>("alpha_t")),
+    _trapping_energy(getParam<Real>("trapping_energy")),
+    _N(getParam<Real>("N")),
+    _Ct0(getFunction("Ct0")),
+    _mobile_concentration(adCoupledValue("mobile_concentration")),
+    _trap_per_free(getParam<Real>("trap_per_free")),
+    _temperature(adCoupledValue("temperature"))
+{
+  _n_other_concs = coupledComponents("other_trapped_concentration_variables");
+
+  // Resize to n_other_concs plus the concentration corresponding to this Kernel's variable
+  _trapped_concentrations.resize(1 + _n_other_concs);
+
+  // Keep pointers to references of the trapped species concentrations
+  for (MooseIndex(_n_other_concs) i = 0; i < _n_other_concs; ++i)
+    _trapped_concentrations[i] =
+        &adCoupledValue("other_trapped_concentration_variables", /*comp*/ 0);
+  _trapped_concentrations[_n_other_concs] =
+      isParamValid("trapped_concentration") ? &adCoupledValue("trapped_concentration") : &_u;
+}
+
+ADReal
+TrappingKernel::computeQpResidual()
+{
+  // Remove filled trapping sites from total number of trapping sites
+  ADReal empty_trapping_sites = _Ct0.value(_t, _q_point[_qp]) * _N;
+  for (const auto & trap_conc : _trapped_concentrations)
+    empty_trapping_sites -= (*trap_conc)[_qp] * _trap_per_free;
+
+  return -_alpha_t * std::exp(-_trapping_energy / _temperature[_qp]) * empty_trapping_sites *
+         _mobile_concentration[_qp] / (_N * _trap_per_free) * _test[_i][_qp];
+}