Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Ver 1kb #201

Open
wants to merge 8 commits into
base: devel
Choose a base branch
from
Open

Ver 1kb #201

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
8 commits
Select commit Hold shift + click to select a range
06038bc
(Ref. #12)
Lee01Atom Oct 25, 2024
bda1e65
delete trailing whitespaces
Lee01Atom Oct 25, 2024
7365d5e
delete trailing whitespaces
Lee01Atom Oct 25, 2024
ce50a86
add preliminary changes
Lee01Atom Oct 27, 2024
381e41d
Apply suggestions from code review
Lee01Atom Oct 30, 2024
72bfe63
add mass conservation and sorption law verification
Lee01Atom Oct 30, 2024
913e354
Apply suggestions from code review
Lee01Atom Oct 31, 2024
bd6bdc6
update of the documentation
Lee01Atom Nov 1, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions doc/content/verification_and_validation/figures/comparison_ver-1kb.py
View file
Open in desktop
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This one is not fixed yet.

Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
/Users/lee/projects/TMAP8/test/tests/ver-1kb/comparison_ver-1kb.py
1 change: 1 addition & 0 deletions doc/content/verification_and_validation/index.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
83 changes: 83 additions & 0 deletions doc/content/verification_and_validation/ver-1kb.md
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,83 @@
# 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
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You should document the values used in this case, like the system dimensions, the initial pressures, the diffusivities and solubilities used, etc. With all the corresponding units, of course.

Make sure to also add the equation for diffusion.


This verification problem is taken from [!cite](ambrosek2008verification).

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 the following equations:

\begin{equation}
\frac{\partial C_1}{\partial t} = D \nabla^2 C_1
Copy link
Collaborator

@simopier simopier Nov 4, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
\frac{\partial C_1}{\partial t} = D \nabla^2 C_1
\frac{\partial C_1}{\partial t} = \nabla D \nabla C_1

Copy link
Collaborator

@simopier simopier Nov 4, 2024 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

D \nabla^2 C_1 is only valid if the diffusion coefficient does not depend on the position. The expression I am proposing is more general.

\end{equation}

\begin{equation}
\frac{\partial C_1}{\partial t} = D \nabla^2 C_1
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
\frac{\partial C_1}{\partial t} = D \nabla^2 C_1
\frac{\partial C_1}{\partial t} = - \nabla D \nabla C_1

\end{equation}

where $C_1$, $C_2$ represent the concentration fields in enclosures 1 and 2 respectively, 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( \frac{C_2 RT}{n} \right)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
C_1 = K P_2^n = K \left( \frac{C_2 RT}{n} \right)
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 the Henry’s law, $n=1$.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
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 the Henry’s law, $n=1$.
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, it is assumed that the pressures in the two enclosures are in equilibrium, implying $K = 1/RT$.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Two subcases are considered. In the first subcase, it is assumed that the pressures in the two enclosures are in equilibrium, implying $K = 1/RT$.
Two subcases are considered. In the first subcase, we assume that $K = 1/RT$.

Consistent with the results from TMAP7, the pressure evolution in both enclosures is shown in [ver-1kb_comparison_time] as a function of time, confirming the equilibrium between the pressures in enclosures 1 and 2. The linear regression in [ver-1kb_comparison_concentration] demonstrates that the concentration values at the interface comply with the sorption law, with a proportionality coefficient consistent with the solubility value $K = \approx 2.4 \times 10^{-4}$. As shown in [ver-1kb_mass_conservation], mass is conserved between the two enclosures over time.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
Consistent with the results from TMAP7, the pressure evolution in both enclosures is shown in [ver-1kb_comparison_time] as a function of time, confirming the equilibrium between the pressures in enclosures 1 and 2. The linear regression in [ver-1kb_comparison_concentration] demonstrates that the concentration values at the interface comply with the sorption law, with a proportionality coefficient consistent with the solubility value $K = \approx 2.4 \times 10^{-4}$. As shown in [ver-1kb_mass_conservation], mass is conserved between the two enclosures over time.
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 pressure find equilibrium and become equal, which is consistent with $C_1 = K *RT * C_2^n$ for $K=1/RT$ and $n=1$. The linear regression in [ver-1kb_comparison_concentration] demonstrates that the concentration values at the interface comply with the sorption law, with a proportionality coefficient consistent with the solubility value $K=1/RT$. As shown in [ver-1kb_mass_conservation], mass is conserved between the two enclosures over time.

Since you know the exact value of K, it's better to use it rather than use an approximate value.


!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_comparison_concentration.png
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
id=ver-1kb_comparison_concentration
caption=Concentration in enclosure 1 as a function of pressure in enclosure 2 at the interface for $K=1/RT$. Validation of the sorption law across the interface.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
caption=Concentration in enclosure 1 as a function of pressure in enclosure 2 at the interface for $K=1/RT$. Validation of the sorption law across the interface.
caption=Concentration in enclosure 1 as a function of pressure in enclosure 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.

Validation requires comparison against experimental data, verification is comparison against analytical solution. This here is a verification, not a validation.


!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$ prevents the pressures in the two enclosures from reaching equilibrium. 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$ for $K=10/RT$. The linear regression in [ver-1kb_comparison_concentration_k10] confirms that the concentration values at the interface adhere to the sorption law, with a proportionality coefficient aligned with the solubility value $K = \approx 2.4 \times 10^{-5}$. Additionally, [ver-1kb_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
In the second subcase, the sorption law with $K = 10/RT$ prevents the pressures in the two enclosures from reaching equilibrium. 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$ for $K=10/RT$. The linear regression in [ver-1kb_comparison_concentration_k10] confirms that the concentration values at the interface adhere to the sorption law, with a proportionality coefficient aligned with the solubility value $K = \approx 2.4 \times 10^{-5}$. Additionally, [ver-1kb_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time.
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 linear regression in [ver-1kb_comparison_concentration_k10] confirms that the concentration values at the interface adhere to the sorption law, with a proportionality coefficient aligned with the solubility value $K=10/RT$. Additionally, [ver-1kb_mass_conservation_k10] verifies that mass is conserved between the two enclosures over time.

The pressure ARE at equilibrium, even if they are not equal.


!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 pressures under Henry’s law for $K=10/RT$.

!media comparison_ver-1kb.py
image_name=ver-1kb_comparison_concentration_k10.png
style=width:50%;margin-bottom:2%;margin-left:auto;margin-right:auto
id=ver-1kb_comparison_concentration_k10
caption=Concentration in enclosure 1 as a function of pressure in enclosure 2 at the interface for $K=10/RT$. Validation of the sorption law across the interface.
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
caption=Concentration in enclosure 1 as a function of pressure in enclosure 2 at the interface for $K=10/RT$. Validation of the sorption law across the interface.
caption=Concentration in enclosure 1 as a function of pressure in enclosure 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$.

## 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
111 changes: 111 additions & 0 deletions test/tests/ver-1kb/comparison_ver-1kb.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,111 @@
import numpy as np
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The new figures created here need to be added to the documentation and commented there.

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.csv"
else: # if in test folder
csv_folder = "./gold/ver-1kb_out.csv"
Comment on lines +13 to +15
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
csv_folder = "../../../../test/tests/ver-1kb/gold/ver-1kb_out.csv"
else: # if in test folder
csv_folder = "./gold/ver-1kb_out.csv"
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

# 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

# Subplot 1: Time vs Pressure
fig1 = plt.figure(figsize=[6, 5.5])
ax1 = fig1.add_subplot(111)
ax1.plot(TMAP8_time / 3600, TMAP8_pressure_enclosure_1, label=r"H2 Encl 1", c='tab:red', linestyle='dotted')
ax1.plot(TMAP8_time / 3600, TMAP8_pressure_enclosure_2, label=r"H2 Encl 2", c='tab:blue', linestyle='dotted')
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')

# Subplot 2: Pressure Encl 2 vs Concentration Encl 1 at interface
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make sure to also calculate the difference between the expected final results and the calculated final results, and discuss how accurate TMAP8 is in the documentation.

k = np.sum(concentration_enclosure_1_at_interface * pressure_enclosure_2_at_interface) / np.sum(pressure_enclosure_2_at_interface ** 2)
line_fit = k * pressure_enclosure_2_at_interface[1:]
fig2 = plt.figure(figsize=[6, 5.5])
ax2 = fig2.add_subplot(111)
ax2.plot(pressure_enclosure_2_at_interface[1:], concentration_enclosure_1_at_interface[1:], marker='o', linestyle='None', c='tab:blue')
ax2.plot(pressure_enclosure_2_at_interface[1:], line_fit, label=f'Fit: y = {k:.2e}x\n', c='tab:red', linestyle='--')
ax2.set_xlabel(r"Pressure Encl 2 at interface (Pa)")
ax2.set_ylabel(r"Concentration Encl 1 at interface (mol/m^3)")
ax2.legend(loc="best")
ax2.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
fig2.savefig('ver-1kb_comparison_concentration.png', bbox_inches='tight')

# 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, 'o', color='tab:blue', markersize=0.7)
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')

# Repeat the same for K=10/RT
# Subplot 1 : Time vs Pressure
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"H2 Encl 1", c='tab:red', linestyle='-')
ax1_k10.plot(TMAP8_time_k10 / 3600, TMAP8_pressure_enclosure_2_k10, label=r"H2 Encl 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')

# Subplot 2 : Pressure Encl 2 vs Concentration Encl 1 at interface
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Make sure to also calculate the difference between the expected final results and the calculated final results, and discuss how accurate TMAP8 is in the documentation.

k_k10 = np.sum(concentration_enclosure_1_at_interface_k10 * pressure_enclosure_2_at_interface_k10) / np.sum(pressure_enclosure_2_at_interface_k10 ** 2)
line_fit_k10 = k_k10 * pressure_enclosure_2_at_interface_k10[1:]
fig2_k10 = plt.figure(figsize=[6, 5.5])
ax2_k10 = fig2_k10.add_subplot(111)
ax2_k10.plot(pressure_enclosure_2_at_interface_k10[1:], concentration_enclosure_1_at_interface_k10[1:], marker='o', linestyle='None', c='tab:blue')
ax2_k10.plot(pressure_enclosure_2_at_interface_k10[1:], line_fit_k10, label=f'Fit: y = {k_k10:.2e}x\n', c='tab:red', linestyle='--')
ax2_k10.set_xlabel(r"Pressure Encl 2 at interface (Pa)")
ax2_k10.set_ylabel(r"Concentration Encl 1 at interface (mol/m^3)")
ax2_k10.legend(loc="best")
ax2_k10.grid(which='major', color='0.65', linestyle='--', alpha=0.3)
fig2_k10.savefig('ver-1kb_comparison_concentration_k10.png', bbox_inches='tight')

# 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, 'o', color='tab:blue', markersize=0.7)
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')

Loading