Ver 1kb

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

@@ -10,6 +10,8 @@ Two enclosures are separated by a membrane that allows diffusion according to He
 
 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 length of this domain is defined by the number of segments (20 segments) and the node size of $ 1.25\cdot 10^{-5}$ m, yielding a total system length of $ 2.5\cdot 10^{-4}$ m. The system operates at a constant temperature of $ 500$ Kelvin, with an ideal gas constant $ R=8.314$ J/(mol K). 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
@@ -35,7 +37,7 @@ where $R$ is the ideal gas constant in J/mol/K, $T$ is the temperature in K, $K$
 ## 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 xxx \%.
+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 $2 \cdot 10^{—6}$ \%.
 
 !media comparison_ver-1kb.py
        image_name=ver-1kb_comparison_time.png
@@ -55,7 +57,7 @@ In the first subcase, consistent with the results from TMAP7, the pressure evolu
        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 xxx \%..
+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 $8 \cdot 10^{—7}$ \%.
 
 !media comparison_ver-1kb.py
        image_name=ver-1kb_comparison_time_k10.png
@@ -75,6 +77,9 @@ In the second subcase, the sorption law with $K=10/RT$ does not lead to equal pr
        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

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

@@ -22,7 +22,60 @@
 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 / 3600, TMAP8_pressure_enclosure_1, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
+ax.plot(TMAP8_time / 3600, 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 (hr)')
+ax.set_ylabel('Pressure (Pa)')
+ax.set_xlim(-0.001, TMAP8_time.max() / 3600)
+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:] / 3600, concentration_ratio[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time[1:] / 3600, 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() / 3600)
+ax.set_xlabel('Time (hr)')
+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 / 3600, 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 (hr)')
+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)
+std_mass_conservation = np.std(mass_conservation_sum_encl1_encl2)
+# print("Standard deviation of mass conservation sum: ", std_mass_conservation)
+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
@@ -36,78 +89,57 @@
 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
-fig1 = plt.figure(figsize=[6, 5.5])
-ax1 = fig1.add_subplot(111)
-ax1.plot(TMAP8_time / 3600, TMAP8_pressure_enclosure_1, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
-ax1.plot(TMAP8_time / 3600, TMAP8_pressure_enclosure_2, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
-ax1.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
-ax1.set_xlabel('Time (hr)')
-ax1.set_ylabel('Pressure (Pa)')
-ax1.set_xlim(0, 3)
-ax1.set_ylim(bottom=0)
-ax1.legend(loc="best")
-ax1.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig1.savefig('ver-1kb_comparison_time.png', bbox_inches='tight', dpi=300)
+# Subplot 1 : Pressure vs time
 
-# Subplot 2: Solubility and concentration ratios vs time
-fig4 = plt.figure(figsize=[6, 5.5])
-ax4 = fig4.add_subplot(111)
-ax4.plot(TMAP8_time[1:] / 3600, concentration_ratio[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
-ax4.plot(TMAP8_time[1:] / 3600, [1] * len(TMAP8_time[1:]), label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
-ax4.set_yticks(np.arange(0, 3, 1))
-ax4.set_xlim(0,TMAP8_time.max() / 3600)
-ax4.set_xlabel('Time (hr)')
-ax4.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
-ax4.legend(loc="best")
-ax4.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig4.savefig('ver-1kb_concentration_ratio.png', bbox_inches='tight', dpi=300)
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
 
-# Subplot 3: Mass Conservation Sum Encl 1 and 2 vs Time
-fig3 = plt.figure(figsize=[6, 5.5])
-ax3 = fig3.add_subplot(111)
-ax3.plot(TMAP8_time / 3600, mass_conservation_sum_encl1_encl2, c='tab:blue')
-ax3.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
-ax3.set_xlabel('Time (hr)')
-ax3.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
-ax3.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig3.savefig('ver-1kb_mass_conservation.png', bbox_inches='tight', dpi=300)
-
-# Repeat the same for K=10/RT
-# Subplot 1 : Pressure vs time
-fig1_k10 = plt.figure(figsize=[6, 5.5])
-ax1_k10 = fig1_k10.add_subplot(111)
-ax1_k10.plot(TMAP8_time_k10 / 3600, TMAP8_pressure_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
-ax1_k10.plot(TMAP8_time_k10 / 3600, TMAP8_pressure_enclosure_2_k10, label=r"T$_2$ Enclosure 2", c='tab:blue', linestyle='-')
-ax1_k10.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.1e}'.format(val)))
-ax1_k10.set_xlabel('Time (hr)')
-ax1_k10.set_ylabel('Pressure (Pa)')
-ax1_k10.set_xlim(0, 3)
-ax1_k10.set_ylim(bottom=0)
-ax1_k10.legend(loc="best")
-ax1_k10.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig1_k10.savefig('ver-1kb_comparison_time_k10.png', bbox_inches='tight', dpi=300)
+ax.plot(TMAP8_time_k10 / 3600, TMAP8_pressure_enclosure_1_k10, label=r"T$_2$ Enclosure 1", c='tab:red', linestyle='-')
+ax.plot(TMAP8_time_k10 / 3600, 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 (hr)')
+ax.set_ylabel('Pressure (Pa)')
+ax.set_xlim(-0.001, TMAP8_time.max() / 3600)
+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
-fig4_k10 = plt.figure(figsize=[6, 5.5])
-ax4_k10 = fig4_k10.add_subplot(111)
-ax4_k10.plot(TMAP8_time[1:] / 3600, concentration_ratio_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
-ax4_k10.plot(TMAP8_time[1:] / 3600, [10] * len(TMAP8_time[1:]), label=r"Solubility Ratio (Analytical)", color='tab:red', linestyle='--')
-ax4_k10.set_yticks(np.arange(0, 21, 10))
-ax4_k10.set_xlim(0,TMAP8_time.max() / 3600)
-ax4_k10.set_xlabel('Time (hr)')
-ax4_k10.set_ylabel(r"Concentrations ratio $C_{\text{encl1}} / C_{\text{encl2}}$")
-ax4_k10.legend(loc="best")
-ax4_k10.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig4_k10.savefig('ver-1kb_concentration_ratio_k10.png', bbox_inches='tight', dpi=300)
+
+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:] / 3600, concentration_ratio_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-')
+ax.plot(TMAP8_time[1:] / 3600, 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() / 3600)
+ax.set_xlabel('Time (hr)')
+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
-fig3_k10 = plt.figure(figsize=[6, 5.5])
-ax3_k10 = fig3_k10.add_subplot(111)
-ax3_k10.plot(TMAP8_time_k10 / 3600, mass_conservation_sum_encl1_encl2_k10, c='tab:blue')
-ax3.yaxis.set_major_formatter(plt.FuncFormatter(lambda val, pos: '{:.3e}'.format(val)))
-ax3_k10.set_xlabel('Time (hr)')
-ax3_k10.set_ylabel(r"Mass Conservation Sum Encl 1 and 2 (mol/m$^3$)")
-ax3_k10.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
-fig3_k10.savefig('ver-1kb_mass_conservation_k10.png', bbox_inches='tight', dpi=300)
+
+fig = plt.figure(figsize=[6.5,5.5])
+gs = gridspec.GridSpec(1,1)
+ax = fig.add_subplot(gs[0])
+
+ax.plot(TMAP8_time_k10 / 3600, 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 (hr)')
+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)
+std_mass_conservation = np.std(mass_conservation_sum_encl1_encl2_k10)
+# print("Standard deviation of mass conservation sum: ", std_mass_conservation)
+fig.savefig('ver-1kb_mass_conservation_k10.png', bbox_inches='tight', dpi=300)