Add MMS to ver-1dc

doc/content/verification_and_validation/figures/spatial_mms.py

@@ -0,0 +1 @@
+../../../../test/tests/ver-1dc/spatial_mms.py

doc/content/verification_and_validation/ver-1dc.md

@@ -90,6 +90,100 @@ The analytical solution for the permeation transient is compared with TMAP8 resu
        id=ver-1dc_comparison_diffusion
        caption=Comparison of TMAP8 calculation with the analytical solution for the permeation transient from [!cite](ambrosek2008verification).
 
+## Further verification using the method of manufactured solutions (MMS)
+
+Although the flux and breakthrough time can be verified using analytical solutions as presented above, the [MMS](mms.md) approach is a powerful method to verify complex system of PDEs such as the one studied in this case.
+Here, we apply the [MMS](mms.md) approach available as a Python-based utility in the [MOOSE framework](https://mooseframework.inl.gov) to verify TMAP8's predictions for ver-1dc.
+
+Below, we (1) describe how to derive the weak form of the equation, which is necessary to apply the MMS, (2) detail how we apply the MMS approach to this case, and (3) discuss the results.
+
+### Derivation of the weak forms of the equations
+
+#### Step 1: Define and rearrange the strong form of the equations
+
+The strong form of the governing equations is provided in [eqn:diffusion_mobile] and [eqn:trapped_rate], and can be slightly rearranged as:
+
+\begin{equation} \label{eqn:diffusion_mobile_step1}
+    \frac{\partial C_M}{\partial t} - \nabla \cdot \left( D \nabla C_M \right) + \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{\partial C_{T_i}}{\partial t} = 0,
+\end{equation}
+and, for $i=1$, $i=2$, and $i=3$:
+\begin{equation} \label{eqn:trapped_rate_step1}
+    \frac{\partial C_{T_i}}{\partial t} - \alpha_t^i \frac{C_{T_i}^{\text{empty}} C_M}{N \cdot \text{trap\_per\_free}} + \alpha_r^i C_{T_i} = 0,
+\end{equation}
+respectively.
+
+#### Step 2: Multiply every term by a test function
+
+Multiply each term of the equations by an appropriate test function $\psi$ (or $\psi_i$ for $C_{T_i}$):
+
+[eqn:diffusion_mobile_step1] becomes
+\begin{equation} \label{eqn:diffusion_mobile_step2}
+    \psi \frac{\partial C_M}{\partial t} - \psi \nabla \cdot \left( D \nabla C_M \right) + \psi \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{\partial C_{T_i}}{\partial t} = 0,
+\end{equation}
+and, for $i=1$, $i=2$, and $i=3$, [eqn:trapped_rate_step1] becomes
+\begin{equation} \label{eqn:trapped_rate_step2}
+    \psi_i \left( \frac{\partial C_{T_i}}{\partial t} - \alpha_t^i \frac{C_{T_i}^{\text{empty}} C_M}{N \cdot \text{trap\_per\_free}} + \alpha_r^i C_{T_i} \right) = 0.
+\end{equation}
+
+#### Step 3: Integrate over the domain $\Omega$
+
+After integration over the domain, we obtain:
+
+\begin{equation} \label{eqn:diffusion_mobile_step3}
+    \int_\Omega \psi \frac{\partial C_M}{\partial t} \, d\Omega - \int_\Omega \psi \nabla \cdot \left( D \nabla C_M \right) \, d\Omega + \int_\Omega \psi \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{\partial C_{T_i}}{\partial t} \, d\Omega = 0,
+\end{equation}
+
+and, for $i=1$, $i=2$, and $i=3$,
+\begin{equation} \label{eqn:trapped_rate_step3}
+    \int_\Omega \psi_i \frac{\partial C_{T_i}}{\partial t} \, d\Omega - \int_\Omega \psi_i \alpha_t^i \frac{C_{T_i}^{\text{empty}} C_M}{N \cdot \text{trap\_per\_free}} \, d\Omega + \int_\Omega \psi_i \alpha_r^i C_{T_i} \, d\Omega = 0.
+\end{equation}
+
+#### Step 4: Integrate by parts and apply the divergence theorem
+
+For the divergence term in [eqn:diffusion_mobile_step3], applying integration by parts and the divergence theorem provides
+\begin{equation}
+    \int_\Omega \psi \nabla \cdot \left( D \nabla C_M \right) \, d\Omega = - \int_\Omega \nabla \psi \cdot \left( D \nabla C_M \right) \, d\Omega + \oint\limits_{\partial \Omega} \psi \left( D \nabla C_M \right) \cdot \mathbf{n} \, d\partial \Omega,
+\end{equation}
+where $\partial \Omega$ is the boundary of the domain and $\mathbf{n}$ is the outward-facing normal vector.
+
+This update term is then substituted back into [eqn:diffusion_mobile_step3], which leads to:
+\begin{equation} \label{eqn:diffusion_mobile_step4}
+    \int_\Omega \psi \frac{\partial C_M}{\partial t} \, d\Omega + \int_\Omega \nabla \psi \cdot \left( D \nabla C_M \right) \, d\Omega - \oint\limits_{\partial \Omega} \psi \left( D \nabla C_M \right) \cdot \mathbf{n} \, d\partial \Omega + \int_\Omega \psi \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{\partial C_{T_i}}{\partial t} \, d\Omega = 0.
+\end{equation}
+
+Since no divergence terms exist in [eqn:trapped_rate_step3], no integration by parts is needed.
+
+#### Step 5: Derive the final weak form in inner product notation
+
+[eqn:diffusion_mobile_step4] and [eqn:trapped_rate_step3] can be expressed in inner product notation, where $\langle a, b \rangle = \int_\Omega a b \, d\Omega$.
+
+[eqn:diffusion_mobile_step4] becomes
+\begin{equation}
+    \langle \psi, \frac{\partial C_M}{\partial t} \rangle + \langle \nabla \psi, D \nabla C_M \rangle - \langle \psi, \text{trap\_per\_free} \cdot \sum_{i=1}^{3} \frac{\partial C_{T_i}}{\partial t} \rangle + \int_{\partial \Omega} \psi \left( D \nabla C_M \right) \cdot \mathbf{n} \, d\Gamma = 0,
+\end{equation}
+and, for $i=1$, $i=2$, and $i=3$, [eqn:trapped_rate_step3] becomes
+\begin{equation}
+    \langle \psi_i, \frac{\partial C_{T_i}}{\partial t} \rangle - \langle \psi_i, \alpha_t^i \frac{C_{T_i}^{\text{empty}} C_M}{N \cdot \text{trap\_per\_free}} \rangle + \langle \psi_i, \alpha_r^i C_{T_i} \rangle = 0.
+\end{equation}
+These weak forms of the equations can now be used for finite element analysis and MMS.
+
+### Application of the MMS
+
+To perform MMS we select a smoothly varying sinusoidal spatial solution for the mobile and trapped species. In order to prevent temporal error from polluting the spatial error we pick temporal dependence that can be exactly represented by the time integrator we choose. In this MMS case we use implicit or backwards Euler which is first order accurate; consequently we choose a linear dependence on time. For our 1D case this leads us to select the following MMS solution for the mobile species:
+
+\begin{equation}
+\cos(x) t
+\end{equation}
+
+The trapping species concentrations are represented similarly. With the manufactured solutions selected we generate forcing functions by substituting the exact solutions into the strong-form PDEs, [eqn:diffusion_mobile_step1] and [eqn:trapped_rate_step1]. These forcing functions are then imposed in the TMAP8 input file using [BodyForce.md] and [UserForcingFunctionNodalKernel.md] respectively. Dirichlet boundary conditions for the mobile species are imposed using the selected exact/MMS solution. For the spatial discretization, we select first order Lagrange basis functions for the mobile concentration and for projecting the trapped species concentrations from nodes into element interiors for coupling in the trapped specie time derivative term in the mobile specie governing equation.
+
+### Results
+
+!media spatial_mms.py
+       image_name=ver-1dc-mms-spatial.png
+       id=ver-1dc_mms
+       caption=Spatial convergence for a diffusion-trapping-release problem modeled after the physics of ver-1dc using first order Lagrange basis functions. The expected quadratic convergence rate of the $L^2$ error is observed.
+
 
 ## Input files
 

doc/moosedocs.py

@@ -12,10 +12,10 @@
 import os
 
 # Locate MOOSE directory
-MOOSE_DIR = os.getenv('MOOSE_DIR', os.path.abspath(os.path.join(os.path.dirname(__name__), '..', 'moose')))
-if not os.path.exists(MOOSE_DIR):
-    MOOSE_DIR = os.path.abspath(os.path.join(os.path.dirname(__name__), '..', '..', 'moose'))
-if not os.path.exists(MOOSE_DIR):
+MOOSE_DIR = os.getenv('MOOSE_DIR', os.path.abspath(os.path.join(os.path.dirname(__file__), '..', 'moose')))
+if not os.path.exists(os.path.join(MOOSE_DIR, 'python')):
+    MOOSE_DIR = os.path.abspath(os.path.join(os.path.dirname(__file__), '..', '..', 'moose'))
+if not os.path.exists(os.path.join(MOOSE_DIR, 'python')):
     raise Exception('Failed to locate MOOSE, specify the MOOSE_DIR environment variable.')
 os.environ['MOOSE_DIR'] = MOOSE_DIR
 

test/tests/val-2a/comparison_val-2a.py

@@ -29,7 +29,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
     return new_tmap_output
 
 # Read simulation data
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2a/gold/val-2a_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/val-2a_out.csv"
@@ -39,7 +39,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
 simulation_recom_flux_right_TMAP4 = simulation_TMAP4_data['scaled_recombination_flux_right']
 
 # Read experiment data
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2a/gold/experiment_data_paper.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/experiment_data_paper.csv"

test/tests/val-2b/comparison_val-2b.py

@@ -48,7 +48,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
 #===============================================================================
 # Extract TMAP8 predictions
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2b/gold/val-2b_heavy_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/val-2b_heavy_out.csv"
@@ -77,7 +77,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
 #===============================================================================
 # Extract experimental data
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2b/gold/experimental_data.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/experimental_data.csv"

test/tests/val-2c/comparison_val-2c.py

@@ -81,7 +81,7 @@ def sample_solution(time_stamps, time_solution, value_solution):
 
 #===============================================================================
 # Extract TMAP8 predictions
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2c/gold/val-2c_immediate_injection_csv.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/val-2c_immediate_injection_csv.csv"
@@ -91,7 +91,7 @@ def sample_solution(time_stamps, time_solution, value_solution):
 tmap8_hto_enclosure_concentration = tmap8_prediction['hto_enclosure_edge_concentration']
 
 # Extract TMAP8 predictions - delay
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2c/gold/val-2c_delay_csv.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/val-2c_delay_csv.csv"
@@ -102,15 +102,15 @@ def sample_solution(time_stamps, time_solution, value_solution):
 
 #===============================================================================
 # Extract Experimental data
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2c/gold/Experimental_data_HTO_concentration.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/Experimental_data_HTO_concentration.csv"
 experimental_data_hto = pd.read_csv(csv_folder)
 experimental_time_hto = experimental_data_hto['time (hr)']
 experimental_hto_enclosure_concentration = experimental_data_hto['Concentration (Ci/m3)']
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2c/gold/Experimental_data_T2_concentration.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/Experimental_data_T2_concentration.csv"

test/tests/val-2d/comparison_val-2d.py

@@ -37,7 +37,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
 
 
 # Read simulation data
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2d/gold/val-2d_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/val-2d_out.csv"
@@ -52,7 +52,7 @@ def numerical_solution_on_experiment_input(experiment_input, tmap_input, tmap_ou
 
 
 # Read experiment data
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/val-2d/gold/experiment_data_paper.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/experiment_data_paper.csv"

test/tests/ver-1a/comparison_ver-1a.py

@@ -152,7 +152,7 @@ def analytical_expression_flux(t, P_0, D, S, V, T, A, l):
     return flux
 
 # Extract data from 'gold' TMAP8 run
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1a/gold/ver-1a_csv.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1a_csv.csv"

test/tests/ver-1b/comparison_ver-1b.py

@@ -15,7 +15,7 @@
 gs = gridspec.GridSpec(1,1)
 ax = fig.add_subplot(gs[0])
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1b/gold/ver-1b_csv.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1b_csv.csv"
@@ -54,7 +54,7 @@
 gs = gridspec.GridSpec(1,1)
 ax = fig.add_subplot(gs[0])
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1b/gold/ver-1b_vector_postproc_line_0250.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1b_vector_postproc_line_0250.csv"
@@ -90,7 +90,7 @@
 gs = gridspec.GridSpec(1,1)
 ax = fig.add_subplot(gs[0])
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1b/gold/ver-1b_csv.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1b_csv.csv"

test/tests/ver-1c/comparison_ver-1c.py

@@ -10,7 +10,7 @@
 script_folder = os.path.dirname(__file__)
 os.chdir(script_folder)
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder4 = "../../../../test/tests/ver-1c/gold/ver-1c_tmap4.csv"
     csv_folder7 = "../../../../test/tests/ver-1c/gold/ver-1c_tmap7.csv"
 else:                                  # if in test folder

test/tests/ver-1d/comparison_ver-1d.py

@@ -55,7 +55,7 @@ def summation_term(num_terms, time):
         sum += (-1)**m * np.exp(-1 * m**2 * time / (2*tau_be))
     return sum
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1d/gold/ver-1d-diffusion_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1d-diffusion_out.csv"
@@ -101,7 +101,7 @@ def summation_term(num_terms, time):
 # breakthrough time.
 analytical_sol = [0, 3.2e18]
 
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1d/gold//ver-1d-trapping_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold//ver-1d-trapping_out.csv"

test/tests/ver-1dc/comparison_ver-1dc.py

@@ -62,7 +62,7 @@ def summation_term(num_terms, time):
     return sum
 
 # Extract data from 'gold' TMAP8 run
-if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+if "/tmap8/doc/" in script_folder.lower():     # if in documentation folder
     csv_folder = "../../../../test/tests/ver-1dc/gold/ver-1dc_out.csv"
 else:                                  # if in test folder
     csv_folder = "./gold/ver-1dc_out.csv"

test/tests/ver-1dc/functions.i

@@ -0,0 +1,46 @@
+[exact_u]
+  type = ParsedFunction
+  expression = 't*cos(x)'
+[]
+[forcing_u]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac1*cos(x) + (1/2)*N*frac2*cos(x) + (1/2)*N*frac3*cos(x) + t*cos(x) + cos(x)'
+  symbol_names = 'N frac1 frac2 frac3'
+  symbol_values = '${N} ${frac1} ${frac2} ${frac3}'
+[]
+[exact_t1]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac1*t*cos(x) + (1/2)*N*frac1'
+  symbol_names = 'frac1 N'
+  symbol_values = '${frac1} ${N}'
+[]
+[forcing_t1]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac1*cos(x) + alphar*((1/2)*N*frac1*t*cos(x) + (1/2)*N*frac1)*exp(-epsilon_1/temperature) - alphat*t*(-1/2*N*frac1*t*cos(x) + (1/2)*N*frac1)*cos(x)/N'
+  symbol_names = 'alphar alphat N frac1 temperature epsilon_1'
+  symbol_values = '${alphar} ${alphat} ${N} ${frac1} ${temperature} ${epsilon_1}'
+[]
+[exact_t2]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac2*t*cos(x) + (1/2)*N*frac2'
+  symbol_names = 'frac2 N'
+  symbol_values = '${frac2} ${N}'
+[]
+[forcing_t2]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac2*cos(x) + alphar*((1/2)*N*frac2*t*cos(x) + (1/2)*N*frac2)*exp(-epsilon_2/temperature) - alphat*t*(-1/2*N*frac2*t*cos(x) + (1/2)*N*frac2)*cos(x)/N'
+  symbol_names = 'alphar alphat N frac2 temperature epsilon_2'
+  symbol_values = '${alphar} ${alphat} ${N} ${frac2} ${temperature} ${epsilon_2}'
+[]
+[exact_t3]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac3*t*cos(x) + (1/2)*N*frac3'
+  symbol_names = 'frac3 N'
+  symbol_values = '${frac3} ${N}'
+[]
+[forcing_t3]
+  type = ParsedFunction
+  expression = '(1/2)*N*frac3*cos(x) + alphar*((1/2)*N*frac3*t*cos(x) + (1/2)*N*frac3)*exp(-epsilon_3/temperature) - alphat*t*(-1/2*N*frac3*t*cos(x) + (1/2)*N*frac3)*cos(x)/N'
+  symbol_names = 'alphar alphat N frac3 temperature epsilon_3'
+  symbol_values = '${alphar} ${alphat} ${N} ${frac3} ${temperature} ${epsilon_3}'
+[]