Decay-chain problem

import os
import time

import ogs6py
import vtuIO

import numpy as np
from scipy import special

import h5py

import matplotlib.pyplot as plt
from matplotlib.pyplot import cm
from IPython.display import display, Image

1. Decay-chain problem

1.1 Problem description

This benchmark is meant to model the migration of radionuclides in the Curium-247 decay chain through a semi-infinite porous column. The diagram below maps the Curium-247 decay chain which contains 6 radionuclides before ending with Actinium-227.

data_dir = os.environ.get('OGS_DATA_DIR', '../../../..')
notebook_dir = os.path.abspath(f"{data_dir}/Parabolic/ComponentTransport/ReactiveTransport/DecayChain")
display(Image(filename=f"{notebook_dir}/chains.png", width=600))

png

According to the mass conservative law, the diffusive transport of each member of the Curium-247 decay chain can be described by

\[ \frac{\partial \left( \phi c_{i} \right)}{\partial t} = \nabla \cdot \left(\phi \mathrm{D_p} \nabla c_{i} \right) + \phi k_{i-1} c_{i-1} - \phi k_{i} c_{i}, \quad \forall i = 1,2,…,n, \]

where \(c_{i}\) [mol/m\(^3\)] is the radionuclide concentration, \(\phi\) [-] is the porosity, \(\mathrm{D_p}\) [m\(^2\)/s] is the pore diffusion coefficient, and \(k_{i}\) [1/s] is the first-order decay constant of a radionuclide, i.e.:

\[ k = \mathrm{ln} 2 \, / \, t_{1/2}, \]

where \(t_{1/2}\) [s] is its half life.

1.2 Model setups

In the 1-D model, the computational domain is set to be 20 meters long, which is sufficient to ensure that the concentration profile does not reach the right-hand-side boundary during simulation. The simulated time is 100 thousand years. Initially, the entire domain is assumed to be solute free, i.e., \(c_{\mathrm{ini}}(x, 0) = 0\) mol/m\(^3\). The inlet concentration for all the radionuclides is held at 10\(^{-3}\) mol/m\(^3\) throughout the simulation, i.e., \(c(0, t) = 1\) mol/L. The spatial domain is discretized by linear line elements and the near-inlet grids are refined in order to achieve sufficiently high accuracy. A fixed time step size of one hundred years is used in the simulation. The linearized governing equation system is iteratively solved using the standard Picard iteration method.

The table below summarizes the transport parameters and spatio-temporal discretization parameters used in the simulation.

Parameter Value Unit
Porosity \(\phi\) 0.12 -
Pore diffusion coefficient \(\mathrm{D_p}\) 1e-11 m\(^2\)/s
Time step size \(\Delta t\) 1e2 year
Grid size \(\Delta x\) 0.00512 - 1.985 m

The following table lists the dataset for the half-life of the radionuclides

Radionuclide Half-life [year] First-order decay constant [1/year]
Cm-247 1.56\(\times\)10\(^{7}\) 4.44\(\times\)10\(^{-8}\)
Am-243 7.37\(\times\)10\(^{3}\) 9.40\(\times\)10\(^{-5}\)
Pu-239 2.41\(\times\)10\(^{4}\) 2.88\(\times\)10\(^{-5}\)
U-235 7.04\(\times\)10\(^{8}\) 9.84\(\times\)10\(^{-10}\)
Pa-231 3.28\(\times\)10\(^{4}\) 2.11\(\times\)10\(^{-5}\)
Ac-227 21.773 3.18\(\times\)10\(^{-5}\)

Notes: The half-life dataset is sourced from Report GRS-289.

1.3 Results

Analytical solution

The analytical solution for the diffusive transport of radionuclides that are degrading sequentially with a series of first-order reactions is given in Sun et al. (1999):

\[ c_{i} = a_{i} - \sum_{j=1}^{i-1} \prod_{l=j}^{i=1} \frac{k_{l}}{k_{l} - k_{i}} c_{j}, \quad \forall i = 1,2,…,n, \]

Below a python script is presented, through which the concentration profiles of radionuclides in the Curium-247 decay chain are analytically calculated over a period of 10\(^5\) years.

def computeProduct(j, i, k, c_inlet):
    value = 1
    for l in range(j, i):
        value *= k[l] / (k[l] - k[i]) * c_inlet[j]
    
    return value

