{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "349e87a2",
   "metadata": {},
   "source": [
    "# Using class EquilibriumSolver\n",
    "\n",
    "This tutorial demonstrates the use of class\n",
    "[EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html)\n",
    "for repeated equilibrium calculations.\n",
    "\n",
    "For this, we model the dissolution of mineral calcite (CaCO<sub>3</sub>) in a\n",
    "saline aqueous solution as we increase the amount of CO<sub>2</sub> into the system.\n",
    "\n",
    "> **Note:** Prefer class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html)\n",
    "instead of method [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n",
    "when many chemical equilibrium calculations are needed.\n",
    "\n",
    "We follow the usual initialization steps, such as importing the `reaktoro`\n",
    "package,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fd29c0df",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reaktoro import *\n",
    "import numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "097e80c9",
   "metadata": {},
   "source": [
    "creating an object of class\n",
    "[Database](https://reaktoro.org/cpp/classReaktoro_1_1Database.html),"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc94a013",
   "metadata": {},
   "outputs": [],
   "source": [
    "db = Database(\"supcrt98.xml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "445440f8",
   "metadata": {},
   "source": [
    "and using this database object to initialize the object of class\n",
    "[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html)\n",
    "to specify the phases of interest in the chemical system:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7505e70",
   "metadata": {},
   "outputs": [],
   "source": [
    "editor = ChemicalEditor(db)\n",
    "editor.addAqueousPhaseWithElements(\"H O Na Cl C Ca\")\n",
    "editor.addGaseousPhase([\"H2O(g)\", \"CO2(g)\", \"O2(g)\"])\n",
    "editor.addMineralPhase(\"Halite\")\n",
    "editor.addMineralPhase(\"Calcite\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7cd34699",
   "metadata": {},
   "source": [
    "We now use the chemical system configuration stored in `editor` to create an\n",
    "instance of the class\n",
    "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html),\n",
    "which will be constructed containing four phases: an aqueous phase, a gaseous\n",
    "phase, and two pure mineral phases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0c45dc9a",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = ChemicalSystem(editor)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "313b0616",
   "metadata": {},
   "source": [
    "We want to perform a sequence of equilibrium calculations. For this, we create\n",
    "first a numpy array containing evenly spaced numbers ranging from 0 to 4 (a\n",
    "total of 101 values). This array represents the amount of CO<sub>2</sub> in\n",
    "each equilibrium problem we shall solve later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d61aaf96",
   "metadata": {},
   "outputs": [],
   "source": [
    "co2_amounts = numpy.linspace(0.0, 4.0, num=101)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e074374b",
   "metadata": {},
   "source": [
    "We now create the list of equilibrium problems stored in the list `problems`, one for each `amount` in array\n",
    "`x`. In each problem, we have 1 mol of calcite, 1 kg of H<sub>2</sub>O, 1 mol\n",
    "of NaCl, and `amount` of CO<sub>2</sub>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4331d3f6",
   "metadata": {},
   "outputs": [],
   "source": [
    "problems = []\n",
    "\n",
    "for amount in co2_amounts:\n",
    "    problem = EquilibriumProblem(system)\n",
    "    problem.setTemperature(70.0, \"celsius\")\n",
    "    problem.setPressure(100.0, \"bar\")\n",
    "    problem.add(\"Calcite\", 1.0, \"mol\")\n",
    "    problem.add(\"H2O\", 1.0, \"kg\")\n",
    "    problem.add(\"NaCl\", 1.0, \"mol\")\n",
    "    problem.add(\"CO2\", amount, \"mol\")\n",
    "\n",
    "    problems.append(problem)  # append the new problem into problems"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "810477ea",
   "metadata": {},
   "source": [
    "Solving this sequence of equilibrium problems can be done with a single\n",
    "instance of the class\n",
    "[EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9422b3e7",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = EquilibriumSolver(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f206238",
   "metadata": {},
   "source": [
    "Each chemical equilibrium calculation results in an individual chemical\n",
    "state. Below, we create a list of\n",
    "[ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html)\n",
    "objects that will contain the solution to each problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1a1eac90",
   "metadata": {},
   "outputs": [],
   "source": [
    "states = [ChemicalState(system) for _ in range(len(problems))]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0af2cf52",
   "metadata": {},
   "source": [
    "We now solve the sequence of chemical equilibrium problems:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d981a33",
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(len(problems)):\n",
    "    solver.solve(states[i], problems[i])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98281971",
   "metadata": {},
   "source": [
    "We want to plot the amount of calcite versus the amount of CO<sub>2</sub>\n",
    "added into the system. We collect below the amount of calcite in each computed\n",
    "equilibrium state in a list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7527264",
   "metadata": {},
   "outputs": [],
   "source": [
    "calcite_amounts = [state.speciesAmount(\"Calcite\") for state in states]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2705884",
   "metadata": {},
   "source": [
    "We now use [bokeh](https://docs.bokeh.org/en/latest/docs/gallery.html#standalone-examples)\n",
    "to do the plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0ed0bc7",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "from bokeh.plotting import figure, show\n",
    "from bokeh.io import output_notebook\n",
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "827b53df",
   "metadata": {},
   "source": [
    "Besides, we define a custom function that would generate figure of a size 600 x 300 with a label\n",
    "`Amount of CO2 [mol]` on the x-axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f82f4b59",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def custom_figure(title, y_axis_label): return figure(plot_width=600, plot_height=300,\n",
    "                                                      title=title,\n",
    "                                                      x_axis_label='Amount of CO2 added [mol]',\n",
    "                                                      y_axis_label=y_axis_label)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19188da2",
   "metadata": {},
   "source": [
    "Create a new plot with the given title and y-axis label, where the amount of calcite in `calcite_amounts` versus\n",
    "the amount of CO<sub>2</sub> in `co2_amounts`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7a3cf9f",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig1 = custom_figure(title=\"Calcite dissolution\", y_axis_label='Amount of Calcite [mol]')\n",
    "fig1.line(co2_amounts, calcite_amounts, line_width=4)\n",
    "show(fig1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "687df78b",
   "metadata": {},
   "source": [
    "To understand why calcite does not dissolve further after\n",
    "about 0.76 mol of CO<sub>2</sub> is added into the system,\n",
    "we plot the amount of gaseous species CO<sub>2</sub>(g):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "88ac12d6",
   "metadata": {},
   "outputs": [],
   "source": [
    "co2gas_amount = [state.speciesAmount(\"CO2(g)\") for state in states]\n",
    "fig2 = custom_figure(title=\"CO2(g)\", y_axis_label='Amount of CO2(g) [mol]')\n",
    "fig2.line(co2_amounts, co2gas_amount, line_width=4)\n",
    "show(fig2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2339d0f",
   "metadata": {},
   "source": [
    "Note above that after about 0.76 mol of CO<sub>2</sub> is added, the aqueous\n",
    "solution becomes CO<sub>2</sub>-saturated, and the amount of\n",
    "CO<sub>2</sub>(g), initially zero, starts to increase."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}