Increase val-2f/val-2c reproducibility

.gitignore

@@ -168,6 +168,7 @@ framework/contrib/asio/
 *.cpr
 *.cpa-*
 *.cpr-*
+*_cp/
 
 # Ignore petsc arch
 petsc/arch-*/*
@@ -245,3 +246,9 @@ doc/content/verification_and_validation/figures/val-2g_environment_history.png
 doc/content/examples/figures/fuel_cycle_*.png
 doc/content/source/interfacekernels/figures/YHx_PCT_*.png
 doc/content/source/interfacekernels/figures/ZrCoHx_PCT_*.png
+
+# temporary testharness_time files
+*.testharness_time
+
+# val-2c PSS results
+test/tests/val-2c/val-2c_pss_results_short

AGENTS.md

@@ -1,7 +1,7 @@
 # AGENTS.md
 
 ## Environment
-- Use the MOOSE conda environment (`conda activate moose`)
+- Ask the user if you need to load a MOOSE conda environment, and if so, what the name of the conda environment is
 - If `./tmap8-opt` or `run_tests` throws an error, stop and ask the user
 - Pre-commit hooks handle all formatting and input file validation on every `git commit` — do not invoke `clang-format`, `black`, or `hit format` manually
 - `moose/framework/contrib/hit/hit check <file>.i` is available to validate input files during development — use it before committing

doc/content/source/kernels/ScaledCoupledTimeDerivative.md

@@ -5,8 +5,9 @@
 ## Overview
 
 This class inherits from [CoupledTimeDerivative.md] and simply scales the
-residual and Jacobian by the user provided `factor`. The default value for
-`factor` is 1.
+residual and Jacobian by the user provided `factor`. It is also used for the
+dimensionless trapping formulation, where `factor` is set to a concentration
+reference ratio. The default value for `factor` is 1.
 
 !syntax parameters /Kernels/ScaledCoupledTimeDerivative
 

doc/content/source/nodal_kernels/ReleasingNodalKernelDimensionless.md

@@ -0,0 +1,30 @@
+# ReleasingNodalKernelDimensionless
+
+!syntax description /NodalKernels/ReleasingNodalKernelDimensionless
+
+## Overview
+
+`ReleasingNodalKernelDimensionless` is the dimensionless analogue of
+[ReleasingNodalKernel](ReleasingNodalKernel.md). It applies the release source term
+to a dimensionless trapped-species variable, $\hat{C}_t = C_t / C_{t,\mathrm{ref}}$
+with $C_t$ the trapped-species variable and $C_{t,\mathrm{ref}}$ the characteristic trapped-species scale in the same unit as $C_t$.
+
+This object implements the residual
+
+\begin{equation}
+\hat{k}_r \exp(-\epsilon / T) \hat{C}_t
+\end{equation}
+
+where $\hat{C}_t$ is the dimensionless trapped-species concentration, $\hat{k}_r$ is
+the dimensionless release rate, $\epsilon$ is the detrapping activation energy in
+Kelvin, and $T$ is the temperature in Kelvin.
+
+Compared with [ReleasingNodalKernel](ReleasingNodalKernel.md), the form is the same
+except that the trapped concentration variable is nondimensionalized and the rate
+coefficient is provided directly as a dimensionless quantity.
+
+!syntax parameters /NodalKernels/ReleasingNodalKernelDimensionless
+
+!syntax inputs /NodalKernels/ReleasingNodalKernelDimensionless
+
+!syntax children /NodalKernels/ReleasingNodalKernelDimensionless

doc/content/source/nodal_kernels/TrappingNodalKernelDimensionless.md

@@ -0,0 +1,49 @@
+# TrappingNodalKernelDimensionless
+
+!syntax description /NodalKernels/TrappingNodalKernelDimensionless
+
+## Overview
+
+`TrappingNodalKernelDimensionless` is the dimensionless analogue of
+[TrappingNodalKernel](TrappingNodalKernel.md). It applies the trapping source term
+to a dimensionless trapped-species variable, $\hat{C}_t = C_t / C_{t,\mathrm{ref}}$
+with $C_t$ the trapped-species variable and $C_{t,\mathrm{ref}}$ the characteristic trapped-species scale in the same unit as $C_t$.
+
+This object implements the residual
+
+\begin{equation}
+\hat{k}_t \exp(-\epsilon / T)
+\left(
+\frac{N C_{t0} - C_{t,\mathrm{ref}} \hat{C}_t - \sum_j C_{t,\mathrm{ref},j} \hat{C}_{t,j}}
+{C_{t,\mathrm{ref}}}
+\right)
+\hat{C}_m
+\end{equation}
+
+where $\hat{k}_t$ is the dimensionless trapping rate, $\epsilon$ is the trapping
+activation energy in Kelvin, $T$ is the temperature in Kelvin, $N$ is the atomic
+number density of the host material, $C_{t0}$ is the fraction of host sites that can
+trap the species, $\hat{C}_t$ is the dimensionless concentration of this trapped
+species, $\hat{C}_{t,j}$ are optional competing trapped-species variables, and
+$\hat{C}_m$ is the coupled mobile concentration variable.
+
+The implementation computes the available trapping sites in physical units,
+
+\begin{equation}
+N C_{t0} - C_{t,\mathrm{ref}} \hat{C}_t - \sum_j C_{t,\mathrm{ref},j} \hat{C}_{t,j},
+\end{equation}
+
+and then divides by $C_{t,\mathrm{ref}}$ so that the residual is written for the
+dimensionless unknown $\hat{C}_t$. Unlike [TrappingNodalKernel](TrappingNodalKernel.md),
+this object does not use the `trap_per_free` scaling factor that the dimensional
+kernel uses to convert between trapped-species and mobile-species concentration
+scales. Here, the conversion is handled explicitly through
+`trap_concentration_reference` ($C_{t,\mathrm{ref}}$) and `other_trap_concentration_references` ($C_{t,\mathrm{ref},j}$), so the
+residual is expected to remain $O(\hat{k}_t)$ when the variables are scaled with
+appropriate reference values.
+
+!syntax parameters /NodalKernels/TrappingNodalKernelDimensionless
+
+!syntax inputs /NodalKernels/TrappingNodalKernelDimensionless
+
+!syntax children /NodalKernels/TrappingNodalKernelDimensionless

doc/content/verification_and_validation/val-2c.md

@@ -159,8 +159,8 @@ The results presented here are updated results from those presented in [!cite](S
 
 [val-2c_comparison_T2] and [val-2c_comparison_HTO] show the comparison of the TMAP8 calculations (both with immediately injected and delayed injected T$_2$) against the experimental data for T$_2$ and HTO concentration in the enclosure over time.
 There is reasonable agreement between the TMAP8 predictions and the experimental data.
-In the case of immediate T$_2$ injection, the root mean square percentage errors (RMSPE) are equal to RMSPE = 58.68% for T$_2$ and RMSPE = 146.23% for HTO, respectively.
-When accounting for a delay in T$_2$ injection, the TMAP8 predictions best match the experimental data, in particular the position of the peak HTO concentration. The RMSPE values decrease to RMSPE = 89.50% for T$_2$ and RMSPE = 75.66% for HTO, respectively.
+In the case of immediate T$_2$ injection, the root mean square percentage errors (RMSPE) are equal to RMSPE = 58.68% for T$_2$ and RMSPE = 146.25% for HTO, respectively.
+When accounting for a delay in T$_2$ injection, the TMAP8 predictions best match the experimental data, in particular the position of the peak HTO concentration. The RMSPE values decrease to RMSPE = 89.50% for T$_2$ and RMSPE = 75.71% for HTO, respectively.
 Note that the model parameters listed in [val-2c_parameters] are somewhat different from [!cite](Holland1986,longhurst1992verification,ambrosek2008verification) to better match the experimental data.
 In particular, [!cite](longhurst1992verification,ambrosek2008verification) did not validate the TMAP predictions against T$_2$ concentration, which we do here in [val-2c_comparison_T2] and in [!cite](Simon2025).
 This affects some of the model parameters.
@@ -178,7 +178,7 @@ This affects some of the model parameters.
        caption=Comparison of TMAP8 calculations against the experimental data for HTO concentration in the enclosure over time. Calibration of the delayed injection model delivers further improvements, with more accurate simulation results when T$_2$ is injected over a given period rather than immediately.
 
 As shown in the red curve in [val-2c_comparison_T2] and [val-2c_comparison_HTO], using [MOOSE's stochastic tools module](modules/stochastic_tools/index.md) notably increased the agreement between the modeling predictions and experimental data for both the T$_2$ and HTO concentrations.
-The RMSPE for T$_2$ decreases from 89.50% to 30.18% and the RMSPE for HTO decreases from 75.66% to 67.07%.
+The RMSPE for T$_2$ decreases from 89.50% to 30.18% and the RMSPE for HTO decreases from 75.71% to 67.18%.
 Note that although the calibration approach is similar to the one presented in [!cite](Simon2025), the results presented here include more simulations and the quality of the calibration is increased here (RMSPE values are further decreased here).
 
 [val-2c_calibration_input] and [val-2c_calibration_output] show the evolution of the model parameter values and of the optimization metric (time integral of $g$ defined in [eq:optimization_metric]) as a function of the number of simulation. The calibrated model corresponds to the highest value.
@@ -237,4 +237,6 @@ Modifying the PSS parameters can reduce the computational cost.
 
 Although a very short PSS study is simulated as a test in [/val-2c/tests] to ensure these files run properly, the full calibration study is not performed regularly in tests to limit computational costs within the TMAP8 testing suite. The gold files [/gold/val-2c_pss_results/val-2c_pss_main_out.json] and [/gold/calibrated_parameter_values.txt] are, therefore, not continuously tested, and the calibrated model parameters used in [/val-2c/tests] are not continuously updated.
 
+In addition, this validation case includes tests using two time integration schemes. The `implicit-euler` scheme (first-order, single-step) provides higher stability and consistency of the solution across different operating systems and processor core counts, making it better suited for use in PSS optimization workflows where numerical stability is important. The `bdf2` scheme (second-order, multi-step) provides higher accuracy and is used for final results calculation. Both schemes are tested for this validation case in [/val-2f/tests].
+
 !bibtex bibliography

doc/content/verification_and_validation/val-2f.md

@@ -264,4 +264,4 @@ To learn more about the `!include` feature, refer to the [application_usage/inpu
 
 Note that both surface conditions can be modeled using [/val-2f.i]. By running it as is, the recombination condition with the updated recombination rate is utilized. [/val-2f/tests] used `cli_args` to modify [/val-2f.i] into using the effectively infinite recombination rate at the surface and impose a null concentration at the surfaces.
 
-To limit the computational costs of the test case, the test runs a version of the file with a smaller and coarser mesh, and fewer time steps. More information about the changes can be found in the test specification file for this case, namely [/val-2f/tests].
+To limit the computational costs of the test case, the test runs a version of the file with a smaller and coarser mesh, and fewer time steps. In addition, this validation case includes tests using two time integration schemes. The `implicit-euler` scheme (first-order, single-step) provides higher stability and consistency of the solution across different operating systems and processor core counts, making it better suited for use in Bayesian optimization workflows where numerical stability is important. The `bdf2` scheme (second-order, multi-step) provides higher accuracy and is used for final results calculation. Both schemes are tested for this validation case in [/val-2f/tests]. More information about the changes can be found in the test specification file for this case, namely [/val-2f/tests].

include/nodal_kernels/ReleasingNodalKernel.h

@@ -8,20 +8,15 @@
 
 #pragma once
 
-#include "NodalKernel.h"
+#include "ReleasingNodalKernelBase.h"
 
-class ReleasingNodalKernel : public NodalKernel
+/**
+ * Releasing NodalKernel for trapped-species concentrations in physical units.
+ */
+class ReleasingNodalKernel : public ReleasingNodalKernelBase
 {
 public:
   ReleasingNodalKernel(const InputParameters & parameters);
 
   static InputParameters validParams();
-
-protected:
-  Real computeQpResidual() override;
-  Real computeQpJacobian() override;
-
-  const Real _alpha_r;
-  const Real _detrapping_energy;
-  const VariableValue & _temperature;
 };

