Skip to content
Snippets Groups Projects
Verified Commit e3ba41cd authored by Lars Bilke's avatar Lars Bilke
Browse files

[nb] Added placeholder image for gas (solute) diffusion.

parent 1ecc1e02
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Tests/Data/TH2M/H/diffusion/diffusion.ipynb
+ 10
− 4
View file @ e3ba41cd
%% Cell type:raw id:94e3d277 tags: %% Cell type:raw id:94e3d277 tags:
title = "Solute Diffusion" title = "Gas Diffusion"
date = "2022-10-19" date = "2022-10-19"
author = "Norbert Grunwald" author = "Norbert Grunwald"
notebook = "TH2M/H2/diffusion/diffusion.ipynb" notebook = "TH2M/H2/diffusion/diffusion.ipynb"
image = "" image = "figures/placeholder_diffusion.png"
web_subsection = "th2m" web_subsection = "th2m"
coupling = "h"
<!--eofm--> <!--eofm-->
%% Cell type:markdown id:40f1eb67 tags: %% Cell type:markdown id:40f1eb67 tags:
|<div style="width:300px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:300px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg" width="300"/></div>| |<div style="width:300px"><img src="https://www.ufz.de/static/custom/weblayout/DefaultInternetLayout/img/logos/ufz_transparent_de_blue.png" width="300"/></div>|<div style="width:300px"><img src="https://discourse.opengeosys.org/uploads/default/original/1X/a288c27cc8f73e6830ad98b8729637a260ce3490.png" width="300"/></div>|<div style="width:330px"><img src="https://tu-freiberg.de/sites/default/files/media/freiberger-alumni-netzwerk-6127/wbm_orig_rgb_0.jpg" width="300"/></div>|
|---|---|--:| |---|---|--:|
%% Cell type:markdown id:404dd348 tags: %% Cell type:markdown id:404dd348 tags:
# 1D Linear Diffusion # 1D Linear Diffusion
## Analytical Solution ## Analytical Solution
The solution of the diffusion equation for a point x in time t reads: The solution of the diffusion equation for a point x in time t reads:
$ c(x,t) = \left(c_0-c_i\right)\operatorname{erfc}\left(\frac{x}{\sqrt{4Dt}}\right)+c_i$ $ c(x,t) = \left(c_0-c_i\right)\operatorname{erfc}\left(\frac{x}{\sqrt{4Dt}}\right)+c_i$
where $c$ is the concentration of a solute in mol$\cdot$m$^{-3}$, $c_i$ and $c_b$ are initial and boundary concentrations; $D$ is the diffusion coefficient of the solute in water, x and t are location and time of solution. where $c$ is the concentration of a solute in mol$\cdot$m$^{-3}$, $c_i$ and $c_b$ are initial and boundary concentrations; $D$ is the diffusion coefficient of the solute in water, x and t are location and time of solution.
%% Cell type:code id:6becccb7 tags: %% Cell type:code id:6becccb7 tags:
``` python ``` python
import numpy as np import numpy as np
from scipy.special import erfc from scipy.special import erfc
# Analytical solution of the diffusion equation # Analytical solution of the diffusion equation
def Diffusion(x,t): def Diffusion(x,t):
if type(t) != int and type(t) != np.float64 and type(t) != float: if type(t) != int and type(t) != np.float64 and type(t) != float:
# In order to avoid a division by zero, the time field is increased # In order to avoid a division by zero, the time field is increased
# by a small time unit at the start time (t=0). This should have no # by a small time unit at the start time (t=0). This should have no
# effect on the result. # effect on the result.
tiny = np.finfo(np.float64).tiny tiny = np.finfo(np.float64).tiny
t[t < tiny] = tiny t[t < tiny] = tiny
d = np.sqrt(4*D*t) d = np.sqrt(4*D*t)
e = (c_b - c_i)*erfc(x/d)+c_i e = (c_b - c_i)*erfc(x/d)+c_i
return e return e
# Utility-function transforming mass fraction into conctration # Utility-function transforming mass fraction into conctration
def concentration(xm_WL): def concentration(xm_WL):
xm_CL = 1. - xm_WL xm_CL = 1. - xm_WL
return xm_CL / beta_c return xm_CL / beta_c
``` ```
%% Cell type:markdown id:e30add18 tags: %% Cell type:markdown id:e30add18 tags:
### Material properties and problem specification ### Material properties and problem specification
%% Cell type:code id:eb9e5db4 tags: %% Cell type:code id:eb9e5db4 tags:
``` python ``` python
# Henry-coefficient and compressibility of solution # Henry-coefficient and compressibility of solution
H = 7.65e-6 H = 7.65e-6
beta_c = 2.e-6 beta_c = 2.e-6
# Diffusion coefficient # Diffusion coefficient
D = 1.0e-9 D = 1.0e-9
# Boundary and initial gas pressures # Boundary and initial gas pressures
pGR_b = 9e5 pGR_b = 9e5
pGR_i = 1e5 pGR_i = 1e5
# Boundary and initial concentration # Boundary and initial concentration
c_b = concentration(1. - (beta_c*H*pGR_b)) c_b = concentration(1. - (beta_c*H*pGR_b))
c_i = concentration(1. - (beta_c*H*pGR_i)) c_i = concentration(1. - (beta_c*H*pGR_i))
``` ```
%% Cell type:markdown id:beb845b8 tags: %% Cell type:markdown id:beb845b8 tags:
## Numerical Solution ## Numerical Solution
%% Cell type:code id:57201254 tags: %% Cell type:code id:57201254 tags:
``` python ``` python
import os import os
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out') out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')
if not os.path.exists(out_dir): if not os.path.exists(out_dir):
os.makedirs(out_dir) os.makedirs(out_dir)
``` ```
%% Cell type:code id:e4daf901 tags: %% Cell type:code id:e4daf901 tags:
``` python ``` python
from ogs6py.ogs import OGS from ogs6py.ogs import OGS
model=OGS(INPUT_FILE="diffusion.prj", PROJECT_FILE=f"{out_dir}/modified.prj") model=OGS(INPUT_FILE="diffusion.prj", PROJECT_FILE=f"{out_dir}/modified.prj")
model.replace_text(1e7, xpath="./time_loop/processes/process/time_stepping/t_end") model.replace_text(1e7, xpath="./time_loop/processes/process/time_stepping/t_end")
model.replace_text(5e4, xpath="./time_loop/processes/process/time_stepping/timesteps/pair/delta_t") model.replace_text(5e4, xpath="./time_loop/processes/process/time_stepping/timesteps/pair/delta_t")
# Write every timestep # Write every timestep
model.replace_text(1, xpath="./time_loop/output/timesteps/pair/each_steps") model.replace_text(1, xpath="./time_loop/output/timesteps/pair/each_steps")
model.write_input() model.write_input()
# Run OGS # Run OGS
model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .") model.run_model(logfile=f"{out_dir}/out.txt", args=f"-o {out_dir} -m .")
``` ```
%% Output %% Output
OGS finished with project file modified.prj. OGS finished with project file modified.prj.
Execution took 44.86126470565796 s Execution took 44.86126470565796 s
%% Cell type:code id:d4d6efe6 tags: %% Cell type:code id:d4d6efe6 tags:
``` python ``` python
#Colors #Colors
cls1 = ['#4a001e', '#731331', '#9f2945', '#cc415a', '#e06e85', '#ed9ab0'] cls1 = ['#4a001e', '#731331', '#9f2945', '#cc415a', '#e06e85', '#ed9ab0']
cls2 = ['#0b194c', '#163670', '#265191', '#2f74b3', '#5d94cb', '#92b2de'] cls2 = ['#0b194c', '#163670', '#265191', '#2f74b3', '#5d94cb', '#92b2de']
``` ```
%% Cell type:code id:a794deb5 tags: %% Cell type:code id:a794deb5 tags:
``` python ``` python
import vtuIO import vtuIO
pvdfile = vtuIO.PVDIO(f"{out_dir}/result_diffusion.pvd", dim=2) pvdfile = vtuIO.PVDIO(f"{out_dir}/result_diffusion.pvd", dim=2)
# Get all written timesteps # Get all written timesteps
time = pvdfile.timesteps time = pvdfile.timesteps
# Select individual timesteps for c vs. x plots for plotting # Select individual timesteps for c vs. x plots for plotting
time_steps = [1e6, 2e6, 4e6, 6e6, 8e6, 1e7] time_steps = [1e6, 2e6, 4e6, 6e6, 8e6, 1e7]
# 'Continuous' space axis for c vs. x plots for plotting # 'Continuous' space axis for c vs. x plots for plotting
length = np.linspace(0,1.0,101) length = np.linspace(0,1.0,101)
# Draws a line through the domain for sampling results # Draws a line through the domain for sampling results
x_axis=[(i,0,0) for i in length] x_axis=[(i,0,0) for i in length]
# Discrete locations for c vs. t plots # Discrete locations for c vs. t plots
location = [0.01,0.05,0.1,0.2,0.5,1.0] location = [0.01,0.05,0.1,0.2,0.5,1.0]
``` ```
%% Cell type:code id:9bee423e tags: %% Cell type:code id:9bee423e tags:
``` python ``` python
# The sample locations have to be converted into a 'dict' for vtuIO # The sample locations have to be converted into a 'dict' for vtuIO
observation_points = dict(('x='+str(x),(x,0.0,0.0)) for x in location) observation_points = dict(('x='+str(x),(x,0.0,0.0)) for x in location)
# Samples concentration field at the observation points for all timesteps # Samples concentration field at the observation points for all timesteps
c_over_t_at_x = pvdfile.read_time_series('xmWL', observation_points) c_over_t_at_x = pvdfile.read_time_series('xmWL', observation_points)
for key in c_over_t_at_x: for key in c_over_t_at_x:
x = c_over_t_at_x[key]; x = c_over_t_at_x[key];
c_over_t_at_x[key] = concentration(x) c_over_t_at_x[key] = concentration(x)
# Samples concentration field along the domain at certain timesteps # Samples concentration field along the domain at certain timesteps
c_over_x_at_t = [] c_over_x_at_t = []
for t in range(0,len(time_steps)): for t in range(0,len(time_steps)):
c_over_x_at_t.append(concentration(pvdfile.read_set_data(time_steps[t], "xmWL", pointsetarray=x_axis))) c_over_x_at_t.append(concentration(pvdfile.read_set_data(time_steps[t], "xmWL", pointsetarray=x_axis)))
``` ```
%% Cell type:code id:c93c9252 tags: %% Cell type:code id:c93c9252 tags:
``` python ``` python
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] = (14, 4) plt.rcParams['figure.figsize'] = (14, 4)
# Plot of concentration vs. time at different locations # Plot of concentration vs. time at different locations
fig1, (ax1, ax2) = plt.subplots(1, 2) fig1, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_xlabel("$t$ / s", fontsize=12) ax1.set_xlabel("$t$ / s", fontsize=12)
ax1.set_ylabel("$c$ / mol m$^{-3}$", fontsize=12) ax1.set_ylabel("$c$ / mol m$^{-3}$", fontsize=12)
ax2.set_xlabel("$t$ / s", fontsize=12) ax2.set_xlabel("$t$ / s", fontsize=12)
ax2.set_ylabel("$\epsilon_\mathrm{abs}$ / mol m$^{-3}$", fontsize=12) ax2.set_ylabel("$\epsilon_\mathrm{abs}$ / mol m$^{-3}$", fontsize=12)
label_x = [] label_x = []
for key, c in c_over_t_at_x.items(): for key, c in c_over_t_at_x.items():
x = observation_points[key][0] x = observation_points[key][0]
label_x.append(key+r" m") label_x.append(key+r" m")
# numerical solution # numerical solution
ax1.plot(time, c_over_t_at_x[key], color=cls1[location.index(x)], linewidth=3, ax1.plot(time, c_over_t_at_x[key], color=cls1[location.index(x)], linewidth=3,
linestyle="--") linestyle="--")
# analytical solution # analytical solution
ax1.plot(time, Diffusion(x,time), color=cls1[location.index(x)], linewidth=2, ax1.plot(time, Diffusion(x,time), color=cls1[location.index(x)], linewidth=2,
linestyle="-") linestyle="-")
# absolute error # absolute error
err_abs = Diffusion(x,time) - c_over_t_at_x[key] err_abs = Diffusion(x,time) - c_over_t_at_x[key]
ax2.plot(time, err_abs, color=cls1[location.index(x)], linewidth=1, linestyle="-", label=key+r" m") ax2.plot(time, err_abs, color=cls1[location.index(x)], linewidth=1, linestyle="-", label=key+r" m")
# Hack to force a custom legend: # Hack to force a custom legend:
from matplotlib.lines import Line2D from matplotlib.lines import Line2D
custom_lines = [] custom_lines = []
for i in range(0,6): for i in range(0,6):
custom_lines.append(Line2D([0], [0], color=cls1[i], lw=4)) custom_lines.append(Line2D([0], [0], color=cls1[i], lw=4))
custom_lines.append(Line2D([0], [0], color='black', lw=3, linestyle="--")) custom_lines.append(Line2D([0], [0], color='black', lw=3, linestyle="--"))
custom_lines.append(Line2D([0], [0], color='black', lw=2, linestyle="-")) custom_lines.append(Line2D([0], [0], color='black', lw=2, linestyle="-"))
label_x.append('OGS-TH2M') label_x.append('OGS-TH2M')
label_x.append('analytical') label_x.append('analytical')
ax1.legend(custom_lines, label_x, loc='right') ax1.legend(custom_lines, label_x, loc='right')
ax2.legend() ax2.legend()
fig1.savefig(f"{out_dir}/diffusion_c_vs_t.pdf") fig1.savefig(f"{out_dir}/diffusion_c_vs_t.pdf")
# Plot of concentration vs. location at different times # Plot of concentration vs. location at different times
fig1, (ax1, ax2) = plt.subplots(1, 2) fig1, (ax1, ax2) = plt.subplots(1, 2)
ax1.set_xlabel("$x$ / m", fontsize=12) ax1.set_xlabel("$x$ / m", fontsize=12)
ax1.set_ylabel("$c$ / mol m$^{-3}$", fontsize=12) ax1.set_ylabel("$c$ / mol m$^{-3}$", fontsize=12)
ax1.set_xlim(0,0.4) ax1.set_xlim(0,0.4)
ax2.set_xlabel("$x$ / m", fontsize=12) ax2.set_xlabel("$x$ / m", fontsize=12)
ax2.set_ylabel("$\epsilon$ / mol $m^{-3}$", fontsize=12) ax2.set_ylabel("$\epsilon$ / mol $m^{-3}$", fontsize=12)
ax2.set_xlim(0,0.4) ax2.set_xlim(0,0.4)
# Plot concentration over domain at five moments # Plot concentration over domain at five moments
label_t = [] label_t = []
for t in range(0,len(time_steps)): for t in range(0,len(time_steps)):
s = r"$t=$"+str(time_steps[t]/1e6)+"$\,$Ms" s = r"$t=$"+str(time_steps[t]/1e6)+"$\,$Ms"
label_t.append(s) label_t.append(s)
# numerical solution # numerical solution
ax1.plot(length, c_over_x_at_t[t], color=cls2[t], linewidth=3, linestyle="--") ax1.plot(length, c_over_x_at_t[t], color=cls2[t], linewidth=3, linestyle="--")
# analytical solution # analytical solution
ax1.plot(length, Diffusion(length,time_steps[t]),color=cls2[t], linewidth=2, \ ax1.plot(length, Diffusion(length,time_steps[t]),color=cls2[t], linewidth=2, \
linestyle="-") linestyle="-")
# absolute error # absolute error
err_abs = Diffusion(length,time_steps[t]) - c_over_x_at_t[t] err_abs = Diffusion(length,time_steps[t]) - c_over_x_at_t[t]
ax2.plot(length, err_abs, color=cls2[t], linewidth=1, linestyle="-", label=s) ax2.plot(length, err_abs, color=cls2[t], linewidth=1, linestyle="-", label=s)
custom_lines = [] custom_lines = []
for i in range(0,6): for i in range(0,6):
custom_lines.append(Line2D([0], [0], color=cls2[i], lw=4)) custom_lines.append(Line2D([0], [0], color=cls2[i], lw=4))
custom_lines.append(Line2D([0], [0], color='black', lw=3, linestyle="--")) custom_lines.append(Line2D([0], [0], color='black', lw=3, linestyle="--"))
custom_lines.append(Line2D([0], [0], color='black', lw=2, linestyle="-")) custom_lines.append(Line2D([0], [0], color='black', lw=2, linestyle="-"))
label_t.append('OGS-TH2M') label_t.append('OGS-TH2M')
label_t.append('analytical') label_t.append('analytical')
ax1.legend(custom_lines, label_t, loc='right') ax1.legend(custom_lines, label_t, loc='right')
ax2.legend() ax2.legend()
fig1.savefig(f"{out_dir}/diffusion_c_vs_x.pdf") fig1.savefig(f"{out_dir}/diffusion_c_vs_x.pdf")
``` ```
%% Output %% Output
%% Cell type:markdown id:041ffcef tags: %% Cell type:markdown id:041ffcef tags:
The numerical approximation approaches the exact solution quite well. Deviations can be reduced if the resolution of the temporal discretisation is increased. The numerical approximation approaches the exact solution quite well. Deviations can be reduced if the resolution of the temporal discretisation is increased.
......
web/content/docs/benchmarks/th2m/th2m_1d-diffusion/index.md deleted 100644 → 0
+ 0
− 15
View file @ 1ecc1e02
+++
author = "Norbert Grunwald"
date = "2018-06-20T14:37:58+01:00"
title = "Gas Diffusion"
image = "placeholder_diffusion.png"
coupling = "h"
+++
{{< data-link >}}
## Problem description
Todo
## Results and evaluation
web/content/docs/benchmarks/th2m/th2m_1d-diffusion/test.png deleted 100644 → 0
+ 0
− 0
View file @ 1ecc1e02
web/content/docs/benchmarks/th2m/th2m_1d-diffusion/test.png

94.1 KiB

0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment