Fix unit scaling factors in InterfaceSorption and add transient tests

doc/content/source/interfacekernels/InterfaceSorption.md

@@ -19,7 +19,7 @@ The value of $n$ determines the sorption law. For example:
 - $n = 1$ corresponds to Henry's law
 - $n = 1/2$ corresponds to Sievert's law.
 
-A penalty is applied to enforce the sorption law at the interface. Note that solubilities in the literature may only be compatible with concentrations expressed in certain units, and those units may differ from those used elsewhere in the simulation, e.g., $mmol/kg$ versus $mol/m^3$. The model accepts concentration unit conversion factors to accommodate data from different sources. For example, unit conversion factors of $1000$ can be supplied to accommodate sets of solubilities for concentrations in $mmol/m^3$. Note that for multi-atomic molecules, it is important to specify if the moles correspond to the moles of gas molecules, or the moles of atoms. The model assumes that the units associated with the partial pressures are equal and that all constants are compatible with absolute temperature.
+A penalty is applied to enforce the sorption law at the interface. Note that solubilities in the literature may only be compatible with concentrations expressed in certain units, and those units may differ from those used elsewhere in the simulation, e.g., $mmol/kg$ versus $mol/m^3$. The model accepts concentration unit conversion factors (`unit_scale` and `unit_scale_neighbor`) to accommodate data from different sources. For example, unit conversion factors of $1000$ can be supplied to accommodate sets of solubilities for concentrations in $mmol/m^3$. Note that for multi-atomic molecules, it is important to specify if the moles correspond to the moles of gas molecules, or the moles of atoms. The model assumes that the units associated with the partial pressures are equal and that all constants are compatible with absolute temperature.
 
 To balance the mass flux across the gap, a second interfacial condition is given as:
 
@@ -34,6 +34,8 @@ When using the default, non-penalty method to enforce flux conditions at the int
 
 Optionally, a penalty can be applied within `InterfaceSorption` to directly enforce the full flux condition at the interface. The penalty method adds the correct residual to the neighbor side of the interface without requiring a second kernel, but it may over-constrain the problem under certain conditions.
 
+Note that the unit conversion factors (`unit_scale` and `unit_scale_neighbor`) do not apply to the flux equation to ensure mass conservation.
+
 ## Input File Usage Example
 
 !listing test/tests/interfacekernels/InterfaceSorption/interface_sorption.i block=InterfaceKernels

src/interfacekernels/InterfaceSorption.C

@@ -32,7 +32,7 @@ InterfaceSorptionTempl<is_ad>::validParams()
                                 "The pre-exponential factor for the Arrhenius law for the "
                                 "solubility $K$ for the relationship $C_i = KP_i^n$");
   params.addParam<Real>("Ea",
-                        0,
+                        0.,
                         "The activation energy for the Arrhenius law for the solubility $K$ for "
                         "the relationship $C_i = KP_i^n$");
   params.addRequiredParam<Real>("n_sorption",
@@ -96,7 +96,10 @@ InterfaceSorptionTempl<is_ad>::computeQpResidual(Moose::DGResidualType type)
   const GenericReal<is_ad> temperature_limited = std::max(small, _T[_qp]);
   const auto R = PhysicalConstants::ideal_gas_constant; // ideal gas constant (J/K/mol)
 
-  GenericReal<is_ad> r = 0; // residual
+  GenericReal<is_ad> r = 0.; // residual
+
+  // The unit scaling affects the residual of the sorption equation, but not the flux equation.
+  // Otherwise mass is not conserved.
 
   switch (type)
   {
@@ -132,7 +135,10 @@ InterfaceSorptionTempl<false>::computeQpJacobian(Moose::DGJacobianType type)
   const Real temperature_limited = std::max(small, _T[_qp]);
   const auto R = PhysicalConstants::ideal_gas_constant; // ideal gas constant (J/K/mol)
 
-  Real jac = 0; // jacobian
+  Real jac = 0.; // jacobian
+
+  // The unit scaling affects the jacobian of the sorption equation, but not the flux equation.
+  // Otherwise mass is not conserved.
 
   switch (type)
   {
@@ -141,10 +147,11 @@ InterfaceSorptionTempl<false>::computeQpJacobian(Moose::DGJacobianType type)
       break;
 
     case Moose::ElementNeighbor:
-      jac = -_test[_i][_qp] * _sorption_penalty * _phi_neighbor[_j][_qp] *
+      jac = -_test[_i][_qp] * _sorption_penalty * _unit_scale_neighbor / _unit_scale *
+            _phi_neighbor[_j][_qp] *
             (_K0 * std::exp(-_Ea / R / temperature_limited) *
              std::pow(R * temperature_limited, _n_sorption) * _n_sorption *
-             std::pow(u_neighbor, _n_sorption - 1));
+             std::pow(u_neighbor, _n_sorption - 1.));
       break;
 
     case Moose::NeighborElement:

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Henry_ad_out.e

@@ -0,0 +1 @@
+interface_sorption_Henry_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Henry_ad_penalty_out.e

@@ -0,0 +1 @@
+interface_sorption_Henry_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Henry_non_ad_penalty_out.e

@@ -0,0 +1 @@
+interface_sorption_Henry_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Sievert_ad_out.e

@@ -0,0 +1 @@
+interface_sorption_Sievert_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Sievert_ad_penalty_out.e

@@ -0,0 +1 @@
+interface_sorption_Sievert_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_Sievert_non_ad_penalty_out.e

@@ -0,0 +1 @@
+interface_sorption_Sievert_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Henry_ad_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Henry_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Henry_ad_penalty_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Henry_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Henry_ad_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Henry_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Henry_non_ad_penalty_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Henry_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Sievert_ad_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Sievert_non_ad_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Sievert_ad_penalty_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Sievert_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Sievert_ad_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Sievert_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/gold/interface_sorption_transient_Sievert_non_ad_penalty_scaling_out.e

@@ -0,0 +1 @@
+interface_sorption_transient_Sievert_non_ad_scaling_out.e

test/tests/interfacekernels/InterfaceSorption/interface_sorption.i

@@ -8,7 +8,13 @@
 # In this input file, we apply the Sievert law with n_sorption=1/2.
 
 # Physical Constants
-R = 8.31446261815324 # Based on PhysicalConstants
+R = 8.31446261815324 # J/mol/K Based on PhysicalConstants
+
+solubility = 1.e-2
+n_sorption = 0.5
+
+unit_scale = 1
+unit_scale_neighbor = 1
 
 
 [GlobalParams]
@@ -20,8 +26,7 @@ R = 8.31446261815324 # Based on PhysicalConstants
   [gen]
     type = GeneratedMeshGenerator
     nx = 10
-    ny = 10
-    dim = 2
+    dim = 1
   []
   [block1]
     type = SubdomainBoundingBoxGenerator
@@ -118,12 +123,12 @@ R = 8.31446261815324 # Based on PhysicalConstants
 [InterfaceKernels]
   [interface]
     type = InterfaceSorption
-    K0 = 1.e-2
+    K0 = ${solubility}
     Ea = 0
-    n_sorption = 0.5
+    n_sorption = ${n_sorption}
     diffusivity = diffusivity
-    unit_scale = 1
-    unit_scale_neighbor = 1
+    unit_scale = ${unit_scale}
+    unit_scale_neighbor = ${unit_scale_neighbor}
     temperature = temperature
     variable = u2
     neighbor_var = u1
@@ -142,7 +147,7 @@ R = 8.31446261815324 # Based on PhysicalConstants
   [properties_2]
     type = GenericConstantMaterial
     prop_names = 'thermal_conductivity diffusivity solubility'
-    prop_values = '2 2 1e-2'
+    prop_values = '2 2 ${solubility}'
     block = 2
   []
 []
@@ -156,9 +161,9 @@ R = 8.31446261815324 # Based on PhysicalConstants
   []
   [residual_concentration]
     type = ParsedFunction
-    symbol_names = 'u_mid_inner u_mid_outer T R solubility'
-    symbol_values = 'u_mid_inner u_mid_outer temperature_mid_inner ${R} 1e-2'
-    expression = 'u_mid_outer - solubility*sqrt(u_mid_inner*R*T)'
+    symbol_names = 'u_mid_inner u_mid_outer T'
+    symbol_values = 'u_mid_inner u_mid_outer temperature_mid_inner'
+    expression = 'u_mid_outer*${unit_scale} - ${solubility}*(u_mid_inner*${unit_scale_neighbor}*${R}*T)^${n_sorption}'
   []
   [flux_error]
     type = ParsedFunction
@@ -178,8 +183,6 @@ R = 8.31446261815324 # Based on PhysicalConstants
 [Executioner]
   type = Steady
   solve_type = NEWTON
-  # automatic_scaling = true
-  # scaling_group_variables = 'u1; u2; temperature'
 
   petsc_options_iname = '-pc_type -pc_factor_mat_solver_package'
   petsc_options_value = 'lu superlu_dist'
@@ -188,21 +191,19 @@ R = 8.31446261815324 # Based on PhysicalConstants
   nl_rel_tol = 1e-15
   nl_abs_tol = 1e-9
   l_tol = 1e-3
-  l_max_its = 50
-  nl_max_its = 50
 []
 
 [Postprocessors]
   [u_mid_inner]
     type = PointValue
     variable = u1
-    point = '0.49999 0.5 0'
+    point = '0.49999 0 0'
     outputs = 'csv console'
   []
   [u_mid_outer]
     type = PointValue
     variable = u2
-    point = '0.50001 0.5 0'
+    point = '0.50001 0 0'
     outputs = 'csv console'
   []
   [u_mid_diff]
@@ -213,13 +214,13 @@ R = 8.31446261815324 # Based on PhysicalConstants
   [temperature_mid_inner]
     type = PointValue
     variable = temperature
-    point = '0.49999 0.5 0'
+    point = '0.49999 0 0'
     outputs = csv
   []
   [temperature_mid_outer]
     type = PointValue
     variable = temperature
-    point = '0.50001 0.5 0'
+    point = '0.50001 0 0'
     outputs = csv
   []
   [residual_concentration]