{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# 6: Subcell limiting with the IDP 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": [
    "In the previous tutorial, the element-wise limiting with `IndicatorHennemannGassner`\n",
    "and `VolumeIntegralShockCapturingHG` was explained. This tutorial contains a short\n",
    "introduction to the idea and implementation of subcell shock capturing approaches in Trixi.jl,\n",
    "which is also based on the DGSEM scheme in flux differencing formulation.\n",
    "Trixi.jl contains the a-posteriori invariant domain-preserving (IDP) limiter which was\n",
    "introduced by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) and\n",
    "[Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.1016/j.compfluid.2022.105627).\n",
    "It is a flux-corrected transport-type (FCT) limiter and is implemented using `SubcellLimiterIDP`\n",
    "and `VolumeIntegralSubcellLimiting`.\n",
    "Since it is an a-posteriori limiter you have to apply a correction stage after each Runge-Kutta\n",
    "stage. This is done by passing the stage callback `SubcellLimiterIDPCorrection` to the\n",
    "time integration method."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Time integration method\n",
    "As mentioned before, the IDP limiting is an a-posteriori limiter. Its limiting process\n",
    "guarantees the target bounds for an explicit (forward) Euler time step. To still achieve a\n",
    "high-order approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta\n",
    "methods, which can be written as convex combinations of forward Euler steps.\n",
    "As such, they preserve the convexity of convex functions and functionals, such as the TVD\n",
    "semi-norm and the maximum principle in 1D, for instance."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since IDP/FCT limiting procedure operates on independent forward Euler steps, its\n",
    "a-posteriori correction stage is implemented as a stage callback that is triggered after each\n",
    "forward Euler step in an SSP Runge-Kutta method. Unfortunately, the `solve(...)` routines in\n",
    "[OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), typically employed for time\n",
    "integration in Trixi.jl, do not support this type of stage callback."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern\n",
    "time integration SSPRK method called with\n",
    "````julia\n",
    "Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)\n",
    "````"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Right now, only the canonical three-stage, third-order SSPRK method (Shu-Osher)\n",
    "`Trixi.SimpleSSPRK33` is implemented."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "# IDP Limiting\n",
    "The implementation of the invariant domain preserving (IDP) limiting approach (`SubcellLimiterIDP`)\n",
    "is based on [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) and\n",
    "[Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.101/j.compfluid.2022.105627).\n",
    "It supports several types of limiting which are enabled by passing parameters individually."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Global bounds\n",
    "If enabled, the global bounds enforce physical admissibility conditions, such as non-negativity\n",
    "of variables. This can be done for conservative variables, where the limiter is of a one-sided\n",
    "Zalesak-type ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2)), and general\n",
    "non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The Newton-bisection algorithm is an iterative method and requires some parameters.\n",
    "It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and\n",
    "relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are\n",
    "sufficient in most cases and therefore used as default. Additionally, there is the parameter\n",
    "`gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation\n",
    "of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`,\n",
    "as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.]\n",
    "that this value guarantees the fulfillment of bounds for a forward-Euler increment."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Very small non-negative values can be an issue as well. That's why we use an additional\n",
    "correction factor in the calculation of the global bounds,\n",
    "$$\n",
    "u^{new} \\geq \\beta * u^{FV}.\n",
    "$$\n",
    "By default, $\\beta$ (named `positivity_correction_factor`) is set to `0.1` which works properly\n",
    "in most of the tested setups."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Conservative variables\n",
    "The procedure to enforce global bounds for a conservative variables is as follows:\n",
    "If you want to guarantee non-negativity for the density of the compressible Euler equations,\n",
    "you pass the specific quantity name of the conservative variable."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Trixi\n",
    "equations = CompressibleEulerEquations2D(1.4)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The quantity name of the density is `rho` which is how we enable its limiting."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "positivity_variables_cons = [\"rho\"]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The quantity names are passed as a vector to allow several quantities.\n",
    "This is used, for instance, if you want to limit the density of two different components using\n",
    "the multicomponent compressible Euler equations."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),\n",
    "                                                       gas_constants = (0.287, 1.578))"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, we just pass both quantity names."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "positivity_variables_cons = [\"rho1\", \"rho2\"]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Alternatively, it is possible to all limit all density variables with a general command using"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "positivity_variables_cons = [\"rho\" * string(i) for i in eachcomponent(equations)]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "#### Non-linear variables\n",
    "To allow limitation for all possible non-linear variables, including variables defined\n",
    "on-the-fly, you can directly pass the function that computes the quantity for which you want\n",
    "to enforce positivity. For instance, if you want to enforce non-negativity for the pressure,\n",
    "do as follows."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "positivity_variables_nonlinear = [pressure]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "### Local bounds\n",
    "Second, Trixi.jl supports the limiting with local bounds for conservative variables using a\n",
    "two-sided Zalesak-type limiter ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2))\n",
    "and for general non-linear variables using a one-sided Newton-bisection algorithm.\n",
    "They allow to avoid spurious oscillations within the global bounds and to improve the\n",
    "shock-capturing capabilities of the method. The corresponding numerical admissibility conditions\n",
    "are frequently formulated as local maximum or minimum principles. The local bounds are computed\n",
    "using the maximum and minimum values of all local neighboring nodes. Within this calculation we\n",
    "use the low-order FV solution values for each node."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "As for the limiting with global bounds you are passing the quantity names of the conservative\n",
    "variables you want to limit. So, to limit the density with lower and upper local bounds pass\n",
    "the following."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "local_twosided_variables_cons = [\"rho\"]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To limit non-linear variables locally, pass the variable function combined with the requested\n",
    "bound (`min` or `max`) as a tuple. For instance, to impose a lower local bound on the modified\n",
    "specific entropy `Trixi.entropy_guermond_etal`, use"
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal, min)]"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Exemplary simulation\n",
    "How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary\n",
    "setup. This will be a simplified version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`.\n",
    "Since the setup is mostly very similar to a pure DGSEM setup as in\n",
    "`tree_2d_dgsem/elixir_euler_blast_wave.jl`, the equivalent parts are used without any explanation\n",
    "here."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using OrdinaryDiffEq\n",
    "using Trixi\n",
    "\n",
    "equations = CompressibleEulerEquations2D(1.4)\n",
    "\n",
    "function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D)\n",
    "    # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> \"medium blast wave\"\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",
    "    phi = atan(y_norm, x_norm)\n",
    "    sin_phi, cos_phi = sincos(phi)\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",
    "    v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi\n",
    "    p = r > 0.5 ? 1.0E-3 : 1.245\n",
    "\n",
    "    return prim2cons(SVector(rho, v1, v2, p), equations)\n",
    "end\n",
    "initial_condition = initial_condition_blast_wave;"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Since the surface integral is equal for both the DG and the subcell FV method, the limiting is\n",
    "applied only in the volume integral."
   ],
   "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`\n",
    "or `flux_kennedy_gruber`, for the DG volume flux."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "surface_flux = flux_lax_friedrichs\n",
    "volume_flux = flux_ranocha"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The limiter is implemented within `SubcellLimiterIDP`. It always requires the\n",
    "parameters `equations` and `basis`. With additional parameters (described above\n",
    "or listed in the docstring) you can specify and enable additional limiting options.\n",
    "Here, the simulation should contain local limiting for the density using lower and upper bounds."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "basis = LobattoLegendreBasis(3)\n",
    "limiter_idp = SubcellLimiterIDP(equations, basis;\n",
    "                                local_twosided_variables_cons = [\"rho\"])"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "The initialized limiter is passed to `VolumeIntegralSubcellLimiting` in addition to the volume\n",
    "fluxes of the low-order and high-order scheme."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;\n",
    "                                                volume_flux_dg = volume_flux,\n",
    "                                                volume_flux_fv = surface_flux)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Then, the volume integral is passed to `solver` as it is done for the standard flux-differencing\n",
    "DG scheme or the element-wise limiting."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "solver = DGSEM(basis, surface_flux, volume_integral)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "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",
    "\n",
    "semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)\n",
    "\n",
    "tspan = (0.0, 2.0)\n",
    "ode = semidiscretize(semi, tspan)\n",
    "\n",
    "summary_callback = SummaryCallback()\n",
    "\n",
    "analysis_interval = 1000\n",
    "analysis_callback = AnalysisCallback(semi, interval = analysis_interval)\n",
    "\n",
    "alive_callback = AliveCallback(analysis_interval = analysis_interval)\n",
    "\n",
    "save_solution = SaveSolutionCallback(interval = 1000,\n",
    "                                     save_initial_solution = true,\n",
    "                                     save_final_solution = true,\n",
    "                                     solution_variables = cons2prim)\n",
    "\n",
    "stepsize_callback = StepsizeCallback(cfl = 0.3)\n",
    "\n",
    "callbacks = CallbackSet(summary_callback,\n",
    "                        analysis_callback, alive_callback,\n",
    "                        save_solution,\n",
    "                        stepsize_callback);"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "As explained above, the IDP limiter works a-posteriori and requires the additional use of a\n",
    "correction stage implemented with the stage callback `SubcellLimiterIDPCorrection`.\n",
    "This callback is passed within a tuple to the time integration method."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "stage_callbacks = (SubcellLimiterIDPCorrection(),)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "Moreover, as mentioned before as well, simulations with subcell limiting require a Trixi-intern\n",
    "SSPRK time integration methods with passed stage callbacks and a Trixi-intern `Trixi.solve(...)`\n",
    "routine."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks);\n",
    "                  dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback\n",
    "                  callback = callbacks);\n",
    "summary_callback() # print the timer summary"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Visualization\n",
    "As for a standard simulation in Trixi.jl, it is possible to visualize the solution using the\n",
    "`plot` routine from Plots.jl."
   ],
   "metadata": {}
  },
  {
   "outputs": [],
   "cell_type": "code",
   "source": [
    "using Plots\n",
    "plot(sol)"
   ],
   "metadata": {},
   "execution_count": null
  },
  {
   "cell_type": "markdown",
   "source": [
    "To get an additional look at the amount of limiting that is used, you can use the visualization\n",
    "approach using the `SaveSolutionCallback`, [`Trixi2Vtk`](https://github.com/trixi-framework/Trixi2Vtk.jl)\n",
    "and [ParaView](https://www.paraview.org/download/). More details about this procedure\n",
    "can be found in the visualization documentation.\n",
    "Unfortunately, the support for subcell limiting data is not yet merged into the main branch\n",
    "of Trixi2Vtk but lies in the branch [`bennibolm/node-variables`](https://github.com/bennibolm/Trixi2Vtk.jl/tree/node-variables)."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "With that implementation and the standard procedure used for Trixi2Vtk you get the following\n",
    "dropdown menu in ParaView."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "![ParaView_Dropdownmenu](https://github.com/trixi-framework/Trixi.jl/assets/74359358/70d15f6a-059b-4349-8291-68d9ab3af43e)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "The resulting visualization of the density and the limiting parameter then looks like this.\n",
    "![blast_wave_paraview](https://github.com/trixi-framework/Trixi.jl/assets/74359358/e5808bed-c8ab-43bf-af7a-050fe43dd630)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "You can see that the limiting coefficient does not lie in the interval [0,1] because Trixi2Vtk\n",
    "interpolates all quantities to regular nodes by default.\n",
    "You can disable this functionality with `reinterpolate=false` within the call of `trixi2vtk(...)`\n",
    "and get the following visualization.\n",
    "![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116)"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "## Bounds checking\n",
    "Subcell limiting is based on the fulfillment of target bounds - either global or local.\n",
    "Although the implementation works and has been thoroughly tested, there are some cases where\n",
    "these bounds are not met.\n",
    "For instance, the deviations could be in machine precision, which is not problematic.\n",
    "Larger deviations can be cause by too large time-step sizes (which can be easily fixed by\n",
    "reducing the CFL number), specific boundary conditions or source terms. Insufficient parameters\n",
    "for the Newton-bisection algorithm can also be a reason when limiting non-linear variables.\n",
    "There are described above."
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "In many cases, it is reasonable to monitor the bounds deviations.\n",
    "Because of that, Trixi.jl supports a bounds checking routine implemented using the stage\n",
    "callback `BoundsCheckCallback`. It checks all target bounds for fulfillment\n",
    "in every RK stage. If added to the tuple of stage callbacks like\n",
    "````julia\n",
    "stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())\n",
    "````\n",
    "and passed to the time integration method, a summary is added to the final console output.\n",
    "For the given example, this summary shows that all bounds are met at all times.\n",
    "````\n",
    "────────────────────────────────────────────────────────────────────────────────────────────────────\n",
    "Maximum deviation from bounds:\n",
    "────────────────────────────────────────────────────────────────────────────────────────────────────\n",
    "rho:\n",
    "- lower bound: 0.0\n",
    "- upper bound: 0.0\n",
    "────────────────────────────────────────────────────────────────────────────────────────────────────\n",
    "````"
   ],
   "metadata": {}
  },
  {
   "cell_type": "markdown",
   "source": [
    "Moreover, it is also possible to monitor the bounds deviations incurred during the simulations.\n",
    "To do that use the parameter `save_errors = true`, such that the instant deviations are written\n",
    "to `deviations.txt` in `output_directory` every `interval` time steps.\n",
    "````julia\n",
    "BoundsCheckCallback(save_errors = true, output_directory = \"out\", interval = 100)\n",
    "````\n",
    "Then, for the given example the deviations file contains all daviations for the current\n",
    "timestep and simulation time.\n",
    "````\n",
    "iter, simu_time, rho_min, rho_max\n",
    "100, 0.29103427131404924, 0.0, 0.0\n",
    "200, 0.5980281923063808, 0.0, 0.0\n",
    "300, 0.9520853560765293, 0.0, 0.0\n",
    "400, 1.3630295622683186, 0.0, 0.0\n",
    "500, 1.8344999624013498, 0.0, 0.0\n",
    "532, 1.9974179806990118, 0.0, 0.0\n",
    "````"
   ],
   "metadata": {}
  }
 ],
 "nbformat_minor": 3,
 "metadata": {
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.5"
  },
  "kernelspec": {
   "name": "julia-1.10",
   "display_name": "Julia 1.10.5",
   "language": "julia"
  }
 },
 "nbformat": 4
}