def computeInitialAuxiliaryVariable(c_inlet, k):
    a_inlet = np.empty((0))
    
    for i in range(len(c_inlet)):
        value = c_inlet[i]
        if (i > 0):
            for j in range(0, i):
                value += computeProduct(j, i, k, c_inlet)
        a_inlet = np.append(a_inlet, value)
    
    return a_inlet

def computeAnalyticalSolution(x, t, c_0, k, v, D):
    t *= 3.1536e7 #unit conversion from year to second
    
    beta = (v**2/4/D**2 + k/D)**0.5
    c_t = c_0/2 * np.exp(v*x/2/D) * (np.exp(-beta*x) \
        * special.erfc((x-(v**2+4*k*D)**0.5*t)/2/(D*t)**0.5) \
        + np.exp(beta*x) * special.erfc((x+(v**2+4*k*D)**0.5*t)/2/(D*t)**0.5))
    
    return c_t

###Model parameters###
#Diffusion coefficient [m2/s]
D = 1e-11
#Pore water velocity [m/s]
v = 0
#Half life [year]
radionuclides = np.array(["[Cm-247]", "[Am-243]", "[Pu-239]", "[U-235]", "[Pa-231]", "[Ac-227]"])
half_lifes = np.array([1.56e7, 7.37e3, 2.41e4, 7.04e8, 3.28e4, 21.773])
#First-order decay constant [1/s]
k = dict([(radionuclide, np.log(2)/half_life/3.1536e7) for radionuclide, half_life in zip(radionuclides, half_lifes)])

###Initial and boundary conditions###
c_inlet = np.ones(6)
a_inlet = dict([(radionuclide, a) for radionuclide, a in zip(radionuclides, computeInitialAuxiliaryVariable(c_inlet, list(k.values())))])

###Spatial and temporal discretization###
#Distance [m]
x = np.linspace(0, 20, num=2001)
#Time [year]
t = 1e5

###Analytical solution###
c = {}
a = {}

c["[Cm-247]"] = computeAnalyticalSolution(x, t, a_inlet["[Cm-247]"], k["[Cm-247]"], v, D)

a["[Am-243]"] = computeAnalyticalSolution(x, t, a_inlet["[Am-243]"], k["[Am-243]"], v, D)
c["[Am-243]"] = a["[Am-243]"] - k["[Cm-247]"] / (k["[Cm-247]"] - k["[Am-243]"]) * c["[Cm-247]"]

a["[Pu-239]"] = computeAnalyticalSolution(x, t, a_inlet["[Pu-239]"], k["[Pu-239]"], v, D)
c["[Pu-239]"] = a["[Pu-239]"] - k["[Cm-247]"] / (k["[Cm-247]"] - k["[Pu-239]"]) * \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Pu-239]"]) * c["[Cm-247]"] - \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Pu-239]"]) * c["[Am-243]"]

a["[U-235]"] = computeAnalyticalSolution(x, t, a_inlet["[U-235]"], k["[U-235]"], v, D)
c["[U-235]"] = a["[U-235]"] - k["[Cm-247]"] / (k["[Cm-247]"] - k["[U-235]"]) * \
        k["[Am-243]"] / (k["[Am-243]"] - k["[U-235]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[U-235]"]) * c["[Cm-247]"] - \
        k["[Am-243]"] / (k["[Am-243]"] - k["[U-235]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[U-235]"]) * c["[Am-243]"] - \
        k["[Pu-239]"] / (k["[Pu-239]"]- k["[U-235]"]) * c["[Pu-239]"]

a["[Pa-231]"] = computeAnalyticalSolution(x, t, a_inlet["[Pa-231]"], k["[Pa-231]"], v, D)
c["[Pa-231]"] = a["[Pa-231]"] - k["[Cm-247]"] / (k["[Cm-247]"] - k["[Pa-231]"]) * \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Pa-231]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Pa-231]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Pa-231]"]) * c["[Cm-247]"] - \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Pa-231]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Pa-231]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Pa-231]"]) * c["[Am-243]"] - \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Pa-231]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Pa-231]"]) * c["[Pu-239]"] - \
        k["[U-235]"] / (k["[U-235]"] - k["[Pa-231]"]) * c["[U-235]"]

