{ "cells": [ { "cell_type": "markdown", "id": "1e067f1e", "metadata": {}, "source": [ "# Reactive transport of CO2-saturated brine along a porous rock column (using transport solver embedded into Reaktoro)\n", "\n", "In this tutorial, we show how Reaktoro can be used for one-dimensional reactive transport calculations for modeling\n", "the geochemical reactions that occur along a porous rock column as an aqueous fluid is continuously injected on its\n", "left side.\n", "\n", "The injected fluid is a brine with 0.9 molal NaCl, 0.05 molal MgCl2, 0.01 molal CaCl2\n", "and almost CO2-saturated, with 0.75 molal of CO2 dissolved.\n", "The porous rock is initially composed of minerals quartz SiO2 and calcite CaCO3. The\n", "initial porosity is 10 %, and the initial volume percentages of the minerals are 98 %vol of quartz and\n", "2 %vol calcite. The initial conditions for the fluid in the rock is a 0.7 molal NaCl brine in\n", "equilibrium with the existing rock minerals calcite and quartz. These initial fluid and rock composition conditions\n", "are uniform throughout the rock core. We assume a rock column length of 100 m at temperature 60 °C and 100 bar\n", "throughout.\n", "\n", "**Assumptions**: To simplify this tutorial, the following assumptions are made:\n", "* Chemical equilibrium is used for modeling the chemical reactions in this problem, not only for reactions between\n", "aqueous-aqueous species but also for those between mineral and aqueous species.\n", "* A uniform constant velocity field is imposed and it is not updated by solving, for example, Darcy equation.\n", "* Both temperature and pressure are also kept constant along the rock.\n", "\n", "## Import the reaktoro Python package (and other packages)\n", "First, we import the **reaktoro** Python package so that we can use its classes\n", "and methods for performing the chemical reaction calculations." ] }, { "cell_type": "code", "execution_count": null, "id": "984bee7b", "metadata": {}, "outputs": [], "source": [ "print('============================================================')\n", "print('Make sure you have the following Python packages installed: ')\n", "print(' numpy, matplotlib, natsort')\n", "print('These can be installed with pip:')\n", "print(' pip install numpy matplotlib natsort')\n", "print('============================================================')\n", "from reaktoro import *\n", "import numpy as np\n", "import os\n", "import matplotlib.pyplot as plt\n", "from natsort import natsorted" ] }, { "cell_type": "markdown", "id": "121568ab", "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, **os** to provide a portable way of using operating system\n", "dependent functionality, **matplotlib** for plotting capabilities, and **natsort** for sorting the lists.\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", "## Defining auxiliary time-related constants\n", "In this step, we initialize auxiliary time-related constants from seconds to years. This is only done for\n", "convenience, so that we can specify later, for example, fluid velocity as 1 m/week." ] }, { "cell_type": "code", "execution_count": null, "id": "8a0e1411", "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": "8f977522", "metadata": {}, "source": [ "## Defining parameters for the reactive transport simulation\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/week (1.16 · 10-5 m/s) and the same\n", "diffusion coefficient of 10-9 m2/s for all fluid species (without dispersivity). The size of the\n", "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": "c5ff70eb", "metadata": {}, "outputs": [], "source": [ "# Discretisation parameters\n", "xl = 0.0 # the x-coordinate of the left boundary\n", "xr = 1.0 # the x-coordinate of the right boundary\n", "ncells = 100 # the number of cells in the discretization\n", "nsteps = 500 # the number of steps in the reactive transport simulation\n", "dt = 10*minute # the time step (30 minutes in units of s)\n", "dx = (xr - xl)/ncells # length of the mesh cells (in units of m)\n", "\n", "# Physical and chemical parameters\n", "D = 1.0e-9 # the diffusion coefficient (in units of m2/s)\n", "v = 1.0/week # the fluid pore velocity (1 m/week in units of m/s)\n", "T = 60.0 # the temperature (in units of degC)\n", "P = 100 # the pressure (in units of bar)\n", "phi = 0.1 # the porosity" ] }, { "cell_type": "markdown", "id": "9303a339", "metadata": {}, "source": [ "Next, we generate the coordinates of the mesh nodes (array `x`) by equally dividing the interval *[xr, xl]* with\n", "the number of cells `ncells`. The length between each consecutive mesh nodes is computed and stored in `dx` (the\n", "length of the mesh cells)." ] }, { "cell_type": "code", "execution_count": null, "id": "41358068", "metadata": {}, "outputs": [], "source": [ "xcells = np.linspace(xl, xr, ncells) # interval [xl, xr] split into ncells" ] }, { "cell_type": "markdown", "id": "5ba3476f", "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": "495a66a9", "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": "1398b864", "metadata": {}, "source": [ "Using **os** package, we create required folders for outputting the obtained results and folders to save video files." ] }, { "cell_type": "code", "execution_count": null, "id": "102373c0", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "folder_results = 'results-rt-calcite-brine'\n", "folder_videos = 'videos-rt-calcite-brine'\n", "def make_results_folders():\n", " os.system('mkdir -p ' + folder_results)\n", " os.system('mkdir -p ' + folder_videos)" ] }, { "cell_type": "markdown", "id": "26cb12d2", "metadata": {}, "source": [ "## Reactive transport simulations\n", "\n", "### Defining the chemical system\n", "\n", "We need to define a chemical system that can represent both our fluid and rock. We use class\n", "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) below to define a system with an\n", "aqueous phase and three mineral phases: quartz, calcite, and dolomite. Initially, our rock has no dolomite\n", "(CaMg(CO3)2), but since this is a mineral that could potentially precipitate given the fluid\n", "composition injected ( containing CaCl2 and MgCl2 dissolved), we add it here in the\n", "chemical system to ensure that the calculations are able to model dolomite precipitation." ] }, { "cell_type": "code", "execution_count": null, "id": "09c9f2b2", "metadata": {}, "outputs": [], "source": [ "db = Database('supcrt98.xml')\n", "editor = ChemicalEditor(db)\n", "editor.addAqueousPhase('H2O(l) H+ OH- Na+ Cl- Ca++ Mg++ HCO3- CO2(aq) CO3--')\n", "editor.addMineralPhase('Quartz')\n", "editor.addMineralPhase('Calcite')\n", "editor.addMineralPhase('Dolomite')" ] }, { "cell_type": "markdown", "id": "557649fb", "metadata": {}, "source": [ "> **Note**: The aqueous phase is defined above by using a list of compounds, which is then broken automatically by\n", "> Reaktoro into a list of element names. These element names are then used to find in the database all the aqueous\n", "> species that could be formed out of them.\n", "\n", "### Constructing the chemical system\n", "\n", "This step is where we create an object of class\n", "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) using the\n", "chemical system definition details stored in the object ``editor``." ] }, { "cell_type": "code", "execution_count": null, "id": "7ef280b2", "metadata": {}, "outputs": [], "source": [ "system = ChemicalSystem(editor)" ] }, { "cell_type": "markdown", "id": "756392a5", "metadata": {}, "source": [ "### Initial condition for the fluid composition\n", "\n", "Below, we define the **initial condition** for the fluid composition in the rock. We want an aqueous fluid that is\n", "0.7 molal of NaCl and in equilibrium with calcite and quartz. To achieve this, we mix 1 kg of water, 0.7 mol of\n", "NaCl, and plenty of calcite and quartz (10 mol each) to ensure that the aqueous solution is saturated with respect\n", "to these minerals." ] }, { "cell_type": "code", "execution_count": null, "id": "1b36e3b4", "metadata": {}, "outputs": [], "source": [ "problem_ic = EquilibriumProblem(system)\n", "problem_ic.setTemperature(T, 'celsius')\n", "problem_ic.setPressure(P, 'bar')\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')" ] }, { "cell_type": "markdown", "id": "8490a0d6", "metadata": {}, "source": [ "### Boundary condition for the fluid composition\n", "\n", "Next, we define the **boundary condition** for the fluid composition on the left side of the rock, which should be\n", "the one that represents the fluid being continuously injected: 0.9 molal NaCl, 0.05 molal MgCl2,\n", "0.01 molal CaCl2 and almost CO2-saturated, with 0.75 molal of CO2 dissolved. To\n", "provide that, we mix 1 kg of HO2 with 0.9 mol of NaCl, 0.05 mol of MgCl2, 0.01 mol\n", "of CaCl2, and 0.75 mol of CO2." ] }, { "cell_type": "code", "execution_count": null, "id": "2d383f57", "metadata": {}, "outputs": [], "source": [ "problem_bc = EquilibriumProblem(system)\n", "problem_bc.setTemperature(T, 'celsius')\n", "problem_bc.setPressure(P, 'bar')\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')" ] }, { "cell_type": "markdown", "id": "0de80e4f", "metadata": {}, "source": [ "### Calculating the IC and BC fluid compositions\n", "In this step, we use the [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n", "function to calculate the chemical equilibrium state of the system with the given initial and boundary equilibrium\n", "conditions stored in the object `problem_ic` and `problem_bc`, respectively. The result is stored in the\n", "corresponding instances of the class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html),\n", "i.e., `state_ic` and `state_bc`." ] }, { "cell_type": "code", "execution_count": null, "id": "742f97d3", "metadata": {}, "outputs": [], "source": [ "state_ic = equilibrate(problem_ic)\n", "state_bc = equilibrate(problem_bc)" ] }, { "cell_type": "markdown", "id": "b376b093", "metadata": {}, "source": [ "### Scaling the phases in the initial condition\n", "\n", "The initial chemical state `state_ic` computed before has, at this point, phases with volumes that do not\n", "correspond to our desired porosity of 10 % and rock mineral composition of 98 %vol of quartz and\n", "2 %vol of calcite.\n", "\n", "To obtain this, we scale the volumes of the aqueous and mineral phases as shown\n", "below:" ] }, { "cell_type": "code", "execution_count": null, "id": "4c1f474c", "metadata": {}, "outputs": [], "source": [ "# Scale the volumes of the phases in the initial condition\n", "state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3') # corresponds to the initial porosity of 10%.\n", "state_ic.scalePhaseVolume('Quartz', 0.882, 'm3') # 0.882 = 0.9 * 0.98\n", "state_ic.scalePhaseVolume('Calcite', 0.018, 'm3') # 0.018 = 0.9 * 0.02" ] }, { "cell_type": "markdown", "id": "06730698", "metadata": {}, "source": [ "> **Note**: After this scaling step, the sum of the phase volumes in ``state_ic`` is 1 m3. This also\n", "> ensures that the amounts of the species in the chemical system are normalized by m3, and thus they can\n", "> be regarded as concentrations in a unit of mol/m3 (*bulk volume, not fluid volume!*).\n", "\n", "### Scaling the boundary condition state\n", "\n", "Next, we scale the boundary condition state to 1 m3, so that we have the amounts of fluid species in\n", "`state_bc` also normalized by m3.\n", "\n", "> **Note**: The chemical state represented by `state_bc` has no other stable phase than the aqueous phase (i.e.,\n", "> all mineral phases have zero or negligible amounts such as 10-21 mol)." ] }, { "cell_type": "code", "execution_count": null, "id": "8467ff57", "metadata": {}, "outputs": [], "source": [ "state_bc.scaleVolume(1.0, 'm3')" ] }, { "cell_type": "markdown", "id": "5dfb12d3", "metadata": {}, "source": [ "### Creating the mesh\n", "\n", "We define the mesh with the class [Mesh](https://reaktoro.org/cpp/classReaktoro_1_1Mesh.html) to use in\n", "the initialization of class [ReactiveTransportSolver](\n", "https://reaktoro.org/cpp/classReaktoro_1_1ReactiveTransportSolver.html) later. Here, we specify the number of cells\n", "in the mesh and the x-coordinates of the left and right boundaries (in m)." ] }, { "cell_type": "code", "execution_count": null, "id": "3d325ace", "metadata": {}, "outputs": [], "source": [ "mesh = Mesh(ncells, xl, xr)" ] }, { "cell_type": "markdown", "id": "9f68c7df", "metadata": {}, "source": [ "### Creating a chemical field object\n", "\n", "We have been using class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) to represent\n", "an individual chemical state. We will now use class [ChemicalField](\n", "https://reaktoro.org/cpp/classReaktoro_1_1ChemicalField.html) to represent a collection of chemical states: one for\n", "each mesh cell.\n", "\n", "> **Note**: Different initial conditions across the mesh cells are possible by assigning different chemical states to\n", "> each mesh cell. Here, the same chemical state in `state_ic` is used for all cells." ] }, { "cell_type": "code", "execution_count": null, "id": "84df0576", "metadata": {}, "outputs": [], "source": [ "field = ChemicalField(mesh.numCells(), state_ic)" ] }, { "cell_type": "markdown", "id": "6b92a39d", "metadata": {}, "source": [ "### Initializing the reactive transport solver\n", "\n", "At last, we define the object responsible for solving the reactive transport problem, which is handled by the\n", "class [ReactiveTransportSolver](https://reaktoro.org/cpp/classReaktoro_1_1ReactiveTransportSolver.html).\n", "Here, we set the mesh and problem parameters such as velocity, diffusion coefficient, the chemical state\n", "representing the boundary condition, and the time step. We also initialize the reactive solver object with the\n", "chemical field object specified on the previous step, at this point containing the initial condition for the\n", "chemical state of each mesh cell." ] }, { "cell_type": "code", "execution_count": null, "id": "396f4265", "metadata": {}, "outputs": [], "source": [ "rt = ReactiveTransportSolver(system)\n", "rt.setMesh(mesh)\n", "rt.setVelocity(v)\n", "rt.setDiffusionCoeff(D)\n", "rt.setBoundaryState(state_bc)\n", "rt.setTimeStep(dt)\n", "rt.initialize(field)" ] }, { "cell_type": "markdown", "id": "c7ef041c", "metadata": {}, "source": [ "### Defining the output quantities\n", "\n", "Before starting the reactive transport calculations, we define the quantities that will be output for every mesh\n", "cell, at every time step. For this, we use an object of the class\n", "[ChemicalOutput](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html).\n", "The name of the output file is to `reactive-transport.txt`. We specify the parameters that 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": "ecae3b52", "metadata": {}, "outputs": [], "source": [ "output = rt.output()\n", "output.add(\"pH\")\n", "output.add(\"speciesMolality(H+)\")\n", "output.add(\"speciesMolality(Ca++)\")\n", "output.add(\"speciesMolality(Mg++)\")\n", "output.add(\"speciesMolality(HCO3-)\")\n", "output.add(\"speciesMolality(CO2(aq))\")\n", "output.add(\"phaseVolume(Calcite)\")\n", "output.add(\"phaseVolume(Dolomite)\")\n", "output.filename(folder_results + '/state.txt') # Set the name of the output files" ] }, { "cell_type": "code", "execution_count": null, "id": "faae451b", "metadata": {}, "outputs": [], "source": [ "# Make auxiliary folders to save generated results, their plots, or videos\n", "make_results_folders()" ] }, { "cell_type": "markdown", "id": "5f0ca913", "metadata": {}, "source": [ "### Running the reactive transport simulation\n", "\n", "As shown below, we perform a sequence of reactive transport calculations, one for each time step, during which the\n", "chemical state of each mesh cell is updated. The iterations continue until the maximum number of steps is\n", "achieved.\n", "\n", "Using **tqdm** we track the progress of simulations using the progress bar. For that, we wrap the while-loop with\n", "the function 'tqdm()'. We then use the method `step` of class\n", "[ReactiveTransportSolver](https://reaktoro.org/cpp/classReaktoro_1_1ReactiveTransportSolver.html) to perform a\n", "single reactive transport time-stepping. This method also produces a new output file containing the requested\n", "output properties for every mesh cell. In each such file, rows correspond to cells, whereas the columns correspond\n", "to the requested output properties, i.e., pH, molality of `H+`, `Ca++`, `Mg++`, `HCO3-`, `CO2(aq)`, as\n", "well as the phase volume of calcite and dolomite." ] }, { "cell_type": "code", "execution_count": null, "id": "76f0f7ce", "metadata": {}, "outputs": [], "source": [ "# Step 15: Perform given number of reactive tranport steps\n", "t = 0.0 # current time variable\n", "step = 0 # current number of steps\n", "\n", "from tqdm.notebook import tqdm\n", "with tqdm(total=nsteps, desc=\"Reactive transport simulations\") as pbar:\n", " while step <= nsteps: # step until the number of steps are achieved\n", "\n", " # Perform one reactive transport time step\n", " rt.step(field)\n", "\n", " # Increment time step and number of time steps\n", " t += dt\n", " step += 1\n", "\n", " # Update a progress bar\n", " pbar.update(1)" ] }, { "cell_type": "markdown", "id": "7359578b", "metadata": {}, "source": [ "## Plotting of the obtained results\n", "The last block of the main routine is dedicated to plotting of the results and generating a video from the plots to\n", "illustrate the time-dependent behavior of the chemical properties. It uses parallel pthread to run `plotfile`\n", "function for each file from the list `files`.\n", "\n", "First, we collect files with results using `listdir` function, which returns the list containing the names of\n", "the entries in the directory given by path `folder_results`:" ] }, { "cell_type": "code", "execution_count": null, "id": "4cf93597", "metadata": {}, "outputs": [], "source": [ "files = [file for file in natsorted( os.listdir(folder_results) ) ]" ] }, { "cell_type": "markdown", "id": "14500e75", "metadata": {}, "source": [ "To generate animations, we exploit **animation** module of the library **matplotlib**, which provides the framework\n", "to build videos from the plots. **Note**: **ffmpeg** must be installed for handling video files and streams." ] }, { "cell_type": "code", "execution_count": null, "id": "4d0a71d7", "metadata": {}, "outputs": [], "source": [ "from matplotlib import animation\n", "from IPython.display import Video\n", "\n", "animation_starts_at_frame = 0 # the first frame index to be considered\n", "animation_ends_at_frame = 10 * 30 # the last frame index to be considered\n", "animation_num_frames_to_jump = 1 # the number of frames to jump between current and next\n", "# Check for the correct end frame number\n", "assert animation_ends_at_frame <= nsteps, \"WARNING: The number of the end frame must be smaller then number of steps! \"" ] }, { "cell_type": "markdown", "id": "64e7dde8", "metadata": {}, "source": [ "Provide the number of frames per second and the time (in milliseconds) to wait between each frame:" ] }, { "cell_type": "code", "execution_count": null, "id": "1bdde0f4", "metadata": {}, "outputs": [], "source": [ "animation_fps = 30 # the number of frames per second\n", "animation_interval_wait = 200 # the time (in milliseconds) to wait between each frame\n", "# Auxiliary animation options\n", "animation_frame_range = range(animation_starts_at_frame, animation_ends_at_frame, animation_num_frames_to_jump)" ] }, { "cell_type": "markdown", "id": "96de4bbb", "metadata": {}, "source": [ "For plotting of the data saved under the `folder_results` folder, we provide the indices corresponding to the columns\n", "written to the `state.txt` files." ] }, { "cell_type": "code", "execution_count": null, "id": "83d60b2b", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "indx_ph = 0\n", "indx_Hcation = 1\n", "indx_Cacation = 2\n", "indx_Mgcation = 3\n", "indx_HCO3anion = 4\n", "indx_CO2aq = 5\n", "indx_calcite = 6\n", "indx_dolomite = 7" ] }, { "cell_type": "markdown", "id": "9294bf0e", "metadata": {}, "source": [ "Routines `plot_animation_ph()`, `plot_animation_calcite_dolomite()`, and `plot_animation_aqueous_species()`\n", "are dedicated to animating the time-dependent behavior of the chemical properties." ] }, { "cell_type": "code", "execution_count": null, "id": "79b364aa", "metadata": { "lines_to_end_of_cell_marker": 0, "lines_to_next_cell": 1 }, "outputs": [], "source": [ "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", "\n", " return 'Time: %2dh %2dm' % (h, m)\n", "\n", "def line(color):\n", " return {'linestyle': '-', 'color': color, 'zorder': 1, 'linewidth': 2}\n", "\n", "def plot_animation_ph():\n", "\n", " # Plot of mineral's volume the space coordinates\n", " fig = plt.figure()\n", " ax = plt.axes(xlim=(-0.01, 0.501), ylim=(4.0, 12.0))\n", " ax.set_xlabel('Distance [m]')\n", " ax.set_ylabel('pH')\n", " ax.set_title(titlestr(0.0))\n", " objects = [\n", " ax.plot([], [], label='pH', **line('teal'))[0],\n", " ]\n", " ax.legend(loc='lower right')\n", "\n", " def init():\n", " return tuple(objects)\n", "\n", " def animate(i):\n", " t = i * dt\n", " filearray = np.loadtxt(folder_results + '/' + files[i], skiprows=1)\n", " data = filearray.T\n", " data_ph = data[indx_ph]\n", " objects[0].set_data(xcells, data_ph)\n", " ax.set_title(titlestr(t))\n", " return tuple(objects)\n", "\n", " print(\"Generating the animation of pH behaviour ...\")\n", " anim = animation.FuncAnimation(fig, animate, init_func=init, frames=animation_frame_range, interval=animation_interval_wait, blit=True)\n", " anim.save(folder_videos + '/pH.mp4', fps=animation_fps, extra_args=['-vcodec', 'libx264'])\n", " print(\"Finished!\")\n", "\n", "def plot_animation_calcite_dolomite():\n", "\n", " # Plot of mineral's volume the space coordinates\n", " fig = plt.figure()\n", " ax = plt.axes(xlim=(-0.01, 0.501), ylim=(-0.1, 2.1))\n", " ax.set_xlabel('Distance [m]')\n", " ax.set_ylabel('Mineral Volume [%$_{\\mathsf{vol}}$]')\n", " ax.set_title(titlestr(0.0))\n", " objects = [\n", " ax.plot([], [], label='Calcite', **line('C0'))[0],\n", " ax.plot([], [], label='Dolomite', **line('C1'))[0],\n", " ]\n", " ax.legend(loc='center right')\n", "\n", "\n", " def init():\n", " return tuple(objects)\n", "\n", "\n", " def animate(i):\n", " t = i * dt\n", " filearray = np.loadtxt(folder_results + '/' + files[i], skiprows=1)\n", " data = filearray.T\n", " data_calcite, data_dolomite = data[indx_calcite], data[indx_dolomite]\n", " objects[0].set_data(xcells, data_calcite * 100/(1 - phi))\n", " objects[1].set_data(xcells, data_dolomite * 100/(1 - phi))\n", " ax.set_title(titlestr(t))\n", " return tuple(objects)\n", "\n", " print(\"Generating the animation of calcite-dolomite behaviour ...\")\n", " anim = animation.FuncAnimation(fig, animate, init_func=init, frames=animation_frame_range, interval=animation_interval_wait, blit=True)\n", " anim.save(folder_videos + '/calcite-dolomite.mp4', fps=animation_fps, extra_args=['-vcodec', 'libx264'])\n", " print(\"Finished!\")\n", "\n", "def plot_animation_aqueous_species():\n", "\n", " # Plot of mineral's volume the space coordinates\n", " fig = plt.figure()\n", " ax = plt.axes(xlim=(-0.01, 0.501), ylim=(0.5e-5, 2))\n", " ax.set_xlabel('Distance [m]')\n", " ax.set_ylabel('Concentration [molal]')\n", " ax.set_yscale('log')\n", " ax.set_title(titlestr(0.0))\n", " objects = [\n", " ax.plot([], [], label=r'$\\mathrm{Ca^{2+}}$', **line('C0'))[0],\n", " ax.plot([], [], label=r'$\\mathrm{Mg^{2+}}$', **line('C1'))[0],\n", " ax.plot([], [], label=r'$\\mathrm{HCO_3^{-}}$',**line('C2'))[0],\n", " ax.plot([], [], label=r'$\\mathrm{CO_2(aq)}$',**line('red'))[0],\n", " ax.plot([], [], label=r'$\\mathrm{H^+}$', **line('darkviolet'))[0],\n", " ]\n", " ax.legend(loc='upper right')\n", "\n", " def init():\n", " return tuple(objects)\n", "\n", " def animate(i):\n", " t = i * dt\n", " filearray = np.loadtxt(folder_results + '/' + files[i], skiprows=1)\n", " data = filearray.T\n", "\n", " data_cacation = data[indx_Cacation]\n", " data_mgcation = data[indx_Mgcation]\n", " data_hco3anion = data[indx_HCO3anion]\n", " data_co2aq = data[indx_CO2aq]\n", " data_hcation = data[indx_Hcation]\n", "\n", " objects[0].set_data(xcells, data_cacation)\n", " objects[1].set_data(xcells, data_mgcation)\n", " objects[2].set_data(xcells, data_hco3anion)\n", " objects[3].set_data(xcells, data_co2aq)\n", " objects[4].set_data(xcells, data_hcation)\n", " ax.set_title(titlestr(t))\n", " return tuple(objects)\n", "\n", " print(\"Generating the animation of aqueous species behaviour ...\")\n", " anim = animation.FuncAnimation(fig, animate, init_func=init, frames=animation_frame_range, interval=animation_interval_wait, blit=True)\n", " anim.save(folder_videos + '/aqueous-species.mp4', fps=animation_fps, extra_args=['-vcodec', 'libx264'])\n", " print(\"Finished!\")" ] }, { "cell_type": "markdown", "id": "6fd55204", "metadata": {}, "source": [ "Generate animation with the ph behaviour:" ] }, { "cell_type": "code", "execution_count": null, "id": "b7fbfbae", "metadata": {}, "outputs": [], "source": [ "plot_animation_ph()" ] }, { "cell_type": "markdown", "id": "ae7832a7", "metadata": {}, "source": [ "Show the resulting video:" ] }, { "cell_type": "code", "execution_count": null, "id": "00aa6c06", "metadata": {}, "outputs": [], "source": [ "Video(folder_videos + '/pH.mp4')" ] }, { "cell_type": "markdown", "id": "20d5ab3e", "metadata": {}, "source": [ "Generate animation with calcite and dolomite dynamics:" ] }, { "cell_type": "code", "execution_count": null, "id": "fcc85034", "metadata": {}, "outputs": [], "source": [ "plot_animation_calcite_dolomite()" ] }, { "cell_type": "markdown", "id": "94df5b12", "metadata": {}, "source": [ "Show the video with precipitating dolomite and dissolving calcite:" ] }, { "cell_type": "code", "execution_count": null, "id": "0e12a40d", "metadata": {}, "outputs": [], "source": [ "Video(folder_videos + '/calcite-dolomite.mp4')" ] }, { "cell_type": "markdown", "id": "b39c2350", "metadata": {}, "source": [ "Generate an animation with aqueous species:" ] }, { "cell_type": "code", "execution_count": null, "id": "31a3d207", "metadata": {}, "outputs": [], "source": [ "plot_animation_aqueous_species()" ] }, { "cell_type": "markdown", "id": "6181202f", "metadata": {}, "source": [ "Show corresponding video:" ] }, { "cell_type": "code", "execution_count": null, "id": "8b0fd001", "metadata": {}, "outputs": [], "source": [ "Video(folder_videos + '/aqueous-species.mp4')" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" } }, "nbformat": 4, "nbformat_minor": 5 }