{
 "cells": [
  {
   "cell_type": "raw",
   "id": "c45a2da1",
   "metadata": {
    "vscode": {
     "languageId": "raw"
    }
   },
   "source": [
    "---\n",
    "title: \"Lecture 06 - Finite volumes\"\n",
    "subtitle: \"Derivation, MUSCL, Riemann Solvers\"\n",
    "author: \"Georgii Oblapenko\"\n",
    "jupyter: fenicsx\n",
    "# execute:\n",
    "#   freeze: auto\n",
    "format:\n",
    "  html:\n",
    "    output-file: 06-notes.html\n",
    "    embed-resources: true\n",
    "  pdf: default\n",
    "  revealjs:\n",
    "    output-file: 06-slides.html\n",
    "    theme: [default, ../styles/theme.scss]\n",
    "    logo: ../images/rwth_acom_en_rgb.png\n",
    "    center: false\n",
    "    smaller: true\n",
    "    slide-number: c/t\n",
    "    show-slide-number: all\n",
    "    sc-sb-title: true\n",
    "    toc: true\n",
    "    toc-depth: 2\n",
    "    toc-title: Contents\n",
    "    slide-level: 3\n",
    "    embed-resources: true\n",
    "    # navigation-mode: vertical\n",
    "format-links:\n",
    "  - pdf\n",
    "  - revealjs\n",
    "  - href: 06.ipynb\n",
    "    text: Jupyter\n",
    "    icon: journal-code\n",
    "execute:\n",
    "  echo: true\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3e67ba4",
   "metadata": {},
   "source": [
    "## Objectives\n",
    "\n",
    "In this lecture, we will learn (or revisit) \"classical\" finite volume methods for hyperbolic PDEs. We will focus on the following parts of the theory:\n",
    "\n",
    "1. Motivation, general formulation\n",
    "2. Exact and approximate Riemann solvers\n",
    "3. Higher-order reconstruction, slope limiters"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b973ac5",
   "metadata": {},
   "source": [
    "## Motivation, general formulation\n",
    "\n",
    "We want to solve hyperbolic PDEs of the form\n",
    "\n",
    "$$\n",
    "\\partial_t \\mathbf{u} + \\nabla \\cdot \\mathbf{f}(\\mathbf{u}) = \\mathbf{s}(\\mathbf{u}),\n",
    "$$\n",
    "\n",
    "where $\\mathbf{u}$ is the vector of conserved variables, $\\mathbf{f}$ is the **flux function**, and $\\mathbf{s}$ is the source term."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93c36e2d",
   "metadata": {},
   "source": [
    "### Examples: linear advection\n",
    "\n",
    "Simplest hyperbolic PDE:\n",
    "\n",
    "$$\n",
    "\\partial_t u + \\mathbf{a} \\cdot \\nabla u = 0,\n",
    "$$\n",
    "\n",
    "i.e. $\\mathbf{u} = u$, $\\mathbf{f}(\\mathbf{u}) = \\mathbf{a} u$, $\\mathbf{s}(\\mathbf{u}) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38598afc",
   "metadata": {},
   "source": [
    "### Examples: inviscid Burgers' equation\n",
    "\n",
    "Non-linear equation that can exhibit formation of shocks:\n",
    "\n",
    "$$\n",
    "\\partial_t u + \\partial_x \\left( \\frac{u^2}{2} \\right) = 0,\n",
    "$$\n",
    "i.e. $\\mathbf{u} = u$, $\\mathbf{f}(\\mathbf{u}) = \\frac{u^2}{2}$, $\\mathbf{s}(\\mathbf{u}) = 0$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6abc231a",
   "metadata": {},
   "source": [
    "### Examples: Compressible Euler equations in 2D\n",
    "\n",
    "$$\n",
    "\\partial_t \\rho + \\nabla \\cdot (\\rho \\mathbf{v}) = 0, \\\\\n",
    "\\partial_t (\\rho \\mathbf{v}) + \\nabla \\cdot (\\rho \\mathbf{v} \\otimes \\mathbf{v}) + \\nabla p = 0, \\\\\n",
    "\\partial_t E + \\nabla \\cdot \\left( (E + p) \\mathbf{v} \\right) = 0,\n",
    "$$\n",
    "\n",
    "i.e. $\\mathbf{u} = (\\rho, \\rho \\mathbf{v}, E)^T$, $\\mathbf{f}(\\mathbf{u}) = (f_x, f_y)$, $\\mathbf{s}(\\mathbf{u}) = 0$, flux components $f_x, f_y$ are given by\n",
    "\n",
    "$$\n",
    "f_x = (\\rho v_x, \\rho v_x^2 + p, \\rho v_x v_y, (E + p) v_x)^T, \\\\\n",
    "f_y = (\\rho v_y, \\rho v_x v_y, \\rho v_y^2 + p, (E + p) v_y)^T\n",
    "$$\n",
    "\n",
    "Closure:\n",
    "\n",
    "$$\n",
    "    E = \\rho \\varepsilon_{int} + \\rho \\frac{v^2}{2},\\quad p=p(\\rho, \\varepsilon_{int})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8d05f579",
   "metadata": {},
   "source": [
    "### Why methods for unsteady equations?\n",
    "\n",
    "\n",
    "::: {.content-visible when-format=\"revealjs\"}\n",
    "In steady-state $\\partial_t \\equiv 0$, the **type** of the problem depends on the regime:\n",
    "\n",
    "- subsonic: elliptic\n",
    "- supersonic: hyperbolic\n",
    "\n",
    "Unsteady approaches more general.\n",
    ":::\n",
    "\n",
    "::: {.content-visible when-format=\"html\" unless-format=\"revealjs\"}\n",
    "Even if we are interested in a steady-state solution of Euler equations, i.e. $\\partial_t \\equiv 0$, the **type** of the problem depends on the regime:\n",
    "\n",
    "- subsonic: elliptic (propagation of perturbations happens infinitely fast, boundary conditions required everywhere)\n",
    "- supersonic: hyperbolic (boundary conditions not required on supersonic outflow boundaries)\n",
    "\n",
    "Therefore, timestepping until a steady-state solution is a common approach that allows for a unified solution methodology.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18d9505d",
   "metadata": {},
   "source": [
    "### Finite volumes\n",
    "\n",
    "Compared to finite differences, FVM allows for\n",
    "\n",
    "- complex geometries and meshes\n",
    "- easier handling of discontinuities\n",
    "- more robust conservation properties"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fc3c7516",
   "metadata": {},
   "source": [
    "## Finite volume method\n",
    "\n",
    "\n",
    "$$\n",
    "\\partial_t \\mathbf{u} + \\nabla \\cdot \\mathbf{f}(\\mathbf{u}) = \\mathbf{s}(\\mathbf{u})\n",
    "$$\n",
    "\n",
    "defined on domain $V$. Decompose $V$ into control volumes:\n",
    "\n",
    "$$\n",
    "V = \\bigcup_{i=1}^N V_i,\\:\\: V_i \\cap V_j = \\emptyset,\\: i \\neq j.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1aa101c",
   "metadata": {},
   "source": [
    "### Integral conservation law\n",
    "\n",
    "Integrate PDE over control volume $V_i$:\n",
    "\n",
    "$$\n",
    "\\int_{V_i} \\partial_t \\mathbf{u} \\, \\mathrm{d}V + \\int_{V_i} \\nabla \\cdot \\mathbf{f}(\\mathbf{u}) = \\int_{V_i} \\mathbf{s}(\\mathbf{u}) \\, \\mathrm{d}V.\n",
    "$$\n",
    "\n",
    "::: {.fragment}\n",
    "Transform the flux term using the divergence theorem:\n",
    "\n",
    "$$\n",
    "\\int_{V_i} \\partial_t \\mathbf{u} \\, \\mathrm{d}V + \\int_{\\partial V_i} \\mathbf{f}(\\mathbf{u}) \\cdot \\mathbf{n} \\, \\mathrm{d}S = \\int_{V_i} \\mathbf{s}(\\mathbf{u}) \\, \\mathrm{d}V,\n",
    "$$\n",
    "\n",
    "or\n",
    "\n",
    "$$\n",
    "\\partial_t  \\int_{V_i} \\mathbf{u} \\, \\mathrm{d}V + \\int_{\\partial V_i} \\mathbf{f}(\\mathbf{u}) \\cdot \\mathbf{n} \\, \\mathrm{d}S = \\int_{V_i} \\mathbf{s}(\\mathbf{u}) \\, \\mathrm{d}V,\n",
    "$$\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "748a9a3c",
   "metadata": {},
   "source": [
    "### Fluxes at cell interfaces\n",
    "\n",
    "Split the cell boundary $\\partial V_i$ into cell faces: $\\partial V_i = \\cup_j f_{ij},\\, j\\in \\mathcal{N}(i)$, $\\mathcal{N}(i)$ is the set of neighbors of cell $i$, $f_{ij} = \\partial V_i \\cap \\partial V_j$.\n",
    "\n",
    "Then we have\n",
    "$$\n",
    "\\int_{\\partial V_i} \\mathbf{f}(\\mathbf{u}) \\cdot \\mathbf{n} = \\sum_{j\\in \\mathcal{N}(i)} \\int_{f_{ij}} \\mathbf{f}(\\mathbf{u})\\cdot \\mathbf{n} \\, \\mathrm{d}S.\n",
    "$$\n",
    "\n",
    "::: {.center}\n",
    "![Cells and cell faces](../images/lecture6/cell_interfaces.png){width=300px}\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1837801",
   "metadata": {},
   "source": [
    "### Fluxes at cell interfaces\n",
    "We can approximate the flux integral over a cell face $f_{ij}$ as\n",
    "$$\n",
    "\\int_{f_{ij}} \\mathbf{f}(\\mathbf{u})\\cdot \\mathbf{n} \\, \\mathrm{d}S \\approx g_{ij}(\\mathbf{u}_{i},\\mathbf{u}_{j}),\n",
    "$$\n",
    "\n",
    "$g_{ij}$ is a numerical approximation of the flux **from cell** $i$ **into cell** $j$.\n",
    "\n",
    "\n",
    "::: {.callout-tip .content-visible when-format=\"revealjs\"}\n",
    "## Approximations made?\n",
    "We made two approximations here, which ones?\n",
    ":::\n",
    "\n",
    "\n",
    "::: {.callout-note .content-visible when-format=\"html\" unless-format=\"revealjs\"}\n",
    "## Approximations made\n",
    "We made two approximations here: 1) we approximate integration over the interface $f_{ij}$ with a single-point quadrature 2) we are also (potentially) approximating the function (flux) we are integrating!\n",
    ":::\n",
    "\n",
    "::: {.fragment .callout-note}\n",
    "## Constraints on numerical flux\n",
    "$g_{ij}$ should satisfy two requirements:\n",
    "\n",
    "1. **Flux conservation**: $g_{ij}(\\mathbf{u}_{i},\\mathbf{u}_{j}) = -g_{ji}(\\mathbf{u}_{j},\\mathbf{u}_{i})$.\n",
    "1. **Consistency**: $g_{ij}(\\mathbf{u},\\mathbf{u}) = \\int_{f_{ij}} \\mathbf{f}(\\mathbf{u})\\cdot \\mathbf{n} \\, \\mathrm{d}S$.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f0383ee",
   "metadata": {},
   "source": [
    "### Final semi-discrete form\n",
    "\n",
    "Define\n",
    "\n",
    "$$\n",
    "\\mathbf{u}_i = \\frac{1}{V_i} \\int_{V_i} \\mathbf{u} \\, \\mathrm{d}V,\n",
    "$$\n",
    "\n",
    "then we have the **semi-discrete** form of the conservation law:\n",
    "\n",
    "$$\n",
    "\\partial_t \\mathbf{u}_i + \\frac{1}{V_i} \\sum_{j \\in \\mathcal{N}(i)} g_{ij}(\\mathbf{u}_i, \\mathbf{u}_j) = \\mathbf{s}_i.\n",
    "$$\n",
    "\n",
    "::: {.callout-note}\n",
    "## Source term treatment\n",
    "Simplest way to treat source term is to simply evaluate it from $\\mathbf{u}_i$ and incorporate it independently of the fluxes (fluxes are not directly affected by source terms).\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cfbe05b",
   "metadata": {},
   "source": [
    "### Cell-centered vs node/vertex-centered\n",
    "\n",
    "Assume we have a grid we would like to compute our solution on. What does that mean exactly?\n",
    "\n",
    "::: {.content-visible when-format=\"html\" unless-format=\"revealjs\"}\n",
    "Solution is defined at the center of the cell, but one could also reasonably expect that solution is defined at the vertices of the cell. For example, if we want to have the value at the boundary.\n",
    ":::\n",
    "\n",
    "\n",
    "::: {.center}\n",
    "![Cells and cell faces](../images/lecture6/init_grid.png){width=300px}\n",
    ":::\n",
    "\n",
    "There are two ways to take a grid and define a solution on it: **cell-centered** and **vertex-centered**."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "409a8024",
   "metadata": {},
   "source": [
    "### Cell-centered\n",
    "\n",
    "In a cell-centered approach, we define the solution at the center of the grid elements.\n",
    "\n",
    "\n",
    "::: {.center}\n",
    "![Cells and cell faces](../images/lecture6/cell_centered.png){width=500px}\n",
    ":::\n",
    "\n",
    "::: {.callout-info}\n",
    "## Pros and cons\n",
    "Direct access to cell values, easier to implement, less memory. Harder to handle boundary nodes, less robust on certain meshes.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "836bc6a5",
   "metadata": {},
   "source": [
    "### Node-centered\n",
    "\n",
    "In a vertex-centered approach, we first construct a **Voronoi**-based **dual grid**.\n",
    "\n",
    "\n",
    "::: {.center}\n",
    "![Cells and cell faces](../images/lecture6/dual_grid.png){width=400px}\n",
    ":::\n",
    "\n",
    "::: {.callout-note}\n",
    "## Pros and cons\n",
    "More robust on skewed meshes, provides more stable computation of gradients, boundary nodes are easier to handle. Harder to implement, requires more memory.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2246432",
   "metadata": {},
   "source": [
    "## Solver overview\n",
    "\n",
    "| **Cell-centered** | **Vertex-centered** |\n",
    "|--------------|-----------------|\n",
    "| OpenFOAM     | SU2             |\n",
    "| Ansys Fluent | DLR Tau         |\n",
    "| Eilmer       |                 |\n",
    "| DLR Coda     |                 |\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f693f7eb",
   "metadata": {},
   "source": [
    "## Riemann solvers\n",
    "\n",
    "We still haven't discussed how we actually compute the flux at the interface."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ee176b14",
   "metadata": {},
   "source": [
    "### Riemann problem\n",
    "\n",
    "For a 1-D hyperbolic conservation law\n",
    "\n",
    "$$\n",
    "    \\partial_t \\mathbf{u} + \\partial_x \\mathbf{f}(\\mathbf{u}) = 0\n",
    "$$\n",
    "\n",
    "define initial conditions\n",
    "\n",
    "$$\n",
    "    \\mathbf{u}(x, t=0) = \\begin{cases}\n",
    "       \\mathbf{u}_L & x < 0 \\\\\n",
    "       \\mathbf{u}_R & x > 0\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "This is called the Riemann problem, we denote the solution as $\\mathbf{u}^{\\mathcal{RP}}(x,t)$.\n",
    "\n",
    "::: {.callout-important}\n",
    "## Self-similarity\n",
    "Solutions to the Riemann problem are self-similar, i.e $\\mathbf{u}^{\\mathcal{RP}}(x,t) = const$ along $x/t = const$.\n",
    "So solution is only dependent on $x/t$; this also means that **for any** $t>0$, $\\mathbf{u}^{\\mathcal{RP}}(0,t) \\equiv const$.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daf9c438",
   "metadata": {},
   "source": [
    "### Example: linear advection\n",
    "\n",
    "Simplest textbook example:\n",
    "\n",
    "$$\n",
    "    \\partial_t u + a \\partial_x u = 0.\n",
    "$$\n",
    "\n",
    "Solution to Riemann problem is given by\n",
    "\n",
    "$$\n",
    "    u^{\\mathcal{RP}}(x,t) = \\begin{cases}\n",
    "       u_L & x/t < a, \\\\\n",
    "       u_R & x/t > a.\n",
    "    \\end{cases}\n",
    "$$\n",
    "\n",
    "In fact, any system of linear transport equations\n",
    "\n",
    "$$\n",
    "    \\partial_t \\mathbf{u} + \\mathbf{A} \\partial_x \\mathbf{u} = 0.\n",
    "$$\n",
    "\n",
    "has a simple linear solution related to the eigenvalues and eigenvectors of $\\mathbf{A}$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "423106f6",
   "metadata": {},
   "source": [
    "### Godunov's method\n",
    "\n",
    "Take our conservation law (in 1D for simplicity)\n",
    "\n",
    "$$\n",
    "\\partial_t  \\int_{x_{i-1/2}}^{x_{i+1/2}} \\mathbf{u} \\, \\mathrm{d}x + \\int_{x_{i-1/2}}^{x_{i+1/2}} \\partial_x \\mathbf{f}(\\mathbf{u}) \\mathrm{d}x = 0,\n",
    "$$\n",
    "\n",
    "and integrate from time $t^n$ to $t^{n+1}$:\n",
    "\n",
    "$$\n",
    "\\int_{t^n}^{t^{n+1}} \\partial_t  \\int_{x_{i-1/2}}^{x_{i+1/2}} \\mathbf{u} \\, \\mathrm{d}x \\mathrm{d}t + \\int_{t^n}^{t^{n+1}} \\int_{x_{i-1/2}}^{x_{i+1/2}} \\partial_x \\mathbf{f}(\\mathbf{u}) \\mathrm{d}x  \\mathrm{d}t = 0 \\implies\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\frac{\\mathbf{u}_{i}^{n+1} - \\mathbf{u}_{i}^{n}}{\\Delta t} + \\frac{1}{\\Delta x} \\int_{t^n}^{t^{n+1}} (\\mathbf{f}(\\mathbf{u}_{i+1/2}) - \\mathbf{f}(\\mathbf{u}_{i-1/2}))\\mathrm{d}t   = 0.\n",
    "$$\n",
    "\n",
    "::: {.callout-tip}\n",
    "## Question\n",
    "How do we 1) compute flux at $x=i \\pm 1/2$ (solution is discontinuous across cell interface) and 2) compute time integral of flux difference?\n",
    ":::\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e07d41d7",
   "metadata": {},
   "source": [
    "### Godunov's method\n",
    "\n",
    "If we take $\\mathbf{u} \\equiv const$ in each cell, then we have\n",
    "\n",
    "1. Riemann problem at each interface\n",
    "2. Solution of Riemann problem at $x/t = 0$ stays constant $\\implies$ flux of Riemann problem solution also stays constant in time\n",
    "\n",
    "So we can write\n",
    "$$\n",
    "\\int_{t^n}^{t^{n+1}} (\\mathbf{f}(\\mathbf{u}_{i+1/2}) - \\mathbf{f}(\\mathbf{u}_{i-1/2}))\\mathrm{d}t = \\Delta t (\\mathbf{f}(u^{\\mathcal{RP}}_{i+1/2}) - \\mathbf{f}(u^{\\mathcal{RP}}_{i-1/2}))\n",
    "$$\n",
    "\n",
    "::: {.callout-important}\n",
    "## CFL restriction\n",
    "Immediately see that we have a restriction $\\lambda \\Delta t / \\Delta x < \\frac{1}{2}$, where $\\lambda$ is fastest wavespeed in Riemann problem: i.e. Riemann problems from neighbouring cells do not interact\n",
    "within the timestep.\n",
    ":::\n",
    "\n",
    "\n",
    "::: {.callout-tip}\n",
    "## Question\n",
    "What are the issues with Godunov's method?\n",
    ":::\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c50b6296",
   "metadata": {},
   "source": [
    "\n",
    "### The birth of modern(-er) CFD\n",
    "\n",
    "\n",
    "::: {.callout-important .content-visible when-format=\"html\" unless-format=\"revealjs\"}\n",
    "## Loss of self-similarity\n",
    "If we use non-constant values within the cell, the solution to the Riemann problem is no longer self-similar.\n",
    ":::\n",
    "\n",
    "::: {.fragment}\n",
    "- Bram van Leer (1979): we can use a linear reconstruction within the cell to achieve 2-nd order accuracy\n",
    "- Philip Roe (1981): we don't need exact solution of the Riemann problem\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6cbbf9a4",
   "metadata": {},
   "source": [
    "### Roe's method\n",
    "\n",
    "The key insight is: better to solve Riemann problem approximately and gain solution accuracy **elsewhere**. Which approximation? Local linearization (at cell interface):\n",
    "\n",
    "$$\n",
    "\\partial_t \\mathbf{u} + \\nabla \\mathbf{f}(\\mathbf{u}) = 0 \\implies \\partial_t \\mathbf{u} + \\mathbf{A}(\\mathbf{u}_L, \\mathbf{u}_R)  \\nabla \\mathbf{u} = 0.\n",
    "$$\n",
    "\n",
    "Conditions on $\\mathbf{A}$ (**Roe matrix**):\n",
    "\n",
    "1. $\\mathbf{A}$ is diagonalizable with real eigenvalues (remain hyperbolic)\n",
    "2. $\\mathbf{A}$ is consistent: $\\mathbf{A}(\\mathbf{u}, \\mathbf{u}) = \\mathbf{A}(\\mathbf{u})$\n",
    "3. $\\mathbf{f}(\\mathbf{u}_R) - \\mathbf{f}(\\mathbf{u}_L) = \\mathbf{A}(\\mathbf{u}_L, \\mathbf{u}_R)(\\mathbf{u}_R - \\mathbf{u}_L)$\n",
    "\n",
    "One way to define this matrix is as $\\mathbf{A} = \\frac{\\partial \\mathbf{f}}{\\partial \\mathbf{u}} \\left( \\mathbf{u}_{avg} \\right)$, where $\\mathbf{u}_{avg}$ is **some** average of the states $\\mathbf{u}_L$ and $\\mathbf{u}_R$."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a4483af",
   "metadata": {},
   "source": [
    "### Roe's method\n",
    "\n",
    "Recipe:\n",
    "\n",
    "1. Compute $\\mathbf{u}_{avg}$, $\\mathbf{A}$\n",
    "2. $\\mathbf{A}$ is diagonalizable, i.e. $\\mathbf{A} = \\mathbf{R} \\mathbf{\\Lambda} \\mathbf{R}^{-1}$\n",
    "3. Compute $|\\mathbf{A}| := \\mathbf{R} |\\mathbf{\\Lambda}| \\mathbf{R}^{-1}$\n",
    "\n",
    "For Euler (1D):\n",
    "\n",
    "- $\\rho_{Roe} = \\sqrt{\\rho_L \\rho_R}$\n",
    "- $u_{Roe} = \\dfrac{\\sqrt{\\rho_L} u_L + \\sqrt{\\rho_R} u_R}{\\sqrt{\\rho_L} + \\sqrt{\\rho_R}}$\n",
    "- $h_{Roe} = \\dfrac{\\sqrt{\\rho_L} h_L + \\sqrt{\\rho_R} h_R}{\\sqrt{\\rho_L} + \\sqrt{\\rho_R}}$ ($h = e + p/\\rho$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d359de56",
   "metadata": {},
   "source": [
    "### Roe's method: summary\n",
    "\n",
    "- Computationally expensive\n",
    "- Not always stable (hypersonic flows):\n",
    "    - if $f'(u_L) < 0 < f'(u_R)$ (rarefaction) leads to unphysical solution\n",
    "    - \"Carbuncle phenomena\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b56c17f",
   "metadata": {},
   "source": [
    "### Other Riemann solvers\n",
    "\n",
    "- HLL (Harten, Lax, Leer): assumes we know two fastest wave speeds of Riemann problem, resolves only them: need *wavespeed estimates*, i.e. min/max eigenvalues of $\\mathbf{A}$\n",
    "    - cannot resolve contact waves\n",
    "- HLLC\n",
    "    - add contact wave (need speed estimate, etc.)\n",
    "- Hybrid methods (Roe + HLL)\n",
    "- AUSM, AUSM+\n",
    "    - high-speed flows, splitting of transport due to pressure and convection\n",
    "    - efficient (no eigendecomposition required)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "392317a6",
   "metadata": {},
   "source": [
    "## Reconstruction and limiting"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87aecce3",
   "metadata": {},
   "source": [
    "### TVD and Godunov's theorem\n",
    "\n",
    "\n",
    "::: {.callout-note}\n",
    "## Godunov's theorem\n",
    "A linear, monotonicity-preserving method can be at most first-order accurate.\n",
    ":::\n",
    "\n",
    "\n",
    "::: {.callout-tip .content-visible when-format=\"revealjs\"}\n",
    "## Question\n",
    "Is all lost?\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d733e441",
   "metadata": {},
   "source": [
    "### MUSCL\n",
    "\n",
    "- Reconstruct piecewise linear solution from neighbouring cell averages\n",
    "- Use values at cell interfaces to compute fluxes\n",
    "\n",
    "$$\n",
    "\\mathbf{u}_L = \\mathbf{u}_i + \\nabla \\mathbf{u}_i \\cdot (\\mathbf{r}_m - \\mathbf{r}_i),\\quad \\mathbf{u}_R = \\mathbf{u}_j + \\nabla \\mathbf{u}_j \\cdot (\\mathbf{r}_m - \\mathbf{r}_j)\n",
    "$$\n",
    "\n",
    "- Gradient computed either based on linear interpolation or least-squares\n",
    "\n",
    "::: {.fragment}\n",
    "::: {.callout-tip}\n",
    "## Question\n",
    "We use a Riemann solver and don't account for the gradient of the solution in the cell. Why do we still gain accuracy?\n",
    ":::\n",
    ":::\n",
    "\n",
    "\n",
    "::: {.fragment}\n",
    "::: {.callout-note}\n",
    "## Taylor series expansion analysis\n",
    "For a smooth scalar field $u$:\n",
    "$$\n",
    "u(x_{i+1/2}) = u_i + \\frac{\\partial u}{\\partial x}\\bigg|_i \\Delta x + \\frac{1}{2} \\frac{\\partial^2 u}{\\partial x^2}\\bigg|_i \\Delta x^2 + \\dots\n",
    "$$\n",
    "Cell-wise constant values disregard the gradient and higher-order terms; Riemann solver adds diffusion proportional to jump at interface - so even though we omit the gradient information within the Riemann solver, we gain accuracy in the other components of the method.\n",
    ":::\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cda6d0c6",
   "metadata": {},
   "source": [
    "### Slope limiting\n",
    "\n",
    "Replace $\\mathbf{u}_i + \\nabla \\mathbf{u}_i \\cdot (\\mathbf{r}_m - \\mathbf{r}_i)$ with $\\mathbf{u}_i + \\Psi \\nabla \\mathbf{u}_i \\cdot (\\mathbf{r}_m - \\mathbf{r}_i)$.\n",
    "\n",
    "$\\Psi = 0$ recovers 1st order scheme, $\\Psi$ - **slope limiter**. Restricts local solution slope to avoid oscillations.\n",
    "\n",
    ":::{.callout-note}\n",
    "### Structured vs unstructured limiters\n",
    "\n",
    "On structured grids (SG), limiter is function of local slope in specific direction. On unstructured grids (UG), we need to consider all directions and apply **strongest** limiter. Many more limiters for SG (minmod, Superbee, van Leer, van Albada) than UG (Barth-Jespersen, Venkatakrishnan, Nishikawa).\n",
    ":::\n",
    "\n",
    "::: {.callout-note}\n",
    "## Flux limiting\n",
    "Alternative approach is to limit values of fluxes $\\mathbf{g}_{ij}$ (and not values of $\\mathbf{u}$ used to compute fluxes). Flux limiting is harder to construct in multidimensional cases/on unstructured grids,\n",
    "but can be shown to be equivalent to slope limiting.\n",
    ":::"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "549e3b2b",
   "metadata": {},
   "source": [
    "## Examples\n",
    "\n",
    "[Online 1D Euler solver](https://euler-1d-solver.streamlit.app/)\n",
    "\n",
    "[Github](https://github.com/fhermet/euler-1d-solver/tree/main)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "487f2fa7",
   "metadata": {},
   "source": [
    "## OpenFOAM\n",
    "\n",
    "`rhoCentralFoam` - compressible Euler/Navier--Stokes equations.\n",
    "\n",
    "Some options:\n",
    "\n",
    "- `ddtSchemes`: `Euler`, `backward`, `CrankNicolson`\n",
    "- `gradSchemes`, `divSchemes`, `interpolationSchemes`: `Gauss`/`leastSquares`, `faceMDLimited`/`cellMDLimited`, `linear`/`vanLeer`/`limitedLinear` (minmod limiter):\n",
    "- `fluxScheme`: `Kurganov`, `HLLC`, `HLL`\n",
    "\n",
    "\n",
    "```\n",
    "gradSchemes\n",
    "{\n",
    "    default cellMDLimited Gauss linear 1;\n",
    "}\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8daa188",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "- [R. LeVeque, Numerical Methods for Conservation Laws](https://link.springer.com/book/10.1007/978-3-0348-8629-1)\n",
    "- [J. Quirk, A contribution to the great Riemann solver debate](https://doi.org/10.1002/fld.1650180603)\n",
    "- [H. Nishikawa et al., Very simple, carbuncle-free, boundary-layer-resolving, rotated-hybrid Riemann solvers, J. Comp. Phys.](https://doi.org/10.1016/j.jcp.2007.11.003)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0794a145",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}