Merge pull request #183 from Lee01Atom/ver-1ka

Ver 1ka

doc/content/verification_and_validation/figures/comparison_ver-1ka.py

@@ -0,0 +1 @@
+../../../../test/tests/ver-1ka/comparison_ver-1ka.py

doc/content/verification_and_validation/index.md

@@ -25,7 +25,7 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
 | ver-1hb | [Convective Gas Outflow Problem in Equilibrating Enclosures](ver-1hb.md)        |
 | ver-1ja | [Radioactive Decay of Mobile Tritium in a Slab](ver-1ja.md)                                       |
 | ver-1jb | [Radioactive Decay of Mobile Tritium in a Slab with a Distributed Trap Concentration](ver-1jb.md) |
-
+| ver-1ka | [Simple Volumetric Source](ver-1ka.md)                                                            |
 
 # List of benchmarking cases
 

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

@@ -0,0 +1,45 @@
+# ver-1ka
+
+# Simple Volumetric Source
+
+## General Case Description
+
+This problem involves two enclosures connected by a diffusive membrane that follows Sieverts law for diffusion. Both enclosures contain hydrogen (H$_2$), tritium (T$_2$), and tritium gas (HT). In the first enclosure, there is an initial inventory of only hydrogen (H$_2$) along with a constant volumetric source rate of tritium (T$_2$). The second enclosure starts out empty.
+
+## Case Set Up
+
+This verification problem is taken from [!cite](ambrosek2008verification). 
+The rise in pressure of T$_2$ molecules in the first enclosure can be monitored by not enabling T$_2$ molecules to traverse the membrane between the two enclosures (no tritium flux). Consequently, the rate of pressure increase in the first enclosure can be expressed as:
+
+\begin{equation} \label{eq:source_term}
+\frac{dP_{T_2}}{dt} = \frac{S}{V} kT,
+\end{equation}
+
+where $S$ represents the volumetric T$_2$ source rate, $V$ is the volume of the enclosure, $k$ is the Boltzmann constant, and $T$ is the temperature of the enclosure.
+In this case, $S$ is set to 10$^{20}$ molecules/m$^{-3}$/s, $V = 1$ m$^3$, and the temperature of the enclosure is constant at $T = 500$ K. 
+
+## Analytical Solution
+
+The analytical solution for [eq:source_term] is simply
+
+\begin{equation}
+P_{T_2}(t) = \frac{S}{V} kT t.
+\end{equation}
+
+## Results
+
+Comparison of the TMAP8 results and the analytical solution is shown in
+[ver-1ka_comparison_time] as a function of time. The TMAP8 code predictions match very well with the analytical solution with a root mean squared percentage error of RMSPE $= 0.00$ %.
+
+!media comparison_ver-1ka.py 
+       image_name=ver-1ka_comparison_time.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1ka_comparison_time
+       caption=Comparison of T$_2$ partial pressure in an enclosure with no loss pathways as function of time calculated through TMAP8 and analytically. TMAP8 matches the analytical solution.
+
+## Input files
+
+!style halign=left
+The input file for this case can be found at [/ver-1ka.i], which is also used as tests in TMAP8 at [/ver-1ka/tests].
+
+!bibtex bibliography

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

@@ -0,0 +1,64 @@
+import numpy as np
+import pandas as pd
+import os
+from matplotlib import gridspec
+import matplotlib.pyplot as plt
+
+# Changes working directory to script directory (for consistent MooseDocs usage)
+script_folder = os.path.dirname(__file__)
+os.chdir(script_folder)
+
+# Constants
+S = 1e20 # Source term (m^-3 * s^-1)
+V = 1 # Volume (m^3)
+kb = 1.380649e-23 # Boltzmann constant (J/K)
+T = 500 # Temperature (K)
+
+# Extract time and pressure data
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder = "../../../../test/tests/ver-1ka/gold/ver-1ka_out.csv"
+else:                                  # if in test folder
+    csv_folder = "./gold/ver-1ka_out.csv"
+expt_data = pd.read_csv(csv_folder)
+TMAP8_time = expt_data['time']
+TMAP8_pressure = expt_data['pressure']
+
+# Calculate the theoretical expression for pressure
+analytical_pressure = (S / V) * kb * T * TMAP8_time
+
+fig = plt.figure(figsize=[6.5, 5.5])
+gs = gridspec.GridSpec(1, 1)
+ax = fig.add_subplot(gs[0])
+
+# Plot the experimental data
+ax.plot(TMAP8_time/3600, TMAP8_pressure, label=r"TMAP8", c='tab:gray',linestyle='-')
+
+# Plot the selected theoretical data
+ax.plot(TMAP8_time/3600, analytical_pressure, label=r"Analytical",c='k', linestyle='--')
+
+RMSE_pressure = np.linalg.norm(TMAP8_pressure-analytical_pressure)
+err_percent_pressure = RMSE_pressure*100/np.mean(analytical_pressure)
+
+# Add text annotation for RMSPE on the plot
+ax.text(0.05, 0.95, '$P_{T_2}$ RMSPE = %.2f ' % err_percent_pressure + '%', fontweight='bold')
+
+# Format the y-axis to use scientific notation
+plt.gca().yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+
+# Label the axes
+ax.set_xlabel('Time (hr)')
+ax.set_ylabel('Pressure (Pa)')
+
+# define axis range
+ax.set_xlim(left=0)
+ax.set_xlim(right=3)
+ax.set_ylim(bottom=0)
+
+# Add a legend
+ax.legend(loc="best")
+
+# Add a grid
+plt.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+
+# Save the plot as a PNG file
+plt.savefig('ver-1ka_comparison_time.png', bbox_inches='tight')

test/tests/ver-1ka/gold/ver-1ka_out.csv

@@ -0,0 +1,22 @@
+time,pressure
+0,0
+500,345.16225000001
+1000,690.32450270146
+1500,1035.4867536019
+2000,1380.6490039021
+2500,1725.8112540022
+3000,2070.9735040355
+3500,2416.1357540466
+4000,2761.2980040503
+4500,3106.4602540515
+5000,3451.6225040519
+5500,3796.784754052
+6000,4141.9470040521
+6500,4487.1092540521
+7000,4832.2715040521
+7500,5177.4337540521
+8000,5522.5960040519
+8500,5867.7582540518
+9000,6212.9205040517
+9500,6558.0827540518
+10000,6903.2450040518

test/tests/ver-1ka/tests

@@ -0,0 +1,17 @@
+[Tests]
+  design = 'ODETimeDerivative.md ParsedODEKernel.md'
+  issues = '#12'
+  verification = 'ver-1ka.md'
+  [ver-1ka_csv]
+    type = CSVDiff
+    input = ver-1ka.i
+    csvdiff = ver-1ka_out.csv
+    requirement = 'The system shall be able to model a tritium volumetric source in one enclosure'
+  []
+  [ver-1ka_comparison]
+    type = RunCommand
+    command = 'python3 comparison_ver-1ka.py'
+    requirement = 'The system shall be able to generate comparison plots between the analytical solution and simulated solution of verification case 1ka, modeling a tritium volumetric source in one enclosure.'
+    required_python_packages = 'matplotlib numpy pandas os git'
+  []
+[]

test/tests/ver-1ka/ver-1ka.i

@@ -0,0 +1,48 @@
+initial_pressure = ${units 0 Pa} # initial internal pressure
+kb = ${units 1.380649e-23 J/K} # Boltzmann constant J/K - from PhysicalConstants.h
+T = ${units 500 K} # Temperature
+S = ${units 1e20 1/m^3/s} # Source term
+V = ${units 1 m^3} # Volume
+end_time = ${units 10000 s}
+time_step = ${units 500 s}
+
+[Mesh]
+  type = GeneratedMesh
+  dim = 2
+  nx = 10
+  ny = 10
+  xmax = 1
+  ymax = 1
+[]
+
+[Variables]
+  [pressure]
+    family = SCALAR
+    order = FIRST
+    initial_condition = '${fparse initial_pressure}'
+  []
+[]
+
+[ScalarKernels]
+  [time]
+    type = ODETimeDerivative
+    variable = pressure
+  []
+  [source]
+    type = ParsedODEKernel
+    variable = pressure
+    expression = '${fparse - S/V * kb * T}'
+  []
+[]
+
+[Executioner]
+  type = Transient
+  dt = ${time_step}
+  end_time = ${end_time}
+  scheme = 'bdf2'
+[]
+
+[Outputs]
+  file_base = 'ver-1ka_out'
+  csv = true
+[]