{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 5: Shock capturing with flux differencing and stage limiter"
   ],
   "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 contains a short summary of the idea of shock capturing for DGSEM with flux differencing\n",
    "and its implementation in [Trixi.jl](https://github.com/trixi-framework/Trixi.jl).\n",
    "In the second part, an implementation of a positivity preserving limiter is added to the simulation."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Shock capturing with flux differencing"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The following rough explanation is on a very basic level. More information about an entropy stable\n",
    "shock-capturing strategy for DGSEM discretizations of advection dominated problems, such as the\n",
    "compressible Euler equations or the compressible Navier-Stokes equations, can be found in\n",
    "[Hennemann et al. (2021)](https://doi.org/10.1016/j.jcp.2020.109935). In\n",
    "[Rueda-Ramírez et al. (2021)](https://doi.org/10.1016/j.jcp.2021.110580) you find the extension to\n",
    "the systems with non-conservative terms, such as the compressible magnetohydrodynamics (MHD) equations."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The strategy for a shock-capturing method presented by Hennemann et al. is based on a hybrid blending\n",
    "of a high-order DG method with a low-order variant. The low-order subcell finite volume (FV) method is created\n",
    "directly with the Legendre-Gauss-Lobatto (LGL) nodes already used for the high-order DGSEM.\n",
    "Then, the final method is a convex combination with regulating indicator $\\alpha$ of these two methods."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since the surface integral is equal for both the DG and the subcell FV method, only the volume integral divides\n",
    "between the two methods."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This strategy for the volume integral is implemented in Trixi.jl under the name of\n",
    "`VolumeIntegralShockCapturingHG` with the three parameters of the indicator and the volume fluxes for\n",
    "the DG and the subcell FV method."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Note, that the DG method is based on the flux differencing formulation. Hence, you have to use a\n",
    "two-point flux, such as `flux_ranocha`, `flux_shima_etal`, `flux_chandrashekar` or `flux_kennedy_gruber`,\n",
    "for the DG volume flux. We would recommend to use the entropy conserving flux `flux_ranocha` by\n",
    "[Ranocha (2018)](https://cuvillier.de/en/shop/publications/7743) for the compressible Euler equations.\n",
    "````julia\n",
    "volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;\n",
    "                                                 volume_flux_dg=volume_flux_dg,\n",
    "                                                 volume_flux_fv=volume_flux_fv)\n",
    "````"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now focus on a choice of the shock capturing indicator `indicator_sc`.\n",
    "A possible indicator is $\\alpha_{HG}$ presented by Hennemann et al. (p.10), which depends on the\n",
    "current approximation with modal coefficients $\\{m_j\\}_{j=0}^N$ of a given `variable`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The indicator is calculated for every DG element by itself. First, we calculate a smooth $\\alpha$ by\n",
    "$$\n",
    "\\alpha = \\frac{1}{1+\\exp(-\\frac{-s}{\\mathbb{T}}(\\mathbb{E}-\\mathbb{T}))}\n",
    "$$\n",
    "with the total energy $\\mathbb{E}=\\max\\big(\\frac{m_N^2}{\\sum_{j=0}^N m_j^2}, \\frac{m_{N-1}^2}{\\sum_{j=0}^{N-1} m_j^2}\\big)$,\n",
    "threshold $\\mathbb{T}= 0.5 * 10^{-1.8*(N+1)^{1/4}}$ and parameter $s=ln\\big(\\frac{1-0.0001}{0.0001}\\big)\\approx 9.21024$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "For computational efficiency, $\\alpha_{min}$ is introduced and used for\n",
    "$$\n",
    "\\tilde{\\alpha} = \\begin{cases}\n",
    "0, & \\text{if } \\alpha<\\alpha_{min}\\\\\n",
    "\\alpha, & \\text{if } \\alpha_{min}\\leq \\alpha \\leq 1- \\alpha_{min}\\\\\n",
    "1, & \\text{if } 1-\\alpha_{min}<\\alpha.\n",
    "\\end{cases}\n",
    "$$"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Moreover, the parameter $\\alpha_{max}$ sets a maximal value for $\\alpha$ by\n",
    "$$\n",
    "\\alpha = \\min\\{\\tilde{\\alpha}, \\alpha_{max}\\}.\n",
    "$$\n",
    "This allows to control the maximal dissipation."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "To remove numerical artifact the final indicator is smoothed with all the neighboring elements'\n",
    "indicators. This is activated with `alpha_smooth=true`.\n",
    "$$\n",
    "\\alpha_{HG} = \\max_E \\{ \\alpha, 0.5 * \\alpha_E\\},\n",
    "$$\n",
    "where $E$ are all elements sharing a face with the current element."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Furthermore, you can specify the variable used for the calculation. For instance you can choose\n",
    "`density`, `pressure` or both with `density_pressure` for the compressible Euler equations.\n",
    "For every equation there is also the option to use the first conservation variable with `first`."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "This indicator is implemented in Trixi.jl and called `IndicatorHennemannGassner` with the parameters\n",
    "`equations`, `basis`, `alpha_max`, `alpha_min`, `alpha_smooth` and `variable`.\n",
    "````julia\n",
    "indicator_sc = IndicatorHennemannGassner(equations, basis,\n",
    "                                         alpha_max=0.5,\n",
    "                                         alpha_min=0.001,\n",
    "                                         alpha_smooth=true,\n",
    "                                         variable=variable)\n",
    "````"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Positivity preserving limiter"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Some numerical solutions are physically meaningless, for instance negative values of pressure\n",
    "or density for the compressible Euler equations. This often results in crashed simulations since\n",
    "the calculation of numerical fluxes or stable time steps uses mathematical operations like roots or\n",
    "logarithms. One option to avoid these cases are a-posteriori positivity preserving limiters.\n",
    "Trixi.jl provides the fully-discrete positivity-preserving limiter of\n",
    "[Zhang, Shu (2011)](https://doi.org/10.1098/rspa.2011.0153)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "It works the following way. For every passed (scalar) variable and for every DG element we calculate\n",
    "the minimal value $value_{min}$. If this value falls below the given threshold $\\varepsilon$,\n",
    "the approximation is slightly adapted such that the minimal value of the relevant variable lies\n",
    "now above the threshold.\n",
    "$$\n",
    "\\underline{u}^{new} = \\theta * \\underline{u} + (1-\\theta) * u_{mean}\n",
    "$$\n",
    "where $\\underline{u}$ are the collected pointwise evaluation coefficients in element $e$ and\n",
    "$u_{mean}$ the integral mean of the quantity in $e$. The new coefficients are a convex combination\n",
    "of these two values with factor\n",
    "$$\n",
    "\\theta = \\frac{value_{mean} - \\varepsilon}{value_{mean} - value_{min}},\n",
    "$$\n",
    "where $value_{mean}$ is the relevant variable evaluated for the mean value $u_{mean}$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The adapted approximation keeps the exact same mean value, but the relevant variable is now greater\n",
    "or equal the threshold $\\varepsilon$ at every node in every element."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We specify the variables the way we did before for the shock capturing variables. For the\n",
    "compressible Euler equations `density`, `pressure` or the combined variable `density_pressure`\n",
    "are a reasonable choice."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "You can implement the limiter in Trixi.jl using `PositivityPreservingLimiterZhangShu` with parameters\n",
    "`threshold` and `variables`.\n",
    "````julia\n",
    "stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=thresholds,\n",
    "                                                     variables=variables)\n",
    "````\n",
    "Then, the limiter is added to the time integration method in the `solve` function. For instance, like\n",
    "````julia\n",
    "CarpenterKennedy2N54(stage_limiter!, williamson_condition=false)\n",
    "````\n",
    "or\n",
    "````julia\n",
    "SSPRK43(stage_limiter!).\n",
    "````"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Simulation with shock capturing and positivity preserving"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we can run a simulation using the described methods of shock capturing and positivity\n",
    "preserving limiters. We want to give an example for the 2D compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq, Trixi\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As our initial condition we use the Sedov blast wave setup."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)\n",
    "  # Set up polar coordinates\n",
    "  inicenter = SVector(0.0, 0.0)\n",
    "  x_norm = x[1] - inicenter[1]\n",
    "  y_norm = x[2] - inicenter[2]\n",
    "  r = sqrt(x_norm^2 + y_norm^2)\n",
    "\n",
    "  r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)\n",
    "  # r0 = 0.5 # = more reasonable setup\n",
    "  E = 1.0\n",
    "  p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)\n",
    "  p0_outer = 1.0e-5 # = true Sedov setup\n",
    "  # p0_outer = 1.0e-3 # = more reasonable setup\n",
    "\n",
    "  # Calculate primitive variables\n",
    "  rho = 1.0\n",
    "  v1  = 0.0\n",
    "  v2  = 0.0\n",
    "  p   = r > r0 ? p0_outer : p0_inner\n",
    "\n",
    "  return prim2cons(SVector(rho, v1, v2, p), equations)\n",
    "end\n",
    "initial_condition = initial_condition_sedov_blast_wave"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "basis = LobattoLegendreBasis(3)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We set the numerical fluxes and divide between the surface flux and the two volume fluxes for the DG\n",
    "and FV method. Here, we are using `flux_lax_friedrichs` and `flux_ranocha`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "surface_flux = flux_lax_friedrichs\n",
    "volume_flux  = flux_ranocha"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we specify the shock capturing indicator $\\alpha$."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "We implement the described indicator of Hennemann, Gassner as explained above with parameters\n",
    "`equations`, `basis`, `alpha_max`, `alpha_min`, `alpha_smooth` and `variable`.\n",
    "Since density and pressure are the critical variables in this example, we use\n",
    "`density_pressure = density * pressure = rho * p` as indicator variable."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "indicator_sc = IndicatorHennemannGassner(equations, basis,\n",
    "                                         alpha_max=0.5,\n",
    "                                         alpha_min=0.001,\n",
    "                                         alpha_smooth=true,\n",
    "                                         variable=density_pressure)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Now, we can use the defined fluxes and the indicator to implement the volume integral using shock\n",
    "capturing."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;\n",
    "                                                 volume_flux_dg=volume_flux,\n",
    "                                                 volume_flux_fv=surface_flux)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We finalize the discretization by implementing Trixi.jl's `solver`, `mesh`, `semi` and `ode`,\n",
    "while `solver` now has the extra parameter `volume_integral`."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(basis, surface_flux, volume_integral)\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=6,\n",
    "                n_cells_max=10_000)\n",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)\n",
    "\n",
    "tspan = (0.0, 1.0)\n",
    "ode = semidiscretize(semi, tspan);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We add some callbacks to get an solution analysis and use a CFL-based time step size calculation."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "analysis_callback = AnalysisCallback(semi, interval=100)\n",
    "\n",
    "stepsize_callback = StepsizeCallback(cfl=0.8)\n",
    "\n",
    "callbacks = CallbackSet(analysis_callback, stepsize_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We now run the simulation using the positivity preserving limiter of Zhang and Shu for the variables\n",
    "density and pressure."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds=(5.0e-6, 5.0e-6),\n",
    "                                                     variables=(Trixi.density, pressure))\n",
    "\n",
    "sol = solve(ode, CarpenterKennedy2N54(stage_limiter!, williamson_condition=false),\n",
    "            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback\n",
    "            save_everystep=false, callback=callbacks);\n",
    "\n",
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "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
}