{
"cells": [
{
"cell_type": "markdown",
"id": "b7fbfc84",
"metadata": {},
"source": [
"# Coupling Reaktoro into other reactive transport codes\n",
"\n",
"In this tutorial, we show how Reaktoro can be used in other codes for reactive transport modeling. Here,\n",
"you will find that we have split the mass transport and chemical reaction calculations so that you can see\n",
"how a dedicated and advanced transport solver could be combined with Reaktoro's solvers for chemical reaction\n",
"calculations. A basic transport solver is used here for the sake of the software coupling demonstration.\n",
"\n",
"The reactive transport problem we consider here is also relatively simple. We proceed with a step-by-step\n",
"explanation.\n",
"\n",
"## Importing python packages\n",
"\n",
"First, we need to import a few Python packages to enable us to perform the numerical calculations and plotting."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d5f109cb",
"metadata": {},
"outputs": [],
"source": [
"print('============================================================')\n",
"print('Make sure you have the following Python packages installed: ')\n",
"print(' numpy, natsort, bokeh')\n",
"print('These can be installed with pip:')\n",
"print(' pip install numpy natsort bokeh')\n",
"print('============================================================')\n",
"from reaktoro import *\n",
"import numpy as np\n",
"from natsort import natsorted\n",
"from tqdm.notebook import tqdm\n",
"import os\n",
"\n",
"# Import components of bokeh library\n",
"from bokeh.io import show, output_notebook\n",
"from bokeh.layouts import column\n",
"from bokeh.plotting import figure\n",
"from bokeh.models import Range1d, ColumnDataSource\n",
"from bokeh.models.widgets import Slider\n"
]
},
{
"cell_type": "markdown",
"id": "121f47ee",
"metadata": {},
"source": [
"We import the **reaktoro** Python package so that we can use its classes and methods for performing chemical\n",
"reaction calculations, **numpy** for working with arrays, **tqdm** for the progress bar functionality and **os**,\n",
"to provide a portable way of using operating system dependent functionality. For plotting capabilities of obtained\n",
"results, we use **bokeh** library.\n",
"\n",
"> **Note**: To make sure that all the widgets are working correctly, make sure to run:\n",
"> `$ jupyter nbextension enable --py widgetsnbextension` and\n",
"> `$jupyter labextension install @jupyter-widgets/jupyterlab-manager`.\n",
"\n",
"## Initializing auxiliary time-related constants\n",
"In this step, we initialize auxiliary time-related constants from seconds up to years used in the rest of the code."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2ca7e55c",
"metadata": {},
"outputs": [],
"source": [
"second = 1\n",
"minute = 60\n",
"hour = 60 * minute\n",
"day = 24 * hour\n",
"week = 7 * day\n",
"year = 365 * day"
]
},
{
"cell_type": "markdown",
"id": "96047bc1",
"metadata": {},
"source": [
"## Defining parameters for the reactive transport simulation\n",
"\n",
"Next, we define reactive transport and numerical discretization parameters. In particular, we specify the considered\n",
"rock domain by setting coordinates of its left and right boundaries to 0.0 m and 100.0 m, respectively. The\n",
"discretization parameters, i.e., the number of cells and steps in time, are both set to 100. The reactive\n",
"transport modeling procedure assumes a constant fluid velocity of 1 m/day (1.16 · 10-5 m/s) and the same\n",
"diffusion coefficient of 10-9 m2/s for all fluid species (without dispersivity). The\n",
"size of the time-step is set to 30 minutes. Temperature and pressure are set to 60 °C and 100 bar, respectively,\n",
"throughout the whole tutorial. "
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4a8e85ff",
"metadata": {},
"outputs": [],
"source": [
"# Discretization parameters\n",
"xl = 0.0 # x-coordinate of the left boundary\n",
"xr = 1.0 # x-coordinate of the right boundary\n",
"ncells = 100 # number of cells in the discretization\n",
"nsteps = 100 # number of steps in the reactive transport simulation\n",
"dx = (xr - xl) / ncells # length of the mesh cells (in units of m)\n",
"dt = 10 * minute # time step\n",
"\n",
"# Physical parameters\n",
"D = 1.0e-9 # diffusion coefficient (in units of m2/s)\n",
"v = 1.0 / day # fluid pore velocity (in units of m/s)\n",
"T = 60.0 + 273.15 # temperature (in units of K)\n",
"P = 100 * 1e5 # pressure (in units of Pa)\n",
"phi = 0.1 # the porosity"
]
},