{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical diffusion\n",
    "\n",
    "Numerical diffusion is a common problem with advection of field data on an Eulerian particle grid (one that is fixed in space). The code below can be used to demonstrate numerical diffusion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting started\n",
    "\n",
    "Here we consider advection of a rock mass with density $\\rho$.\n",
    "\n",
    "\\begin{equation}\n",
    "  \\frac{\\partial \\rho}{\\partial t} = -v_{x} \\left( \\frac{\\partial \\rho}{\\partial x} \\right)\n",
    "\\end{equation}\n",
    "\n",
    "where $x$ is the spatial coordinate and $t$ is time.\n",
    "\n",
    "The upwind finite-difference solution to this advection equation is\n",
    "\n",
    "\\begin{equation}\n",
    "  \\rho_{n}^{i} = \\rho_{n}^{i-1} -v_{x} \\Delta t \\frac{\\rho_{n}^{i-1} - \\rho_{n-1}^{i-1}}{\\Delta x}.\n",
    "\\end{equation}\n",
    "\n",
    "The code below calculates this finite-difference solution and produces a plot of the output."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Setting up the problem\n",
    "\n",
    "First we can import the necessary libraries and configure the problem geometry.  We start by setting the spatial scale `xmax`, number of grid points `nx`, advection velocity `vx`, and time step `dt`. Then we can create the spatial coordinate (`x`) and density (`density`, `densitynew`) arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xmax = 0.25                      # Length scale of problem\n",
    "nx = 30                          # Number of points to consider along x axis\n",
    "vx = 1.0                         # Advection velocity\n",
    "dt = 0.001                       # Time step\n",
    "nsteps = 21                      # Number of time steps\n",
    "\n",
    "x = np.linspace(0.0, xmax, nx)   # Spatial array\n",
    "dx = x[1] - x[0]                 # Grid spacing\n",
    "density = np.zeros(nx) + 3200    # Density array 1\n",
    "densitynew = np.copy(density)    # Density array 2\n",
    "\n",
    "# We initialize the density to 3300 between x = 0.05 and x = 0.10\n",
    "density[(x >= 0.05) & (x <= 0.10)] = 3300"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the solution and plotting the results\n",
    "\n",
    "Now we can calculate the finite difference solution and plot the results. Note that we plot only every 5th time step for the output. Note that if you make changes to the geometry of the problem (cell above) you will need to run both that cell and the one below to create a new plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(12,6))\n",
    "for t in range(nsteps):\n",
    "    for i in range(1, len(density)-1):\n",
    "        densitynew[i] = density[i] - vx * dt * (density[i] - density[i-1]) / dx\n",
    "    # Plot every 5th time step\n",
    "    if t%5 == 0:\n",
    "        plt.plot(x, density, '-o', label=\"Step {0:d}\".format(t))\n",
    "    density = densitynew\n",
    "\n",
    "#plt.plot(x, densityOrig)\n",
    "plt.ylim([3180, 3320])\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What to do?\n",
    "\n",
    "Feel free to play with the time step, number of time steps, advection velocity, etc.\n",
    "You should find that numerical diffusion is problematic no matter what you do (other than `vx = 0`)."
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}