05-advection-ETBS-dissipation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5af539f-b704-4e1f-8ad7-4cfd71bebd76",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%config InlineBackend.figure_format = 'retina'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e798996c-ac5f-4391-98f9-f4e6e6b29196",
   "metadata": {},
   "source": [
    "### Problem setup\n",
    "\n",
    "We set up a problem with constant advection.  There are a few things to adjust in this code block:\n",
    "* `T` is the period.  You can observe more dissipation with a longer run.\n",
    "* `nx` is the spatial resolution.\n",
    "* `gamma` is the factor $a h_t / h_x$ and will be a central component of $s(kh_x)$\n",
    "* `kmode` to observe the dissipation directly, start with a specific Fourier mode $\\sin(k_{\\text{mode}} \\pi x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a7241d6f-ee52-4961-b394-959175d34941",
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 1.0 # advection speed\n",
    "T = 1.0 # end time\n",
    "\n",
    "nx = 100# spacial size\n",
    "\n",
    "x = np.linspace(0, 1, nx, endpoint=False)\n",
    "hx = x[1] - x[0]\n",
    "xx = np.linspace(0, 1, 1000, endpoint=False)\n",
    "\n",
    "gamma = 0.95 # gamma parameter\n",
    "ht = hx * gamma / c\n",
    "nt = int(T/ht)\n",
    "\n",
    "print(f'     T = {T}')\n",
    "print(f'tsteps = {nt}')\n",
    "print(f'    hx = {hx}')\n",
    "print(f'    ht = {ht}')\n",
    "print(f' gamma = {gamma}')\n",
    "\n",
    "def pulse(x):\n",
    "    return np.exp(-100*(x-0.5)**2)\n",
    "\n",
    "kmode = 6 # which mode?\n",
    "def sinusoidal(x):\n",
    "    return np.sin(kmode * np.pi * x)\n",
    "\n",
    "def hump(x):\n",
    "    u = np.zeros_like(x)\n",
    "    I = np.where((x>=0.2) & (x<=0.4))[0]\n",
    "    u[I] = -100 * (x[I]-0.2) * (x[I]-0.4)\n",
    "    return u\n",
    "\n",
    "uinit = sinusoidal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40979b43-68ca-4706-9dbe-23df90331866",
   "metadata": {},
   "source": [
    "#### Plot the initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08abb6e4-4e55-42c2-9434-0b2eeae8bd91",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(3, 3))\n",
    "ax.plot(xx, uinit(xx), lw=2, clip_on=False)\n",
    "ax.plot(x, uinit(x), lw=2, clip_on=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "63a7f943-61fa-4c70-9775-7c6d25e5fb69",
   "metadata": {},
   "source": [
    "### Run `nt` steps\n",
    "\n",
    "Here we can run ETBS for `nt` steps (set by `gamma` above).\n",
    "\n",
    "Of note -- at the end of `T`, plot the exact solution (blue) and the numerical solution (organge).  Then calculate\n",
    "$$\n",
    "s(k h_x) = 1 - \\gamma + \\gamma e^{- \\imath \\pi k h_x}\n",
    "$$\n",
    "for $k=k_{\\text{mode}}$.  Plot this with a dashed line.\n",
    "\n",
    "Notice the reduction in mode $k$ after `nt` steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f2759ed-3268-427b-b87b-77add8a958f8",
   "metadata": {},
   "outputs": [],
   "source": [
    "J = np.arange(0, nx)  # all vertices\n",
    "Jm1 = np.roll(J, 1)\n",
    "Jp1 = np.roll(J, -1)\n",
    "\n",
    "u = uinit(x)\n",
    "for tstep in range(1, nt+1):\n",
    "    \n",
    "    u[J] = u[J] - gamma * (u[J] - u[Jm1])\n",
    "    \n",
    "    #u[J] = u[J] - (gamma / 2.0) * (u[Jp1] - u[Jm1]) \\\n",
    "    #     + (gamma**2 / 2.0) * (u[Jp1] - 2 * u[J] + u[Jm1])\n",
    "    \n",
    "    uexact = uinit((xx - c * tstep * ht) % 1.0)\n",
    "\n",
    "plt.plot(xx, uexact, 'tab:blue')\n",
    "plt.plot(x, u, 'tab:orange')\n",
    "\n",
    "s = 1 - gamma + gamma * np.exp(- 1j * np.pi * kmode * hx)\n",
    "amp = np.abs(s)**nt\n",
    "plt.hlines(amp, 0, 1, 'k', '--')\n",
    "print(amp)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6fe3f0af-1968-417a-8a36-bbad5bb6b04a",
   "metadata": {},
   "source": [
    "#### Why?\n",
    "\n",
    "Let's take a look at the amplification factor which controls dissipation, namely $|s|$ in the term\n",
    "$$\n",
    "z_{j,\\ell} = z_0 |s|^{\\ell} e^{-i \\pi k j h_x - (\\frac{\\phi}{h_t})\\ell h_t}\n",
    "$$\n",
    "Notice that this is for *any* mode $k$ since $s = s(kh_x)$.\n",
    "\n",
    "Let's look at a few things.  Let's plot *every* value of $s(kh_x)$ (for a fixed $h_x$) and for all modes $k$.  Also note in this form we're using a modification with $\\pi$ in the exponent, in order to map to the periodic $[0,1]$ domain.  Here is $s(kh_x)$:\n",
    "$$\n",
    "s(k h_x) = 1 - \\gamma + \\gamma e^{- \\imath \\pi k h_x}\n",
    "$$\n",
    "\n",
    "This unit circle is in black.  The circle formed by $s()$ is dashed.  Then the $k$ modes from $0\\ldots n_x-1$ are plotted on top of that.  In these cases, it appears that the highest mode $n_x-1$ has the lowest value of $|s|$.  The k-mode above is highlight with a black/white circle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0b60c2d9-bad5-484a-aea6-c1b0ac2392ae",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots()\n",
    "\n",
    "# unit circle\n",
    "t = np.linspace(0, np.pi * 2.0, 1000)\n",
    "ax.plot(np.cos(t), np.sin(t), 'k.', ms=1)\n",
    "\n",
    "# All k\n",
    "K = np.arange(nx)\n",
    "s = 1 - gamma + gamma * np.exp(- 1j * np.pi * K * hx)\n",
    "l=ax.scatter(s.real, s.imag,\n",
    "             c=K, cmap='rainbow')\n",
    "\n",
    "# The k-mode above\n",
    "s0 = 1 - gamma + gamma * np.exp(- 1j * np.pi * kmode * hx)\n",
    "ax.plot(s0.real, s0.imag, 'ko', markerfacecolor='w', ms=4)\n",
    "\n",
    "c=plt.colorbar(l)\n",
    "c.set_label('mode #, $k$')\n",
    "ax.set_aspect(1)\n",
    "ax.spines['left'].set_position('center')\n",
    "ax.spines['bottom'].set_position('center')\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.set_title(f'$\\gamma$={gamma}');"
   ]
  }
 ],
 "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
}