1D_Inviscid_Burgers_nodal_DG.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b8d6c733",
   "metadata": {},
   "source": [
    "#### Ignoring Runtime Warnings! This could be dangerous"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca9269e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\", category=RuntimeWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "986c5451",
   "metadata": {},
   "source": [
    "# Discretizing the inviscid Burgers' equation with SBP Operators"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9bcb6c52",
   "metadata": {},
   "source": [
    "Let us solve the one-dimensional advection equation with periodic boundary conditions\n",
    "$$\n",
    "u_t + \\left(\\frac12 u^2\\right)_x = 0\\;\\; \\text{on }[a,b], \\\\\n",
    "u(x,0)=u_0(x), \\\\\n",
    "u(a,t)=g(t).\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dfe1c0bd",
   "metadata": {},
   "source": [
    "Importing some libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "661b6d4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.optimize import fsolve\n",
    "import math\n",
    "\n",
    "#Local libraries\n",
    "from basis_nodes import generate_lagrange_poly, generate_LGL_points\n",
    "from SBP_matrices import first_order_D, first_order_P_Q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a41cbfa2",
   "metadata": {},
   "source": [
    "Defining parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5af6bc4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "## Computational domain\n",
    "\n",
    "#Each element i will be enclosed by x_{i} and x_{i+1}\n",
    "n_elements = 40\n",
    "xlower = 0.\n",
    "xupper = 2.\n",
    "x = np.linspace(xlower,xupper, n_elements+1)\n",
    "\n",
    "#Advection speed\n",
    "c = 1.\n",
    "\n",
    "#Spatial formal order of accuracy\n",
    "order_x = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05cc0ce",
   "metadata": {},
   "source": [
    "Computing nodes at reference element"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "44cfbca7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#We'll need a Lagrange basis of degree order_x-1\n",
    "#Thus, we will need order_x nodes in each element\n",
    "degree_basis = order_x-1\n",
    "if degree_basis >0:\n",
    "    (_, _,\n",
    "    _, _,\n",
    "    _, _,\n",
    "    xi_LGL_unsorted, w_LGL_unsorted) = generate_LGL_points(degree_basis)\n",
    "else:\n",
    "    xi_LGL_unsorted, w_LGL_unsorted = np.array([0.]), np.array([1.])\n",
    "\n",
    "#######################\n",
    "#######################\n",
    "#The root finder returns unsorted roots,\n",
    "#We fix that in these lines\n",
    "\n",
    "#Indexes that would sort x_LG using value\n",
    "sort_idxs_LGL = np.argsort(xi_LGL_unsorted)\n",
    "#Sorting x_LG and w_LG using those indexes\n",
    "xi_LGL = xi_LGL_unsorted[sort_idxs_LGL]\n",
    "w_LGL = w_LGL_unsorted[sort_idxs_LGL]\n",
    "#######################\n",
    "#######################\n",
    "\n",
    "n_nodes = len(xi_LGL)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6399b06a",
   "metadata": {},
   "source": [
    "Mapping reference element to whole grid. The interface points are stored twice for simplicity. This will be just used for plotting purposes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a43269de",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "#We map the nodes from the reference element to one of the elements \n",
    "#from our computational domain (this suffices since we are using an uniform grid)\n",
    "x_min = x[0]\n",
    "x_max = x[1]\n",
    "\n",
    "#We just need to map the LGL nodes once since we are using a uniform grid\n",
    "x_element = (xi_LGL*(x_max-x_min)+(x_min+x_max))/2\n",
    "len_element = (x_max-x_min)/2. #Constant in this case\n",
    "w_element = (x_max-x_min)*w_LGL\n",
    "list_elements = [x_element]\n",
    "\n",
    "#We define our grid (interface points will be repeated)\n",
    "x_grid = np.zeros(n_elements*n_nodes)\n",
    "\n",
    "#print(f\"Numerical grid with double interface nodes\")\n",
    "for i in range(n_elements): \n",
    "    j = i*n_nodes\n",
    "    x_grid[j:j+n_nodes] = x_element+x[i]\n",
    "    list_elements.append(x_grid[j:j+n_nodes])\n",
    "#    print(f\"Element {i+1:}\\t {x_grid[j:j+n_nodes]}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "964a75bd",
   "metadata": {},
   "source": [
    "## SBP Operators\n",
    "Obtaining local relevant SBP operators and resizing them over the whole grid with Kronecker products"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "188d1d46",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "#Local opretaros P and Q\n",
    "P_LGL, Q_LGL = first_order_P_Q(x_Lagrange_nodes=xi_LGL, x_abcissae=xi_LGL, w_abcissae=w_LGL)\n",
    "\n",
    "#Global operators P and Q\n",
    "P = np.kron(np.eye(n_elements), P_LGL)\n",
    "Q = np.kron(np.eye(n_elements), Q_LGL)\n",
    "\n",
    "diagP=np.diag(P)\n",
    "diagPinv=1./diagP\n",
    "\n",
    "#Restriction operators R, B\n",
    "#Local\n",
    "R_LGL = np.zeros((2,n_nodes))\n",
    "R_LGL[0,0] = R_LGL[-1,-1] = 1\n",
    "B_LGL = np.zeros((2,2))\n",
    "B_LGL[0,0] = -1; B_LGL[-1,-1]=1\n",
    "#Global\n",
    "R = np.kron(np.eye(n_elements), R_LGL)\n",
    "B = np.kron(np.eye(n_elements), B_LGL)\n",
    "\n",
    "\n",
    "#Differential operator D\n",
    "#Local\n",
    "D_LGL = first_order_D(x_nodes=xi_LGL)\n",
    "#Global\n",
    "D = np.kron(np.eye(n_elements), D_LGL)\n",
    "\n",
    "\n",
    "#Vector and matrix with advection speed\n",
    "c_vec = c*np.ones_like(x_grid)\n",
    "c_mat = c*np.eye(len(x_grid))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7eb9008f",
   "metadata": {},
   "source": [
    "## Numerical Flux"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e15d20c3",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Flux function\n",
    "flux_fun = lambda u: 0.5*(u)**2\n",
    "\n",
    "#Vectorized flux function (maybe not needed at this point)\n",
    "flux_fun_vec = np.vectorize(flux_fun)\n",
    "\n",
    "#Numerical flux\n",
    "#Local\n",
    "#f_num_loc = lambda uL, uR : (flux_fun(uL)+flux_fun(uR))/2 #Centered flux\n",
    "f_num_loc = lambda uL, uR : flux_fun(uL) #Upwind flux\n",
    "#Global\n",
    "def f_num(u,n_elements=n_elements,n_nodes=n_nodes):\n",
    "    #Computing numerical flux just at the interface of elements\n",
    "    #Fill numerical flux vector for first element\n",
    "    \n",
    "    if n_elements == 1:\n",
    "        f = [flux_fun(1.), flux_fun(u[0])]    #The first component should be c*source at the left\n",
    "                            #Here the source is g(t)=1\n",
    "        return np.array(f)\n",
    "    \n",
    "    #For general periodic BCs\n",
    "    f = [flux_fun(1.),         #Left interface\n",
    "        flux_fun(u[n_nodes-1])]  #Right interface\n",
    "    for idx_elem in range(1,n_elements-1):\n",
    "        #Left interface\n",
    "        idx_R = idx_elem*n_nodes\n",
    "        idx_L = idx_R-1\n",
    "        f.append(f_num_loc(u[idx_L],u[idx_R]))\n",
    "        #print(f\"Element {idx_elem+1}\")\n",
    "        #print(f\"Left interface\")\n",
    "        #print(f\"{idx_L, idx_R}\")\n",
    "        #Right interface\n",
    "        idx_R = (idx_elem+1)*n_nodes\n",
    "        idx_L = idx_R-1\n",
    "        f.append(f_num_loc(u[idx_L],u[idx_R]))\n",
    "        #print(f\"Right interface\")\n",
    "        #print(f\"{idx_L, idx_R}\")\n",
    "        \n",
    "    #Fill numerical flux vector for last element\n",
    "    f.append(flux_fun(u[-n_nodes-1])) #Left interface\n",
    "    f.append(flux_fun(u[-1]))       #Right interface\n",
    "        #print(f\"Last element\")\n",
    "        #print(f\"Left interface\")\n",
    "        #print(f\"{-n_nodes-1, -n_nodes}\")\n",
    "        #print(f\"Right interface\")\n",
    "        #print(f\"{n_nodes*n_elements-1, 0}\")\n",
    "    return np.array(f)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f09eb90",
   "metadata": {},
   "source": [
    "## Semi-discrete RHS operator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "be06c23e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def RHS(u):\n",
    "    #split_form_interior = -0.5*D@c_mat@u -0.5*c_mat@D@u -0.5*np.diag(u)@D@c_vec\n",
    "    non_split_form_interior = -D@flux_fun(u)\n",
    "    elem_boundary_terms = -diagPinv*(R.T@B@(f_num(u=u)-R@flux_fun(u)))\n",
    "    return (1./len_element)*(non_split_form_interior+elem_boundary_terms)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09e01370",
   "metadata": {},
   "source": [
    "## Plotting function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "91022b26",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from matplotlib.animation import FuncAnimation\n",
    "\n",
    "def plotting_sol(u_frames,x_grid, t_eval, critical_time, plotexact):\n",
    "    # Set up the figure and axis\n",
    "    fig, ax = plt.subplots(dpi=200, figsize=[8,4])\n",
    "    ax.set_xlim(x_grid.min(), x_grid.max())\n",
    "    minu = u_frames[0].min()\n",
    "    maxu = u_frames[0].max()\n",
    "    margin = 0.5\n",
    "    ax.set_ylim(minu-margin*max(np.abs(minu),1.), maxu+margin*max(np.abs(maxu),1.))\n",
    "    line_approx, = ax.plot(x_grid,u_frames[0], label= \"SBP solution\")#ax.plot([], [])\n",
    "\n",
    "     \n",
    "    dx = x_grid[1]-x_grid[0]\n",
    "    u0 = u_frames[0]\n",
    "    if plotexact:\n",
    "        u_exact = np.array([exact_sol(x_grid,time) for time in t_eval])\n",
    "        line_exact, = ax.plot(x_grid,u_exact[0], linestyle='dotted', color='r', label=\"exact solution\")\n",
    "    else:\n",
    "        u_exact = np.zeros_like(u_frames[0])\n",
    "    \n",
    "\n",
    "    def animate(i):\n",
    "        u_approx = u_frames[i]\n",
    "        line_approx.set_data(x_grid, u_approx)\n",
    "        mass = np.round(np.linalg.norm(u_approx,ord=2)*dx,2)\n",
    "        error = np.round(np.linalg.norm(u_exact[i]-u_approx, ord=2)*np.sqrt(dx),4)\n",
    "        ax.set_title(f\"Time: {np.round(t_eval[i])}    ||   Mass={mass}    ||   $L^2$ Error={error}\")\n",
    "        #line_approx.set_label(\"SBP solution\")\n",
    "        \n",
    "        #if critical_time > t_eval[i]:\n",
    "        #    line_exact.remove()  # remove line_exact from the plot\n",
    "        #    plotexact = False  # set plotexact to False to avoid calling set_data on line_exact\n",
    "        \n",
    "        if plotexact:\n",
    "            if critical_time <= t_eval[i]:\n",
    "                line_exact.set_data(x_grid, 1e20*np.ones_like(x_grid))  # remove line_exact from the plot\n",
    "            else: \n",
    "                line_exact.set_data(x_grid, u_exact[i])\n",
    "            line_exact.set_label(\"exact solution\")\n",
    "            return line_approx, line_exact\n",
    "        else:\n",
    "            return line_approx,\n",
    "\n",
    "    import matplotlib\n",
    "    from matplotlib import animation\n",
    "    from IPython.display import HTML\n",
    "    anim = FuncAnimation(fig, animate,\n",
    "                                       frames=len(u_frames),\n",
    "                                       interval=200,\n",
    "                                       repeat=False)\n",
    "    plt.legend()\n",
    "    plt.close()\n",
    "    return HTML(anim.to_jshtml())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33b1eaf8",
   "metadata": {},
   "source": [
    "## Time Integration\n",
    "We use Strong Stability Preserving RK schemes as described in Gottlieb et al. (2011)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da268651",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Local library\n",
    "import SSPRK"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2b89bbf9",
   "metadata": {},
   "source": [
    "## Testing implementation with different initial conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d5598d9",
   "metadata": {},
   "source": [
    "### 1.- Gaussian function\n",
    "$$\n",
    "u(x,0) = 1+\\exp\\left({-25\\, (x-\\frac12)^2}\\right)%\\sin(x).\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9922416d",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "initial_c = lambda x: np.exp(-25*(x-0.5)**2)+1\n",
    "\n",
    "\n",
    "#To find shock formation for exact solution\n",
    "initial_cprime = lambda x: np.exp(-25*(x-0.5)**2)*(-50*(x-0.5))\n",
    "critical_time = -1/np.min(initial_cprime(x_grid))\n",
    "\n",
    "\n",
    "u0 = initial_c(x_grid)\n",
    "u = np.copy(u0)\n",
    "n_grid = len(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95496b99",
   "metadata": {},
   "source": [
    "Defining exact solution for plotting purposes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4d1bca3f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_sol(x,t):\n",
    "    def nonlinear_char(xi):\n",
    "        return xi+initial_c(xi)*t-x\n",
    "    xi = fsolve(nonlinear_char, x)\n",
    "    return initial_c(xi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28cc253d",
   "metadata": {},
   "source": [
    "Solving problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aacd9d05",
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = 0.\n",
    "tfinal = 1.\n",
    "dx = min(np.min(np.diff(xi_LGL)),np.min(np.diff(x_grid[:n_nodes])))\n",
    "CFL = 0.1\n",
    "dt = CFL*dx/np.max(u0+0.5)\n",
    "nframes = 50\n",
    "u_frames, t_eval = SSPRK.SSPRK33(fun=RHS, u0=u0, dt=dt, t0=t0, tfinal=tfinal, nframes=nframes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3af9ac15",
   "metadata": {},
   "source": [
    "Plotting solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "80bce7dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "plotting_sol(u_frames,x_grid, t_eval, critical_time=critical_time, plotexact=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0df79c7f",
   "metadata": {},
   "source": [
    "### 2.- Trigonometric functions\n",
    "$$\n",
    "u(x,0) = -\\frac12 \\tanh(3\\,x-3)+\\frac12\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2a677691",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "initial_c = lambda x: -0.5*np.tanh(3*x-3)+0.5\n",
    "\n",
    "#To find shock formation for exact solution\n",
    "initial_cprime = lambda x: -3*0.5*(1./np.cosh(3*x-3)**2)\n",
    "critical_time = -1/np.min(initial_cprime(x_grid))\n",
    "\n",
    "u0 = initial_c(x_grid)\n",
    "u = np.copy(u0)\n",
    "n_grid = len(u)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1d3e2bf4",
   "metadata": {},
   "source": [
    "Defining exact solution for plotting purposes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0846b6ef",
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_sol(x,t):\n",
    "    def nonlinear_char(xi):\n",
    "        return xi+initial_c(xi)*t-x\n",
    "    xi = fsolve(nonlinear_char, x)\n",
    "    return initial_c(xi)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fc98ffb",
   "metadata": {},
   "source": [
    "Solving problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "21a6f299",
   "metadata": {},
   "outputs": [],
   "source": [
    "t0 = 0.\n",
    "tfinal = 1.\n",
    "dx = min(np.min(np.diff(xi_LGL)),np.min(np.diff(x_grid[:n_nodes])))\n",
    "CFL = 0.1\n",
    "dt = CFL*dx/10.#np.max(u0+0.5)\n",
    "nframes = 50\n",
    "u_frames, t_eval = SSPRK.SSPRK33(fun=RHS, u0=u0, dt=dt, t0=t0, tfinal=tfinal, nframes=nframes)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "534ef9d7",
   "metadata": {},
   "source": [
    "Plotting solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55834582",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plotting_sol(u_frames,x_grid, t_eval, critical_time=critical_time, plotexact=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "23c0ab45",
   "metadata": {},
   "source": [
    "Some debugging lines for the time integrator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0fa5e05b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from nodepy import rk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "12bc8de7",
   "metadata": {},
   "outputs": [],
   "source": [
    "#rkm = rk.ExplicitRungeKuttaMethod(A_Heun,b_Heun)\n",
    "#rkm.plot_stability_region();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca2a300e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#rkm.imaginary_stability_interval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99a17088",
   "metadata": {
    "hideOutput": true
   },
   "outputs": [],
   "source": [
    "#rkm.order()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "3.8.9",
   "language": "python",
   "name": "pyenv_3.8.9"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}