
{
"cell_type": "markdown",
"id": "f063f707",
"metadata": {},
"source": [
"Next, we generate the coordinates of the mesh nodes (array `xcells`) by equally dividing the interval *[xr, xl]* with\n",
"the number of cells `ncells`. The length between each consecutive mesh node is computed and stored in `dx` (the\n",
"length of the mesh cells)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fec3e428",
"metadata": {},
"outputs": [],
"source": [
"xcells = np.linspace(xl, xr, ncells + 1) # interval [xl, xr] split into ncells"
]
},
{
"cell_type": "markdown",
"id": "c44a3501",
"metadata": {},
"source": [
"The boolean variable `dirichlet` is set to `True` or `False` depending on which boundary condition is considered in\n",
"the numerical calculation. `False` corresponds to imposing the flux of the injected fluid, otherwise, `True` means\n",
"imposing the composition of the fluid on the left boundary."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8037aa41",
"metadata": {},
"outputs": [],
"source": [
"dirichlet = False # parameter that determines whether Dirichlet BC must be used"
]
},
{
"cell_type": "markdown",
"id": "adb3c9cc",
"metadata": {},
"source": [
"To make sure that the applied finite-volume scheme is stable, we need to keep track of Courant–Friedrichs–Lewy (CFL)\n",
"number, which should be less than 1.0."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d2be91e7",
"metadata": {},
"outputs": [],
"source": [
"CFL = v * dt / dx\n",
"assert CFL <= 1.0, f\"Make sure that CFL = {CFL} is less that 1.0\""
]
},
{
"cell_type": "markdown",
"id": "615d73c3",
"metadata": {},
"source": [
"## Specifying the quantities and properties to be outputted\n",
"Before running the reactive transport simulations, we specify the list of parameters we are interested in\n",
"outputting. In this case, it is pH, molality of `H+`, `Ca++`, `Mg++`, `HCO3-`, `CO2(aq)`, as well as a phase volume\n",
"of calcite and dolomite."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d764af9d",
"metadata": {},
"outputs": [],
"source": [
"output_quantities = \"\"\"\n",
" pH\n",
" speciesMolality(H+)\n",
" speciesMolality(Ca++)\n",
" speciesMolality(Mg++)\n",
" speciesMolality(HCO3-)\n",
" speciesMolality(CO2(aq))\n",
" phaseVolume(Calcite)\n",
" phaseVolume(Dolomite)\n",
"\"\"\".split()"
]
},
{
"cell_type": "markdown",
"id": "ff214d4a",
"metadata": {},
"source": [
"Then, we define the list of name for the DataFrame columns. Note, that they must correspond\n",
"to the order of the properties defined in the `output_quantities` list:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "532824ac",
"metadata": {},
"outputs": [],
"source": [
"column_quantities = \"\"\"\n",
" pH\n",
" Hcation\n",
" Cacation\n",
" Mgcation\n",
" HCO3anion\n",
" CO2aq\n",
" calcite\n",
" dolomite\n",
"\"\"\".split()"
]
},
{
"cell_type": "markdown",
"id": "c9963d0d",
"metadata": {},
"source": [
"We create a DataFrame storage `df` using python library **pandas**. The columns in it match include the number of\n",
"reactive transport step, the space coordinates of the discretization cells, and the name of the properties we\n",
"assigned to be output (see `output_quantities` list above)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1354b098",
"metadata": {},
"outputs": [],
"source": [
"# Create the list of columns stored in dataframes\n",
"columns = ['step', 'x'] + column_quantities\n",
"import pandas as pd"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "00eb12be",
"metadata": {},
"outputs": [],
"source": [
"# Initialize dataframes with above defined columns\n",
"df = pd.DataFrame(columns=columns)"
]
},
{
"cell_type": "markdown",
"id": "682eaa44",
"metadata": {},
"source": [
"## Organization of the program\n",
"\n",
"The main part of the program (at the bottom of this tutorial) consists of three parts, each represented by a Python\n",
"function and documented in the following sections:\n",
"* creation of folders for the results (function `make_results_folders()`),\n",
"* simulation of reactive transport problem (method `simulate()`), and\n",
"* plotting of the obtained results.\n",
"\n",
"## Creating folders for the outputted results\n",
"Using **os** package, we create required folders for outputting the obtained results and for the plot and video\n",
"files later."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c1d0d6f6",
"metadata": {},
"outputs": [],
"source": [
"folder_results = 'results-rt-coupling'\n",
"def make_results_folders():\n",
" os.system('mkdir -p ' + folder_results)"
]
},
{
"cell_type": "markdown",
"id": "56e6f6c1",
"metadata": {},
"source": [
"## Performing the reactive transport simulation\n",
"The reactive transport simulation is performed in the function `simulate()`, which consists of several building\n",
"blocks:\n",
"* initialization of the reactive transport problem and\n",
"* performing the reactive transport simulation along defined time interval.\n",
"\n",
"The *preparatory initialization* step consists of the following sub-steps:\n",
"* definition of chemical system with its phases and species using `define_chemical_system()`,\n",
"* definition of the initial condition of the reactive transport problem in `define_initial_condition()`,\n",
"* definition of the boundary condition of the reactive transport problem in `define_initial_condition()`,\n",
"* generation of auxiliary indices to partition elements using `partition_indices()` and elements' partitioning \n",
"corresponding to fluid and solid species with function `partition_elements_in_mesh_cell()`, and finally\n",
"* definition of instance of [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html).\n",
"\n",
"The *simulation of the reactive transport problem* is represented by the loop over discretized time interval until\n",
"the final time is reached. On each step of this loop, the following functionality of performed:\n",
"* transport calculations using method `transport()`,\n",
"* reactive chemical calculations with `reactive_chemistry()` function, and\n",
"* saving the obtained results by means of `outputstate()`.\n",
"\n",
"Performing of the transport and reactive chemistry sequentially is possible due to the *operator splitting\n",
"procedure*, in which we first update the amounts of elements `b`. These updated amounts of\n",
"elements in each cell are used to evaluate its new chemical equilibrium state, thus producing new amounts of the\n",
"species in both the fluid and solid phases (available in the list `states`, the instances of\n",
"[ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) classes). This chemical reaction\n",
"equilibrium calculation step, at each mesh cell, permits, for example, aqueous species and minerals to react,\n",
"and thus causes mineral dissolution or precipitation, depending on how much the amount of mineral species changes.\n",
"This can then be used, for example, to compute new porosity value for the cell.\n",
"\n",
"The structure of the function `simulate`, performing the simulation of reactive transport problem, can be inspected\n",
"below:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5b5a06b0",
"metadata": {},
"outputs": [],
"source": [
"def simulate():\n",
" # Construct the chemical system with its phases and species\n",
" system = define_chemical_system()\n",
"\n",
" # Define the initial condition of the reactive transport modeling problem\n",
" state_ic = define_initial_condition(system)\n",
"\n",
" # Define the boundary condition of the reactive transport modeling problem\n",
" state_bc = define_boundary_condition(system)\n",
"\n",
" # Generate indices of partitioning fluid and solid species\n",
" nelems, ifluid_species, isolid_species = partition_indices(system)\n",
"\n",
" # Partitioning fluid and solid species\n",
" b, bfluid, bsolid, b_bc = partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc)\n",
"\n",
" # Create a list of chemical states for the mesh cells (one for each cell, initialized to state_ic)\n",
" states = [state_ic.clone() for _ in range(ncells + 1)]\n",
"\n",
" # Create the equilibrium solver object for the repeated equilibrium calculation\n",
" solver = EquilibriumSolver(system)\n",
"\n",
" # Running the reactive transport simulation loop\n",
" step = 0 # the current step number\n",
" t = 0.0 # the current time (in seconds)\n",
" # Save the initial state of the reactive transport calculation to the dataframe\n",
" outputstate_df(step, system, states)\n",
"\n",
" with tqdm(total=nsteps, desc=\"Reactive transport simulations\") as pbar:\n",
" while step <= nsteps:\n",
" # Perform transport calculations\n",
" bfluid, bsolid, b = transport(states, bfluid, bsolid, b, b_bc, nelems, ifluid_species, isolid_species)\n",
"\n",
" # Perform reactive chemical calculations\n",
" states = reactive_chemistry(solver, states, b)\n",
"\n",
" # Increment time step and number of time steps\n",
" t += dt\n",
" step += 1\n",
"\n",
" # Save the current state of the reactive transport calculation to the dataframe\n",
" outputstate_df(step, system, states)\n",
"\n",
" # Update a progress bar\n",
" pbar.update(1)\n",
"\n",
" print(\"Finished!\")"
]
},