<commit_before><commit_msg>problem with dispersion example [ci skip]<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Predicting dispersion coefficient\n",
    "This example will show how to find the dispersion coefficient of a porous medium."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:28.086712Z",
     "iopub.status.busy": "2021-06-24T11:30:28.084191Z",
     "iopub.status.idle": "2021-06-24T11:30:28.810748Z",
     "shell.execute_reply": "2021-06-24T11:30:28.809367Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import openpnm as op\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import special\n",
    "from scipy.optimize import curve_fit\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_formats = ['svg']\n",
    "np.random.seed(10)\n",
    "np.set_printoptions(precision=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating Network\n",
    "A 2D 40 X 40 ``Cubic`` network is generated with a spacing of $10^{-4}$m, but a 3D network would work as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:28.818282Z",
     "iopub.status.busy": "2021-06-24T11:30:28.816859Z",
     "iopub.status.idle": "2021-06-24T11:30:29.000881Z",
     "shell.execute_reply": "2021-06-24T11:30:28.999636Z"
    }
   },
   "outputs": [],
   "source": [
    "shape = [40, 40, 1]\n",
    "pn = op.network.Cubic(shape=shape, spacing=1e-4)\n",
    "pn.add_model_collection(op.models.collections.geometry.spheres_and_cylinders)\n",
    "pn.regenerate_models()\n",
    "water = op.phase.Water(network=pn)\n",
    "water.add_model_collection(op.models.collections.physics.basic)\n",
    "water.regenerate_models()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining Effective Pore Volume\n",
    "\n",
    "The accumulation of mass in the network occurs only in the pores, where the concentration is solved.  In order for mass  accumulate properly, it is necessary to assign the throat volumes to their surrounded throats.  This creates an effective pore volume.  We can define this in a custom pore-scale model, making use of the ``numpy.add.at`` function, to add 1/2 the volume of each throat to its neighboring pores.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.009162Z",
     "iopub.status.busy": "2021-06-24T11:30:29.007566Z",
     "iopub.status.idle": "2021-06-24T11:30:29.011686Z",
     "shell.execute_reply": "2021-06-24T11:30:29.010456Z"
    }
   },
   "outputs": [],
   "source": [
    "def effective_pore_volume(target, throat_volume='throat.volume', pore_volume='pore.volume'):\n",
    "    Pvol = pn['pore.volume']\n",
    "    Tvol = pn['throat.volume']\n",
    "    Vtot = Pvol.sum() + Tvol.sum()\n",
    "    np.add.at(Pvol, pn.conns[:, 0], pn['throat.volume']/2)\n",
    "    np.add.at(Pvol, pn.conns[:, 1], pn['throat.volume']/2)\n",
    "    assert np.isclose(Pvol.sum(), Vtot)  # Ensure total volume has been added to Pvol\n",
    "    return Pvol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.018563Z",
     "iopub.status.busy": "2021-06-24T11:30:29.017158Z",
     "iopub.status.idle": "2021-06-24T11:30:29.022326Z",
     "shell.execute_reply": "2021-06-24T11:30:29.021036Z"
    }
   },
   "outputs": [],
   "source": [
    "pn.add_model(propname='pore.effective_volume', model=effective_pore_volume)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform Stokes flow\n",
    "The advection diffusion algorithm assumes a velocity field. Therefore, Stokes flow in the pore netwok is solved. The ``StokesFlow`` algorthm is solved prior to running the ``AdvectionDiffusion`` algorthim. For more information there is a seperate tutorial on Stokes Flow."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.031313Z",
     "iopub.status.busy": "2021-06-24T11:30:29.030041Z",
     "iopub.status.idle": "2021-06-24T11:30:29.128824Z",
     "shell.execute_reply": "2021-06-24T11:30:29.128178Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                                                                                                                       \r"
     ]
    }
   ],
   "source": [
    "sf = op.algorithms.StokesFlow(network=pn, phase=water)\n",
    "sf.set_value_BC(pores=pn.pores('back'), values=50.0)\n",
    "sf.set_value_BC(pores=pn.pores('front'), values=0)\n",
    "sf.run();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results obtained from the StokesFlow algorthim must be attached to the water phase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.132649Z",
     "iopub.status.busy": "2021-06-24T11:30:29.131998Z",
     "iopub.status.idle": "2021-06-24T11:30:29.134363Z",
     "shell.execute_reply": "2021-06-24T11:30:29.133789Z"
    }
   },
   "outputs": [],
   "source": [
    "water.update({'pore.pressure':sf['pore.pressure']})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add Diffusive Conductance Model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.139009Z",
     "iopub.status.busy": "2021-06-24T11:30:29.138420Z",
     "iopub.status.idle": "2021-06-24T11:30:29.143880Z",
     "shell.execute_reply": "2021-06-24T11:30:29.144352Z"
    }
   },
   "outputs": [],
   "source": [
    "mod = op.models.physics.ad_dif_conductance.ad_dif\n",
    "water.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Transient Advection Diffusion\n",
    "An algorthim for transient advection diffusion is defined here. It is assigned to the network and the phase, and will be able to retrieve all information that will be needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.149060Z",
     "iopub.status.busy": "2021-06-24T11:30:29.148495Z",
     "iopub.status.idle": "2021-06-24T11:30:29.150670Z",
     "shell.execute_reply": "2021-06-24T11:30:29.151146Z"
    }
   },
   "outputs": [],
   "source": [
    "ad = op.algorithms.TransientAdvectionDiffusion(network=pn, phase=water) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Dirichlet boundary conditions and the inital conditions are next defined as follows. If the inital condition is not defined then it is assumed to be zero, so it is redundunt in this case. The boundary conditions can be defined as ``value``, ``outflow``, and ``rate``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.155386Z",
     "iopub.status.busy": "2021-06-24T11:30:29.154800Z",
     "iopub.status.idle": "2021-06-24T11:30:29.158625Z",
     "shell.execute_reply": "2021-06-24T11:30:29.159095Z"
    }
   },
   "outputs": [],
   "source": [
    "inlet  = pn.pores('back') \n",
    "outlet = pn.pores('front')\n",
    "ad.set_value_BC(pores=inlet, values=1.0)\n",
    "ad.set_outflow_BC(pores=outlet)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setup the Transient Algorithim\n",
    "The settings of the transient algorthim can be updated here. We first define the time span:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.162863Z",
     "iopub.status.busy": "2021-06-24T11:30:29.162300Z",
     "iopub.status.idle": "2021-06-24T11:30:29.164061Z",
     "shell.execute_reply": "2021-06-24T11:30:29.164524Z"
    }
   },
   "outputs": [],
   "source": [
    "tspan = (0, 100)\n",
    "saveat = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We must also tell the algorithm to use the effective pore volume rather than the default which is just 'pore.volume'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.168035Z",
     "iopub.status.busy": "2021-06-24T11:30:29.167478Z",
     "iopub.status.idle": "2021-06-24T11:30:29.169248Z",
     "shell.execute_reply": "2021-06-24T11:30:29.169705Z"
    }
   },
   "outputs": [],
   "source": [
    "ad.settings['pore_volume'] = 'pore.effective_volume'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "openpnm.algorithms.TransientAdvectionDiffusion : trans_ad_01\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "#     Properties                                    Valid Values\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "1     pore.bc.outflow                                  40 / 1600 \n",
      "2     pore.bc.rate                                      0 / 1600 \n",
      "3     pore.bc.value                                    40 / 1600 \n",
      "4     pore.concentration                             1600 / 1600 \n",
      "5     pore.ic                                        1600 / 1600 \n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "#     Labels                                        Assigned Locations\n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n",
      "1     pore.all                                      1600      \n",
      "2     throat.all                                    3120      \n",
      "――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――――\n"
     ]
    }
   ],
   "source": [
    "print(ad)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The algorthim than can be run, but we must pass the initial conditions (could be a scalar or an array), time span, and optionally the intervals at which the solution is desired to be stored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:29.173152Z",
     "iopub.status.busy": "2021-06-24T11:30:29.172593Z",
     "iopub.status.idle": "2021-06-24T11:30:29.348307Z",
     "shell.execute_reply": "2021-06-24T11:30:29.348820Z"
    }
   },
   "outputs": [],
   "source": [
    "soln = ad.run(x0=0, tspan=tspan, saveat=saveat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving for the Predicti\n",
    "\n",
    "The following equation given by Fried (1971) is used to solve the longitudinal dispersion coefficient:\n",
    "\n",
    "$$\\frac{C}{C_{0}} = \\frac{1}{2}erfc\\Bigl(\\frac{x-Ut}{2(D_{L}t)^{\\frac{1}{2}}}\\Bigr)+\\frac{1}{2}exp\\Bigl(\\frac{Ux}{D_{L}}\\Bigl)erfc\\Bigr(\\frac{x+Ut}{2(D_{L}t)^{\\frac{1}{2}}}\\Bigr)$$\n",
    "\n",
    "Where $x$ is the length between the inlet and the outlet, $t$ is the time, $D_{L}$ is the longitudinal dispersion coefficient, $U$ is the average pore velocity, $C_{0}$ is the inlet concentration, and $C$ is the concentration at the given time. Since we defined the inlet concentration as being equal to 1, solving for $C$ is effictivly equal to solving for $\\frac{C}{C_{0}}$. ``erfc`` is the complementary error function, which is imported from ``scipy``.\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:30.180580Z",
     "iopub.status.busy": "2021-06-24T11:30:30.179193Z",
     "iopub.status.idle": "2021-06-24T11:30:30.183072Z",
     "shell.execute_reply": "2021-06-24T11:30:30.184182Z"
    }
   },
   "outputs": [],
   "source": [
    "q_throat = sf.rate(throats=pn.Ts, mode='single')\n",
    "A_throat = pn['throat.cross_sectional_area']\n",
    "v_throat = q_throat/A_throat\n",
    "v_pred = sum(q_throat*v_throat)/sum(q_throat)\n",
    "\n",
    "def elution(step,v,DL):\n",
    "    x = 40*1e-4\n",
    "    el1 = 0.5*(special.erfc((x-step*v)/(2*(DL*step)**(1/2))))\n",
    "    el2 = 0.5*np.exp(v*x/DL)\n",
    "    el3 = special.erfc((x+step*v)/(2*(DL*step)**(1/2)))\n",
    "    return el1+el2*el3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:30.191876Z",
     "iopub.status.busy": "2021-06-24T11:30:30.190439Z",
     "iopub.status.idle": "2021-06-24T11:30:30.224446Z",
     "shell.execute_reply": "2021-06-24T11:30:30.225572Z"
    }
   },
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'SolutionContainer' object has no attribute 't'",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[1;32m<ipython-input-22-8bc900e84543>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m      1\u001b[0m \u001b[0mg\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mv_pred\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;36m1e-3\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m----> 2\u001b[1;33m \u001b[0mxdata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mfloat\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mx\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mx\u001b[0m \u001b[1;32min\u001b[0m \u001b[0msoln\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m      3\u001b[0m \u001b[0mydata\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mc_avg\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      4\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m      5\u001b[0m \u001b[0mpopt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpcov\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mcurve_fit\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0melution\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mxdata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mydata\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mp0\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mg\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
      "\u001b[1;31mAttributeError\u001b[0m: 'SolutionContainer' object has no attribute 't'"
     ]
    }
   ],
   "source": [
    "g = [v_pred, 1e-3]\n",
    "xdata = [float(x) for x in soln.t]\n",
    "ydata = c_avg\n",
    "\n",
    "popt, pcov = curve_fit(elution, xdata, ydata, p0=g)\n",
    "disp_coeff = popt[1]\n",
    "v_fit = popt[0]\n",
    "\n",
    "print('Dispersion Coefficient = ', \"{0:.4E}\".format(disp_coeff), ' m^2/s')\n",
    "print('v_pred = ', \"{0:.4E}\".format(v_pred), ' m/s')\n",
    "print('v_fitted = ', \"{0:.4E}\".format(v_fit), ' m/s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2021-06-24T11:30:30.250496Z",
     "iopub.status.busy": "2021-06-24T11:30:30.245161Z",
     "iopub.status.idle": "2021-06-24T11:30:30.351452Z",
     "shell.execute_reply": "2021-06-24T11:30:30.352594Z"
    }
   },
   "outputs": [],
   "source": [
    "el = np.zeros(len(ydata))\n",
    "for i in range(len(ydata)):\n",
    "    el[i] = elution(xdata[i], popt[0], popt[1])\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(xdata, ydata, label=\"simulation\")\n",
    "ax.plot(xdata, el, \".\", label=\"fitted\")\n",
    "ax.legend()\n",
    "ax.set_xlabel('time (s)')\n",
    "ax.set_ylabel('concentration');"
   ]
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
{'fraction_non_alphanumeric': 0.22107050831576483, 'fraction_numerical': 0.1158936519091122, 'mean_word_length': 2.8776112624886467, 'pattern_counts': {'":': 219, '<': 5, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6192862', 'n_tokens_mistral': 8521, 'n_tokens_neox': 7375, 'n_words': 1063}