05-advection-1d-template.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem description"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you will set up the problem for\n",
    "$$ u_t + c u_x = 0$$\n",
    "with periodic BC on the interval [0,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c = 1.0\n",
    "T = 1.0 / c # end time"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Set up the grid"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`hx` is the grid spacing in the `x`-direction\n",
    "\n",
    "`x` are the grid coordinates\n",
    "\n",
    "`xx` are really fine grid coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nx = 32\n",
    "x = np.linspace(0, 1, nx, endpoint=False)\n",
    "hx = x[1] - x[0]\n",
    "xx = np.linspace(0, 1, 100, endpoint=False)\n",
    "print(hx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now define an initial condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def square(x):\n",
    "    u = np.zeros(x.shape)\n",
    "    u[np.intersect1d(np.where(x>0.4), np.where(x<0.6))] = 1.0\n",
    "    return u\n",
    "f = square\n",
    "plt.plot(x, f(x), lw=3, clip_on=False)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Setting time step $h_t$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we need a time step.  Let\n",
    "$$ \\lambda = a \\frac{h_t}{h_x}$$\n",
    "and\n",
    "$$ h_t = h_x \\frac{\\lambda}{c}$$\n",
    "\n",
    "So we first set the parameter $\\lambda$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lmbda = 0.95\n",
    "ht = hx * lmbda / c\n",
    "nt = int(T/ht)\n",
    "\n",
    "print('Setting...')\n",
    "print(f'lambda = {lmbda}')\n",
    "print(f'     T = {T}')\n",
    "print(f'    hx = {hx}')\n",
    "print('...results in...')\n",
    "print(f'tsteps = {nt}')\n",
    "print(f'    ht = {ht}')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now make an index list, called $J$, so that we can access $J+1$ and $J-1$ easily"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "J = np.arange(0, nx)  # all vertices\n",
    "Jm1 = np.roll(J, 1)\n",
    "Jp1 = np.roll(J, -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "err = np.zeros((nt+1, nx))\n",
    "u = f(x)\n",
    "\n",
    "for tstep in range(nt+1):\n",
    "    u[J] = u[J] - lmbda * (u[J] - u[Jm1])           \n",
    "\n",
    "    err[tstep, :] = f((x - c * tstep * ht) % 1.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Check the error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "uexact = f((x - c * (nt+1) * ht) % 1.0)\n",
    "error = u - uexact\n",
    "l2err = np.sqrt(hx * np.sum(error**2))\n",
    "print(l2err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "error = np.ones(len(u))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "l2err = np.sqrt(hx * np.sum(error**2))\n",
    "print(l2err)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": 4
}