
{
"cell_type": "markdown",
"id": "c824676c",
"metadata": {},
"source": [
"Subsections below correspond to the methods responsible for each of the functional parts of `simulate()` method.\n",
"\n",
"### Construction of the chemical system with its phases and species\n",
"Method `define_chemical_system()` begins by defining the chemical system using [ChemicalEditor](\n",
"https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) class. In particular, we specify aqueous and mineral\n",
"phases that should be considered in the chemical system. For performance reasons, the aqueous phase is defined by\n",
"manually specifying the chemical species, instead of automatic collection of species from the database. There are\n",
"three pure mineral phase considered: quartz (SiO2), calcite (CaCO3), and dolomite\n",
"(CaMg(CO3)2). Finally, using class\n",
"[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html),\n",
"we define the chemical system."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4bbdde11",
"metadata": {},
"outputs": [],
"source": [
"def define_chemical_system():\n",
" # Define thermodynamical database to be used for chemical simulations\n",
" db = Database('supcrt98.xml')\n",
"\n",
" # Define chemical editor to work with the chemical system\n",
" editor = ChemicalEditor(db)\n",
"\n",
" # Define phases and corresponding to these phases species\n",
" editor.addAqueousPhaseWithElements('H O Na Cl Ca Mg C Si Ca')\n",
" editor.addMineralPhase('Quartz')\n",
" editor.addMineralPhase('Calcite')\n",
" editor.addMineralPhase('Dolomite')\n",
"\n",
" # Create the chemical system using the configured editor\n",
" system = ChemicalSystem(editor)\n",
"\n",
" return system"
]
},
{
"cell_type": "markdown",
"id": "46c4d51b",
"metadata": {},
"source": [
"### Initial condition of the reactive transport problem\n",
"After constructing the chemical system of interest, we can proceed to the definition of a chemical equilibrium\n",
"problem to set the *initial condition* for the fluid and rock composition that is consistent with our intention of\n",
"modeling reactive transport in a porous rock column composed of quartz and calcite, at 60 °C and 100 bar,\n",
"with a 0.7 NaCl molal brine in equilibrium with the rock minerals. To ensure equilibrium with both quartz and\n",
"calcite, the equilibrium problem setup below considers a relatively large amount for these minerals (10 moles each).\n",
"\n",
"Next, we use method [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n",
"to calculate the chemical equilibrium state of the system with the given initial conditions stored in the objects\n",
"`problem_ic`. The numerical solution of each problem results in the objects `state_ic` of class\n",
"[ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html), which stores the temperature,\n",
"pressure, and the amounts of every species in the system.\n",
"\n",
"The function ends with scaling the volumes of the aqueous and mineral phases so that they are consistent with a 10 %\n",
"porosity and the required volume percentages of the rock minerals (98 %vol of quartz and 2 %vol\n",
"of calcite)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3dac4745",
"metadata": {},
"outputs": [],
"source": [
"def define_initial_condition(system):\n",
" problem_ic = EquilibriumProblem(system)\n",
" problem_ic.setTemperature(T)\n",
" problem_ic.setPressure(P)\n",
" problem_ic.add('H2O', 1.0, 'kg')\n",
" problem_ic.add('NaCl', 0.7, 'mol')\n",
" problem_ic.add('CaCO3', 10, 'mol')\n",
" problem_ic.add('SiO2', 10, 'mol')\n",
"\n",
" # Calculate the equilibrium states for the initial conditions\n",
" state_ic = equilibrate(problem_ic)\n",
"\n",
" # Scale the volumes of the phases in the initial condition\n",
" state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')\n",
" state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')\n",
" state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')\n",
"\n",
" return state_ic"
]
},
{
"cell_type": "markdown",
"id": "96dbcf8a",
"metadata": {},
"source": [
"### Boundary condition of the reactive transport problem\n",
"For the boundary condition, we need to specify the composition of the fluid that is injected into the rock. This is\n",
"done below, by defining an equilibrium problem that will later be solved to produce an aqueous fluid with 0.9 molal\n",
"NaCl, 0.05 molal MgCl2, 0.01 CaCl2, and 0.75 molal CO2, in a state\n",
"very close to CO2 saturation.\n",
"\n",
"After equilibration, the obtained chemical state representing the boundary condition for the injected fluid\n",
"composition, we scale its volume to 1 m3. This is done so that the amounts of the species in the fluid are\n",
"consistent with a mol/m3 scale."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d0fe673a",
"metadata": {},
"outputs": [],
"source": [
"def define_boundary_condition(system):\n",
" # Define the boundary condition of the reactive transport modeling problem\n",
" problem_bc = EquilibriumProblem(system)\n",
" problem_bc.setTemperature(T)\n",
" problem_bc.setPressure(P)\n",
" problem_bc.add('H2O', 1.0, 'kg')\n",
" problem_bc.add('NaCl', 0.90, 'mol')\n",
" problem_bc.add('MgCl2', 0.05, 'mol')\n",
" problem_bc.add('CaCl2', 0.01, 'mol')\n",
" problem_bc.add('CO2', 0.75, 'mol')\n",
"\n",
" # Calculate the equilibrium states for the boundary conditions\n",
" state_bc = equilibrate(problem_bc)\n",
" # Scale the boundary condition state to 1 m3\n",
" state_bc.scaleVolume(1.0, 'm3')\n",
"\n",
" return state_bc"
]
},