06-advection-general.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be9d6bd0-37ab-46e6-b3cd-b401b3343956",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "u_t + (c(x) u)_x = 0 on [a,b]\n",
    "with u(b) = u(a)\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import animation\n",
    "from matplotlib import rc\n",
    "rc('animation', html='jshtml')\n",
    "from tqdm.notebook import tqdm\n",
    "%config InlineBackend.figure_format='retina'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e93f244f-4bf2-4ee5-bea7-6b4bdc7b39ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x):\n",
    "    u = np.exp(-100 * (x - 0.25)**2)\n",
    "    return u\n",
    "\n",
    "def step(x):\n",
    "    u = np.zeros(x.shape)\n",
    "    for j in range(len(x)):\n",
    "        if (x[j] >= 0.6) and (x[j] <= 0.8):\n",
    "            u[j] = 1.0\n",
    "\n",
    "    return u\n",
    "\n",
    "def c(x):\n",
    "    return 1.0 + 0.0 * x\n",
    "    #return 2.0 + np.sin(2 * np.pi * x)\n",
    "\n",
    "T = 1.0\n",
    "gamma = 0.95\n",
    "nx = 128\n",
    "\n",
    "x, hx = np.linspace(0, 1, nx, endpoint=False, retstep=True)\n",
    "xx = np.linspace(0, 1-hx, 1000)\n",
    "\n",
    "ht = hx * gamma / c(x).max()\n",
    "nt = int(np.ceil(T/ht))\n",
    "ht = T/nt\n",
    "\n",
    "print('T = %g' % T)\n",
    "print('tsteps = %d' % nt)\n",
    "print('    hx = %g' % hx)\n",
    "print('    ht = %g' % ht)\n",
    "print('lambda = %g' % gamma)\n",
    "\n",
    "K = np.arange(0, nx)  # all vertices\n",
    "Km1 = np.roll(K, 1)\n",
    "Kp1 = np.roll(K, -1)\n",
    "\n",
    "u = gaussian(x) + step(x)\n",
    "u0 = u.copy()\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(x, u0, 'r-', linewidth=1)\n",
    "uline, = ax.plot([], [], '-', lw=3, color='tab:blue')\n",
    "txt = ax.text(1.0 / 3.0, 1.1, 't=0, i=0', fontsize=16)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "96978b67-759b-4cb2-bdb2-c7ec04894cb5",
   "metadata": {},
   "outputs": [],
   "source": [
    "def minmod(a, b):\n",
    "    c = 0 * a\n",
    "    for i in range(len(a)):\n",
    "        if a[i] * b[i] > 0:\n",
    "            if abs(a[i]) <= abs(b[i]):\n",
    "                c[i] = a[i]\n",
    "            else:\n",
    "                c[i] = b[i]\n",
    "        else:\n",
    "            c[i] = 0.0\n",
    "    return c\n",
    "\n",
    "\n",
    "def maxmod(a, b):\n",
    "    c = 0 * a\n",
    "    for i in range(len(a)):\n",
    "        if a[i] * b[i] > 0:\n",
    "            if abs(a[i]) >= abs(b[i]):\n",
    "                c[i] = a[i]\n",
    "            else:\n",
    "                c[i] = b[i]\n",
    "        else:\n",
    "            c[i] = 0.0\n",
    "    return c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7012f8c-f5a2-4046-a1bd-ac1b194cb389",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_tv(u):\n",
    "    return np.sum(np.abs(u[K]-u[Km1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "266e2d73-8dd3-4f51-aea4-f278c3d3c5ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "TV = []\n",
    "method = 'linear'\n",
    "sigtype = 'superbee'\n",
    "u = u0.copy()\n",
    "def init():\n",
    "    return (uline,)  \n",
    "def timestep(n):\n",
    "    if n==0:\n",
    "        uline.set_data(x, u)\n",
    "        txt.set_text('t=0, i=0')\n",
    "        ax.axis([0, 1, -0.25, 1.25])\n",
    "        TV.append(calculate_tv(u))\n",
    "        return (uline,) \n",
    "\n",
    "    if method == 'constant':\n",
    "        ci = c(x[Kp1])  # at x = i + 1/2\n",
    "        fi = 0 * ci   # at x = i + 1/2\n",
    "        pos = np.where(ci >= 0)[0]\n",
    "        fi[pos] = ci[pos] * u[K[pos]]\n",
    "        neg = np.where(ci < 0)[0]\n",
    "        fi[neg] = ci[neg] * u[Kp1[neg]]\n",
    "\n",
    "        u[K] = u[K] - (ht / hx) * (fi[K] - fi[Km1])\n",
    "\n",
    "    if method == 'linear':\n",
    "\n",
    "        if sigtype == 'LaxW':\n",
    "            sig = (u[Kp1] - u[K]) / hx\n",
    "\n",
    "        if sigtype == 'BW':\n",
    "            sig = (u[K] - u[Km1]) / hx\n",
    "\n",
    "        if sigtype == 'Fromm':\n",
    "            sig = (u[Kp1] - u[Km1]) / (2.0 * hx)\n",
    "\n",
    "        if sigtype == 'minmod' or sigtype == 'superbee':\n",
    "            sigm = (u[K] - u[Km1]) / hx\n",
    "            sigp = (u[Kp1] - u[K]) / hx\n",
    "\n",
    "            if sigtype == 'minmod':\n",
    "                sig = minmod(sigm, sigp)\n",
    "\n",
    "            if sigtype == 'superbee':\n",
    "                sig1 = minmod(sigm, 2 * sigp)\n",
    "                sig2 = minmod(2 * sigm, sigp)\n",
    "                sig = maxmod(sig1, sig2)\n",
    "\n",
    "        u[K] = u[K] - (c(x[K]) * ht / hx) * (u[K] - u[Km1])\\\n",
    "                    - (c(x[K]) * ht / (2.0 * hx)) * (sig[K] - sig[Km1])\\\n",
    "                                                  * (hx - c(x[K]) * ht)\n",
    "    TV.append(calculate_tv(u))\n",
    "    uline.set_data(x, u)\n",
    "    txt.set_text(f't={n*ht}, i={n}')\n",
    "    ax.axis([0, 1, -0.25, 1.25])\n",
    "    return (uline,)\n",
    "\n",
    "\n",
    "anim = animation.FuncAnimation(fig, timestep,\n",
    "                               init_func=init,\n",
    "                               frames=tqdm(range(0,nt+1)), interval=20,\n",
    "                               blit=True, repeat=False);\n",
    "\n",
    "anim"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "493a865c-bb0a-4f59-84e4-873986e7ef6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(TV, 'ko')\n",
    "plt.xlabel('timestep')\n",
    "plt.ylabel('total variation')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5ca931a-799b-41c6-b7c3-94c069b99265",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}