pgd-small-hat.ipynb.fixme

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipywidgets import interact, IntSlider, FloatSlider\n",
    "%matplotlib notebook\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Introduction\n",
    "\n",
    "The purpose of this hands-on session is two-fold: \n",
    "\n",
    " * Explaining how to evaluate a quasi-analytical solution, that is singular of the pressureless gas dynamics\n",
    " \n",
    " * Show that once the solution is becoming singular, the dynamics of the singularity is not the one of the original"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# General setting and smooth solution\n",
    "\n",
    "We consider the following Cauchy problem:\n",
    "\n",
    "$\\rho(t, \\mathbf{x})$, $\\mathbf{u}(t, \\mathbf{x})$ satisfying\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left\\{\n",
    "    \\begin{aligned}\n",
    "    \\frac{\\partial \\rho}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{u}) = 0 \\\\\n",
    "    \\frac{\\partial \\rho \\mathbf{u}}{\\partial t} + \\nabla \\cdot (\\rho \\mathbf{u} \\mathbf{u} ) = 0\n",
    "    \\end{aligned}\n",
    "    \\right.\n",
    "\\end{equation}\n",
    "\n",
    "with  $\\rho(0, \\mathbf{x}) = \\rho_0(x) = 1$\n",
    "\n",
    "\\begin{equation}\n",
    "    \\left\\{\n",
    "        \\begin{aligned}\n",
    "            \\mathbf{u}(0, x) = 0, \\; x \\in (-\\infty, -1] \\cup [-1, +\\infty) \\\\\n",
    "            \\mathbf{u}(0, x) = 1 + x, \\; x \\in [-1, 0] \\\\\n",
    "            \\mathbf{u}(0, x) = 1 - x, \\; x \\in [0, 1]           \n",
    "        \\end{aligned}\n",
    "    \\right.    \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Burgers equation\n",
    "\n",
    "The method of characteristics, provides a convenient way to calculate the evolution of both $\\rho\\left(t,\\mathbf{x}\\right)$ and $\\mathbf{u}\\left(t,\\mathbf{x}\\right)$.\\\\\n",
    "The Burgers equation can be written as:\n",
    "\n",
    "$$\n",
    "\\frac{\\partial u}{\\partial t}+ \\frac{\\partial}{\\partial x}\\left(\\frac{u^2}{2}\\right)=0\n",
    "$$\n",
    "\n",
    "That is the Burgers equation for which we know the solution. The density field becomes discontinuous but can be evaluated for all $t \\in [0,1)$.\\\\\n",
    "Anywhere on $\\left(\\left.-\\infty,-1\\right]\\right. \\cup \\left[\\left.1,+\\infty\\right)\\right.$, we still have $\\rho(t,\\mathbf{x})=1 \\forall t \\in \\left[\\left. 0,1 \\right) \\right.$.\\\\\n",
    "In the expansion wave, we have:\n",
    "\n",
    "$$\n",
    "\\frac{D}{Dt}\\log \\rho = -\\frac{\\partial u}{\\partial x} =\\frac{-1}{1+t}\n",
    "$$so that:\n",
    "\n",
    "$$\n",
    "    \\rho\\left(t,x\\right)=\\frac{1}{1+t}\n",
    "$$\n",
    "\n",
    "with $t\\in \\left[0,1\\right]$ and $x\\in\\left[t,1\\right]$.\\\\\n",
    "Considering the mass conservation equation, in the compression wave, we have:\n",
    "\n",
    "$$\n",
    "    \\frac{D}{Dt}\\log\\rho = -\\frac{\\partial u}{\\partial x}=\\frac{1}{1-t}\n",
    "$$\n",
    "\n",
    "so that:\n",
    "\n",
    "$$\n",
    "    \\rho(t,x)=\\frac{1}{1-t}\n",
    "$$with $t\\in\\left[0,1\\right]$ and $x\\in\\left[t,1\\right]$. It clearely blows-up at $t=1$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Rankine-Hugoniot generalized problem\n",
    "\n",
    "We consider the typical singularity which can appear is a Delta shock, a shock in velocity, accompained with a Dirac delta function in density.\\\\\n",
    "At time $t=1$, conserving mass and momentum allows to deduce that the amount of mass, which concentrate in the delta shock is:\n",
    "$$\n",
    "    \\int_0^1\\rho^0(x)\\mathrm{d}x=1\n",
    "$$and the momentum:\n",
    "$$\n",
    "    \\int_0^1\\rho^0(x)u^0(x)\\mathrm{d}x=\\frac{1}{2}\n",
    "$$Thus we will start at $t=1$ by a solution given by the theory of characteristics anywhere in $x\\in R\\setminus\\{1\\}$ to which one \"add\" the following function:\n",
    "$$\n",
    "    m^1\\delta\\left(x-p^1\\right)\n",
    "$$\n",
    "with:\n",
    "$$\n",
    "\\begin{cases}\n",
    "    m(1)=m^1=1\\\\\n",
    "    p(1)=p^1=1\\\\\n",
    "    v(1)=v^1=1\n",
    "\\end{cases}\n",
    "$$$m(t)$ represents the mass in the delta shock, $p(t)$ its position, $v(t)$ its velocity.\\\\In order to resolve the dynamics of this singularity, we need generalized Rankine-Hugoniot jump conditions:\\begin{align}\n",
    "\\begin{cases}\n",
    "    \\dfrac{d\\rho(t)}{dt} &= v(t)\\\\\n",
    "    \\dfrac{d m(t)}{dt} &= -\\left(\\rho^{-}-\\rho^{+}\\right)v(t)+\\rho^{-}u^{-}-\\rho^{+}u^{+}\\\\\n",
    "    \\dfrac{d\\left(m(t)v(t)\\right)}{dt} &= -\\left(\\rho^{-}u^{-} -\\rho^{+}u^{+}\\right)v(t)+\\rho^{-}\\left(u^{-}\\right)^2-\\rho^{+}\\left(u^{+}\\right)^2\n",
    "\\end{cases}\n",
    "\\end{align}We will test two conditions:\n",
    "\n",
    "## Condition 1\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{cases}\n",
    "    \\rho^{+} &= 1\\\\\n",
    "    \\rho^{-} &=-\\frac{1}{1+t}\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{cases}\n",
    "    u^{+} &= 0\\\\\n",
    "    u^{-} &=\\frac{1+p(t)}{1+t}\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$The system above will become:\\begin{align}\n",
    "\\begin{cases}\n",
    "    \\dfrac{d\\rho(t)}{dt} &= v(t)\\\\\n",
    "    \\dfrac{d m(t)}{dt} &= \\frac{1}{1+t}\\left(tv(t)+\\frac{1+p(t)}{1+t}\\right)\\\\\n",
    "    \\dfrac{d\\left(m(t)v(t)\\right)}{dt} &=\\frac{1}{1+t}\\left[\\frac{1+p(t)}{1+t}-v(t)\\right]^2-v(t)^2\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "\n",
    "## Condition 2\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{cases}\n",
    "    \\rho^{+} &= 2\\\\\n",
    "    \\rho^{-} &=\\frac{1}{1+t}\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$\n",
    "$$\n",
    "\\begin{align}\n",
    "\\begin{cases}\n",
    "    u^{+} &= 0\\\\\n",
    "    u^{-} &=\\frac{1+p(t)}{1+t}\n",
    "\\end{cases}\n",
    "\\end{align}\n",
    "$$The system above will become:\\begin{align}\n",
    "\\begin{cases}\n",
    "    \\dfrac{d\\rho(t)}{dt} &= v(t)\\\\\n",
    "    \\dfrac{d m(t)}{dt} &= \\frac{1}{1+t}\\left(v(t)-\\frac{1+p(t)}{1-t}\\right)+2v(t)\\\\\n",
    "    \\dfrac{d\\left(m(t)v(t)\\right)}{dt} &=\\left(\\frac{1}{1+t}\\right)\\left(\\frac{1+p(t)}{1+t}\\right)\\left[v-\\frac{1+p(t)}{1+t}\\right]\n",
    "\\end{cases}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "\n",
    "def rh_gen(mvp, t, rho_plus):\n",
    "    m = mvp[0]\n",
    "    v = mvp[1]\n",
    "    p = mvp[2]\n",
    "    \n",
    "    rho_minus = 1/(1+t)\n",
    "    u_plus = 0\n",
    "    u_minus = (1+p)/(1+t)\n",
    "    mdot = -(rho_minus-rho_plus)*v+rho_minus*u_minus-rho_plus*u_plus\n",
    "    vdot = 1/m*(-(rho_minus*u_minus-rho_plus*u_plus)*v+rho_minus*u_minus**2-rho_plus*u_plus**2-mdot*v)\n",
    "    pdot = v\n",
    "    \n",
    "    return [mdot, vdot, pdot]\n",
    "\n",
    "t = np.linspace(1, 15, 1000)\n",
    "\n",
    "def burgers_sol(t):\n",
    "    v = 1/np.sqrt(2*(1+t))\n",
    "    p = np.sqrt(2*(1+t)) - 1\n",
    "    \n",
    "    return v, p\n",
    "\n",
    "@interact(rho_plus=FloatSlider(value=1, min=1, max=2, step=0.1))\n",
    "def solve_ode(rho_plus):\n",
    "    sol = odeint(rh_gen, [1, 0.5, 1], t, args=(rho_plus,))\n",
    "    m = sol[:, 0]\n",
    "    v = sol[:, 1]\n",
    "    p = sol[:, 2]\n",
    "    \n",
    "    v_burgers, p_burgers = burgers_sol(t)\n",
    "\n",
    "\n",
    "    fig = plt.figure()\n",
    "    fig.add_subplot(131)\n",
    "    plt.plot(t,m)\n",
    "    plt.xlabel(\"t\")\n",
    "    plt.ylabel(\"Delta-shock mass\")\n",
    "    \n",
    "    fig.add_subplot(132)\n",
    "    plt.plot(t, v)\n",
    "    plt.plot(t, v_burgers)\n",
    "    plt.xlabel(\"t\")\n",
    "    plt.ylabel(\"Delta-shock velocity\")\n",
    "    \n",
    "    \n",
    "    fig.add_subplot(133)\n",
    "    plt.plot(p,t)\n",
    "    plt.plot(p_burgers,t)\n",
    "    plt.ylabel(\"t\")\n",
    "    plt.xlabel(\"Delta-shock position\")\n",
    "    \n",
    "    fig.tight_layout()\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.geom import CircleArc, Line\n",
    "from josie.mesh import Mesh\n",
    "from josie.bc import make_periodic, Direction\n",
    "\n",
    "left = Line([0, 0], [0, 1])\n",
    "bottom = CircleArc([0, 0], [1, 0], [0.2, 0.2])\n",
    "# bottom = Line([0, 0], [1, 0])\n",
    "right = Line([1, 0], [1, 1])\n",
    "top = Line([0, 1], [1, 1])\n",
    "\n",
    "left, right = make_periodic(left, right, Direction.X)\n",
    "bottom, top = make_periodic(bottom, top, Direction.Y)\n",
    "\n",
    "mesh = Mesh(left, bottom, right, top)\n",
    "mesh.interpolate(20, 20)\n",
    "mesh.generate()\n",
    "plt.figure()\n",
    "mesh.plot()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "from josie.solver import Problem, StateTemplate, State, Solver\n",
    "from josie.mesh import Cell\n",
    "\n",
    "class Advection(Problem):\n",
    "    Q = StateTemplate(\"u\")\n",
    "    V = np.array([1, 0])\n",
    "    \n",
    "    @classmethod\n",
    "    def flux(cls, Q: State) -> np.ndarray:\n",
    "        return cls.V*Q\n",
    "    \n",
    "    \n",
    "def init(cell: Cell) -> State: \n",
    "    xc, yc = cell.centroid\n",
    "    \n",
    "    if xc > 0.45:\n",
    "        return Advection.Q(1)\n",
    "    else:\n",
    "        return Advection.Q(0)\n",
    "    \n",
    "\n",
    "    \n",
    "def upwind(cell: Cell) -> State: \n",
    "    Q = Advection.Q(0)\n",
    "\n",
    "    for neigh in cell:\n",
    "        norm = neigh.face.normal\n",
    "        flux = Advection.flux\n",
    "        S = neigh.face.surface\n",
    "\n",
    "        un = Advection.V.dot(norm)\n",
    "\n",
    "        if un >= 0:\n",
    "            Q = Q + flux(cell.value).dot(norm)*S\n",
    "        else:\n",
    "            Q = Q + flux(neigh.value).dot(norm)*S\n",
    "\n",
    "    return Q\n",
    "\n",
    "\n",
    "solver = Solver(mesh, Advection)\n",
    "solver.init(init)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.animation as animation\n",
    "\n",
    "dt = 0.01\n",
    "time = np.arange(0, 2, dt)\n",
    "fig = plt.figure()\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u')\n",
    "ims = []\n",
    "for i, t in enumerate(time):\n",
    "    solver.step(dt, upwind)\n",
    "    \n",
    "    # Pick one line\n",
    "    data = [(np.array(cell.centroid), cell.value) for cell in solver.mesh.cells[:, -1]]\n",
    "    data = np.asarray(data)  \n",
    "    x = np.array([x  for x, y in data[:, 0] ])\n",
    "    val = np.concatenate(data[:, 1])\n",
    "    ims.append(plt.plot(x, val, 'ko-'))\n",
    "    ani = animation.ArtistAnimation(fig, ims, interval=10)\n",
    "ani\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Numerical resolution of PGD\n",
    "We propose the solution of the following Riemann problem:\n",
    "$$\n",
    "\\rho^0(x) =\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\t\\rho_l & & x&<0 \\\\\n",
    "\t\\rho_r & & x&\\geq 0\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "$$\n",
    "u^0(x) =\n",
    "\\left\\{\n",
    "\\begin{aligned}\n",
    "\tu_l & & x&<0 \\\\\n",
    "\tu_r & & x&\\geq 0\n",
    "\\end{aligned}\n",
    "\\right.\n",
    "$$\n",
    "\n",
    "with $u_l>u_r$\n",
    "\n",
    "for the PGD equation.\n",
    "\n",
    "## Numerical method\n",
    "In order to solve the PGD equation a Godunov scheme is considered.  Given:\n",
    "$$\n",
    "\tw(t)=\\sqrt{\\rho_l\\rho_r}\\left(u_l-u_r\\right)t\n",
    "$$\n",
    "$$\n",
    "\tu_{\\delta}=\\frac{\\sqrt{\\rho_l}u_l+\\sqrt{\\rho_r}u_r}{\\sqrt{\\rho_l}+\\sqrt{\\rho_r}}\n",
    "$$\n",
    "\n",
    "We can derive the Godunov's fluxes:\n",
    "\n",
    "\n",
    "* if $u_l\\geq u_r$:\n",
    "\t$$\n",
    "\tF_r =\n",
    "\t\\left\\{\n",
    "\t\\begin{aligned}\n",
    "\t\t\t\\left(\\rho_l u_l, \\rho_l u_l^2\\right) & & u_{\\delta}&>0 \\\\\n",
    "\t\t\t\\left(\\rho_r u_r, \\rho_r u_r^2\\right) & & u_{\\delta}&<0 \\\\\n",
    "              \t\t\\left(\\left(\\rho_r u_r+\\rho_l u_l\\right)/2, \\rho_r u_r^2\\right) & & u_{\\delta}&=0 \t\t\t\\end{aligned}\n",
    "       \\right.\n",
    "\t$$\n",
    "* $u_l< u_r$:\n",
    "\t$$\n",
    "\tF_r =\n",
    "\t\\left\\{\n",
    "\t\\begin{aligned}\n",
    "\t\t\t\\left(\\rho_l u_l, \\rho_l u_l^2\\right) & & u_{l}&>0 \\\\\n",
    "\t\t\t\\left(\\rho_r u_r, \\rho_r u_r^2\\right) & & u_{r}&<0 \\\\\n",
    "\t\t        (0, 0) & & &\\text{otherwise} \t\t\t\\end{aligned}\n",
    "       \\right.\n",
    "\t$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.geom import CircleArc, Line\n",
    "from josie.mesh import Mesh\n",
    "from josie.bc import make_periodic, Direction\n",
    "\n",
    "left = Line([-2, 0], [-2, 1])\n",
    "bottom = Line([-2, 0], [4, 0])\n",
    "right = Line([4, 0], [4, 1])\n",
    "top = Line([-2, 0], [4, 1])\n",
    "\n",
    "left, right = make_periodic(left, right, Direction.X)\n",
    "top.bc = None\n",
    "bottom.bc = None\n",
    "\n",
    "mesh = Mesh(left, bottom, right, top)\n",
    "mesh.interpolate(300, 1)\n",
    "mesh.generate()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from josie.solver import Problem, StateTemplate, State, Solver\n",
    "from josie.mesh import Cell\n",
    "\n",
    "class PGD(Problem):\n",
    "    Q = StateTemplate(\"rho\", \"rhoU\", \"rhoV\")\n",
    "    @classmethod\n",
    "    def flux(cls, Q: State) -> np.ndarray:\n",
    "        U = Q.rhoU/Q.rho\n",
    "        V = Q.rhoV/Q.rho\n",
    "        \n",
    "        return np.array([\n",
    "            [Q.rhoU, Q.rhoV],\n",
    "            [Q.rhoU*U, Q.rhoU*V],\n",
    "            [Q.rhoV*U, Q.rhoV*V]\n",
    "        ])\n",
    "    \n",
    "    \n",
    "def init(cell: Cell) -> State: \n",
    "    xc, yc = cell.centroid\n",
    "    \n",
    "    if xc > 1:\n",
    "        rho = 1\n",
    "    else: \n",
    "        rho = 1\n",
    "            \n",
    "    \n",
    "    if (xc >= -1) and (xc <= 0):\n",
    "        U = 1 + xc\n",
    "        V = 0\n",
    "    elif (xc > 0) and (xc <= 1):\n",
    "        U = 1 - xc\n",
    "        V = 0\n",
    "    else:\n",
    "        U = 0\n",
    "        V = 0\n",
    "\n",
    "    return PGD.Q(rho, rho*U, rho*V)\n",
    "\n",
    "solver = Solver(mesh, PGD)\n",
    "solver.init(init)\n",
    "\n",
    "# Pick one line\n",
    "data = [(np.array(cell.centroid), cell.value) for cell in solver.mesh.cells[:, -1]]\n",
    "data = np.asarray(data)\n",
    "x = np.array([x  for x, y in data[:, 0] ])\n",
    "u = [s.rhoU/s.rho for s in data[:, 1]]\n",
    "rho = [s.rho for s in data[:, 1]]\n",
    "fig = plt.figure()\n",
    "fig.add_subplot(121)\n",
    "plt.plot(x, u, 'ko-')\n",
    "fig.add_subplot(122)\n",
    "plt.plot(x, rho, 'ko-')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def godunov(cell):\n",
    "    Q = solver.problem.Q(0, 0, 0)\n",
    "    \n",
    "    rho_cell = cell.value.rho\n",
    "    rhoU_cell = cell.value.rhoU\n",
    "    rhoV_cell = cell.value.rhoV\n",
    "    \n",
    "    U_cell = rhoU_cell/rho_cell\n",
    "    V_cell = rhoV_cell/rho_cell\n",
    "    \n",
    "    UU_cell = np.array([U_cell, V_cell])\n",
    "\n",
    "    for neigh in cell:        \n",
    "        norm = neigh.face.normal\n",
    "        flux = solver.problem.flux\n",
    "        S = neigh.face.surface\n",
    "        \n",
    "        rho_neigh = neigh.value.rho\n",
    "        rhoU_neigh = neigh.value.rhoU\n",
    "        rhoV_neigh = neigh.value.rhoV\n",
    "    \n",
    "        U_neigh = rhoU_neigh/rho_neigh\n",
    "        V_neigh = rhoV_neigh/rho_neigh\n",
    "        \n",
    "        UU_neigh = np.array([U_neigh, V_neigh])\n",
    "        \n",
    "        un_cell = UU_cell.dot(norm)\n",
    "        un_neigh = UU_neigh.dot(norm)\n",
    "        \n",
    "        sqrtrho_cell = np.sqrt(rho_cell)\n",
    "        sqrtrho_neigh = np.sqrt(rho_neigh)\n",
    "        \n",
    "        num = sqrtrho_cell*un_cell + sqrtrho_neigh*un_neigh\n",
    "        den = sqrtrho_cell + sqrtrho_neigh\n",
    "                \n",
    "        u_delta = num/den\n",
    "\n",
    "        if un_cell >= un_neigh:\n",
    "            if u_delta > 0:\n",
    "                Q = Q + flux(cell.value).dot(norm)*S  # (rhoU(j), rhoU(j)U(j), rhoV(j)U(j))\n",
    "            elif u_delta < 0:\n",
    "                Q = Q + flux(neigh.value).dot(norm)*S\n",
    "            else:\n",
    "                Q_cell = flux(cell.value).dot(norm)*S\n",
    "                Q_neigh = flux(neigh.value).dot(norm)*S\n",
    "                Q_mean = (Q_cell + Q_neigh)/2\n",
    "                \n",
    "                # assert np.array_equal(Q_cell, Q_neigh)\n",
    "                                \n",
    "                Q = Q + np.array([Q_mean[0], Q_cell[1], Q_cell[2]])     \n",
    "                # Q = Q + Q_mean\n",
    "        else:\n",
    "            if un_cell > 0:\n",
    "                Q = Q + flux(cell.value).dot(norm)*S\n",
    "            elif un_neigh < 0:\n",
    "                Q = Q + flux(neigh.value).dot(norm)*S\n",
    "            else:\n",
    "                pass\n",
    "\n",
    "\n",
    "    return Q\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def kinetic(cell):\n",
    "    Q = solver.problem.Q(0, 0, 0)\n",
    "    \n",
    "    rho_cell = cell.value.rho\n",
    "    rhoU_cell = cell.value.rhoU\n",
    "    rhoV_cell = cell.value.rhoV\n",
    "    \n",
    "    U_cell = rhoU_cell/rho_cell\n",
    "    V_cell = rhoV_cell/rho_cell\n",
    "    \n",
    "    UU_cell = np.array([U_cell, V_cell])\n",
    "\n",
    "    for neigh in cell:\n",
    "        # Neighbour state\n",
    "        rho_neigh = neigh.value.rho\n",
    "        rhoU_neigh = neigh.value.rhoU\n",
    "        rhoV_neigh = neigh.value.rhoV\n",
    "        \n",
    "        U_neigh = rhoU_neigh/rho_neigh\n",
    "        V_neigh = rhoV_neigh/rho_neigh\n",
    "        UU_neigh = np.array([U_neigh, V_neigh])\n",
    "        \n",
    "        # Geometry stuff\n",
    "        norm = neigh.face.normal\n",
    "        flux = solver.problem.flux\n",
    "        S = neigh.face.surface\n",
    "        \n",
    "        un_cell = UU_cell.dot(norm)\n",
    "        un_neigh = UU_neigh.dot(norm)\n",
    "        \n",
    "        F_plus = (un_cell > 0)*flux(cell.value)\n",
    "        F_minus = (un_neigh < 0)*flux(neigh.value)\n",
    "        \n",
    "        Q = Q + (F_plus + F_minus).dot(norm)*S\n",
    "            \n",
    "    return Q"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.animation as animation\n",
    "\n",
    "t = 0\n",
    "dt = 0.01\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u')\n",
    "ims = []\n",
    "solver.init(init)\n",
    "\n",
    "def burgers_sol(t):\n",
    "    def pos(t):\n",
    "        return np.sqrt(2*(1+t)) - 1\n",
    "    def deltaU(t):\n",
    "        return np.sqrt(2/(1+t))\n",
    "    \n",
    "    if t<1:\n",
    "        return ((-1,0), (t,1), (1,0))\n",
    "    else:\n",
    "        return ((-1,0), (pos(t),deltaU(t)), (pos(t),0))\n",
    "        \n",
    "    \n",
    "for i in range(1000):\n",
    "    print(f'[{i}] - t:{t}, dt:{dt}')\n",
    "    solver.step(dt, kinetic)\n",
    "\n",
    "    \n",
    "    # Pick one line\n",
    "    data = [(np.array(cell.centroid), cell.value) for cell in solver.mesh.cells[:, -1]]\n",
    "    data = np.asarray(data)\n",
    "    x = np.array([x  for x, y in data[:, 0] ])\n",
    "    \n",
    "    rho = [s.rho for s in data[:, 1]]\n",
    "    im1, = ax1.plot(x, rho, 'k-')\n",
    "    \n",
    "    U = [s.rhoU/s.rho for s in data[:, 1]]\n",
    "    A, B, C = burgers_sol(t)\n",
    "    \n",
    "    im2, = ax2.plot(x, U, 'k-')\n",
    "    im3, = ax2.plot([A[0],B[0],C[0]], [A[1],B[1],C[1]], 'r-')\n",
    "    ims.append([im1, im2, im3])\n",
    "    t = t + dt\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, ims, interval=1)\n",
    "ani"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def init_2(cell: Cell) -> State: \n",
    "    xc, yc = cell.centroid\n",
    "    \n",
    "    if xc > 1:\n",
    "        rho = 2\n",
    "    else: \n",
    "        rho = 1\n",
    "            \n",
    "    \n",
    "    if (xc >= -1) and (xc <= 0):\n",
    "        U = 1 + xc\n",
    "        V = 0\n",
    "    elif (xc > 0) and (xc <= 1):\n",
    "        U = 1 - xc\n",
    "        V = 0\n",
    "    else:\n",
    "        U = 0\n",
    "        V = 0\n",
    "\n",
    "    return PGD.Q(rho, rho*U, rho*V)\n",
    "\n",
    "t = 0\n",
    "dt = 0.01\n",
    "fig = plt.figure()\n",
    "ax1 = fig.add_subplot(121)\n",
    "ax2 = fig.add_subplot(122)\n",
    "plt.xlabel('x')\n",
    "plt.ylabel('u')\n",
    "ims = []\n",
    "solver.init(init_2)\n",
    "\n",
    "def burgers_sol(t):\n",
    "    def pos(t):\n",
    "        return np.sqrt(2*(1+t)) - 1\n",
    "    def deltaU(t):\n",
    "        return np.sqrt(2/(1+t))\n",
    "    \n",
    "    if t<1:\n",
    "        return ((-1,0), (t,1), (1,0))\n",
    "    else:\n",
    "        return ((-1,0), (pos(t),deltaU(t)), (pos(t),0))\n",
    "        \n",
    "    \n",
    "for i in range(1000):\n",
    "    print(f'[{i}] - t:{t}, dt:{dt}')\n",
    "    solver.step(dt, godunov)\n",
    "\n",
    "    \n",
    "    # Pick one line\n",
    "    data = [(np.array(cell.centroid), cell.value) for cell in solver.mesh.cells[:, -1]]\n",
    "    data = np.asarray(data)\n",
    "    x = np.array([x  for x, y in data[:, 0] ])\n",
    "    \n",
    "    rho = [s.rho for s in data[:, 1]]\n",
    "    im1, = ax1.plot(x, rho, 'k-')\n",
    "    \n",
    "    U = [s.rhoU/s.rho for s in data[:, 1]]\n",
    "    A, B, C = burgers_sol(t)\n",
    "    \n",
    "    im2, = ax2.plot(x, U, 'k-')\n",
    "    im3, = ax2.plot([A[0],B[0],C[0]], [A[1],B[1],C[1]], 'r-')\n",
    "    ims.append([im1, im2, im3])\n",
    "    t = t + dt\n",
    "\n",
    "ani = animation.ArtistAnimation(fig, ims, interval=1)\n",
    "ani"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}