
{
"cell_type": "markdown",
"id": "015e4b8f",
"metadata": {},
"source": [
"### Indices of partitioning fluid and solid species\n",
"Only species in fluid phases are mobile and transported by advection and diffusion mechanisms. The solid phases are\n",
"immobile. The code below identifies the indices of the fluid and solid species. We use methods\n",
"[indicesFluidSpecies](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html#ac2a8b713f46f7a66b2731ba63faa95ad)\n",
"and [indicesSolidSpecies](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html#a8b0c237fff1d827f7bf2dbc911fa5bbf)\n",
"of class [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) to get the indices of the\n",
"fluid and solid species, which are stored in the lists `ifluid_species` and `isolid_species`, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "bf9e6bdc",
"metadata": {},
"outputs": [],
"source": [
"def partition_indices(system):\n",
" nelems = system.numElements()\n",
"\n",
" ifluid_species = system.indicesFluidSpecies()\n",
" isolid_species = system.indicesSolidSpecies()\n",
"\n",
" return nelems, ifluid_species, isolid_species"
]
},
{
"cell_type": "markdown",
"id": "767ee132",
"metadata": {},
"source": [
"### Partitioning fluid and solid species\n",
"In this function, we create arrays to keep track of the amounts of elements in the fluid and solid partition\n",
"(i.e., the amounts of elements among all fluid phases, here only an aqueous phase, and the amounts of elements among\n",
"all solid phases, here the mineral phases). For that, we define the arrays `b`, `bfluid`, `bsolid`, that\n",
"will store, respectively, the concentrations (mol/m3) of each element in the system, in the fluid\n",
"partition, and in the solid partition at every time step.\n",
"\n",
"The array `b` is initialized with the concentrations of the elements at the initial chemical state, `state_ic`,\n",
"using method [elementAmounts](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html#a827457e68a90f89920c13f0cc06fda78)\n",
"of class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html). The array `b_bc` stores\n",
"the concentrations of each element on the boundary in mol/m3fluid."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5690dfc9",
"metadata": {},
"outputs": [],
"source": [
"def partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc):\n",
" # The concentrations of each element in each mesh cell (in the current time step)\n",
" b = np.zeros((ncells, nelems))\n",
" # Initialize the concentrations (mol/m3) of the elements in each mesh cell\n",
" b[:] = state_ic.elementAmounts()\n",
"\n",
" # The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell\n",
" bfluid = np.zeros((ncells, nelems))\n",
"\n",
" # The concentrations (mol/m3) of each element in the solid partition, in each mesh cell\n",
" bsolid = np.zeros((ncells, nelems))\n",
"\n",
" # Initialize the concentrations (mol/m3) of each element on the boundary\n",
" b_bc = state_bc.elementAmounts()\n",
"\n",
" return b, bfluid, bsolid, b_bc"
]
},
{
"cell_type": "markdown",
"id": "73bd1da6",
"metadata": {},
"source": [
"### Reactive transport cycle\n",
"\n",
"#### Transport\n",
"This step updates in the fluid partition `bfluid` using the transport equations (without reactions).\n",
"The `transport_fullimplicit()` function below is responsible for solving an advection-diffusion equation, that is\n",
"later applied to transport the concentrations mol/m3 of elements in the fluid partition (*a\n",
"simplification that is possible because of common diffusion coefficients and velocities of the fluid species,\n",
"otherwise the transport of individual fluid species would be needed*).\n",
"\n",
"To match the units of concentrations of the elements in the fluid measured in mol/m3bulk and the\n",
"imposed concentration `b_bc[j]` mol/m3fluid, we need to multiply it by the porosity `phi_bc`\n",
"m3fluid/m3bulk on the boundary cell . We use function\n",
"[properties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html#ad3fa8fd9e1b948da7a698eb020513f3d)\n",
"of the class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) to retrieve fluid volume\n",
"m3fluid and total volume m3bulk in the inflow boundary cell.\n",
"\n",
"The updated amounts of elements in the fluid partition are then summed with the amounts of elements in the solid partition\n",
"`bsolid`, which remained constant during the transport step), and thus updating the amounts of elements in the\n",
"chemical system `b`. Reactive transport calculations involve the solution of a system of\n",
"advection-diffusion-reaction equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "87fd4ea8",
"metadata": {},
"outputs": [],
"source": [
"def transport(states, bfluid, bsolid, b, b_bc, nelems, ifluid_species, isolid_species):\n",
" # Collect the amounts of elements from fluid and solid partitions\n",
" for icell in range(ncells):\n",
" bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)\n",
" bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)\n",
"\n",
" # Get the porosity of the boundary cell\n",
" bc_cell = 0\n",
" phi_bc = states[bc_cell].properties().fluidVolume().val / states[bc_cell].properties().volume().val\n",
"\n",
" # Transport each element in the fluid phase\n",
" for j in range(nelems):\n",
" transport_fullimplicit(bfluid[:, j], dt, dx, v, D, phi_bc * b_bc[j])\n",
"\n",
" # Update the amounts of elements in both fluid and solid partitions\n",
" b[:] = bsolid + bfluid\n",
"\n",
" return bfluid, bsolid, b"
]
},