include/nodal_kernels/ReleasingNodalKernelBase.h

@@ -0,0 +1,36 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2025 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#pragma once
+
+#include "NodalKernel.h"
+
+/**
+ * Shared implementation for trapped-species release nodal kernels.
+ *
+ * This base class factors the common Arrhenius release residual and Jacobian
+ * used by both the dimensional and dimensionless releasing kernels.
+ */
+class ReleasingNodalKernelBase : public NodalKernel
+{
+public:
+  ReleasingNodalKernelBase(const InputParameters & parameters, Real release_rate_coefficient);
+
+  static InputParameters validParams();
+
+protected:
+  Real computeQpResidual() override;
+  Real computeQpJacobian() override;
+
+  /// Effective release-rate coefficient multiplying the Arrhenius factor and variable value.
+  const Real _release_rate_coefficient;
+  /// Detrapping activation energy, expressed in Kelvin.
+  const Real _detrapping_energy;
+  /// Coupled temperature field.
+  const VariableValue & _temperature;
+};

include/nodal_kernels/ReleasingNodalKernelDimensionless.h

@@ -0,0 +1,31 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2025 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#pragma once
+
+#include "ReleasingNodalKernelBase.h"
+
+/**
+ * Releasing NodalKernel for a dimensionless trapped-species variable
+ * Ct_hat = C_t / C_t_ref.
+ *
+ * The residual is:
+ *   R = +k_r_hat * exp(-E_r / T) * Ct_hat
+ *
+ * where k_r_hat = t_ref * alpha_r.
+ *
+ * The quadrature point residual should be dimensionless and O(k_r_hat) because Ct_hat should be
+ * dimensionless and O(1).
+ */
+class ReleasingNodalKernelDimensionless : public ReleasingNodalKernelBase
+{
+public:
+  ReleasingNodalKernelDimensionless(const InputParameters & parameters);
+
+  static InputParameters validParams();
+};

include/nodal_kernels/TrappingNodalKernel.h

@@ -8,43 +8,15 @@
 
 #pragma once
 
-#include "NodalKernel.h"
+#include "TrappingNodalKernelBase.h"
 
-#include "metaphysicl/dualdynamicsparsenumberarray.h"
-
-using MetaPhysicL::DualNumber;
-using MetaPhysicL::DynamicSparseNumberArray;
-
-class Function;
-
-// Forward Declarations
-typedef DualNumber<Real, DynamicSparseNumberArray<Real, unsigned int>> LocalDN;
-
-class TrappingNodalKernel : public NodalKernel
+/**
+ * Trapping NodalKernel for trapped-species concentrations in physical units.
+ */
+class TrappingNodalKernel : public TrappingNodalKernelBase
 {
 public:
   TrappingNodalKernel(const InputParameters & parameters);
 
   static InputParameters validParams();
-
-protected:
-  Real computeQpResidual() override;
-  Real computeQpJacobian() override;
-  Real computeQpOffDiagJacobian(unsigned int jvar) override;
-
-  const Real _alpha_t;
-  const Real _trapping_energy;
-  const Real _N;
-  const Function & _Ct0;
-  const VariableValue & _mobile_concentration;
-  unsigned int _n_other_concs;
-  std::vector<const VariableValue *> _trapped_concentrations;
-  std::vector<unsigned int> _var_numbers;
-  const Node * _last_node;
-  const Real _trap_per_free;
-  const VariableValue & _temperature;
-  LocalDN _jacobian;
-
-private:
-  void ADHelper();
 };

