{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hydrogeophysical modeling with [pyGIMLi](https://www.pygimli.org)\n",
    "Coupled hydrogeophysical modeling example [presented at AGU17](https://zenodo.org/record/1095621) in a Jupyter Notebook. \n",
    "(*Hint: Press `Shift-Enter` to execute a cell and go to the next.*)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Imports and plot settings"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import ImageGrid\n",
    "import seaborn\n",
    "seaborn.set(context=\"notebook\", style=\"white\", font_scale=1.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Model creation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygimli as pg\n",
    "import pygimli.meshtools as mt\n",
    "\n",
    "# Create geometry definition for the modeling domain\n",
    "world = mt.createWorld(start=[-20, 0], end=[20, -16], \n",
    "                       layers=[-2, -8], worldMarker=False)\n",
    "\n",
    "# Create a heterogeneous block\n",
    "block = mt.createRectangle(start=[-6, -3.5], end=[6, -6.0],\n",
    "                           marker=4,  boundaryMarker=10, area=0.1)\n",
    "\n",
    "# Merge geometrical entities\n",
    "geom = mt.mergePLC([world, block])\n",
    "pg.show(geom, boundaryMarker=True)\n",
    "\n",
    "# Create a mesh from the geometry definition\n",
    "mesh = mt.createMesh(geom, quality=32, area=0.2, smooth=[1, 10])\n",
    "\n",
    "pg.show(mesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Modeling groundwater flow"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fluid flow in a porous medium of slow non-viscous and non-frictional hydraulic movement is governed by Darcy's Law according to:\n",
    "\n",
    "\\begin{align}\n",
    "K^{-1}\\mathbf{v} + \\nabla p & = 0 \\\\\n",
    "\\nabla \\cdot \\mathbf{v} & = 0\\\\\n",
    "\\text{leading}\\,\\,\\text{to}\\,\\,\n",
    "\\nabla\\cdot(K \\nabla p) & = 0 \\quad \\text{on} \\quad\\Omega\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We begin by defining isotropic values of hydraulic conductivity $K$ and mapping these to each mesh cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Map regions to hydraulic conductivity in $m/s$\n",
    "kMap = [[1, 1e-8], [2, 5e-3], [3, 1e-4], [4, 8e-4]]\n",
    "\n",
    "# Map conductivity value per region to each cell in the given mesh\n",
    "K = pg.solver.parseMapToCellArray(kMap, mesh)\n",
    "\n",
    "pg.show(mesh, data=K, label='Hydraulic conductivity $K$ in m$/$s',\n",
    "        cMin=1e-5, cMax=1e-2, cmap='viridis', grid=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The problem further boundary conditions of the hydraulic potential. We use $p=p_0=0.75$ m on the left and $p=0$ on the right boundary of the modelling domain, equaling a hydraulic gradient of 1.75%."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Dirichlet conditions for hydraulic potential\n",
    "pBound = [[[1, 2, 3], 0.75], [[5, 6, 7], 0.0]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now call the finite element solver with the generated mesh, hydraulic conductivity and the boundary condition. The sought hydraulic velocity distribution can then be calculated as the\n",
    "gradient field of $\\mathbf{v}=-\\nabla p$ and visualized using the generic `pg.show()` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solve for hydraulic potential\n",
    "p = pg.solver.solveFiniteElements(mesh, a=K, uB=pBound)\n",
    "\n",
    "# Solve velocity as gradient of hydraulic potential\n",
    "vel = -pg.solver.grad(mesh, p) * np.asarray([K, K, K]).T\n",
    "\n",
    "ax, _ = pg.show(mesh, data=pg.abs(vel) * 1000, logScale=False, cMin=0.05, cMax=0.15,\n",
    "                label='Velocity $v$ in mm$/$s', cmap='YlGnBu', hold=True)\n",
    "ax, _ = pg.show(mesh, data=vel, ax=ax, color='k',  linewidth=0.8,\n",
    "                dropTol=1e-5, hold=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Solving transport"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next step, we use this velocity field to simulate the dynamic movement of a particle (e.g., salt) concentration $c(\\mathbf{r}, t)$ in the aquifer by using the advection-diffusion equation:\n",
    "\n",
    "$$\\frac{\\partial c}{\\partial t} = \\underbrace{\\nabla\\cdot(D \\nabla c)}_{\\text{Diffusion / Dispersion}} - \\underbrace{\\nabla \\cdot (\\mathbf{v}\\nabla c)}_{\\text{Advection}} + S$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pg.RVector(mesh.cellCount(), 0.0)\n",
    "# Fill injection source vector for a fixed injection position\n",
    "sourceCell = mesh.findCell([-19.1, -4.6])\n",
    "S[sourceCell.id()] = 1.0 / sourceCell.size()  # g/(l s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pg.show(mesh, S)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a time vector and common velocity-depending dispersion coefficient $D = \\alpha |\\mathbf{v}|$ with a dispersivity $\\alpha = 1\\cdot10^{-2}$ m. We solve the advection-diffsuion equation on the equation level with the finite volume solver, which results in a particle concentration $c(\\mathbf{r},t)$ (in g$/$l) for each cell center and time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Choose 800 time steps for 6 days in seconds\n",
    "t = pg.utils.grange(0, 6 * 24 * 3600, n=800)\n",
    "\n",
    "# Create dispersitivity, depending on the absolute velocity\n",
    "dispersion = pg.abs(vel) * 1e-2\n",
    "\n",
    "# Solve for injection time, but we need velocities on cell nodes\n",
    "veln = mt.cellDataToNodeData(mesh, vel)\n",
    "c1 = pg.solver.solveFiniteVolume(mesh, a=dispersion, f=S, vel=veln,\n",
    "                                 times=t, uB=[1, 0],\n",
    "                                 scheme='PS', verbose=0)\n",
    "\n",
    "# Solve without injection starting with last result\n",
    "c2 = pg.solver.solveFiniteVolume(mesh, a=dispersion, f=0, vel=veln,\n",
    "                                 u0=c1[-1], times=t, uB=[1, 0],\n",
    "                                 scheme='PS', verbose=0)\n",
    "# Stack results together\n",
    "c = np.vstack((c1, c2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now visualize the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import clear_output\n",
    "\n",
    "for ci in c[1:][::25]:\n",
    "    pg.show(mesh, data=ci*0.001, cMin=0, cMax=3, cmap=\"rocket_r\", label=\"Concentration c in $g/l$\")\n",
    "    clear_output(wait=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Simulate time-lapse electrical resistivity measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a dipole-dipole measurement scheme and a suitable mesh for ERT forward simulations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygimli.physics.ert as ert\n",
    "\n",
    "\n",
    "ertScheme = ert.createERTData(pg.utils.grange(-20, 20, dx=1.0),\n",
    "                              schemeName='dd')\n",
    "\n",
    "meshERT = mt.createParaMesh(ertScheme, quality=33, paraMaxCellSize=0.2,\n",
    "                            boundaryMaxCellSize=50, smooth=[1, 2])\n",
    "pg.show(meshERT)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use simulated concentrations to calculate bulk resistivity using Archie's Law and fill matrix with apparent resistivity ratios with respect to a background model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pygimli.physics.petro as petro\n",
    "\n",
    "# Select 10 time frame to simulate ERT data\n",
    "timesERT = pg.IndexArray(np.floor(np.linspace(0, len(c)-1, 10)))\n",
    "\n",
    "# Create conductivity of fluid for salt concentration $c$\n",
    "sigmaFluid = c[timesERT] * 0.1 + 0.01\n",
    "\n",
    "# Calculate bulk resistivity based on Archie's Law\n",
    "resBulk = petro.resistivityArchie(rFluid=1. / sigmaFluid,\n",
    "                                  porosity=0.3, m=1.3,\n",
    "                                  mesh=mesh, meshI=meshERT, fill=1)\n",
    "\n",
    "# apply background resistivity model\n",
    "rho0 = np.zeros(meshERT.cellCount()) + 1000.\n",
    "for c in meshERT.cells():\n",
    "    if c.center()[1] < -8:\n",
    "        rho0[c.id()] = 150.\n",
    "    elif c.center()[1] < -2:\n",
    "        rho0[c.id()] = 500.\n",
    "resis = pg.RMatrix(resBulk)\n",
    "for i, rbI in enumerate(resBulk):\n",
    "    resis[i] = 1. / ((1./rbI) + 1./rho0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initialize and call the ERT manager for electrical simulation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pygimli.physics.ert import ERTManager\n",
    "\n",
    "# Initialize ert method manager\n",
    "ERT = ERTManager(verbose=False)\n",
    "# Run  simulation for  the apparent resistivities\n",
    "rhoa = ERT.simulate(meshERT, resis, ertScheme, verbose=0, returnArray=True)\n",
    "\n",
    "for i in range(4):\n",
    "    ERT.showData(ertScheme, vals=rhoa[i]/rhoa[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The inversion part would exceed the computational resources available for this live demo, but the corresponding Python script can be downloaded [here](http://cg17.pygimli.org). See [pygimli.org](https://www.pygimli.org) for more examples and tutorials."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}