a["[Ac-227]"] = computeAnalyticalSolution(x, t, a_inlet["[Ac-227]"], k["[Ac-227]"], v, D)
c["[Ac-227]"] = a["[Ac-227]"] - k["[Cm-247]"] / (k["[Cm-247]"] - k["[Ac-227]"]) * \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Ac-227]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Ac-227]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Ac-227]"]) * \
        k["[Pa-231]"] / (k["[Pa-231]"] - k["[Ac-227]"]) * c["[Cm-247]"] - \
        k["[Am-243]"] / (k["[Am-243]"] - k["[Ac-227]"]) * \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Ac-227]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Ac-227]"]) * \
        k["[Pa-231]"] / (k["[Pa-231]"] - k["[Ac-227]"]) * c["[Am-243]"] - \
        k["[Pu-239]"] / (k["[Pu-239]"] - k["[Ac-227]"]) * \
        k["[U-235]"] / (k["[U-235]"] - k["[Ac-227]"]) * \
        k["[Pa-231]"] / (k["[Pa-231]"] - k["[Ac-227]"]) * c["[Pu-239]"] - \
        k["[U-235]"] / (k["[U-235]"] - k["[Ac-227]"]) * \
        k["[Pa-231]"] / (k["[Pa-231]"] - k["[Ac-227]"]) * c["[U-235]"] - \
        k["[Pa-231]"] / (k["[Pa-231]"] - k["[Ac-227]"]) * c["[Pa-231]"]    

###Plot figure###
fig, ax = plt.subplots()
    
ax.set(xlim=(0,20), ylim=(0,1.4))
ax.set(xlabel="Distance [m]", ylabel="Concentration [mol/L]")
    
color_map=cm.rainbow(np.linspace(0,1,radionuclides.size))

for radionuclide,color in zip(radionuclides, color_map):
    ax.plot(x, c[radionuclide], linestyle='-', lw=1.5, label=radionuclide, c=color, zorder=10, clip_on=False)

ax.legend(frameon=False, loc='upper right', numpoints=1, fontsize=12, ncol=1)
    
ax.xaxis.grid(color='gray', linestyle='dashed')
ax.yaxis.grid(color='gray', linestyle='dashed')

png

Numerical solution

Two numerical models are presented here, one based on the global implicit approach (GIA) and the other on the operator splitting approach (OS). The OGS input files of these two models can be found here. Then, the obtained numerical solutions both are plotted against the analytical solution for verification.

###Run two numerical models###
prj_name = "1d_decay_chain"
prj_file_GIA = f"{notebook_dir}/GlobalImplicitApproach/{prj_name}_GIA.prj"
prj_file_OS = f"{notebook_dir}/{prj_name}_OS.prj"
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')

if not os.path.exists(out_dir):
    os.makedirs(out_dir)
os.chdir(out_dir)

print(f"ogs {prj_file_GIA} > out.txt")
start_time = time.time()
! ogs {prj_file_GIA} > out.txt
end_time =time.time()
runtime_GIA = round(end_time - start_time, 2)
print("Execution time for the GIA model is ", runtime_GIA, "s")

###Read simulation results###
pvdfile_GIA = vtuIO.PVDIO(f"{notebook_dir}/GlobalImplicitApproach/{prj_name}_GIA.pvd", dim=1)
#Given the fact that the runtime of the OS model (about 1800s) is 
#far longer than the time constraint specified (600s), we decide not
#to use the OS simulation results obtained from automated testing.
#Instead, the pre-prepared reference simulation results are used.
pvdfile_OS = vtuIO.PVDIO(f"{notebook_dir}/{prj_name}_OS.pvd", dim=1)

fig, (ax1, ax2) = plt.subplots(2,1, figsize=(6, 9))

color_map=cm.rainbow(np.linspace(0,1,radionuclides.size))

###Plot subfigure 1###
ax1.set(xlim=(0,20), ylim=(0,1.4))
ax1.set(xlabel="Distance [m]", ylabel="Concentration [mol/L]")

#Analytical solution     
for radionuclide,color in zip(radionuclides, color_map):
    ax1.plot(x, c[radionuclide], linestyle='-', lw=1.5, label="Exact - "+radionuclide, \
             c=color, zorder=10, clip_on=False)      
    
#GIA numerical solution    
for radionuclide,color in zip(radionuclides, color_map):
    c_gia = pvdfile_GIA.read_set_data(t*3.1536e7, radionuclide, data_type="point", pointsetarray=[(i,0,0) for i in x])
    ax1.plot(x, c_gia, label="GIA - " + radionuclide, color=color, marker='o', \
             markevery=50, linestyle="", zorder=10, clip_on=False)
    
#OS numerical solution    
for radionuclide,color in zip(radionuclides, color_map):
    c_os = pvdfile_OS.read_set_data(t*3.1536e7, radionuclide, data_type="point", pointsetarray=[(i,0,0) for i in x])
    ax1.plot(x, c_os, label="OS - " + radionuclide, color=color, marker='^', \
             markevery=50, linestyle="", zorder=10, clip_on=False)

