{
 "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_\\text{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}^\\text{new} = \\theta * \\underline{u} + (1-\\theta) * u_\\text{mean}\n",
    "$$\n",
    "where $\\underline{u}$ are the collected pointwise evaluation coefficients in element $e$ and\n",
    "$u_\\text{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_\\text{mean} - \\varepsilon}{value_\\text{mean} - value_\\text{min}},\n",
    "$$\n",
    "where $value_\\text{mean}$ is the relevant variable evaluated for the mean value $u_\\text{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 OrdinaryDiffEqLowStorageRK\n",
    "using 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",
    "            ode_default_options()..., callback = callbacks);\n",
    "\n",
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Entropy bounded limiter"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As argued in the description of the positivity preserving limiter above it might sometimes be\n",
    "necessary to apply advanced techniques to ensure a physically meaningful solution.\n",
    "Apart from the positivity of pressure and density, the physical entropy of the system should increase\n",
    "over the course of a simulation, see e.g. [this](https://doi.org/10.1016/0168-9274(86)90029-2) paper by Tadmor where this property is\n",
    "shown for the compressible Euler equations.\n",
    "As this is not necessarily the case for the numerical approximation (especially for the high-order, non-diffusive DG discretizations),\n",
    "Lv and Ihme devised an a-posteriori limiter in [this paper](https://doi.org/10.1016/j.jcp.2015.04.026) which can be applied after each Runge-Kutta stage.\n",
    "This limiter enforces a non-decrease in the physical, thermodynamic entropy $S$\n",
    "by bounding the entropy decrease (entropy increase is always tolerated) $\\Delta S$ in each grid cell.\n",
    "\n",
    "This translates into a requirement that the entropy of the limited approximation $S\\Big(\\mathcal{L}\\big[\\boldsymbol u(\\Delta t) \\big] \\Big)$ should be\n",
    "greater or equal than the previous iterates' entropy $S\\big(\\boldsymbol u(0) \\big)$, enforced at each quadrature point:\n",
    "$$\n",
    "S\\Big(\\mathcal{L}\\big[\\boldsymbol u(\\Delta t, \\boldsymbol{x}_i) \\big] \\Big) \\overset{!}{\\geq} S\\big(\\boldsymbol u(0, \\boldsymbol{x}_i) \\big), \\quad i = 1, \\dots, (k+1)^d\n",
    "$$\n",
    "where $k$ denotes the polynomial degree of the element-wise solution and $d$ is the spatial dimension.\n",
    "For an ideal gas, the thermodynamic entropy $S(\\boldsymbol u) = S(p, \\rho)$ is given by\n",
    "$$\n",
    "S(p, \\rho) = \\ln \\left( \\frac{p}{\\rho^\\gamma} \\right) \\: .\n",
    "$$\n",
    "Thus, the non-decrease in entropy can be reformulated as\n",
    "$$\n",
    "p(\\boldsymbol{x}_i) - e^{ S\\big(\\boldsymbol u(0, \\boldsymbol{x}_i) \\big)} \\cdot \\rho(\\boldsymbol{x}_i)^\\gamma \\overset{!}{\\geq} 0, \\quad i = 1, \\dots, (k+1)^d \\: .\n",
    "$$\n",
    "In a practical simulation, we might tolerate a maximum (exponentiated) entropy decrease per element, i.e.,\n",
    "$$\n",
    "\\Delta e^S \\coloneqq \\min_{i} \\left\\{ p(\\boldsymbol{x}_i) - e^{ S\\big(\\boldsymbol u(0, \\boldsymbol{x}_i) \\big)} \\cdot \\rho(\\boldsymbol{x}_i)^\\gamma \\right\\} < c\n",
    "$$\n",
    "with hyper-parameter $c$ which is to be specified by the user.\n",
    "The default value for the corresponding parameter $c=$ `exp_entropy_decrease_max` is set to $-10^{-13}$, i.e., slightly less than zero to\n",
    "avoid spurious limiter actions for cells in which the entropy remains effectively constant.\n",
    "Other values can be specified by setting the `exp_entropy_decrease_max` keyword in the constructor of the limiter:\n",
    "```julia\n",
    "stage_limiter! = EntropyBoundedLimiter(exp_entropy_decrease_max=-1e-9)\n",
    "```\n",
    "Smaller values (larger in absolute value) for `exp_entropy_decrease_max` relax the entropy increase requirement and are thus less diffusive.\n",
    "On the other hand, for larger values (smaller in absolute value) of `exp_entropy_decrease_max` the limiter acts more often and the solution becomes more diffusive.\n",
    "\n",
    "In particular, we compute again a limiting parameter $\\vartheta \\in [0, 1]$ which is then used to blend the\n",
    "unlimited nodal values $\\boldsymbol u$ with the mean value $\\boldsymbol u_{\\text{mean}}$ of the element:\n",
    "$$\n",
    "\\mathcal{L} [\\boldsymbol u](\\vartheta) \\coloneqq (1 - \\vartheta) \\boldsymbol u + \\vartheta \\cdot \\boldsymbol u_{\\text{mean}}\n",
    "$$\n",
    "For the exact definition of $\\vartheta$ the interested user is referred to section 4.4 of the paper by Lv and Ihme.\n",
    "Note that therein the limiting parameter is denoted by $\\epsilon$, which is not to be confused with the threshold $\\varepsilon$ of the Zhang-Shu limiter."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As for the positivity preserving limiter, the entropy bounded limiter may be applied after every Runge-Kutta stage.\n",
    "Both fixed timestep methods such as [`CarpenterKennedy2N54`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) and\n",
    "adaptive timestep methods such as [`SSPRK43`](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/) are supported.\n",
    "We would like to remark that of course every `stage_limiter!` can also be used as a `step_limiter!`, i.e.,\n",
    "acting only after the full time step has been taken."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As an example, we consider a variant of the [1D medium blast wave example](https://github.com/trixi-framework/Trixi.jl/blob/main/examples/tree_1d_dgsem/elixir_euler_blast_wave.jl)\n",
    "wherein the shock capturing method discussed above is employed to handle the shock.\n",
    "Here, we use a standard DG solver with HLLC surface flux:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi\n",
    "\n",
    "solver = DGSEM(polydeg = 3, surface_flux = flux_hllc)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The remaining setup is the same as in the standard example:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEqSSPRK\n",
    "\n",
    "###############################################################################"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "semidiscretization of the compressible Euler equations"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "equations = CompressibleEulerEquations1D(1.4)\n",
    "\n",
    "function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations1D)\n",
    "    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> \"medium blast wave\"\n",
    "    # Set up polar coordinates\n",
    "    inicenter = SVector(0.0)\n",
    "    x_norm = x[1] - inicenter[1]\n",
    "    r = abs(x_norm)\n",
    "    # The following code is equivalent to\n",
    "    # phi = atan(0.0, x_norm)\n",
    "    # cos_phi = cos(phi)\n",
    "    # in 1D but faster\n",
    "    cos_phi = x_norm > 0 ? one(x_norm) : -one(x_norm)\n",
    "\n",
    "    # Calculate primitive variables\n",
    "    rho = r > 0.5 ? 1.0 : 1.1691\n",
    "    v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi\n",
    "    p = r > 0.5 ? 1.0E-3 : 1.245\n",
    "\n",
    "    return prim2cons(SVector(rho, v1, p), equations)\n",
    "end\n",
    "initial_condition = initial_condition_blast_wave\n",
    "\n",
    "coordinates_min = (-2.0,)\n",
    "coordinates_max = (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",
    "###############################################################################"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "ODE solvers, callbacks etc."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "tspan = (0.0, 12.5)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "summary_callback = SummaryCallback()\n",
    "\n",
    "analysis_interval = 100\n",
    "\n",
    "analysis_callback = AnalysisCallback(semi, interval = analysis_interval)\n",
    "\n",
    "alive_callback = AliveCallback(analysis_interval = analysis_interval)\n",
    "\n",
    "stepsize_callback = StepsizeCallback(cfl = 0.5)\n",
    "\n",
    "callbacks = CallbackSet(summary_callback,\n",
    "                        analysis_callback, alive_callback,\n",
    "                        stepsize_callback)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We specify the `stage_limiter!` supplied to the classic SSPRK33 integrator\n",
    "with strict entropy increase enforcement\n",
    "(recall that the tolerated exponentiated entropy decrease is set to -1e-13)."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "stage_limiter! = EntropyBoundedLimiter()"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "We run the simulation with the SSPRK33 method and the entropy bounded limiter:"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = solve(ode, SSPRK33(stage_limiter!);\n",
    "            dt = 1.0,\n",
    "            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\", \"OrdinaryDiffEqLowStorageRK\", \"OrdinaryDiffEqSSPRK\", \"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.10"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.10",
   "language": "julia"
  }
 },
 "nbformat": 4
}