Ver 1kb

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

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

doc/content/verification_and_validation/index.md

@@ -32,6 +32,7 @@ TMAP8 also contains [example cases](examples/tmap_index.md), which showcase how
 | 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)                                                            |
+| ver-1kb | [Henry’s Law Boundaries with No Volumetric Source](ver-1kb.md)                                    |
 
 # List of benchmarking cases
 

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

@@ -0,0 +1,88 @@
+# ver-1kb
+
+# Henry’s Law Boundaries with No Volumetric Source
+
+## General Case Description
+
+Two enclosures are separated by a membrane that allows diffusion according to Henry’s law, with no volumetric source present. Enclosure 2 has twice the volume of Enclosure 1.
+
+## Case Set Up
+
+This verification problem is taken from [!cite](ambrosek2008verification).
+
+This setup describes a diffusion system in which tritium T$_2$ is modeled across a one-dimensional domain split into two enclosures. The total system length is $2.5 \times 10^{-4}$ m, divided into 100 segments. The system operates at a constant temperature of 500 Kelvin. Initial tritium pressures are specified as $10^{5}$ Pa for Enclosure 1 and $10^{-10}$ Pa for Enclosure 2.
+
+Over time, the pressures of T$_2$, which diffuses across the membrane in accordance with Henry’s law, will gradually equilibrate between the two enclosures.
+
+The diffusion process in each of the two enclosures can be described by
+
+\begin{equation}
+\frac{\partial C_1}{\partial t} = \nabla D \nabla C_1,
+\end{equation}
+and
+\begin{equation}
+\frac{\partial C_2}{\partial t} = \nabla D \nabla C_2,
+\end{equation}
+
+where $C_1$ and $C_2$ represent the concentration fields in enclosures 1 and 2 respectively, $t$ is the time, and $D$ denotes the diffusivity.
+
+The concentration in Enclosure 1 is related to the partial pressure and concentration in Enclosure 2 via the interface sorption law:
+
+\begin{equation}
+C_1 = K P_2^n = K \left( C_2 RT \right)^n
+\end{equation}
+
+where $R$ is the ideal gas constant in J/mol/K, $T$ is the temperature in K, $K$ is the solubility, and $n$ is the exponent of the sorption law. For Henry’s law, $n=1$.
+
+## Results
+
+Two subcases are considered. In the first subcase, we assume that $K=1/RT$ as is done in [!cite](ambrosek2008verification), which is expected to lead to $C_1 = C_2$ at equilibrium. In the second, $K=10/RT$, which is expected to lead to $C_1 = 10 C_2$. This second case is added to exercise TMAP8 in a case with a concentration jump.
+In the first subcase, consistent with the results from TMAP7, the pressure evolution in both enclosures is shown in [ver-1kb_comparison_time] as a function of time. Both pressures find equilibrium and become equal, which is consistent with $C_1 = K (RT C_2)^n$ for $K=1/RT$ and $n=1$. The concentration ratio between enclosures 1 and 2 in [ver-1kb_concentration_ratio] shows that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $K R T=1$. As shown in [ver-1kb_mass_conservation], mass is conserved between the two enclosures over time, with a variation in mass of only $1.0$ %. This variation in mass can be further minimized by refining the mesh, i.e., increasing the number of segments in the domain.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_comparison_time.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_comparison_time
+       caption=Equilibration of species pressures under Henry’s law for $K = 1/RT$.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_concentration_ratio.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_concentration_ratio
+       caption=Concentrations ratio between enclosures 1 and 2 at the interface for $K = 1/RT$. This verifies TMAP8's ability to apply the sorption law across the interface without a concentration jump.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_mass_conservation.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_mass_conservation
+       caption=Total mass conservation across both enclosures over time for $K = 1/RT$.
+
+In the second subcase, the sorption law with $K=10/RT$ does not lead to equal pressure in both enclosure. As illustrated in [ver-1kb_comparison_time_k10], the pressure jump maintains a ratio of $C_1/C_2 \approx 10$, which is consistent with the relationship $C_1 = K (RT C_2)^n$ for $K=10/RT$ and $n=1$. The concentration ratio between enclosures 1 and 2 in [ver-1kb_concentration_ratio_k10] shows that the results obtained with TMAP8 are consistent with the analytical results derived from the sorption law for $K RT=10$. Additionally, [ver-1kb_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time, with a variation in mass of only $0.25$ %. As in the previous case, this variation in mass can be further reduced by refining the mesh.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_comparison_time_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_comparison_time_k10
+       caption=Pressures jump of species under Henry’s law for $K = 10/RT$.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_concentration_ratio_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_concentration_ratio_k10
+       caption=Concentrations ratio between enclosures 1 and 2 at the interface for $K = 10/RT$. This verifies TMAP8's ability to apply the sorption law across the interface with a concentration jump.
+
+!media comparison_ver-1kb.py
+       image_name=ver-1kb_mass_conservation_k10.png
+       style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
+       id=ver-1kb_mass_conservation_k10
+       caption=Total mass conservation across both enclosures over time for $K = 10/RT$.
+
+!alert note title=A Comparison with TMAP7 Results: Impact of Diffusivity Variations on Kinetics
+The kinetics observed in our results differ from those presented in TMAP7. We attribute this discrepancy to a variation in the diffusivity value used, which significantly affects the diffusion rate.
+
+## Input files
+
+!style halign=left
+The input file for this case can be found at [/ver-1kb.i], which is also used as tests in TMAP8 at [/ver-1kb/tests].
+
+!bibtex bibliography

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

@@ -0,0 +1,145 @@
+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)
+
+# Extract columns for time, pressures, concentration_enclosure_1_at_interface, and concentration_enclosure_2_at_interface
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder = "../../../../test/tests/ver-1kb/gold/ver-1kb_out_k1.csv"
+else:                                  # if in test folder
+    csv_folder = "./gold/ver-1kb_out_k1.csv"
+expt_data = pd.read_csv(csv_folder)
+TMAP8_time = expt_data['time']
+TMAP8_pressure_enclosure_1 = expt_data['pressure_enclosure_1']
+TMAP8_pressure_enclosure_2 = expt_data['pressure_enclosure_2']
+concentration_enclosure_1_at_interface = expt_data['concentration_enclosure_1_at_interface']
+pressure_enclosure_2_at_interface = expt_data['pressure_enclosure_2_at_interface']
+mass_conservation_sum_encl1_encl2 = expt_data['mass_conservation_sum_encl1_encl2'].values
+concentration_ratio = expt_data['concentration_ratio']
+
+# Subplot 1: Pressure vs time
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time, TMAP8_pressure_enclosure_1, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
+ax.plot(TMAP8_time, TMAP8_pressure_enclosure_2, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel('Pressure (Pa)')
+ax.set_xlim(0, TMAP8_time.max())
+ax.set_ylim(bottom=0)
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+fig.savefig('ver-1kb_comparison_time.png', bbox_inches='tight', dpi=300)
+
+# Subplot 2: Solubility and concentration ratios vs time
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+solubility_ratio = [1] * len(TMAP8_time[1:])
+ax.plot(TMAP8_time[1:], concentration_ratio[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
+ax.set_yticks(np.arange(0, 3, 1))
+ax.set_xlim(0,TMAP8_time.max())
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((concentration_ratio[1:]-solubility_ratio)**2) )
+RMSPE = RMSE*100/np.mean(solubility_ratio)
+x_pos = TMAP8_time.max() / 7200
+y_pos = 0.9 * ax.get_ylim()[1]
+ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
+fig.savefig('ver-1kb_concentration_ratio.png', bbox_inches='tight', dpi=300)
+
+# Subplot 3: Mass Conservation Sum Encl 1 and 2 vs Time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time, mass_conservation_sum_encl1_encl2, c='tab:blue')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+mass_variation_percentage = (np.max(mass_conservation_sum_encl1_encl2)-np.min(mass_conservation_sum_encl1_encl2))/np.min(mass_conservation_sum_encl1_encl2)*100
+print("Percentage of mass variation: ", mass_variation_percentage)
+fig.savefig('ver-1kb_mass_conservation.png', bbox_inches='tight', dpi=300)
+
+# Repeat the same for K=10/RT
+
+if "/TMAP8/doc/" in script_folder:     # if in documentation folder
+    csv_folder_k10 = "../../../../test/tests/ver-1kb/gold/ver-1kb_out_k10.csv"
+else:                                  # if in test folder
+    csv_folder_k10 = "./gold/ver-1kb_out_k10.csv"
+expt_data_k10 = pd.read_csv(csv_folder_k10)
+TMAP8_time_k10 = expt_data_k10['time']
+TMAP8_pressure_enclosure_1_k10 = expt_data_k10['pressure_enclosure_1']
+TMAP8_pressure_enclosure_2_k10 = expt_data_k10['pressure_enclosure_2']
+concentration_enclosure_1_at_interface_k10 = expt_data_k10['concentration_enclosure_1_at_interface']
+pressure_enclosure_2_at_interface_k10 = expt_data_k10['pressure_enclosure_2_at_interface']
+mass_conservation_sum_encl1_encl2_k10 = expt_data_k10['mass_conservation_sum_encl1_encl2'].values
+concentration_ratio_k10 = expt_data_k10['concentration_ratio']
+
+# Subplot 1 : Pressure vs time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time_k10, TMAP8_pressure_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
+ax.plot(TMAP8_time_k10, TMAP8_pressure_enclosure_2_k10, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel('Pressure (Pa)')
+ax.set_xlim(0, TMAP8_time.max())
+ax.set_ylim(bottom=0)
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+fig.savefig('ver-1kb_comparison_time_k10.png', bbox_inches='tight', dpi=300)
+
+# Subplot 2: Solubility and concentration ratios vs time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+solubility_ratio = [10] * len(TMAP8_time[1:])
+ax.plot(TMAP8_time[1:], concentration_ratio_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time[1:], solubility_ratio, label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
+ax.set_yticks(np.arange(0, 21, 10))
+ax.set_xlim(0,TMAP8_time.max())
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
+ax.legend(loc="best")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+RMSE = np.sqrt(np.mean((concentration_ratio_k10[1:]-solubility_ratio)**2))
+RMSPE = RMSE*100/np.mean(solubility_ratio)
+x_pos = TMAP8_time.max() / 7200
+y_pos = 0.9 * ax.get_ylim()[1]
+ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold')
+fig.savefig('ver-1kb_concentration_ratio_k10.png', bbox_inches='tight', dpi=300)
+
+# Subplot 3 : Mass Conservation Sum Encl 1 and 2 vs Time
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time_k10, mass_conservation_sum_encl1_encl2_k10, c='tab:blue')
+ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
+ax.set_xlabel('Time (s)')
+ax.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
+ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
+mass_variation_percentage = (np.max(mass_conservation_sum_encl1_encl2_k10)-np.min(mass_conservation_sum_encl1_encl2_k10))/np.min(mass_conservation_sum_encl1_encl2_k10)*100
+print("Percentage of mass variation: ", mass_variation_percentage)
+fig.savefig('ver-1kb_mass_conservation_k10.png', bbox_inches='tight', dpi=300)
+