03sim/x3p36_wave2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "5fd2afa0",
   "metadata": {},
   "source": [
    "This is the very fast vectorized solver for the explicit, leapfrog scheme. This vectorization will also speed up the code on a parallel, multi-core system. Not only are there no loops over space in the code, but there are no matrices that need to be assembled and stored, not even sparse ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "85bb13e1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import newaxis\n",
    "\n",
    "def solver(U0, U1, f, c, Lx, Ly, Nx, Ny, dt, T, exact_solution):\n",
    "\n",
    "    x = np.linspace(0, Lx, Nx+1)  # Mesh points in x dir\n",
    "    y = np.linspace(0, Ly, Ny+1)  # Mesh points in y dir\n",
    "    # Make sure dx, dy, and dt are compatible with x, y and t\n",
    "    dx = x[1] - x[0]\n",
    "    dy = y[1] - y[0]\n",
    "    #dt = t[1] - t[0]\n",
    "\n",
    "    xv = x[:,newaxis]          # For vectorized function evaluations\n",
    "    yv = y[newaxis,:]\n",
    "\n",
    "    stability_limit = (1/float(c))*(1/np.sqrt(1/dx**2 + 1/dy**2))\n",
    "    if dt <= 0:                # max time step?\n",
    "        safety_factor = -dt    # use negative dt as safety factor\n",
    "        dt = safety_factor*stability_limit\n",
    "    elif dt > stability_limit:\n",
    "        print ('error: dt=%g exceeds the stability limit %g' % \\\n",
    "              (dt, stability_limit))\n",
    "    Nt = int(round(T/float(dt)))\n",
    "    t  = np.linspace(0, Nt*dt, Nt+1)    # mesh points in time\n",
    "    dt = t[1] - t[0]  # to avoid rounding errors\n",
    "    Cx2 = (c*dt/dx)**2;  Cy2 = (c*dt/dy)**2    # helper variables\n",
    "    dt2 = dt**2\n",
    "\n",
    "    u     = np.zeros((Nx+1,Ny+1))  # Solution array\n",
    "    u_n   = np.zeros((Nx+1,Ny+1))  # Solution at t-dt\n",
    "    u_nm1 = np.zeros((Nx+1,Ny+1))  # Solution at t-2*dt\n",
    "    f_a   = np.zeros((Nx+1,Ny+1))  # For compiled loops\n",
    "\n",
    "    Ix = range(0, u.shape[0])\n",
    "    Iy = range(0, u.shape[1])\n",
    "    It = range(0, t.shape[0])\n",
    "\n",
    "    import time; t0 = time.process_time()  # For measuring CPU time\n",
    "\n",
    "    # Load initial condition into u_n\n",
    "    u_n[:,:] = U0(xv, yv)\n",
    "\n",
    "    n = 0\n",
    "    # First step requires a special formula\n",
    "    f_a[:,:] = f(xv, yv, t[n])  # precompute, size as u\n",
    "    U1_a     = U1(xv, yv)\n",
    "    u_xx = u_n[:-2,1:-1] - 2*u_n[1:-1,1:-1] + u_n[2:,1:-1]\n",
    "    u_yy = u_n[1:-1,:-2] - 2*u_n[1:-1,1:-1] + u_n[1:-1,2:]\n",
    "    u[1:-1,1:-1] = u_n[1:-1,1:-1] +  dt*U1_a[1:-1, 1:-1] + \\\n",
    "                   0.5*Cx2*u_xx + 0.5*Cy2*u_yy + 0.5*dt2*f_a[1:-1,1:-1]\n",
    "    # Boundary condition u=0\n",
    "    j = 0 ;             u[:,j] = 0\n",
    "    j = u.shape[1]-1 ;  u[:,j] = 0\n",
    "    i = 0 ;             u[i,:] = 0\n",
    "    i = u.shape[0]-1 ;  u[i,:] = 0\n",
    "    # Update solution for next step\n",
    "    #u_nm1[:] = u_n;  u_n[:] = u  # safe, but slow\n",
    "    u_nm1, u_n, u = u_n, u, u_nm1\n",
    "    # Loop over time\n",
    "    for n in It[1:-1]:\n",
    "        f_a[:,:] = f(xv, yv, t[n])  # precompute, size as u\n",
    "        u_xx = u_n[:-2,1:-1] - 2*u_n[1:-1,1:-1] + u_n[2:,1:-1]\n",
    "        u_yy = u_n[1:-1,:-2] - 2*u_n[1:-1,1:-1] + u_n[1:-1,2:]\n",
    "        u[1:-1,1:-1] = 2*u_n[1:-1,1:-1] - u_nm1[1:-1,1:-1] + \\\n",
    "                       Cx2*u_xx + Cy2*u_yy + dt2*f_a[1:-1,1:-1]\n",
    "        # Boundary condition u=0\n",
    "        j = 0 ;             u[:,j] = 0\n",
    "        j = u.shape[1]-1 ;  u[:,j] = 0\n",
    "        i = 0 ;             u[i,:] = 0\n",
    "        i = u.shape[0]-1 ;  u[i,:] = 0\n",
    "        # Update solution for next step\n",
    "        u_nm1, u_n, u = u_n, u, u_nm1\n",
    "\n",
    "    # Error at t=T\n",
    "    u_e  = exact_solution(xv, yv, t[n])\n",
    "    diff = abs(u_nm1 - u_e).max()\n",
    "    print ('Error at time t = %g is %g' % (T, diff))\n",
    "    # Important to set u = u_n if u is to be returned!\n",
    "    t1 = time.process_time()\n",
    "    # dt might be computed in this function so return the value\n",
    "    return dt, t1 - t0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77138d54",
   "metadata": {},
   "source": [
    "This is the definition of the quadratic test problem, for which the numerical solution should be exact, to within machine precision. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "52b9e9bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def quadratic(Nx, Ny):\n",
    "\n",
    "    def exact_solution(x, y, t):\n",
    "        return x*(Lx - x)*y*(Ly - y)*(1 + 0.5*t)\n",
    "\n",
    "    def U0(x, y):\n",
    "        return exact_solution(x, y, 0)\n",
    "\n",
    "    def U1(x, y):\n",
    "        return 0.5*exact_solution(x, y, 0)\n",
    "\n",
    "    def f(x, y, t):\n",
    "        return 2*c**2*(1 + 0.5*t)*(y*(Ly - y) + x*(Lx - x))\n",
    "\n",
    "    Lx = 5;  Ly = 2\n",
    "    c = 1.5\n",
    "    dt = -1 # use longest possible steps\n",
    "    T = 18\n",
    "\n",
    "    new_dt, cpu = solver(\n",
    "        U0, U1, f, c, Lx, Ly, Nx, Ny, dt, T, exact_solution)\n",
    "    return new_dt, cpu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "45efbdd6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error at time t = 18 is 1.20792e-13\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.06189844605901729, 0.038959999999999884)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quadratic(20,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f181ea00",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error at time t = 18 is 2.13163e-13\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.030949223029508643, 0.09991799999999995)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quadratic(40,40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "46d73800",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "582.5242718446602"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "18/.0309"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "fc13357f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error at time t = 18 is 3.48166e-13\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(0.012379689211803458, 0.51606)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "quadratic(100,100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "0dc70312",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1453.9579967689822"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "18/0.01238"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "cell_metadata_filter": "-all",
   "executable": "/usr/bin/env python",
   "notebook_metadata_filter": "-all",
   "text_representation": {
    "extension": ".py",
    "format_name": "light"
   }
  },
  "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}