include/nodal_kernels/TrappingNodalKernelBase.h

@@ -0,0 +1,84 @@
+/************************************************************/
+/*                DO NOT MODIFY THIS HEADER                 */
+/*   TMAP8: Tritium Migration Analysis Program, Version 8   */
+/*                                                          */
+/*   Copyright 2021 - 2025 Battelle Energy Alliance, LLC    */
+/*                   ALL RIGHTS RESERVED                    */
+/************************************************************/
+
+#pragma once
+
+#include "NodalKernel.h"
+
+#include "metaphysicl/dualdynamicsparsenumberarray.h"
+
+using MetaPhysicL::DualNumber;
+using MetaPhysicL::DynamicSparseNumberArray;
+
+class Function;
+
+/// Local AD number type used to assemble manual Jacobians for trapping kernels.
+typedef DualNumber<Real, DynamicSparseNumberArray<Real, unsigned int>> LocalDN;
+
+/**
+ * Shared implementation for trapped-species trapping nodal kernels.
+ *
+ * This base class owns the common occupancy bookkeeping, empty-site assembly,
+ * and manual Jacobian logic used by both dimensional and dimensionless
+ * trapping kernels.
+ */
+class TrappingNodalKernelBase : public NodalKernel
+{
+public:
+  /**
+   * Build the shared trapping kernel state from the supplied rate and residual scaling.
+   */
+  TrappingNodalKernelBase(const InputParameters & parameters,
+                          Real trapping_rate_coefficient,
+                          Real residual_denominator);
+
+  static InputParameters validParams();
+
+protected:
+  /**
+   * Initialize the trapped-species variables that contribute to site occupancy.
+   * @param other_weights Physical occupancy weight for each coupled trapped species
+   * @param self_weight Physical occupancy weight for this kernel's primary variable
+   */
+  void initializeOccupancyTracking(const std::vector<Real> & other_weights, Real self_weight);
+
+  Real computeQpResidual() override;
+  Real computeQpJacobian() override;
+  Real computeQpOffDiagJacobian(unsigned int jvar) override;
+
+  /// Effective trapping-rate coefficient multiplying the Arrhenius factor.
+  const Real _trapping_rate_coefficient;
+  /// Trapping activation energy, expressed in Kelvin.
+  const Real _trapping_energy;
+  /// Atomic number density of the host material.
+  const Real _N;
+  /// Fraction of host sites available for trapping as a function of position.
+  const Function & _Ct0;
+  /// Coupled mobile-species concentration.
+  const VariableValue & _mobile_concentration;
+  /// Number of additional trapped concentrations coupled into the occupancy calculation.
+  unsigned int _n_other_concs;
+  /// Coupled concentration fields that consume trapping sites, including this kernel's variable.
+  std::vector<const VariableValue *> _occupancy_concentrations;
+  /// Physical occupancy weight associated with each trapped concentration.
+  std::vector<Real> _occupancy_weights;
+  /// Variable numbers corresponding to occupancy concentrations plus the mobile concentration.
+  std::vector<unsigned int> _var_numbers;
+  /// Most recent node for which the cached Jacobian state was assembled.
+  const Node * _last_node;
+  /// Coupled temperature field.
+  const VariableValue & _temperature;
+  /// Denominator used to scale the residual after empty-site assembly.
+  const Real _residual_denominator;
+  /// Cached AD residual used to extract diagonal and off-diagonal Jacobian entries.
+  LocalDN _jacobian;
+
+private:
+  /// Rebuild the cached AD residual when the active node changes.
+  void ADHelper();
+};