#numerical solution by reference code
#added once the bc value is double-checked
porosity = 0.12
with h5py.File(f"{notebook_dir}/solution_reference_code.hdf5","r") as f:
    species_ = f["species"][:]
    x_ = f["x"][:]
    for s_, radionuclide, color in zip(species_, radionuclides, color_map):
        data_ = f[s_][:]
        ax1.plot(x_, data_[:,1]/porosity, label="Reference code - " + radionuclide, \
                 color=color, marker='D', markevery=5, linestyle="", zorder=10, clip_on=False)
        
ax1.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=4)
    
ax1.xaxis.grid(color='gray', linestyle='dashed')
ax1.yaxis.grid(color='gray', linestyle='dashed')

###Plot subfigure 2###
ax2.set(xlim=(0,0.6), ylim=(0,1.2))
ax2.set(xlabel="Distance [m]", ylabel="[Ac-227] Concentration [mol/L]")

#Analytical solution
ax2.plot(x[np.where(x<0.6)], c["[Ac-227]"][np.where(x<0.6)], linestyle='-', lw=1.5, label="Exact solution", \
         c=color_map[-1], zorder=10, clip_on=False)      

#GIA numerical solution
Ac_227_gia = pvdfile_GIA.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", pointsetarray=[(i,0,0) for i in x])
ax2.plot(x[np.where(x<0.6)], Ac_227_gia[np.where(x<0.6)], label="GIA solution", color=color_map[2], \
         linestyle="--", zorder=10, clip_on=False)

#OS numerical solution
Ac_227_os = pvdfile_OS.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", pointsetarray=[(i,0,0) for i in x])
ax2.plot(x[np.where(x<0.6)], Ac_227_os[np.where(x<0.6)], label="OS solution", color=color_map[-2], \
         linestyle="-.", zorder=10, clip_on=False)

#numerical solution by reference code
with h5py.File(f"{notebook_dir}/solution_reference_code.hdf5","r") as f:
    Ac_227_ = f["species"][-1]
    x_ = f["x"][:]
    Ac_227_ = f[Ac_227_][:]
    ax2.plot(f["x"][np.where(x_<0.7)], Ac_227_[np.where(x_<0.7),1][0]/porosity, label="Reference code", \
                 color=color_map[0], linestyle="--", zorder=10, clip_on=False)
        
ax2.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=1)
    
ax2.xaxis.grid(color='gray', linestyle='dashed')
ax2.yaxis.grid(color='gray', linestyle='dashed')
ogs /home/jovyan/work/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/DecayChain/GlobalImplicitApproach/1d_decay_chain_GIA.prj > out.txt
Execution time for the GIA model is  15.15 s

png

As can be seen from Subfigure 1, for all radionuclides except the fastest decaying Actinium-227, the predicted concentration profiles using GIA and OS both agree well with the analytical solutions. For Actinium-227 concentration, the GIA solution perfectly match with the reference result (see Subfigure 2). By contrast, the OS solution near the inlet deviates significantly from the reference result. The observed deviation is due to the fact that the time step size adopted in the simulation (\(\Delta t\) = 100 years) is too large compared to the half-life of Actinium-227 (\(t_{1/2}\) = 21.773 years).

Error analysis

Furthermore, we perform an error analysis of the GIA and OS solutions. As can be observed from Subfigures 1 and 3, the absolute error of the GIA solution for all the radionuclides is suppressed at a low level. The maximum absolute error is only 3.5\(\times\)10\(^{-4}\) mol/L. In contrast, the absolute error of the OS solution is two to three orders of magnitude greater than that of the GIA solution (see Subfigures 2 and 3). And the error becomes more apparent as you approach the inlet.

fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(6, 13.5))

color_map=cm.rainbow(np.linspace(0,1,radionuclides.size))

###Plot subfigure 1###
ax1.set(xlim=(0,20), ylim=(0, 0.0005))
ax1.set(xlabel="Distance [m]", ylabel="Concentration [mol/L]")
    
#GIA numerical solution    
for radionuclide,color in zip(radionuclides[:-1], color_map):
    c_gia = pvdfile_GIA.read_set_data(t*3.1536e7, radionuclide, data_type="point", pointsetarray=[(i,0,0) for i in x])
    ax1.plot(x, abs(c_gia - c[radionuclide]), label="GIA - " + radionuclide, color=color, linestyle='-', \
             lw=1.5, zorder=10, clip_on=False)
    
ax1.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=1)
    
ax1.xaxis.grid(color='gray', linestyle='dashed')
ax1.yaxis.grid(color='gray', linestyle='dashed')

###Plot subfigure 2###
ax2.set(xlim=(0,20), ylim=(0, 0.01))
ax2.set(xlabel="Distance [m]", ylabel="Concentration [mol/L]")

#OS numerical solution    
for radionuclide,color in zip(radionuclides[:-1], color_map):
    c_os = pvdfile_OS.read_set_data(t*3.1536e7, radionuclide, data_type="point", pointsetarray=[(i,0,0) for i in x])
    ax2.plot(x, abs(c_os - c[radionuclide]), label="OS - " + radionuclide, color=color, linestyle='-.', \
             lw=1.5, zorder=10, clip_on=False)
    
ax2.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=1)
    
ax2.xaxis.grid(color='gray', linestyle='dashed')
ax2.yaxis.grid(color='gray', linestyle='dashed')

###Plot subfigure 3###
ax3.set(xlim=(0,0.6), ylim=(0,1.0))
ax3.set(xlabel="Distance [m]", ylabel="[Ac-227] Concentration [mol/L]")

#GIA numerical solution
Ac_227_gia = pvdfile_GIA.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", pointsetarray=[(i,0,0) for i in x])
ax3.plot(x[np.where(x<0.6)], abs(Ac_227_gia[np.where(x<0.6)] - c["[Ac-227]"][np.where(x<0.6)]), \
         label="GIA - [Ac-227]", color=color_map[-1], linestyle='-', lw=1.5, zorder=10, clip_on=False)

#OS numerical solution
Ac_227_os = pvdfile_OS.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", pointsetarray=[(i,0,0) for i in x])
ax3.plot(x[np.where(x<0.6)], abs(Ac_227_os[np.where(x<0.6)] - c["[Ac-227]"][np.where(x<0.6)]), \
         label="OS - [Ac-227]", color=color_map[-1], linestyle='-.', lw=1.5, zorder=10, clip_on=False)

ax3.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=1)
    
ax3.xaxis.grid(color='gray', linestyle='dashed')
ax3.yaxis.grid(color='gray', linestyle='dashed')

png

The effect of time step size

As shown in the figure below, the OS solution is rather sensitive to the time stepping scheme. By reducing the time step size from 100 years to 5 years, the solution accuracy near the inlet is greatly improved. Nevertheless, the OS solution accuracy is yet lower than that of the GIA solution despite a larger time step size.

###Read OS solution with a smaller time step size###
pvdfile_OS_small_ts = vtuIO.PVDIO(f"{notebook_dir}/{prj_name}_OS_5year.pvd", dim=1)

###Plot figure###
fig, ax = plt.subplots()

ax.set(xlim=(0, 0.6), ylim=(0, 1.2))
ax.set(xlabel='Distance [m]', ylabel='[Ac-227] Concentration [mol/L]')

color_map=cm.rainbow(np.linspace(0,1,radionuclides.size))

#Analytical solution
ax.plot(x[np.where(x<0.6)], c["[Ac-227]"][np.where(x<0.6)], linestyle='-', lw=1.5, \
        label="Exact solution", color=color_map[-1], zorder=10, clip_on=False)      

#OS solution with a time step size of 100 years
Ac_227_os = pvdfile_OS.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", pointsetarray=[(i,0,0) for i in x])
ax.plot(x[np.where(x<0.6)], Ac_227_os[np.where(x<0.6)], label="Time step size $\Delta$t = 100 years", \
        color=color_map[1], linestyle="--", zorder=10, clip_on=False)

#OS solution with a time step size of 5 years
Ac_227_os_small_ts = pvdfile_OS_small_ts.read_set_data(t*3.1536e7, "[Ac-227]", data_type="point", \
                                                       pointsetarray=[(i,0,0) for i in x])
ax.plot(x[np.where(x<0.6)], Ac_227_os_small_ts[np.where(x<0.6)], label="Time step size $\Delta$t = 5 years", \
        color=color_map[-2], linestyle="-.", zorder=10, clip_on=False)

ax.legend(bbox_to_anchor=(1.04,1), loc="upper left", numpoints=1, fontsize=12, ncol=1)
    
ax.xaxis.grid(color='gray', linestyle='dashed')
ax.yaxis.grid(color='gray', linestyle='dashed')

png

Computational performance of two solution strategies

The runtime of the GIA model is only 22.13 seconds, approximately two orders of magnitude faster than the OS model in this case (see the figure below). Further time profiling analysis for the OS model indicates that Phreeqc consumes over 90 percent of the runtime for solving the local decay-chain problem.

prj_file_GIA_4 = f"{notebook_dir}/GlobalImplicitApproach/MPI/4Processors/{prj_name}_GIA.prj"
prj_file_GIA_8 = f"{notebook_dir}/GlobalImplicitApproach/MPI/8Processors/{prj_name}_GIA.prj"
out_dir = os.environ.get('OGS_TESTRUNNER_OUT_DIR', '_out')

if not os.path.exists(out_dir):
    os.makedirs(out_dir)
os.chdir(out_dir)

print(f"mpirun -np 4 ogs {prj_file_GIA_4} > out.txt")
start_time = time.time()
! mpirun -np 4 ogs {prj_file_GIA_4} > out.txt
end_time =time.time()
runtime_GIA_4 = round(end_time - start_time, 2)
print("Execution time for the parallelized GIA model with 4 processors is ", runtime_GIA_4, "s")

print(f"mpirun -np 8 ogs {prj_file_GIA_8} > out.txt")
start_time = time.time()
! mpirun -np 8 ogs {prj_file_GIA_8} > out.txt
end_time =time.time()
runtime_GIA_8 = round(end_time - start_time, 2)
print("Execution time for the parallelized GIA model with 8 processors is ", runtime_GIA_8, "s")

#runtime [second]
runtime = {
    "OS - 1 Processor": 1980,
    "GIA - 1 Processor": runtime_GIA,
    "GIA - 4 Processor": runtime_GIA_4,
    "GIA - 8 Processor": runtime_GIA_8}

###Plot figure###
fig, ax = plt.subplots()

ax.set(xlim=(-1, 2), ylim=(0, 2500))

ax.set_ylabel('Runtime [second]')
ax.set_yticks(np.arange(0, 3000, 500))

ax.yaxis.grid(color='gray', linestyle='dashed', zorder=0)

ax.bar(list(runtime.keys())[:2], list(runtime.values())[:2], width=0.5, zorder=3)

for i in range(0, 2):
    ax.annotate(list(runtime.values())[i],(i,list(runtime.values())[i]+50), ha="center")
mpirun -np 4 ogs /home/jovyan/work/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/DecayChain/GlobalImplicitApproach/MPI/4Processors/1d_decay_chain_GIA.prj > out.txt
Execution time for the parallelized GIA model with 4 processors is  5.39 s
mpirun -np 8 ogs /home/jovyan/work/Tests/Data/Parabolic/ComponentTransport/ReactiveTransport/DecayChain/GlobalImplicitApproach/MPI/8Processors/1d_decay_chain_GIA.prj > out.txt
Execution time for the parallelized GIA model with 8 processors is  17.95 s

png

To summarize, the global implicit approach is clearly superior to the operator splitting approach in terms of solution accuracy and computational performance. Moreover, the global implicit approach even allows for a larger time step size.

The effect of parallelization on computational performance

In order to further improve the computational performance of the GIA model, MPI parallelization is employed. As its result, the simulation runtime is drastically reduced down to 4.82 seconds with the number of processors increased to 8. This results in a 4.6X speed boost.

###Plot figure###
fig, ax = plt.subplots()

ax.set(ylim=(0, 30))

ax.set_ylabel('Runtime [second]')
ax.yaxis.grid(color='gray', linestyle='dashed', zorder=0)

ax.bar(list(runtime.keys())[1:], list(runtime.values())[1:], width=0.5, zorder=3)

for i in range(1, 4):
    ax.annotate(list(runtime.values())[i],(i-1,list(runtime.values())[i]+2), ha="center")

png

References

(GRS 289) VSG: Radiologische Konsequenzenanalyse (AP 10). J. Larue, B. Baltes, H. Fischer, G. Frieling, I. Kock, M. Navarro, H. Seher.

Sun, Y., Petersen, J. N., & Clement, T. P. (1999). Analytical solutions for multiple species reactive transport in multiple dimensions. Journal of Contaminant Hydrology, 35(4), 429-440.

Credits:

Renchao Lu, Christoph Behrens, Dmitri Naumov, Haibing Shao


This article was written by Renchao Lu, Christoph Behrens, Dmitri Naumov, Haibing Shao. If you are missing something or you find an error please let us know.
Generated with Hugo 0.101.0 in CI job 237559 | Last revision: August 5, 2022