-
Notifications
You must be signed in to change notification settings - Fork 21
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
- Loading branch information
Showing
9 changed files
with
1,672 additions
and
0 deletions.
There are no files selected for viewing
1 change: 1 addition & 0 deletions
1
doc/content/verification_and_validation/figures/comparison_ver-1kc.py
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1 @@ | ||
| ../../../../test/tests/ver-1kc/comparison_ver-1kc.py |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,83 @@ | ||
| # ver-1kc | ||
|
|
||
| # Sieverts’ Law Boundaries with No Volumetric Source | ||
|
|
||
| ## General Case Description | ||
|
|
||
| Two enclosures are separated by a membrane that allows diffusion according to Sievert'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. | ||
|
|
||
| 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 Sieverts' law, $n=0.5$. | ||
|
|
||
| ## Results | ||
|
|
||
| Two subcases are considered. In the first subcase, we assume that $K=1/RT$, which is expected to lead to $C_1=C_2^{1/2}$ at equilibrium. In the second, $K=10/RT$, which is expected to lead to $C_1 = 10 C_2^{1/2}$. This second case is added to exercise TMAP8 in a case with a concentration jump. | ||
| In the first subcase, the pressure evolution in both enclosures is shown in [ver-1kc_comparison_time] as a function of time. Both pressures find equilibrium, which is consistent with $C_1 = K (RT C_2)^n$ for $K=1/RT$ and $n=0.5$. The concentration ratio between enclosures 1 and 2 in [ver-1kc_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-1kc_mass_conservation], mass is conserved between the two enclosures over time, with a variation in mass of only $1.3$ %. 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-1kc.py | ||
| image_name=ver-1kc_comparison_time.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_comparison_time | ||
| caption=Equilibration of species pressures under Sieverts' law for $K = 1/RT$. | ||
|
|
||
| !media comparison_ver-1kc.py | ||
| image_name=ver-1kc_concentration_ratio.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_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-1kc.py | ||
| image_name=ver-1kc_mass_conservation.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_mass_conservation | ||
| caption=Total mass conservation across both enclosures over time for $K = 1/RT$. | ||
|
|
||
| In the second subcase, as illustrated in [ver-1kc_comparison_time_k10], the pressure jump maintains a ratio of $C_1/C_2^{1/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-1kc_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-1kc_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time, with a variation in mass of only $0.4$ %. As in the previous case, this variation in mass can be further reduced by refining the mesh. | ||
|
|
||
| !media comparison_ver-1kc.py | ||
| image_name=ver-1kc_comparison_time_k10.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_comparison_time_k10 | ||
| caption=Pressures jump of species under Sieverts' law for $K = 10/RT$. | ||
|
|
||
| !media comparison_ver-1kc.py | ||
| image_name=ver-1kc_concentration_ratio_k10.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_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-1kc.py | ||
| image_name=ver-1kc_mass_conservation_k10.png | ||
| style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto | ||
| id=ver-1kc_mass_conservation_k10 | ||
| caption=Total mass conservation across both enclosures over time for $K = 10/RT$. | ||
|
|
||
| ## Input files | ||
|
|
||
| !style halign=left | ||
| The input file for this case can be found at [/ver-1kc.i], which is also used as tests in TMAP8 at [/ver-1kc/tests]. | ||
|
|
||
| !bibtex bibliography |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -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-1kc/gold/ver-1kc_out_k1.csv" | ||
| else: # if in test folder | ||
| csv_folder = "./gold/ver-1kc_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-1kc_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-1kc_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-1kc_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-1kc/gold/ver-1kc_out_k10.csv" | ||
| else: # if in test folder | ||
| csv_folder_k10 = "./gold/ver-1kc_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_k10.max()) | ||
| ax.set_ylim(bottom=0) | ||
| ax.legend(loc="best") | ||
| ax.grid(which='major', color='0.65', linestyle='--', alpha=0.3) | ||
| fig.savefig('ver-1kc_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_k10[1:]) | ||
| ax.plot(TMAP8_time_k10[1:], concentration_ratio_k10[1:], label=r"Concentration Ratio (TMAP8)", color='tab:blue', linestyle='-') | ||
| ax.plot(TMAP8_time_k10[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_k10.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_k10.max() / 7200 | ||
| y_pos = 0.9 * ax.get_ylim()[1] | ||
| ax.text(x_pos, y_pos, 'RMSPE = %.3f ' % RMSPE + '%', fontweight='bold') | ||
| fig.savefig('ver-1kc_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-1kc_mass_conservation_k10.png', bbox_inches='tight', dpi=300) | ||
|
|
Oops, something went wrong.