
{
"cell_type": "markdown",
"id": "1e4b2db3",
"metadata": {},
"source": [
"##### Transport calculation with finite-volume scheme\n",
"The function `transport()` expects a conservative property (argument `u`) (e.g., the concentration mol/m3\n",
"of *j*th element in the fluid given by `bfluid[j]`), the time step (`dt`), the mesh cell length (`dx`),\n",
"the fluid velocity (`v`), the diffusion coefficient (`D`), and the boundary condition of the conservative property\n",
"(`g`) (e.g., the concentration of the *j*th element in the fluid on the left boundary).\n",
"\n",
"The transport equations are solved with a finite volume method, where diffusion and convection are treated implicitly.\n",
"Its discretization in space and time (implicit) results in the constants `alpha` and `beta`. These correspond to\n",
"the diffusion and advection terms in the equation: `D*dt/dx**2` and `v*dt/dx`, respectively.\n",
"\n",
"Arrays `a`, `b`, `c` are the diagonals in the tridiagonal matrix that results by writing all discretized equations\n",
"in a matrix equation. This system of linear equations is solved by the tridiagonal matrix algorithm, also known\n",
"as the Thomas algorithm."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "27d8d6b4",
"metadata": {},
"outputs": [],
"source": [
"def transport_fullimplicit(u, dt, dx, v, D, ul):\n",
" # Number of DOFs\n",
" n = len(u)\n",
" alpha = D * dt / dx ** 2\n",
" beta = v * dt / dx\n",
"\n",
" # Upwind finite volume scheme\n",
" a = np.full(n, -beta - alpha)\n",
" b = np.full(n, 1 + beta + 2 * alpha)\n",
" c = np.full(n, -alpha)\n",
"\n",
" # Set the boundary condition on the left cell\n",
" if dirichlet:\n",
" # Use Dirichlet BC boundary conditions\n",
" b[0] = 1.0\n",
" c[0] = 0.0\n",
" u[0] = ul\n",
"\n",
" else:\n",
" # Flux boundary conditions (implicit scheme for the advection)\n",
" # Left boundary\n",
" b[0] = 1 + alpha + beta\n",
" c[0] = -alpha # stays the same as it is defined -alpha\n",
" u[0] += beta * ul # = dt/dx * v * g, flux that we prescribe is equal v * ul\n",
"\n",
" # Right boundary is free\n",
" a[-1] = - beta\n",
" b[-1] = 1 + beta\n",
"\n",
" # Solve a tridiagonal matrix equation\n",
" thomas(a, b, c, u)"
]
},
{
"cell_type": "markdown",
"id": "e108544f",
"metadata": {},
"source": [
"##### Solving the system of equations obtained from finite volume discretization\n",
"The tridiagonal matrix equation is solved using the Thomas algorithm (or the TriDiagonal Matrix Algorithm (TDMA)).\n",
"It is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4661aa98",
"metadata": {},
"outputs": [],
"source": [
"def thomas(a, b, c, d):\n",
" n = len(d)\n",
" c[0] /= b[0]\n",
" for i in range(1, n - 1):\n",
" c[i] /= b[i] - a[i] * c[i - 1]\n",
" d[0] /= b[0]\n",
" for i in range(1, n):\n",
" d[i] = (d[i] - a[i] * d[i - 1]) / (b[i] - a[i] * c[i - 1])\n",
" x = d\n",
" for i in reversed(range(0, n - 1)):\n",
" x[i] -= c[i] * x[i + 1]\n",
" return x"
]
},
{
"cell_type": "markdown",
"id": "637602fb",
"metadata": {},
"source": [
"#### Reactive chemistry\n",
"The chemical equilibrium calculations performed in each mesh cell, using *Gibbs energy minimisation* algorithm (\n",
"provided by the class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html))."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "3d462ab7",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def reactive_chemistry(solver, states, b):\n",
" # Equilibrating all cells with the updated element amounts\n",
" for icell in range(ncells):\n",
" solver.solve(states[icell], T, P, b[icell])\n",
" return states"
]
},
{
"cell_type": "markdown",
"id": "dba89f71",
"metadata": {},
"source": [
"### Saving and analyzing of the results\n",
"\n",
"Function `outputstate_df` is the auxiliary function to add data to the DataFrame at each time step."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "cbecec66",
"metadata": {},
"outputs": [],
"source": [
"def outputstate_df(step, system, states):\n",
" # Define the instance of ChemicalQuantity class\n",
" quantity = ChemicalQuantity(system)\n",
"\n",
" # Create the list with empty values to populate with chemical properties\n",
" values = [None] * len(columns)\n",
" for state, x in zip(states, xcells):\n",
"\n",
" # Populate values with number of reactive transport step and spacial coordinates\n",
" values[0] = step\n",
" values[1] = x\n",
"\n",
" # Update the\n",
" quantity.update(state)\n",
" for quantity_name, i in zip(output_quantities, range(2, len(states))):\n",
" values[i] = quantity.value(quantity_name) * (100 / (1 - phi) if \"phaseVolume\" in quantity_name else 1)\n",
" df.loc[len(df)] = values"
]
},
{
"cell_type": "markdown",
"id": "e525df24",
"metadata": {},
"source": [
"# Main parts of the tutorial\n",
"\n",
"First, we create folders for the results:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "43647f98",
"metadata": {},
"outputs": [],
"source": [
"make_results_folders()"
]
},
{
"cell_type": "markdown",
"id": "435b9b13",
"metadata": {},
"source": [
"Run the reactive transport simulations:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6175a92",
"metadata": {},
"outputs": [],
"source": [
"simulate()"
]
},
{
"cell_type": "markdown",
"id": "335dbf0b",
"metadata": {},
"source": [
"To inspect the collected data, one can run:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "30c0ad16",
"metadata": {},
"outputs": [],
"source": [
"df"
]
},
{
"cell_type": "markdown",
"id": "d4170199",
"metadata": {},
"source": [
"## Plotting of the obtained results\n",
"\n",
"The last block of the main routine is dedicated to the plotting of the results in a Jupyter app generated by the\n",
"library **bokeh**. It is an interactive visualization library that provides elegant, concise construction of\n",
"versatile graphics, and affords high-performance interactivity over large or streaming datasets.\n",
"\n",
"First, we have to collect files with results corresponding to each step:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6b9b10bf",
"metadata": {},
"outputs": [],
"source": [
"print(\"Collecting files...\")\n",
"files = [file for file in natsorted(os.listdir(folder_results))]"
]
},
{
"cell_type": "markdown",
"id": "f74e3144",
"metadata": {},
"source": [
"To make a Bokeh application that helps us study the results of the reactive transport simulations, a function\n",
"`modify_doc_df(doc)` that creates Bokeh content and adds it to the app."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "0a4dcdac",
"metadata": {},
"outputs": [],
"source": [
"def modify_doc_df(doc):\n",
" # Initialize the data by the initial chemical state\n",
" source = ColumnDataSource(df[df['step'] == 0])\n",
"\n",
" # Auxiliary function that returns a string for the title of a figure in the format Time: #h##m\n",
" def titlestr(t):\n",
" t = t / minute # Convert from seconds to minutes\n",
" h = int(t) / 60 # The number of hours\n",
" m = int(t) % 60 # The number of remaining minutes\n",
" return 'Time: %2dh %2dm' % (h, m)\n",
"\n",
" # Plot for ph\n",
" p1 = figure(plot_width=600, plot_height=250)\n",
" p1.line(x='x', y='pH',\n",
" color='teal', line_width=2, legend_label='pH', source=source)\n",
" p1.x_range = Range1d(-0.001, 1.001)\n",
" p1.y_range = Range1d(4.0, 9.0)\n",
" p1.xaxis.axis_label = 'Distance [m]'\n",
" p1.yaxis.axis_label = 'pH'\n",
" p1.legend.location = 'bottom_right'\n",
" p1.title.text = titlestr(0 * dt)\n",
"\n",
" # Plot for calcite and dolomite\n",
" p2 = figure(plot_width=600, plot_height=250)\n",
" p2.line(x='x', y='calcite', color='blue', line_width=2, legend_label='Calcite',\n",
" muted_color='blue', muted_alpha=0.2,\n",
" source=source)\n",
" p2.line(x='x', y='dolomite', color='orange', line_width=2, legend_label='Dolomite',\n",
" muted_color='orange', muted_alpha=0.2,\n",
" source=source)\n",
" p2.x_range = Range1d(-0.001, 1.001)\n",
" p2.y_range = Range1d(-0.1, 2.1)\n",
" p2.xaxis.axis_label = 'Distance [m]'\n",
" p2.yaxis.axis_label = 'Mineral Volume [%vol]'\n",
" p2.legend.location = 'center_right'\n",
" p2.title.text = titlestr(0 * dt)\n",
" p2.legend.click_policy = 'mute'\n",
"\n",
" # Plot for aqueous species\n",
" p3 = figure(plot_width=600, plot_height=300, y_axis_type='log')\n",
" p3.line(x='x', y='Cacation', color='blue', line_width=2, legend_label='Ca++', source=source)\n",
" p3.line(x='x', y='Mgcation', color='orange', line_width=2, legend_label='Mg++', source=source)\n",
" p3.line(x='x', y='HCO3anion', color='green', line_width=2, legend_label='HCO3-', source=source)\n",
" p3.line(x='x', y='CO2aq', color='red', line_width=2, legend_label='CO2(aq)', source=source)\n",
" p3.line(x='x', y='Hcation', color='darkviolet', line_width=2, legend_label='H+', source=source)\n",
" p3.x_range = Range1d(-0.001, 1.001)\n",
" p3.y_range = Range1d(0.5e-6, 1e0)\n",
" p3.xaxis.axis_label = 'Distance [m]'\n",
" p3.yaxis.axis_label = 'Concentration [molal]'\n",
" p3.legend.location = 'top_right'\n",
" p3.title.text = titlestr(0 * dt)\n",
" p3.legend.click_policy = 'mute'\n",
"\n",
" step_slider = Slider(start=0,\n",
" end=nsteps,\n",
" value=0,\n",
" step=1,\n",
" title='Reactive Transport Step')\n",
"\n",
" # Callback function that changes the data and the title of the plots\n",
" def update_step(attrname, old, new):\n",
" step = step_slider.value\n",
" new_source = ColumnDataSource(df[df['step'] == step])\n",
" source.data = dict(x=new_source.data['x'],\n",
" pH=new_source.data['pH'],\n",
" calcite=new_source.data['calcite'],\n",
" dolomite=new_source.data['dolomite'],\n",
" Hcation=new_source.data['Hcation'],\n",
" Cacation=new_source.data['Cacation'],\n",
" Mgcation=new_source.data['Mgcation'],\n",
" HCO3anion=new_source.data['HCO3anion'],\n",
" CO2aq=new_source.data['CO2aq'])\n",
" p1.title.text = titlestr(step * dt)\n",
" p2.title.text = titlestr(step * dt)\n",
" p3.title.text = titlestr(step * dt)\n",
"\n",
" # Put the update function on the slider\n",
" step_slider.on_change('value', update_step)\n",
"\n",
" layout = column(step_slider,\n",
" column(p1, p2, p3))\n",
"\n",
" doc.add_root(layout)"
]
},
{
"cell_type": "markdown",
"id": "261025f7",
"metadata": {},
"source": [
"Outputting the plots to the notebook requires the call of `output_notebook()` that specifies outputting the plot\n",
"inline in the Jupyter notebook. Finally, the function `modify_doc()` must be passed to `show`, so that the app defined\n",
"by it is displayed inline.\n",
"\n",
"> **Important:** If you run this tutorial in the *localhost*, make sure that number provided to the variable\n",
"`notebook_url` below coincides with the number of the localhost you have in your browser.\n",
"\n",
"In the app below, we provide an interactive *slider* witget that defines the number of the reactive transport time\n",
"step. This automatically updates the plots with ph, volume phases of calcite and dolomite, and mollalities of aqueous\n",
"species (in logarithmic scale)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2a1707d0",
"metadata": {},
"outputs": [],
"source": [
"output_notebook()\n",
"show(modify_doc_df, notebook_url=\"http://localhost:8890\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}