{ "cells": [ { "cell_type": "markdown", "metadata": { "lines_to_next_cell": 0, "toc": true }, "source": [ "

Table of Contents

\n", "
" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "lines_to_next_cell": 0, "nbgrader": { "checksum": "cb877aca98898ed1cd0ab4964947a0a6", "grade": false, "grade_id": "cell-0a4af3607542b85a", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "# Introduction" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "lines_to_next_cell": 0, "nbgrader": { "checksum": "9982245b6ca545953821546f2f2d5969", "grade": false, "grade_id": "cell-6f2c0b04c4f4fdf6", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "Let us recall that $J_{EC}$ and $J_{WC}$ represent the flux of mass [M/T] towards the center. If these are positive, the mass in the center will increase. If there is a source of mass in the center, the mass will increase as well. Therefore, the steady-state equation writes:\n", "\n", "\\begin{equation}\n", "J_{EC}+J_{WC} + Q = 0\n", "\\end{equation}\n", "\n", "With Q being written in terms of mg/s (Q: [M/T]) . Therefore mass balance equation with a source term writes:\n", "\n", "\\begin{equation}\n", "\\left( -2 \\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} + \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} + \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) \\Delta y \\Delta z + q \\Delta x \\Delta y \\Delta z = 0\n", "\\end{equation}\n", "\n", "with $q$ being the volumetric source ([M/T/L$^3$])\n", "\n", "We can rearrange the latter equation as:\n", "\n", "\\begin{equation}\n", "\\left( 2\\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) \\Delta y \\Delta z = q \\Delta x \\Delta y \\Delta z\n", "\\end{equation}\n", "\n", "We can also divide everything by $\\Delta x \\Delta y \\Delta z$:\n", "\n", "\\begin{equation}\n", "\\frac{1}{\\Delta x} \\left( 2\\frac{D\\theta}{\\Delta x} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x} c_{\\left[i+1\\right]} \\right) = q\n", "\\end{equation}\n", "\n", "which gives\n", "\\begin{equation}\n", "2\\frac{D\\theta}{\\Delta x^2} c_{\\left[i\\right]} - \\frac{D\\theta}{\\Delta x^2}c_{\\left[i-1\\right]} - \\frac{D\\theta}{\\Delta x^2} c_{\\left[i+1\\right]} = q\n", "\\end{equation}" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "a6432e612b59976fd2986d4fd92d0d85", "grade": false, "grade_id": "cell-98a865dfd995366b", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "import matplotlib.animation as animation\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "from IPython.display import display\n", "from IPython.display import HTML\n", "from numpy.testing import assert_allclose\n", "\n", "# the next two imports are to show animations and video!" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "8999184f8bef2c50b25590e8f81933dd", "grade": false, "grade_id": "cell-bd0b49f6ae1acef3", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "# this function deals with harmonic averaging when diffusion is not the same everywhere.\n", "# It doesn't change anything when diffusion is homogeneous but \n", "# you can try to see how it affects the behavior.\n", "# To modify a read-only cell, create a new cell below and copy these lines\n", "# into that cell\n", "\n", "\n", "def avg(Di, Dj):\n", " \"\"\"\n", " Computes the harmonic average between two values Di and Dj\n", " Returns 0 if either of them is zero\n", " \"\"\"\n", " if (Di * Dj) == 0:\n", " return 0\n", " else:\n", " return 2 / (1 / Di + 1 / Dj)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "e75ad6e41c6e2024decba30ea27f7183", "grade": false, "grade_id": "cell-54d67b8f8547140a", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "def Build_1D_Inhomo_Matrix_Source(bc_left, bc_right, n, D, Width, poro, Q):\n", " \"\"\"\n", " Constructs a coefficient matrix and an array with varying diffusion coefficient and a source term\n", "\n", " Parameters\n", " ----------\n", " bc_left: (float)\n", " left boundary condition (mg/L)\n", " bc_right: (float)\n", " right boundary condition (mg/L)\n", " n: (int)\n", " number of cells\n", " D: (float vector)\n", " values of the diffusion coefficient ((dm)^2/day)\n", " Width: (float)\n", " Total phyiscal width of the domain (dm)\n", " poro: (float)\n", " porosity value\n", " Q: (float vector)\n", " volumetric source term (mg/L/day)\n", "\n", " Returns\n", " -------\n", "\n", " A: nxn float array\n", " b: vector of length n\n", "\n", " that solve the\n", " discretized 1D diffusion problem Ax = b\n", " \"\"\"\n", " Matrix = np.zeros((n, n))\n", " RHS = np.zeros(n)\n", " dx = Width / (n - 1)\n", " coef = poro / dx / dx\n", " for i in range(n):\n", " if i == 0:\n", " RHS[i] = bc_left\n", " Matrix[i, i] = 1\n", " elif i == n - 1:\n", " RHS[i] = bc_right\n", " Matrix[i, i] = 1\n", " else:\n", " RHS[i] = Q[i]\n", " East = coef * avg(D[i], D[i + 1])\n", " West = coef * avg(D[i], D[i - 1])\n", " Matrix[i, i] = East + West\n", " Matrix[i, i + 1] = -East\n", " Matrix[i, i - 1] = -West\n", " return Matrix, RHS" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Initial conditions" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "096a66315aa17e9064fcbe0285554d69", "grade": false, "grade_id": "cell-b41a29613812db59", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "We will use similar conditions than in the previous assignment. However, it is a good practise to use units which are representative of the current problem, otherwise we have to deal with very large or very small numbers.\n", "\n", "The diffusion coefficient of solutes in pure water is usually 2$\\times$10$^{-9}$ m$^2$/s. Concentrations are usually expressed in mg/L, corresponding to mg/dm$^3$. To avoid future unit problems, let us define\n", "\n", "- one day as the time unit\n", "- one dm as the length unit (leading to Liters being the adequate volume unit)\n", "- mg as the mass unit\n", "\n", "In these units, the diffusion is expressed in dm$^2$/day, the width in dm, and the rate in mg/L/day.\n", "We will here change the values of the parameters used previously.." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "lines_to_next_cell": 2, "nbgrader": { "checksum": "37fb27d49858b607d7f3873589359181", "grade": false, "grade_id": "cell-b9bd7eaf90619a64", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "c_left = 1 # mg/L\n", "c_right = 0 # mg/L\n", "n = 51\n", "Diff = 2e-9 * 100 * 24 * 3600 # dm²/day\n", "D = Diff * np.ones(n)\n", "Q = np.zeros(n)\n", "Q0 = 5e-2 # mg/L/day\n", "Q[int(n / 4) : int(n / 2)] = Q0 # mg/L/day\n", "Width = 2 # dm\n", "poro = 0.4\n", "nTstp = 201\n", "dt = 0.5 # days\n", "\n", "c_init = np.zeros((nTstp, n))\n", "c_init[:, 0] = 1 # Boundary condition" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# If you want to change the parameters value, you are welcome to do so here.\n", "# But know that your results will be tested based on the parameters given above\n", "# We encourage you to try different things\n", "# For example, you change the diffusion at certain places, change the source term, ...\n", "# Trying out these things will help you build physical intuition\n", "# So we leave this cell for you to change these parameters.\n", "# As before, make a writeable copy of that cell to make changes\n", "#\n" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "d92264a16faddb9384a8da5039a65dad", "grade": false, "grade_id": "cell-e983be851ed2d326", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "## Save these into a dict for later use" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "lines_to_next_cell": 2, "nbgrader": { "checksum": "8f998c70ca8002d9e83a48a9427a51fd", "grade": false, "grade_id": "cell-b19ecc741595f838", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "init_dict = dict(\n", " c_left=c_left,\n", " c_right=c_right,\n", " c_init=c_init,\n", " n=n,\n", " Diff=Diff,\n", " D=D,\n", " Q=Q,\n", " Width=Width,\n", " poro=poro,\n", " nTstp=nTstp,\n", " dt=dt,\n", ")" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "cc0f2bb121edca31220a3bd57a32f786", "grade": false, "grade_id": "cell-8ce4321aef65caa0", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "# Steady state" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "82374c24f37ea28b1abfc1053419d339", "grade": false, "grade_id": "cell-1a17c575606df6ff", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "x = np.linspace(0, Width, n)\n", "A, b = Build_1D_Inhomo_Matrix_Source(c_left, c_right, n, D, Width, poro, Q)\n", "\n", "c_final = np.linalg.solve(A, b)\n", "plt.plot(x, c_final, label=\"Concentration\")\n", "plt.xlabel(\"x-axis (dm)\")\n", "plt.ylabel(\"Concentration (mg/L)\")" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "2950f08533e97e81d43a296e14fd3f80", "grade": false, "grade_id": "cell-ec63e59809fdbde1", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "# Transient behavior" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "81d334c918425615e048ef7f40c71739", "grade": false, "grade_id": "cell-7adbd2fb47c2104c", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "The transient behavior that we have discussed in the previous notebook simply adds these terms to the previous equation.\n", "\n", "\\begin{equation}\n", "J_{EC}+J_{WC} + Q = \\frac{dm}{dt}\n", "\\end{equation}\n", "\n", "Using the same developements expressed before, using\n", "\n", "\\begin{equation}\n", "\\frac{d\\theta c}{dt} = \\theta \\frac{c(t+\\Delta t) - c(t)}{\\Delta t}\n", "\\end{equation}\n", "\n", "In the following equation, we will write $c(t+\\Delta t)$ as $c_i$ (i being the subscript referring to the position index), while $c(t)$ is the old concentration at the same point $c_{i}(\\text{old})$.\n", "\n", "\\begin{equation}\n", "\\frac{\\theta c_{i}}{\\Delta t} + 2\\frac{D\\theta}{\\Delta x^2} c_{i} - \\frac{D\\theta}{\\Delta x^2}c_{i-1} - \\frac{D\\theta}{\\Delta x^2} c_{i+1} = q + \\frac{\\theta c_{i}(\\text{old})}{\\Delta t}\n", "\\end{equation}\n", "\n", "So we can solve this problem by using the same matrixes used before, but we need to add extra terms.\n", "\n", "Based on the previous assignment, most of these terms are already incorporated in the matrix A and b. There are only a couple of terms to add to be able to solve the transient problem. We will keep the structure of matrix A and b defined.\n", "\n", "You have got to modify the next cell (keep its original structure), and incorporate these new terms in the defined arrays Abis and Bbis." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "checksum": "9df7c44e742f7d7614bf0c88b7c4ff05", "grade": true, "grade_id": "cell-89ed183cacbcbc6d", "locked": false, "points": 5, "schema_version": 1, "solution": true } }, "outputs": [], "source": [ "# HERE YOU HAVE TO ADD TWO LINES OF CODE (see commented lines). \n", "# WE RECOMMEND THAT YOU DO NOT CHANGE ANYTHING ELSE\n", "# Refer to the equations developed before!\n", "\n", "\n", "def calc_conc(\n", " Diff=None,\n", " D=None,\n", " Width=None,\n", " c_left=None,\n", " c_right=None,\n", " poro=None,\n", " Q=None,\n", " c_init=None,\n", " n=None,\n", " nTstp=None,\n", " dt=None,\n", "):\n", "\n", " x = np.linspace(0, Width, n)\n", " A, b = Build_1D_Inhomo_Matrix_Source(c_left, c_right, n, D, Width, poro, Q)\n", "\n", " c = c_init\n", " Abis = np.zeros((n, n))\n", " Bbis = np.zeros(n)\n", "\n", " for t in range(nTstp - 1):\n", " for i in range(n):\n", " Abis[i, i] = 0 # Here is where you need to change!\n", " Bbis[i] = 0 # Here is where you need to change!\n", " Aa = A + Abis\n", " bb = b + Bbis\n", " c[t + 1, :] = np.linalg.solve(Aa, bb)\n", "\n", " return x, c\n", "\n", "\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "30b9027224a41ad040a3bde5a2ae9157", "grade": false, "grade_id": "cell-41e25b3ba511d765", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "x, c = calc_conc(**init_dict)\n", "\n", "plt.plot(x, c_final, label=\"Asymptotic concentration\")\n", "plt.plot(x, c[0, :], label=\"Initial concentration\")\n", "plt.plot(x, c[10, :], label=\"Concentration after 10 timestep\")\n", "plt.plot(x, c[30, :], label=\"Concentration after 30 timesteps\")\n", "plt.plot(x, c[60, :], label=\"Concentration after 60 timesteps\")\n", "plt.plot(x, c[nTstp - 1, :], label=\"Concentration after 200 timesteps\")\n", "plt.xlabel(\"x-axis (dm)\")\n", "plt.ylabel(\"Concentration (mg/L)\")\n", "plt.legend(bbox_to_anchor=(1.01, 1), loc=\"upper left\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The next cell tests one of the cell concentrations" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "a49f72af411ce3dcdff3e79c4a8edb14", "grade": true, "grade_id": "cell-e3dd3d495d9a1bd8", "locked": true, "points": 3, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "# Hidden test for a cell concentration (our choice of cell)" ] }, { "cell_type": "markdown", "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "518a1b20f861210ef4954707e896bec9", "grade": false, "grade_id": "cell-faad482a85fd293b", "locked": true, "schema_version": 1, "solution": false } }, "source": [ "What is the concentration at 1.6 dm from the left boundary condition after 60 days?\n", "Put this value in the variable answer in the next cell.\n", "That is, replace the line\n", "answer = -1\n", "with your answer" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "nbgrader": { "checksum": "20f7bbef0d84da78a988d17477bea750", "grade": false, "grade_id": "cell-b73f549641823e89", "locked": false, "schema_version": 1, "solution": true } }, "outputs": [], "source": [ "answer = -1\n", "# YOUR CODE HERE\n", "raise NotImplementedError()" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "b3f27d2e979d2a19da24eaded5887467", "grade": true, "grade_id": "cell-fecfbea162708f5a", "locked": true, "points": 2, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "# Hidden test for answer" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "lines_to_next_cell": 2, "nbgrader": { "checksum": "3278b159aabf8ae1c8a4f28690b744b4", "grade": false, "grade_id": "cell-22d5514ace8b854a", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "\n", "\n", "def make_animation(init_dict):\n", " image_list = []\n", " x, c = calc_conc(**init_dict)\n", " fig, ax = plt.subplots()\n", " ax.set_xlim((0, 2.0))\n", " ax.set_ylim((0, 2.5))\n", " nsteps, nvals = c.shape\n", " for index in range(nsteps):\n", " line = ax.plot(x, c[index, :], \"r-\", lw=2)\n", " image_list.append(line)\n", " return fig, image_list\n", "\n", "\n", "fig, image_list = make_animation(init_dict)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# If you want to see the animation, just change movie to True\n", "# (might not work depending on your config and browser, working\n", "# for chrome on macos\n", "\n", "movie = False" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "5798ee1037f5f7e6654ad3530575498a", "grade": false, "grade_id": "cell-e09cb384b6b1f9c5", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "if movie:\n", " anim = animation.ArtistAnimation(\n", " fig, image_list, interval=50, blit=True, repeat_delay=1000\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "deletable": false, "editable": false, "nbgrader": { "checksum": "1d8f2704db5c33523c6ce4f8145f7ba8", "grade": false, "grade_id": "cell-61cb781fde1572ce", "locked": true, "schema_version": 1, "solution": false } }, "outputs": [], "source": [ "if movie:\n", " out = HTML(anim.to_html5_video())\n", " display(out)" ] } ], "metadata": { "jupytext": { "cell_metadata_filter": "all", "formats": "ipynb", "notebook_metadata_filter": "all", "text_representation": { "extension": ".py", "format_name": "percent", "format_version": "1.2", "jupytext_version": "1.0.0-rc5" } }, "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.7.3" }, "nbsphinx": { "execute": "never" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }