{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 4: DGSEM with flux differencing"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "**Note:** To improve responsiveness via caching, the notebooks are updated only once a week. They are only\n",
    "available for the latest stable release of Trixi.jl at the time of caching."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This tutorial starts with a presentation of the weak formulation of the discontinuous Galerkin\n",
    "spectral element method (DGSEM) in order to fix the notation of the used operators.\n",
    "Then, the DGSEM formulation with flux differencing (split form DGSEM) and its implementation in\n",
    "[Trixi.jl](https://github.com/trixi-framework/Trixi.jl) is shown."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We start with the one-dimensional conservation law\n",
    "$$\n",
    "u_t + f(u)_x = 0, \\qquad t\\in \\mathbb{R}^+, x\\in\\Omega\n",
    "$$\n",
    "with the physical flux $f$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We split the domain $\\Omega$ into elements $K$ with center $x_K$ and size $\\Delta x$. With the\n",
    "transformation mapping $x(\\xi)=x_K + \\frac{\\Delta x}{2} \\xi$ we can transform the reference element\n",
    "$[-1,1]$ to every physical element. So, the equation can be restricted to the reference element using the\n",
    "determinant of the Jacobian matrix of the transformation mapping\n",
    "$J=\\frac{\\partial x}{\\partial \\xi}=\\frac{\\Delta x}{2}$.\n",
    "$$\n",
    "J u_t + f(u)_{\\xi} = 0, \\qquad t\\in \\mathbb{R}^+, \\xi\\in [-1,1]\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## The weak form of the DGSEM\n",
    "We consider the so-called discontinuous Galerkin spectral element method (DGSEM) with collocation.\n",
    "It results from choosing a nodal DG ansatz using $N+1$ Gauss-Lobatto nodes $\\xi_i$ in $[-1,1]$\n",
    "with matching interpolation weights $w_i$, which are used for numerical integration and interpolation with\n",
    "the Lagrange polynomial basis $l_i$ of degree $N$. The Lagrange functions are created with those nodes and\n",
    "hence fulfil a Kronecker property at the GL nodes.\n",
    "The weak formulation of the DGSEM for one element is\n",
    "$$\n",
    "J \\underline{\\dot{u}}(t) = - M^{-1} B \\underline{f}^* + M^{-1} D^T M \\underline{f}\n",
    "$$\n",
    "where $\\underline{u}=(u_0, u_1, \\dots, u_N)^T\\in\\mathbb{R}^{N+1}$ is the collected pointwise evaluation\n",
    "of $u$ at the discretization nodes and $\\dot{u} = \\partial u / \\partial t = u_t$ is the temporal derivative.\n",
    "The nodal values of the flux function $f$ results with collocation in $\\underline{f}$, since\n",
    "$\\underline{f}_j=f(\\underline{u}_j)$. Moreover, we got the numerical flux $f^*=f^*(u^-, u^+)$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will now have a short overview over the operators we used."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The **derivative matrix** $D\\in\\mathbb{R}^{(N+1)\\times (N+1)}$ mimics a spatial derivation on a\n",
    "discrete level with $\\underline{f}_x \\approx D \\underline{f}$. It is defined by $D_{ij} = l_j'(\\xi_i)$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The diagonal **mass matrix** $M$ is defined by $M_{ij}=\\langle l_j, l_i\\rangle_N$ with the\n",
    "numerical scalar product $\\langle \\cdot, \\cdot\\rangle_N$ defined for functions $f$ and $g$ by\n",
    "$$\n",
    "\\langle f, g\\rangle_N := \\int_{-1, N}^1 f(\\xi) g(\\xi) d\\xi := \\sum_{k=0}^N f(\\xi_k) g(\\xi_k) w_k.\n",
    "$$\n",
    "The multiplication by $M$ matches a discrete integration\n",
    "$$\n",
    "  \\int_{-1}^1 f(\\xi) \\underline{l}(\\xi) d\\xi \\approx M \\underline{f},\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The **boundary matrix** $B=\\text{diag}([-1, 0,..., 0, 1])$ represents an evaluation of a\n",
    "function at the boundaries $\\xi_0=-1$ and $\\xi_N=1$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For these operators the following property holds:\n",
    "$$\n",
    "  M D + (M D)^T = B.\n",
    "$$\n",
    "This is called the summation-by-parts (SBP) property since it mimics integration by parts on a\n",
    "discrete level ([Gassner (2013)](https://doi.org/10.1137/120890144))."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The explicit definitions of the operators and the construction of the 1D algorithm can be found\n",
    "for instance in the tutorial introduction to DG methods\n",
    "or in more detail in [Kopriva (2009)](https://link.springer.com/book/10.1007/978-90-481-2261-5)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This property shows the equivalence between the weak form and the following strong formulation\n",
    "of the DGSEM.\n",
    "$$\n",
    "\\begin{align*}\n",
    "J \\underline{\\dot{u}}(t)\n",
    "&= - M^{-1} B \\underline{f}^* + M^{-1} D^T M \\underline{f}\\\\[5pt]\n",
    "&= - M^{-1} B \\underline{f}^* + M^{-1} (B - MD) \\underline{f}\\\\[5pt]\n",
    "&= - M^{-1} B (\\underline{f}^* - \\underline{f}) - D \\underline{f}\n",
    "\\end{align*}\n",
    "$$\n",
    "More information about the equivalence you can find in [Kopriva, Gassner (2010)](https://doi.org/10.1007/s10915-010-9372-3)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## DGSEM with flux differencing\n",
    "When using the diagonal SBP property it is possible to rewrite the application of the derivative\n",
    "operator $D$ in the calculation of the volume integral into a subcell based finite volume type\n",
    "differencing formulation ([Fisher, Carpenter (2013)](https://doi.org/10.1016/j.jcp.2013.06.014)).\n",
    "Generalizing\n",
    "$$\n",
    "(D \\underline{f})_i = \\sum_j D_{i,j} \\underline{f}_j\n",
    "= 2\\sum_j \\frac{1}{2} D_{i,j} (\\underline{f}_j + \\underline{f}_i)\n",
    "\\eqqcolon 2\\sum_j  D_{i,j} f_\\text{central}(u_i, u_j),\n",
    "$$\n",
    "we replace $D \\underline{f}$ in the strong form by $2D \\underline{f}_{vol}(u^-, u^+)$ with\n",
    "the consistent two-point volume flux $f_{vol}$ and receive the DGSEM formulation with flux differencing\n",
    "(split form DGSEM) ([Gassner, Winters, Kopriva (2016)](https://doi.org/10.1016/j.jcp.2016.09.013))."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "J \\underline{\\dot{u}}(t) &= - M^{-1} B (\\underline{f}^* - \\underline{f}) - 2D \\underline{f}_{vol}(u^-, u^+)\\\\[5pt]\n",
    "&= - M^{-1} B (\\underline{f}^* - \\underline{f}_{vol}(\\underline{u}, \\underline{u})) - 2D \\underline{f}_{vol}(u^-, u^+)\\\\[5pt]\n",
    "&= - M^{-1} B \\underline{f}_{surface}^* - (2D - M^{-1} B) \\underline{f}_{vol}\\\\[5pt]\n",
    "&= - M^{-1} B \\underline{f}_{surface}^* - D_{split} \\underline{f}_{vol}\n",
    "\\end{align*}\n",
    "$$\n",
    "This formulation is in a weak form type formulation and can be implemented by using the derivative\n",
    "split matrix $D_{split}=(2D-M^{-1}B)$ and two different fluxes. We divide between the surface\n",
    "flux $f=f_{surface}$ used for the numerical flux $f_{surface}^*$ and the already mentioned volume\n",
    "flux $f_{vol}$ especially for this formulation."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This formulation creates a more stable version of DGSEM, because it fulfils entropy stability.\n",
    "Moreover it allows the construction of entropy conserving discretizations without relying on\n",
    "exact integration. This is achieved when using a two-point entropy conserving flux function as\n",
    "volume flux in the volume flux differencing formulation.\n",
    "Then, the numerical surface flux can be used to control the dissipation of the discretization and to\n",
    "guarantee decreasing entropy, i.e. entropy stability."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Implementation in Trixi.jl\n",
    "Now, we have a look at the implementation of DGSEM with flux differencing with [Trixi.jl](https://github.com/trixi-framework/Trixi.jl)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq, Trixi"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We implement a simulation for the compressible Euler equations in 2D\n",
    "$$\n",
    "\\partial_t \\begin{pmatrix} \\rho \\\\ \\rho v_1 \\\\ \\rho v_2 \\\\ \\rho e \\end{pmatrix}\n",
    "+ \\partial_x \\begin{pmatrix} \\rho v_1 \\\\ \\rho v_1^2 + p \\\\ \\rho v_1 v_2 \\\\ (\\rho e +p) v_1 \\end{pmatrix}\n",
    "+ \\partial_y \\begin{pmatrix} \\rho v_2 \\\\ \\rho v_1 v_2 \\\\ \\rho v_2^2 + p \\\\ (\\rho e +p) v_2 \\end{pmatrix}\n",
    "= \\begin{pmatrix} 0 \\\\ 0 \\\\ 0 \\\\ 0 \\end{pmatrix}\n",
    "$$\n",
    "for an ideal gas with ratio of specific heats $\\gamma=1.4$.\n",
    "Here, $\\rho$ is the density, $v_1$, $v_2$ the velocities, $e$ the specific total energy and\n",
    "$$\n",
    "p = (\\gamma - 1) \\left( \\rho e - \\frac{1}{2} \\rho (v_1^2+v_2^2) \\right)\n",
    "$$\n",
    "the pressure."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "gamma = 1.4\n",
    "equations = CompressibleEulerEquations2D(gamma)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As our initial condition we will use a weak blast wave from [Hennemann, Gassner (2020)](https://arxiv.org/abs/2008.12044).\n",
    "The primitive variables are defined by\n",
    "$$\n",
    "\\begin{pmatrix} \\rho \\\\ v_1 \\\\ v_2 \\\\ p \\end{pmatrix}\n",
    "= \\begin{pmatrix} 1.0 \\\\ 0.0 \\\\ 0.0 \\\\ 1.0 \\end{pmatrix} \\text{if } \\|x\\|_2 > 0.5,\\;\n",
    "\\text{and } \\begin{pmatrix} \\rho \\\\ v_1 \\\\ v_2 \\\\ p \\end{pmatrix}\n",
    "= \\begin{pmatrix} 1.1691 \\\\ 0.1882 * \\cos(\\phi) \\\\ 0.1882 * \\sin(\\phi) \\\\ 1.245 \\end{pmatrix} \\text{else}\n",
    "$$\n",
    "with $\\phi = \\tan^{-1}(\\frac{x_2}{x_1})$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This initial condition is implemented in Trixi.jl under the name `initial_condition_weak_blast_wave`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "initial_condition = initial_condition_weak_blast_wave"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "In Trixi.jl, flux differencing for the volume integral can be implemented with\n",
    "`VolumeIntegralFluxDifferencing` using symmetric two-point volume fluxes.\n",
    "First, we set up a simulation with the entropy conserving and kinetic energy preserving\n",
    "flux `flux_ranocha` by [Hendrik Ranocha (2018)](https://cuvillier.de/en/shop/publications/7743)\n",
    "as surface and volume flux."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We will confirm the entropy conservation property numerically."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "volume_flux = flux_ranocha # = f_vol\n",
    "solver = DGSEM(polydeg=3, surface_flux=volume_flux,\n",
    "               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we implement Trixi.jl's `mesh`, `semi` and `ode` in a simple framework. For more information please\n",
    "have a look at the documentation, the basic tutorial introduction to DG methods\n",
    "or some basic elixirs."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "coordinates_min = (-2.0, -2.0)\n",
    "coordinates_max = ( 2.0,  2.0)\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max,\n",
    "                initial_refinement_level=5,\n",
    "                n_cells_max=10_000,\n",
    "                periodicity=true)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                    boundary_conditions=boundary_condition_periodic)\n",
    "\n",
    "# ODE solvers\n",
    "tspan = (0.0, 0.4)\n",
    "ode = semidiscretize(semi, tspan);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To analyse the entropy conservation of the approximation, we will use the analysis calllback\n",
    "implemented in Trixi. It provides some information about the approximation including the entropy change."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "analysis_callback = AnalysisCallback(semi, interval=100);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now run the simulation using `flux_ranocha` for both surface and volume flux."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,\n",
    "            ode_default_options()..., callback=analysis_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "A look at the change in entropy $\\sum \\partial S/\\partial U \\cdot U_t$ in the analysis callback\n",
    "confirms that the flux is entropy conserving since the change is about machine precision."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can plot the approximated solution at the time `t=0.4`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we can use for instance the dissipative flux `flux_lax_friedrichs` as surface flux\n",
    "to get an entropy stable method."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq, Trixi\n",
    "\n",
    "gamma = 1.4\n",
    "equations = CompressibleEulerEquations2D(gamma)\n",
    "\n",
    "initial_condition = initial_condition_weak_blast_wave\n",
    "\n",
    "volume_flux = flux_ranocha # = f_vol\n",
    "solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs,\n",
    "               volume_integral=VolumeIntegralFluxDifferencing(volume_flux))\n",
    "\n",
    "coordinates_min = (-2.0, -2.0)\n",
    "coordinates_max = ( 2.0,  2.0)\n",
    "mesh = TreeMesh(coordinates_min, coordinates_max,\n",
    "                initial_refinement_level=5,\n",
    "                n_cells_max=10_000,\n",
    "                periodicity=true)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,\n",
    "                                    boundary_conditions=boundary_condition_periodic)\n",
    "\n",
    "# ODE solvers\n",
    "tspan = (0.0, 0.4)\n",
    "ode = semidiscretize(semi, tspan);\n",
    "\n",
    "analysis_callback = AnalysisCallback(semi, interval=100);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now run the simulation using the volume flux `flux_ranocha` and surface flux `flux_lax_friedrichs`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, RDPK3SpFSAL49(); abstol=1.0e-6, reltol=1.0e-6,\n",
    "            ode_default_options()..., callback=analysis_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The change in entropy confirms the expected entropy stability."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Of course, you can use more than these two fluxes in Trixi. Here, we will give a short list\n",
    "of possible fluxes for the compressible Euler equations.\n",
    "For the volume flux Trixi.jl provides for example `flux_ranocha`, `flux_shima_etal`,\n",
    "`flux_chandrashekar`, `flux_kennedy_gruber`.\n",
    "As surface flux you can use all volume fluxes and additionally for instance `flux_lax_friedrichs`,\n",
    "`flux_hll`, `flux_hllc`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Package versions"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "These results were obtained using the following versions."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using InteractiveUtils\n",
    "versioninfo()\n",
    "\n",
    "using Pkg\n",
    "Pkg.status([\"Trixi\", \"OrdinaryDiffEq\", \"Plots\"],\n",
    "           mode=PKGMODE_MANIFEST)"
   ],
   "metadata": {},
   "execution_count": null
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.3",
   "language": "julia"
  }
 },
 "nbformat": 4
}