1d_diffusion.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1-D Advection–Diffusion with a Continuous Point Release\n",
    "==\n",
    "\n",
    "This notebook simulates the solution to the unsteady 1-D advection diffusion with a continuous point\n",
    "source at the origin, derives the analytical solution at steady state, and compares the two.\n",
    "\n",
    "The PDE is:\n",
    "\n",
    "$$\\frac{\\partial C}{\\partial t} + u\\frac{\\partial C}{\\partial x} = D \\frac{\\partial^2 C}{\\partial x^2} + \\dot{m}\\delta (x)$$\n",
    "\n",
    "and the domain extends infinitely far, $x\\in [-\\infty,\\infty]$. This is a 1-D problem, and we take the mass source $\\dot{m}$ to be per-unit area in the $y-z$ plane (*Aside: what are the units of the delta function $\\delta(x)$?*)\n",
    "\n",
    "\n",
    "**Shift-Enter** will step through the code block-by-block.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Some setup:\n",
    "# Libraries for computation\n",
    "import numpy as np\n",
    "from scipy import sparse\n",
    "from scipy.sparse.linalg import spsolve\n",
    "# Libraries for interactive plotting\n",
    "import matplotlib.pyplot as plt\n",
    "import ipywidgets as widgets\n",
    "%matplotlib notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%html\n",
    "<!-- fix icon for rewind button -->\n",
    "<style> i.fa-stop:before { content: \"❚◀\"; } </style>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Initial parameters for the simulation\n",
    "--\n",
    "Before tackling the steady-state analytical solution, define and run the simulation. The details of how the simulation are defined are covered in CE200B. For our purposes we can take the simulation as a black box (but do take a look at the code if you're curious).\n",
    "\n",
    "We'll simulate $C(x,t)$ over a subset of the $x$ axis, which we break up or discretize into many small control volumes (which is to say that the code below is a a *finite volume* method). (*Aside: while the original problem extends to $\\pm \\infty$, the simulation does not and we have to define boundary conditions. The code below defines the boundary conditions as $C=0$ on the upstream side, and a constant gradient on the downstream side)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the simulation parameters:\n",
    "N=1001 # number of control volumes\n",
    "unum=0.05 # advection velocity [m/s]\n",
    "\n",
    "# Parameters set below interactively:\n",
    "#  mdot: source strength ~ [g/(m^2*s)]\n",
    "#     D: diffusion coefficient, [m^2/s]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define simulation domain:\n",
    "x=np.linspace(-8,10,N) # centers of the control volumes\n",
    "dx=x[1]-x[0] # width of each control volume\n",
    "xedge=np.r_[x-dx/2,x[-1]+dx/2] # edges of the control volumes\n",
    "# Define initial condition:\n",
    "t=0.0 # time [s]\n",
    "Cnum=np.zeros_like(x) # concentration in each control volume"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# The main attraction: plot and define how C evolves with each timestep.\n",
    "\n",
    "# Plot C vs x\n",
    "fig,ax=plt.subplots(figsize=[10,4])\n",
    "l=ax.plot(x,Cnum,lw=2)[0] # keep a handle to the plot of concentration\n",
    "def t_str():\n",
    "    return f\"time={t:7.2f}s\"\n",
    "t_label=ax.text(0.02,0.95,t_str(),transform=ax.transAxes,va='top')\n",
    "ax.axis(xmin=x[0],xmax=x[-1],ymin=-0.5,ymax=100) # for default values [0,100] is a good y range.\n",
    "ax.set_xlabel('x [m]')\n",
    "ax.set_ylabel('C [g/m$^3$]')\n",
    "\n",
    "# ____BEGIN 200B CONTENT____\n",
    "# Implicit diffusion operator is mat_imp=I+D*dt/dx**2 * F\n",
    "I=sparse.eye(N)\n",
    "# F is the part of operator that scales with Dnum*dt/dx\n",
    "d0=2*np.ones(N)\n",
    "dp1=-np.ones(N)\n",
    "dp1[1]*=2\n",
    "dm1=-np.ones(N)\n",
    "# dm1[-2:]*=2 # zero diffusive flux on downstream BC\n",
    "dm1[-2:]=0\n",
    "d0[-1]=0\n",
    "F=sparse.spdiags( [dm1,d0,dp1], [-1,0,1], N, N).tocsr()\n",
    "\n",
    "last_step=0\n",
    "def plot_update(step,D,dt,src,u):\n",
    "    global Dnum,Cnum,t,mdot,unum,last_step\n",
    "    Dnum=D\n",
    "    mdot=src\n",
    "    unum=u\n",
    "    \n",
    "    if step==0: # reset to initial condition.\n",
    "        Cnum=np.zeros_like(x)\n",
    "        t=0.0\n",
    "    elif step>last_step:\n",
    "        # Step forward in time\n",
    "        t+=dt\n",
    "        # Advection:\n",
    "        Cnum=np.interp(x,x+dt*unum,Cnum,left=0,right=0)\n",
    "        # Source w/advection\n",
    "        # add Csrc=mdot/u for the interval [0,u*dt]\n",
    "        # integrate to get the cumulative distribution, then interpolate\n",
    "        # and diff to get discrete distribution:\n",
    "        # for x<0, cdf is 0.\n",
    "        # over [0,u*dt] this should ramp up to mdot*dt \n",
    "        # Make an x coord array for the edge of the control volumes\n",
    "        # The cumulative mass over [-inf,x_i] for each x_i in xedge,\n",
    "        # for mass released over [0,dt]\n",
    "        cmdf_x=np.interp( xedge, [0,unum*dt], [0,mdot*dt])\n",
    "        mdf_x=np.diff(cmdf_x) # mass per CV\n",
    "        Cnum[:]+=mdf_x/dx # concentratino per CV\n",
    "\n",
    "        # Diffusion:\n",
    "        mat_imp = I+Dnum*dt/dx**2*F\n",
    "        Cnum=sparse.linalg.spsolve(mat_imp.tocsr(),Cnum)\n",
    "    # Update the plot\n",
    "    l.set_ydata(Cnum)\n",
    "    t_label.set_text(t_str())\n",
    "    last_step=step\n",
    "    \n",
    "# ____END 200B content____\n",
    "\n",
    "# Make the plot interactive\n",
    "widgets.interact(plot_update,\n",
    "                 step=widgets.Play(value=0,\n",
    "                                   min=0, max=1000000,\n",
    "                                   step=1, interval=20,\n",
    "                                   show_repeat=False),\n",
    "                 D=widgets.FloatSlider(value=0.05,min=0,max=1.0,step=0.001,\n",
    "                                        description='D:',readout_format='.3f'),\n",
    "                 src=widgets.FloatSlider(value=4.0,min=0,max=10.0, step=0.01,description='mdot'),\n",
    "                 dt=widgets.FloatSlider(value=0.4,min=0.01,max=4.0,step=0.01,\n",
    "                                        description='timestep:', readout_format='.2f'),\n",
    "                 u=widgets.FloatSlider(value=0.05,min=0.00,max=0.5,step=0.01,\n",
    "                                        description='u:', readout_format='.2f')\n",
    "                )\n",
    ";"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analytical solution at steady state\n",
    "--\n",
    "\n",
    "The simulation may or may not be correct!  Solve for the steady state $C(x)$ and compare to the simulation.\n",
    "\n",
    "The equation we're trying to solve is:\n",
    "$$ \\frac{\\partial C}{\\partial t} + u\\frac{\\partial C}{\\partial x} = D \\frac{\\partial^2 C}{\\partial x^2} + \\dot{m}\\delta(x)$$\n",
    "\n",
    "The Peclet number tells us that as we get far from the source, advection dominates (taking $L\\sim|x|$, recall $Pe\\equiv \\frac{UL}{D}$, and $Pe \\gg 1$ implies advection-dominated). In our infinite domain, we can evaluate the behavior as $x$ goes to $\\pm \\infty$ by considering only advection.\n",
    "\n",
    " * For $x\\rightarrow -\\infty$: mass from the source is introduced downstream and flowing in the opposite direction, so $C_{-\\infty}=0$\n",
    " * For $x\\rightarrow \\infty$: the advective flux is $J_A=uC$ and at steady state this must match the mass introduce by the source, \n",
    "   $\\dot{m}$. Thus the concentration far downstream is $C_{+\\infty}=\\dot{m}/u$.\n",
    "\n",
    "We now have suitable boundary conditions. What happens near $x=0$? The point source at the origin introduces a delta function $\\dot{m} \\delta(x)$ into the equation. Drop the unsteadiness term and integrate. Note that $\\int \\delta(x) = H(x)$, the Heaviside step function.\n",
    "\n",
    "$$ uC = D \\frac{dC}{dx} + \\dot{m} H(x) + A $$   \n",
    "\n",
    "with a constant of integration $A$.\n",
    "\n",
    "For $x>0$,\n",
    "\n",
    "$$uC -D \\frac{dC}{dx} = \\dot{m} + A$$\n",
    "\n",
    "Note that the left side is the net flux, with advective and diffusive contributions. All mass introduced at the sources eventually advects towards $+\\infty$, such that the net flux must be $J_{net}|_{x>0} = \\dot{m}$ (remember $\\dot{m}$ is for a 1-D source, per unit area in the $y-z$ plane, same units as $J$), and the constant of integration $A=0$. \n",
    "\n",
    "This is a *nonhomogeneous* ODE, and has a solution:\n",
    "\n",
    "$$ C|_{x>0} = B e^\\frac{ux}{D} + \\dot{m}/u $$\n",
    "\n",
    "Where $B$ is from another constant of integration. Knowing that $C(x)$ goes to $\\dot{m}/u$ as $x\\rightarrow\\infty$, $B=0$, $C=\\dot{m}/u$ for $x>0$.\n",
    "\n",
    "For $x<0$, the step function $H(x<0)=0$, and we can solve to get:\n",
    "\n",
    "$$C|_{x<0}=C|_0 e^{ux/D}$$\n",
    "\n",
    "For $D>0$, at steady state we require that the diffusive flux is finite, and specifically $\\partial C/\\partial x$ is finite. Thus there cannot be a discontinuity across $x=0$, and $C|_0=\\dot{m}/u$. The final solution is then:\n",
    "\n",
    "$$C(x) = \n",
    "\\left\\{\\begin{array}{ll}\n",
    "\\frac{\\dot{m}}{u} e^{ux/D}, & x<0 \\\\\n",
    "\\frac{\\dot{m}}{u}, & x \\geq 0 \\\\\n",
    "\\end{array}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "Plot that on top of the simulation results!\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Overlay the analytical solution at steady state\n",
    "C0=mdot/unum # steady state concentration for x>0\n",
    "left=(x<=0) \n",
    "\n",
    "Csoln=np.where( x<0, \n",
    "                C0*np.exp(unum*x/max(Dnum,1e-6)),\n",
    "                C0 )\n",
    "ax.plot(x,Csoln,'k--',lw=5,alpha=0.5) ;"
   ]
  }
 ],
 "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.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}