Coupling Reaktoro into other reactive transport codes¶
In this tutorial, we show how Reaktoro can be used in other codes for reactive transport modeling. Here, you will find that we have split the mass transport and chemical reaction calculations so that you can clearly see how a dedicated and advanced transport solver could be combined with Reaktoro’s solvers for chemical reaction calculations. A basic transport solver is used here instead for the sake of software coupling demonstration.
The reactive transport problem we consider here is also relatively simple, similar to the one presented in the tutorial Reactive transport of CO2-saturated brine along a porous rock column. Unlike other tutorials, we proceed here first with a step-by-step explanation of the full script found at the end of this tutorial.
Importing python packages¶
First, we need to import a few Python packages to enable us to perform the numerical calculations and plotting.
from reaktoro import *
from numpy import *
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
import os, time
We import the reaktoro Python package so that we can use its classes and methods for performing chemical reaction calculations, numpy for working with arrays, matplotlib for plotting capabilities, joblib for simple parallel computing, and, finally, os, to provide a portable way of using operating system dependent functionality.
Defining parameters for the reactive transport simulation¶
Next, we define reactive transport and numerical discretization parameters.
xl = 0.0 # the x-coordinate of the left boundary
xr = 1.0 # the x-coordinate of the right boundary
nsteps = 100 # the number of steps in the reactive transport simulation
ncells = 100 # the number of cells in the discretization
D = 1.0e-9 # the diffusion coefficient (in units of m2/s)
v = 1.0/day # the fluid pore velocity (in units of m/s)
dt = 10*minute # the time step (in units of s)
T = 60.0 + 273.15 # the temperature (in units of K)
P = 100 * 1e5 # the pressure (in units of Pa)
dirichlet = False # the parameter that determines whether Dirichlet BC must be used
We specify the considered rock domain by setting coordinates of its left and right boundaries to 0.0 m and 100.0 m, respectively. The discretization parameters, i.e., the number of cells and steps in time, are both set to 100. The reactive transport modeling procedure assumes a constant fluid velocity of 1 m/day (1.16 · 10-5 m/s) and the same diffusion coefficient of 10-9 m2/s for all fluid species (without dispersivity). The size of the time-step is set to 10 minutes. Temperature and pressure are set to 60 °C and 100 bar, respectively, throughout the whole tutorial. The chemical equilibrium calculations performed in each mesh cell, at every time step, using Gibbs energy minimisation algorithm (provided by the class EquilibriumSolver).
The boolean variable dirichlet is set to True or False depending on
which boundary condition is considered in the numerical calculation. Set to
False for imposing the flux of the injected fluid, otherwise, set to
True for imposing the composition of the fluid on the left boundary.
Specifying the quantities and properties to be outputted¶
Before running the reactive transport simulations, we specify the list of parameters we are interested in outputting. In this case, it is pH, molality of H+, Ca2+, Mg2+, HCO3-, CO2 (aq), as well as a phase volume of calcite and dolomite.
output_quantities = """
pH
speciesMolality(H+)
speciesMolality(Ca++)
speciesMolality(Mg++)
speciesMolality(HCO3-)
speciesMolality(CO2(aq))
phaseVolume(Calcite)
phaseVolume(Dolomite)
""".split()
Organizing the main function in three parts¶
Below is the main function in our script, which shows three parts, each represented by a Python function documented in the next sections:
if __name__ == '__main__':
# Create folders for the results
make_results_folders()
# Run the reactive transport simulations
simulate()
# Plotting the result
plot()
Here, we first create the required folders for the results, we then run the reactive transport simulation, and, finally, plot the outputted results.
Creating folders for the outputted results¶
def make_results_folders():
os.system('mkdir -p results')
os.system('mkdir -p figures/ph')
os.system('mkdir -p figures/aqueous-species')
os.system('mkdir -p figures/calcite-dolomite')
os.system('mkdir -p videos')
Using os package, we create required folders for outputting the obtained results and for the plot and video files later.
Perform the reactive transport simulation¶
The reactive transport simulation is performed in the function simulate
(you may want to have a quick look in the code below and move to a
detailed explanation in the subsequent sections):
def simulate():
# Step 7.1: Construct the chemical system with its phases and species
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
# Step 7.2: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
# Step 7.3: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
# Step 7.4: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
# Step 7.5: Calculate the equilibrium states for the initial and boundary conditions
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
# Step 7.6: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
# Step 7.7: Partitioning fluid and solid species
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
# Step 7.8: Create a list of chemical states for the mesh cells
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
# Step 7.9: Create the equilibrium solver object for the repeated equilibrium calculation
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
# Step 7.10: Create the equilibrium solver object for the repeated equilibrium calculation
solver = EquilibriumSolver(system)
# Step 7.11: The auxiliary function to create an output file each time step
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
# Step 7.12: Running the reactive transport simulation loop
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
We now explain step-by-step the function simulate. We begin by defining the
chemical system using the ChemicalEditor class:
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
We specify aqueous and mineral phases that should be considered in the chemical system. For performance reasons, the aqueous phase is defined by manually specifying the chemical species, instead of automatic collection of species from database. There are three pure mineral phase considered: quartz (SiO2), calcite (CaCO3), and dolomite (CaMg(CO3)2).
Creating the chemical system¶
This step is where we create an object of class ChemicalSystem using the
chemical system definition details stored in the object editor.
system = ChemicalSystem(editor)
Initial condition (IC) of the reactive transport problem¶
After constructing the chemical system of interest, we can proceed to the definition of a chemical equilibrium problem to set the initial condition for the fluid and rock composition that is consistent with our intention of modeling reactive transport in a porous rock column composed of quartz and calcite, at 60 °C and 100 bar, with a 0.7 NaCl molal brine in equilibrium with the rock minerals. To ensure equilibrium with both quartz and calcite, the equilibrium problem setup below considers a relatively large amount for these minerals (10 moles each).
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
We will later scale the volumes of the aqueous and mineral phases so that they are consistent with a 10 % porosity and the required volume percentages of the rock minerals (98 %vol of quartz and 2 %vol of calcite).
Boundary condition (BC) of the reactive transport problem¶
For the boundary condition, we need to specify the composition of the fluid that is injected into the rock. This is done below, by defining an equilibrium problem that will later be solved to produce an aqueous fluid with 0.9 molal NaCl, 0.05 molal MgCl2, 0.01 CaCl2, and 0.75 molal CO2, in a state very close to CO2 saturation.
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
The composition of the injected aqueous fluid results from the solution of the above equilibrium problem, in which 1 kg of H2O is mixed with 0.90 moles of NaCl, 0.05 moles of MgCl2, 0.01 moles of CaCl2, and 0.75 moles of CO2.
Calculate the equilibrium states for the initial and boundary conditions¶
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
In this step, we use method equilibrate to calculate the chemical
equilibrium state of the system with the given initial and boundary equilibrium
conditions stored in the objects problem_ic and problem_bc. The
numerical solution of each problem results in the objects state_ic
and state_bc, respectively, of class ChemicalState, which stores the
temperature, pressure, and the amounts of every species in the system.
Scaling the volumes of the phases¶
We scale the volumes of the phases in the initial condition according to the condition that porosity is 10 % (ratio of fluid volume to total volume), and the volume fractions of minerals are 98 %vol of quartz (SiO2) and 2 %vol of calcite (ratio of mineral volume to solid volume).
For the chemical state representing the boundary condition for the injected fluid composition, we scale its volume to 1 m3. This is done so that the amounts of the species in the fluid is consistent with a mol/m3 scale.
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
Partitioning fluid and solid species¶
Only species in fluid phases are mobile and transported by advection and diffusion mechanisms. The solid phases are immobile.
The code below identify the indices of the fluid and solid species, and also create arrays to keep track of the amounts of elements in the fluid and solid partition (i.e., the amounts of elements among all fluid phases, here only an aqueous phase, and the amounts of elements among all solid phases, here the mineral phases).
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
We use methods indicesFluidSpecies and indicesSolidSpecies of class
ChemicalSystem to get the indices of the fluid and solid species, which are
then stored in the lists ifluid_species and isolid_species,
respectively. Then, we define the arrays b, bfluid, bsolid, that
will store, respectively, the concentrations (mol/m3) of each element in the
system, in the fluid partition, and in the solid partition at every time step.
The array b is initialized with the concentrations of the elements at the
initial chemical state, state_ic, using method elementAmounts of class
ChemicalState. The array b_bc stores the concentrations of each element
on the boundary.
Create a list of chemical states for the mesh cells¶
Every mesh cell needs a ChemicalState object, which is done below with the
creation of a list states of length ncells, initialized with state_ic
for each cell (i.e., all cells have initially the same fluid and rock composition).
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
Cell coordinates and lengths¶
Next, we generate the coordinates of the mesh nodes (array x) by equally
dividing the interval [xr, xl] into ncells. The length between each
consecutive mesh nodes is computed and stored in dx (the length of the mesh
cells).
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
Creating the equilibrium solver¶
For the repeated equilibrium calculation, we define an equilibrium solver object using EquilibriumSolver class.
solver = EquilibriumSolver(system)
Define an auxiliary function for creating the output file¶
The auxiliary function outputstate uses the class ChemicalOutput to
output quantities and properties at each mesh cell, at every time step.
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
Running the reactive transport simulation¶
Before proceeding to the reactive transport calculations, we set the initial time and a counter for the number of time steps.
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
The reactive transport simulation continues until the maximum number of time
steps is achieved. At each time step, we output the progress of the simulation
using function outputstate. Then, we collect the amounts of elements from
fluid and solid partition. We now apply an operator splitting procedure, in
which we first update the amounts of elements in the fluid partition
(bfluid) using the transport equations (without reactions). These updated
amounts of elements in the fluid partition are now summed with the amounts of
elements in the solid partition (bsolid, which remained constant during the
transport step), and thus updating the amounts of elements in the chemical
system (b). We now use these updated amounts of elements in the cell
(b) to evaluate its new chemical equilibrium state, thus producing new
amounts of the species in both the fluid and solid phases (available in the
list states of ChemicalState objects). This chemical reaction
equilibrium calculation step, at each mesh cell, permits, for example,
aqueous species and minerals to react, and thus causes mineral dissolution or
precipitation, depending on how much the amount of mineral species changes. This
can then be used, for example, to compute new porosity value for the cell.
Solving the transport equations¶
Reactive transport calculations involve the solution of a system of advection-diffusion-reaction equations. Due to the use of operator splitting scheme in this tutorial, we solve the transport equations and then perform the chemical reaction calculations for each cell before proceeding to the next time step.
The transport function below is responsible for solving an
advection-diffusion equation, that is later applied to transport the
concentrations (mol/m3) of elements in the fluid partition (a simplification
that is possible because of common diffusion coefficients and velocities of the
fluid species, otherwise the transport of individual fluid species would be
needed).
def transport(u, dt, dx, v, D, g):
# Number of DOFs
n = len(u)
alpha = D*dt/dx**2
beta = v*dt/dx
a = full(n, -beta - alpha)
b = full(n, 1 + beta + 2*alpha)
c = full(n, -alpha)
# Set the boundary condition
if dirichlet:
# Use Dirichlet BC boundary conditions
b[0] = 1.0
c[0] = 0.0
u[0] = g
else:
# Use flux boundary conditions
u[0] += beta*g
b[0] = 1 + beta + alpha
b[-1] = 1 + beta + alpha
# Solve a tridiagonal matrix equation
thomas(a, b, c, u)
The function transport expects a conservative property (argument u)
(e.g., the concentration (mol/m3) of j-th element in the fluid given by
bfluid[j]), the time step (dt), the mesh cell length (dx), the
fluid velocity (v), the diffusion coefficient (D), and the boundary
condition of the conservative property (g) (e.g., the concentration of the
j-th element in the fluid on the left boundary).
The transport equations are solved with a finite volume method. Its
discretization in space and time (implicit) results in the constants alpha
and beta. These correspond to the diffusion and advection terms in the
equation: D*dt/dx**2 and v*dt/dx, respectively.
Arrays a, b, c are the diagonals in the tridiagonal matrix that
results by writing all discretized equations in a matrix equation. This system
of linear equations is solved by the tridiagonal matrix algorithm, also known
as the Thomas algorithm:
def thomas(a, b, c, d):
n = len(d)
c[0] /= b[0]
for i in range(1, n - 1):
c[i] /= b[i] - a[i]*c[i - 1]
d[0] /= b[0]
for i in range(1, n):
d[i] = (d[i] - a[i]*d[i - 1])/(b[i] - a[i]*c[i - 1])
x = d
for i in reversed(range(0, n - 1)):
x[i] -= c[i]*x[i + 1]
return x
Plotting of the results¶
The last block of the main routine is plotting of the results.
def plot():
# Plot all result files
files = sorted(os.listdir('results'))
Parallel(n_jobs=16)(delayed(plotfile)(file) for file in files)
# Create videos for the figures
ffmpegstr = 'ffmpeg -y -r 30 -i figures/{0}/%03d.png -codec:v mpeg4 -flags:v +qscale -global_quality:v 0 videos/{0}.mp4'
os.system(ffmpegstr.format('calcite-dolomite'))
os.system(ffmpegstr.format('aqueous-species'))
os.system(ffmpegstr.format('ph'))
The full script for the reactive transport calculation is shown below.
from reaktoro import *
from numpy import *
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
import os, time
# Step 2: Auxiliary time related constants
second = 1
minute = 60
hour = 60 * minute
day = 24 * hour
year = 365 * day
# Step 3: Parameters for the reactive transport simulation
xl = 0.0 # the x-coordinate of the left boundary
xr = 1.0 # the x-coordinate of the right boundary
nsteps = 100 # the number of steps in the reactive transport simulation
ncells = 100 # the number of cells in the discretization
D = 1.0e-9 # the diffusion coefficient (in units of m2/s)
v = 1.0/day # the fluid pore velocity (in units of m/s)
dt = 10*minute # the time step (in units of s)
T = 60.0 + 273.15 # the temperature (in units of K)
P = 100 * 1e5 # the pressure (in units of Pa)
dirichlet = False # the parameter that determines whether Dirichlet BC must be used
# Step 4: The list of quantities to be output for each mesh cell, each time step
output_quantities = """
pH
speciesMolality(H+)
speciesMolality(Ca++)
speciesMolality(Mg++)
speciesMolality(HCO3-)
speciesMolality(CO2(aq))
phaseVolume(Calcite)
phaseVolume(Dolomite)
""".split()
# Step 7: Perform the reactive transport simulation
def simulate():
# Step 7.1: Construct the chemical system with its phases and species
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
# Step 7.2: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
# Step 7.3: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
# Step 7.4: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
# Step 7.5: Calculate the equilibrium states for the initial and boundary conditions
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
# Step 7.6: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
# Step 7.7: Partitioning fluid and solid species
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
# Step 7.8: Create a list of chemical states for the mesh cells
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
# Step 7.9: Create the equilibrium solver object for the repeated equilibrium calculation
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
# Step 7.10: Create the equilibrium solver object for the repeated equilibrium calculation
solver = EquilibriumSolver(system)
# Step 7.11: The auxiliary function to create an output file each time step
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
# Step 7.12: Running the reactive transport simulation loop
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
# Step 10: Return a string for the title of a figure in the format Time: #h##m
def titlestr(t):
t = t / minute # Convert from seconds to minutes
h = int(t) / 60 # The number of hours
m = int(t) % 60 # The number of remaining minutes
return 'Time: {:>3}h{:>2}m'.format(h, str(m).zfill(2))
# Step 9: Generate figures for a result file
def plotfile(file):
step = int(file.split('.')[0])
print('Plotting figure', step, '...')
t = step * dt
filearray = loadtxt('results/' + file, skiprows=1)
data = filearray.T
ndigits = len(str(nsteps))
plt.figure()
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=2.5, top=10.5)
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('pH')
plt.plot(data[0], data[1])
plt.tight_layout()
plt.savefig('figures/ph/{}.png'.format(str(step).zfill(ndigits)))
plt.figure()
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=-0.1, top=2.1)
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('Mineral Volume [%$_{\mathsf{vol}}$]')
plt.plot(data[0], data[7] * 100, label='Calcite')
plt.plot(data[0], data[8] * 100, label='Dolomite')
plt.legend(loc='center right')
plt.tight_layout()
plt.savefig('figures/calcite-dolomite/{}.png'.format(str(step).zfill(ndigits)))
plt.figure()
plt.yscale('log')
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=0.5e-5, top=2)
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('Concentration [molal]')
plt.plot(data[0], data[3], label='Ca++')
plt.plot(data[0], data[4], label='Mg++')
plt.plot(data[0], data[5], label='HCO3-')
plt.plot(data[0], data[6], label='CO2(aq)')
plt.plot(data[0], data[2], label='H+')
plt.legend(loc='lower right')
plt.tight_layout()
plt.savefig('figures/aqueous-species/{}.png'.format(str(step).zfill(ndigits)))
plt.close('all')
# Step 8: Plot all result files and generate a video
def plot():
# Plot all result files
files = sorted(os.listdir('results'))
Parallel(n_jobs=16)(delayed(plotfile)(file) for file in files)
# Create videos for the figures
ffmpegstr = 'ffmpeg -y -r 30 -i figures/{0}/%03d.png -codec:v mpeg4 -flags:v +qscale -global_quality:v 0 videos/{0}.mp4'
os.system(ffmpegstr.format('calcite-dolomite'))
os.system(ffmpegstr.format('aqueous-species'))
os.system(ffmpegstr.format('ph'))
# Step 7.10.2: Solve a tridiagonal matrix equation using Thomas algorithm (or TriDiagonal Matrix Algorithm (TDMA))
def thomas(a, b, c, d):
n = len(d)
c[0] /= b[0]
for i in range(1, n - 1):
c[i] /= b[i] - a[i]*c[i - 1]
d[0] /= b[0]
for i in range(1, n):
d[i] = (d[i] - a[i]*d[i - 1])/(b[i] - a[i]*c[i - 1])
x = d
for i in reversed(range(0, n - 1)):
x[i] -= c[i]*x[i + 1]
return x
# Step 7.10.1: Perform a transport step
def transport(u, dt, dx, v, D, g):
# Number of DOFs
n = len(u)
alpha = D*dt/dx**2
beta = v*dt/dx
a = full(n, -beta - alpha)
b = full(n, 1 + beta + 2*alpha)
c = full(n, -alpha)
# Set the boundary condition
if dirichlet:
# Use Dirichlet BC boundary conditions
b[0] = 1.0
c[0] = 0.0
u[0] = g
else:
# Use flux boundary conditions
u[0] += beta*g
b[0] = 1 + beta + alpha
b[-1] = 1 + beta + alpha
# Solve a tridiagonal matrix equation
thomas(a, b, c, u)
# Step 6: Creating folders for the results' output
def make_results_folders():
os.system('mkdir -p results')
os.system('mkdir -p figures/ph')
os.system('mkdir -p figures/aqueous-species')
os.system('mkdir -p figures/calcite-dolomite')
os.system('mkdir -p videos')
# Step 5: Define the main function
if __name__ == '__main__':
# Create folders for the results
make_results_folders()
# Run the reactive transport simulations
simulate()
# Plotting the result
plot()
Todo
Add figures and videos here in this tutorial.
Have you got an issue?¶
Have you found any issue or error in this tutorial? Do you have any recommendations or you think something is not clear enough? Please, let us know by filling a new issue here:
You’ll need a GitHub account - but this is easy to sort out if you don’t have one yet!