From 506315ee47c6220dd48399e5135326413a6424f1 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 20 Mar 2024 16:55:40 +0100 Subject: [PATCH 01/17] First version of tutorial for subcell IDP limiting --- docs/literate/src/files/index.jl | 34 ++++---- .../src/files/subcell_shock_capturing.jl | 78 +++++++++++++++++++ docs/make.jl | 1 + 3 files changed, 98 insertions(+), 15 deletions(-) create mode 100644 docs/literate/src/files/subcell_shock_capturing.jl diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index 6d5800c3b66..9bd75b8aace 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -51,20 +51,24 @@ # explained and added to an exemplary simulation of the Sedov blast wave with the 2D compressible Euler # equations. -# ### [6 Non-periodic boundary conditions](@ref non_periodic_boundaries) +# ### [6 Subcell limiting with the IDP Limiter](@ref subcell_shock_capturing) +#- +# TODO + +# ### [7 Non-periodic boundary conditions](@ref non_periodic_boundaries) #- # Thus far, all examples used periodic boundaries. In Trixi.jl, you can also set up a simulation with # non-periodic boundaries. This tutorial presents the implementation of the classical Dirichlet # boundary condition with a following example. Then, other non-periodic boundaries are mentioned. -# ### [7 DG schemes via `DGMulti` solver](@ref DGMulti_1) +# ### [8 DG schemes via `DGMulti` solver](@ref DGMulti_1) #- # This tutorial is about the more general DG solver [`DGMulti`](@ref), introduced [here](@ref DGMulti). # We are showing some examples for this solver, for instance with discretization nodes by Gauss or # triangular elements. Moreover, we present a simple way to include pre-defined triangulate meshes for # non-Cartesian domains using the package [StartUpDG.jl](https://github.com/jlchan/StartUpDG.jl). -# ### [8 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2) +# ### [9 Other SBP schemes (FD, CGSEM) via `DGMulti` solver](@ref DGMulti_2) #- # Supplementary to the previous tutorial about DG schemes via the `DGMulti` solver we now present # the possibility for `DGMulti` to use other SBP schemes via the package @@ -72,7 +76,7 @@ # For instance, we show how to set up a finite differences (FD) scheme and a continuous Galerkin # (CGSEM) method. -# ### [9 Upwind FD SBP schemes](@ref upwind_fdsbp) +# ### [10 Upwind FD SBP schemes](@ref upwind_fdsbp) #- # General SBP schemes can not only be used via the [`DGMulti`](@ref) solver but # also with a general `DG` solver. In particular, upwind finite difference SBP @@ -80,42 +84,42 @@ # schemes in the `DGMulti` framework, the interface is based on the package # [SummationByPartsOperators.jl](https://github.com/ranocha/SummationByPartsOperators.jl). -# ### [10 Adding a new scalar conservation law](@ref adding_new_scalar_equations) +# ### [11 Adding a new scalar conservation law](@ref adding_new_scalar_equations) #- # This tutorial explains how to add a new physics model using the example of the cubic conservation # law. First, we define the equation using a `struct` `CubicEquation` and the physical flux. Then, # the corresponding standard setup in Trixi.jl (`mesh`, `solver`, `semi` and `ode`) is implemented # and the ODE problem is solved by OrdinaryDiffEq's `solve` method. -# ### [11 Adding a non-conservative equation](@ref adding_nonconservative_equation) +# ### [12 Adding a non-conservative equation](@ref adding_nonconservative_equation) #- # In this part, another physics model is implemented, the nonconservative linear advection equation. # We run two different simulations with different levels of refinement and compare the resulting errors. -# ### [12 Parabolic terms](@ref parabolic_terms) +# ### [13 Parabolic terms](@ref parabolic_terms) #- # This tutorial describes how parabolic terms are implemented in Trixi.jl, e.g., # to solve the advection-diffusion equation. -# ### [13 Adding new parabolic terms](@ref adding_new_parabolic_terms) +# ### [14 Adding new parabolic terms](@ref adding_new_parabolic_terms) #- # This tutorial describes how new parabolic terms can be implemented using Trixi.jl. -# ### [14 Adaptive mesh refinement](@ref adaptive_mesh_refinement) +# ### [15 Adaptive mesh refinement](@ref adaptive_mesh_refinement) #- # Adaptive mesh refinement (AMR) helps to increase the accuracy in sensitive or turbolent regions while # not wasting resources for less interesting parts of the domain. This leads to much more efficient # simulations. This tutorial presents the implementation strategy of AMR in Trixi.jl, including the use of # different indicators and controllers. -# ### [15 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping) +# ### [16 Structured mesh with curvilinear mapping](@ref structured_mesh_mapping) #- # In this tutorial, the use of Trixi.jl's structured curved mesh type [`StructuredMesh`](@ref) is explained. # We present the two basic option to initialize such a mesh. First, the curved domain boundaries # of a circular cylinder are set by explicit boundary functions. Then, a fully curved mesh is # defined by passing the transformation mapping. -# ### [16 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial) +# ### [17 Unstructured meshes with HOHQMesh.jl](@ref hohqmesh_tutorial) #- # The purpose of this tutorial is to demonstrate how to use the [`UnstructuredMesh2D`](@ref) # functionality of Trixi.jl. This begins by running and visualizing an available unstructured @@ -124,26 +128,26 @@ # software in the Trixi.jl ecosystem, and then run a simulation using Trixi.jl on said mesh. # In the end, the tutorial briefly explains how to simulate an example using AMR via `P4estMesh`. -# ### [17 P4est mesh from gmsh](@ref p4est_from_gmsh) +# ### [18 P4est mesh from gmsh](@ref p4est_from_gmsh) #- # This tutorial describes how to obtain a [`P4estMesh`](@ref) from an existing mesh generated # by [`gmsh`](https://gmsh.info/) or any other meshing software that can export to the Abaqus # input `.inp` format. The tutorial demonstrates how edges/faces can be associated with boundary conditions based on the physical nodesets. -# ### [18 Explicit time stepping](@ref time_stepping) +# ### [19 Explicit time stepping](@ref time_stepping) #- # This tutorial is about time integration using [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl). # It explains how to use their algorithms and presents two types of time step choices - with error-based # and CFL-based adaptive step size control. -# ### [19 Differentiable programming](@ref differentiable_programming) +# ### [20 Differentiable programming](@ref differentiable_programming) #- # This part deals with some basic differentiable programming topics. For example, a Jacobian, its # eigenvalues and a curve of total energy (through the simulation) are calculated and plotted for # a few semidiscretizations. Moreover, we calculate an example for propagating errors with Measurement.jl # at the end. -# ### [20 Custom semidiscretization](@ref custom_semidiscretization) +# ### [21 Custom semidiscretization](@ref custom_semidiscretization) #- # This tutorial describes the [semidiscretiations](@ref overview-semidiscretizations) of Trixi.jl # and explains how to extend them for custom tasks. diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl new file mode 100644 index 00000000000..3e611821235 --- /dev/null +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -0,0 +1,78 @@ +#src # Subcell limiting with the IDP Limiter + +# In the previous tutorial, the element-wise limiting with [`IndicatorHennemannGassner`](@ref) +# and [`VolumeIntegralShockCapturingHG`](@ref) was explained. This tutorial contains a short +# introduction to the idea and implementation of subcell shock capturing approaches in Trixi.jl, +# which is also based on the DGSEM scheme in flux differencing formulation. +# Trixi.jl contains the a-posteriori invariant domain-preserving (IDP) limiter which was +# introduced by [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.1016/j.compfluid.2022.105627) +# and [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876). It is a flux-corrected +# transport-type (FCT) limiter and is implemented using [`SubcellLimiterIDP`](@ref) and +# [`VolumeIntegralSubcellLimiting`](@ref). +# Since it is an a-posteriori limiter you have to apply a correction stage after each Runge-Kutta +# stage. This is done by passing the stage callback [`SubcellLimiterIDPCorrection`](@ref) to the +# time integration method. + +# TODO: Some comment about the time integration method. (Don't have to be here) +# - Using Trixi-intern SSPRK methods with `Trixi.solve(ode, algorithm; ...)` + +# TODO: Some comments about +# - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) +# - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds) + +using Trixi + +# # `SubcellLimiterIDP` +# The IDP limiter supports several options of limiting which are passed very flexible as parameters to +# the limiter individually. + +# ## Global bounds +# First, there is the use of global bounds. If enabled, they enforce physical admissibility +# conditions, such as non-negativity of variables. +# This can be done for conservative variables, where the limiter is of a one-sided Zalesak-type, and +# general non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds. + +# ### Conservative variables +# The procedure to enforce global bounds for a conservative variables is as follows: +# If you want to guarantee non-negativity for the density of compressible Euler equations, +# you pass the specific quantity name of the conservative variable. +equations = CompressibleEulerEquations2D(1.4) + +# The quantity name of the density is `rho` shich is how we enable its limiting. +positivity_variables_cons = ["rho"] + +# The quantity names are passed as a vector to allow several quantities. +# This is for instance used if you want to limit the density of two different components using +# the multicomponent compressible Euler equations. +equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), + gas_constants = (0.287, 1.578)) + +# Then, we just pass both quantity names. +positivity_variables_cons = ["rho1", "rho2"] + +# Alternatively, it is possible to all limit all density variables with a general command using +# ````julia +# positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)] +# ```` + +# ### Non-linear variables +# To allow limitation for all possible non-linear variables including on-the-fly defined ones, +# you directly pass function here. +# For instance, if you want to enforce non-negativity for the pressure, do as follows. +positivity_variables_nonlinear = [pressure] + +# ## Local bounds (Shock capturing) +# Second, Trixi.jl supports the limiting with local bounds for conservative variables. They +# allow to avoid spurious oscillations within the global bounds and to improve the +# shock-capturing capabilities of the method. The corresponding numerical admissibility +# conditions are frequently formulated as local maximum or minimum principles. +# There are different option to choose local bounds. We calculate them using the low-order FV +# solution. + +# As for the limiting with global bounds you are passing the quantity names of the conservative +# variables you want to limit. So, to limit the density with lower and upper local bounds pass +# the following. +local_minmax_variables_cons = ["rho"] + +# ## Exemplary simulation +# TODO diff --git a/docs/make.jl b/docs/make.jl index 8427c4049bf..33779d591ad 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -60,6 +60,7 @@ files = [ "Introduction to DG methods" => "scalar_linear_advection_1d.jl", "DGSEM with flux differencing" => "DGSEM_FluxDiff.jl", "Shock capturing with flux differencing and stage limiter" => "shock_capturing.jl", + "Subcell limiting with the IDP Limiter" => "subcell_shock_capturing.jl", "Non-periodic boundaries" => "non_periodic_boundaries.jl", "DG schemes via `DGMulti` solver" => "DGMulti_1.jl", "Other SBP schemes (FD, CGSEM) via `DGMulti` solver" => "DGMulti_2.jl", From b1b3f82ea173b137a5d00635caa83ef86272d9ca Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 20 Mar 2024 19:04:46 +0100 Subject: [PATCH 02/17] Add paragraph about time integration method --- .../src/files/subcell_shock_capturing.jl | 22 +++++++++++++++++-- 1 file changed, 20 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 3e611821235..1d33169bef7 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -13,8 +13,26 @@ # stage. This is done by passing the stage callback [`SubcellLimiterIDPCorrection`](@ref) to the # time integration method. -# TODO: Some comment about the time integration method. (Don't have to be here) -# - Using Trixi-intern SSPRK methods with `Trixi.solve(ode, algorithm; ...)` +# ## Time integration method +# As mentioned before, the IDP limiting is an a-posteriori limiter. Its limiting process +# guarantees the target bounds for a simple Euler evolution. To still achieve a high-order +# approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta method +# which can be written as convex combination of these forward Euler steps. +#- +# Due to this functionality of the limiting procedure the correcting stage and therefore the stage +# callbacks has to be applied to the solution after the forward Euler step and before further +# computation. Unfortunately, the `solve(...)` routines of +# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), which is normally used in +# Trixi.jl for the time integration, does not support calculations via callback at this point +# in the simulation. +#- +# Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern +# time integration SSPRK method called with +# ````julia +# Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)`. +# ```` +#- +# Right now, only the third-order SSPRK method [`SimpleSSPRK33`](@ref) is supported. # TODO: Some comments about # - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) From 65a499a5f9fa6d33ee2ef0d44bc989bcc15dff71 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Wed, 20 Mar 2024 19:07:13 +0100 Subject: [PATCH 03/17] Redo code as julia comments --- .../src/files/subcell_shock_capturing.jl | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 1d33169bef7..cb571e14045 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -57,7 +57,9 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # The quantity name of the density is `rho` shich is how we enable its limiting. -positivity_variables_cons = ["rho"] +# ````julia +# positivity_variables_cons = ["rho"] +# ```` # The quantity names are passed as a vector to allow several quantities. # This is for instance used if you want to limit the density of two different components using @@ -66,7 +68,9 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), gas_constants = (0.287, 1.578)) # Then, we just pass both quantity names. -positivity_variables_cons = ["rho1", "rho2"] +# ````julia +# positivity_variables_cons = ["rho1", "rho2"] +# ```` # Alternatively, it is possible to all limit all density variables with a general command using # ````julia @@ -77,7 +81,9 @@ positivity_variables_cons = ["rho1", "rho2"] # To allow limitation for all possible non-linear variables including on-the-fly defined ones, # you directly pass function here. # For instance, if you want to enforce non-negativity for the pressure, do as follows. -positivity_variables_nonlinear = [pressure] +# ````julia +# positivity_variables_nonlinear = [pressure] +# ```` # ## Local bounds (Shock capturing) # Second, Trixi.jl supports the limiting with local bounds for conservative variables. They @@ -90,7 +96,9 @@ positivity_variables_nonlinear = [pressure] # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass # the following. -local_minmax_variables_cons = ["rho"] +# ````julia +# local_minmax_variables_cons = ["rho"] +# ```` # ## Exemplary simulation # TODO From a0b527506abc7e080eca75cbd031b97cb29e6b2b Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 21 Mar 2024 08:22:49 +0100 Subject: [PATCH 04/17] Hopefully fix reference --- docs/literate/src/files/subcell_shock_capturing.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index cb571e14045..375dfa391a7 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -32,7 +32,7 @@ # Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)`. # ```` #- -# Right now, only the third-order SSPRK method [`SimpleSSPRK33`](@ref) is supported. +# Right now, only the third-order SSPRK method [`Trixi.SimpleSSPRK33`](@ref) is implemented. # TODO: Some comments about # - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) From d5259c5eb584d570d13ce77f9a34682b8ad2edae Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 21 Mar 2024 14:56:21 +0100 Subject: [PATCH 05/17] Add section about Visualizaton --- .../src/files/subcell_shock_capturing.jl | 133 +++++++++++++++++- 1 file changed, 130 insertions(+), 3 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 375dfa391a7..cb0b3782af4 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -32,14 +32,13 @@ # Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)`. # ```` #- -# Right now, only the third-order SSPRK method [`Trixi.SimpleSSPRK33`](@ref) is implemented. +# Right now, only the canonical three-stage, third-order SSPRK method (Shu-Osher) +# [`Trixi.SimpleSSPRK33`](@ref) is implemented. # TODO: Some comments about # - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) # - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds) -using Trixi - # # `SubcellLimiterIDP` # The IDP limiter supports several options of limiting which are passed very flexible as parameters to # the limiter individually. @@ -54,6 +53,7 @@ using Trixi # The procedure to enforce global bounds for a conservative variables is as follows: # If you want to guarantee non-negativity for the density of compressible Euler equations, # you pass the specific quantity name of the conservative variable. +using Trixi equations = CompressibleEulerEquations2D(1.4) # The quantity name of the density is `rho` shich is how we enable its limiting. @@ -101,4 +101,131 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # ```` # ## Exemplary simulation +# How to set up a simulation using the IDP limiting becomes clearer when lokking at a exemplary +# setup. This will be a simplyfied version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`. +# Since the setup is mostly very similar to a pure DGSEM setup as in +# `tree_2d_dgsem/elixir_euler_blast_wave.jl`, the equivalent parts are without any explanation +# here. +using OrdinaryDiffEq +using Trixi + +equations = CompressibleEulerEquations2D(1.4) + +""" + initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + +A medium blast wave taken from +- Sebastian Hennemann, Gregor J. Gassner (2020) + A provably entropy stable subcell shock capturing approach for high order split form DG + [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) +""" +function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) + # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + # Set up polar coordinates + inicenter = SVector(0.0, 0.0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + sin_phi, cos_phi = sincos(phi) + + # Calculate primitive variables + rho = r > 0.5 ? 1.0 : 1.1691 + v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi + v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi + p = r > 0.5 ? 1.0E-3 : 1.245 + + return prim2cons(SVector(rho, v1, v2, p), equations) +end +initial_condition = initial_condition_blast_wave + +############################################################################### +# TODO: Some explanation +surface_flux = flux_lax_friedrichs +volume_flux = flux_ranocha +basis = LobattoLegendreBasis(3) +limiter_idp = SubcellLimiterIDP(equations, basis; + local_minmax_variables_cons = ["rho"]) +volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; + volume_flux_dg = volume_flux, + volume_flux_fv = surface_flux) +solver = DGSEM(basis, surface_flux, volume_integral) + + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 5, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver) + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 500 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 200, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.3) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback) + +############################################################################### +# As explained above, the IDP limiter works a-posteriori and requires the additional use of a +# correction stage implemented with the stage callback [`SubcellLimiterIDPCorrection`](@ref). +# This callback is passed within a tuple to the time integration method. +#- +# Moreover, as mentioned before as well, simulations with subcell limiting require a Trixi-intern +# SSPRK time integration methods with passed stage callbacks and a Trixi-intern `Trixi.solve(...)` +# routine. +stage_callbacks = (SubcellLimiterIDPCorrection(),) + +sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + callback = callbacks); +summary_callback() # print the timer summary + + +# ## Visualizaton +# As for a standard simulation in Trixi.jl, it is possible to visualize the solution using Plots +# `plot` routine. +using Plots +plot(sol) + +# To get an additional look at the amount of limiting that is used, you can use the visualization +# approach using the [`SaveSolutionCallback`](@ref), [`Trixi2Vtk`](https://github.com/trixi-framework/Trixi2Vtk.jl) +# and [ParaView](https://www.paraview.org/download/). More details about this procedure +# can be found in the [visualization documentation](@ref visualization). +# Unfortunately, the support for subcell limiting data is not yet merge into the main branch +# of Trixi2Vtk but lies in the branch `bennibolm/node-variables`. +#- +# With that implementation and the standard procedure used for Trixi2Vtk you get the following +# dropdown menu in ParaView. +# ![ParaView_Dropdownmenu](https://github.com/trixi-framework/Trixi.jl/assets/74359358/70d15f6a-059b-4349-8291-68d9ab3af43e) + +# The resulting visualization of the density and the limiting parameter then looks like this. +# ![blast_wave_paraview](https://github.com/trixi-framework/Trixi.jl/assets/74359358/e5808bed-c8ab-43bf-af7a-050fe43dd630) + +# You can see that the limiting coefficient does not lie in the interval [0,1], what actually was +# expected due to its calculation. +# TODO: Did I write something about this calculation? +# This is due to the reconstruction functionality which is defaultly enabled in Trixi2Vtk. +# You can disabled it with `reinterpolate=false` within the call of `trixi2vtk(...)` and get the +# following visualization. +# ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116) + + +# ## Target bounds checking # TODO From f7132f56cc1dafcf8ce2a5b00f9fe5bceef8c0a6 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 21 Mar 2024 14:57:59 +0100 Subject: [PATCH 06/17] Fix typos --- docs/literate/src/files/subcell_shock_capturing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index cb0b3782af4..4bbb7e47bf1 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -102,7 +102,7 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # ## Exemplary simulation # How to set up a simulation using the IDP limiting becomes clearer when lokking at a exemplary -# setup. This will be a simplyfied version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`. +# setup. This will be a simplified version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`. # Since the setup is mostly very similar to a pure DGSEM setup as in # `tree_2d_dgsem/elixir_euler_blast_wave.jl`, the equivalent parts are without any explanation # here. @@ -198,7 +198,7 @@ sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); summary_callback() # print the timer summary -# ## Visualizaton +# ## Visualization # As for a standard simulation in Trixi.jl, it is possible to visualize the solution using Plots # `plot` routine. using Plots From 6bb1e26a54cbcd9e055c56aab1eeee2421d35e74 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Thu, 21 Mar 2024 17:59:56 +0100 Subject: [PATCH 07/17] Small changes --- .../src/files/subcell_shock_capturing.jl | 41 ++++++++----------- 1 file changed, 16 insertions(+), 25 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 4bbb7e47bf1..946dad6cfc5 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -29,7 +29,7 @@ # Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern # time integration SSPRK method called with # ````julia -# Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...)`. +# Trixi.solve(ode, method(stage_callbacks = stage_callbacks); ...) # ```` #- # Right now, only the canonical three-stage, third-order SSPRK method (Shu-Osher) @@ -43,13 +43,13 @@ # The IDP limiter supports several options of limiting which are passed very flexible as parameters to # the limiter individually. -# ## Global bounds +# ### Global bounds # First, there is the use of global bounds. If enabled, they enforce physical admissibility # conditions, such as non-negativity of variables. # This can be done for conservative variables, where the limiter is of a one-sided Zalesak-type, and # general non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds. -# ### Conservative variables +# #### Conservative variables # The procedure to enforce global bounds for a conservative variables is as follows: # If you want to guarantee non-negativity for the density of compressible Euler equations, # you pass the specific quantity name of the conservative variable. @@ -77,7 +77,7 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)] # ```` -# ### Non-linear variables +# #### Non-linear variables # To allow limitation for all possible non-linear variables including on-the-fly defined ones, # you directly pass function here. # For instance, if you want to enforce non-negativity for the pressure, do as follows. @@ -85,7 +85,7 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # positivity_variables_nonlinear = [pressure] # ```` -# ## Local bounds (Shock capturing) +# ### Local bounds (Shock capturing) # Second, Trixi.jl supports the limiting with local bounds for conservative variables. They # allow to avoid spurious oscillations within the global bounds and to improve the # shock-capturing capabilities of the method. The corresponding numerical admissibility @@ -111,17 +111,9 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) -""" - initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - -A medium blast wave taken from -- Sebastian Hennemann, Gregor J. Gassner (2020) - A provably entropy stable subcell shock capturing approach for high order split form DG - [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044) -""" function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquations2D) - # Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" - # Set up polar coordinates + ## Modified From Hennemann & Gassner JCP paper 2020 (Sec. 6.3) -> "medium blast wave" + ## Set up polar coordinates inicenter = SVector(0.0, 0.0) x_norm = x[1] - inicenter[1] y_norm = x[2] - inicenter[2] @@ -129,7 +121,7 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation phi = atan(y_norm, x_norm) sin_phi, cos_phi = sincos(phi) - # Calculate primitive variables + ## Calculate primitive variables rho = r > 0.5 ? 1.0 : 1.1691 v1 = r > 0.5 ? 0.0 : 0.1882 * cos_phi v2 = r > 0.5 ? 0.0 : 0.1882 * sin_phi @@ -137,9 +129,8 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation return prim2cons(SVector(rho, v1, v2, p), equations) end -initial_condition = initial_condition_blast_wave +initial_condition = initial_condition_blast_wave; -############################################################################### # TODO: Some explanation surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha @@ -149,9 +140,9 @@ limiter_idp = SubcellLimiterIDP(equations, basis; volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral) - +solver = DGSEM(basis, surface_flux, volume_integral); +#- coordinates_min = (-2.0, -2.0) coordinates_max = (2.0, 2.0) mesh = TreeMesh(coordinates_min, coordinates_max, @@ -165,12 +156,12 @@ ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() -analysis_interval = 500 +analysis_interval = 1000 analysis_callback = AnalysisCallback(semi, interval = analysis_interval) alive_callback = AliveCallback(analysis_interval = analysis_interval) -save_solution = SaveSolutionCallback(interval = 200, +save_solution = SaveSolutionCallback(interval = 1000, save_initial_solution = true, save_final_solution = true, solution_variables = cons2prim) @@ -180,9 +171,8 @@ stepsize_callback = StepsizeCallback(cfl = 0.3) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, save_solution, - stepsize_callback) + stepsize_callback); -############################################################################### # As explained above, the IDP limiter works a-posteriori and requires the additional use of a # correction stage implemented with the stage callback [`SubcellLimiterIDPCorrection`](@ref). # This callback is passed within a tuple to the time integration method. @@ -213,6 +203,7 @@ plot(sol) #- # With that implementation and the standard procedure used for Trixi2Vtk you get the following # dropdown menu in ParaView. +#- # ![ParaView_Dropdownmenu](https://github.com/trixi-framework/Trixi.jl/assets/74359358/70d15f6a-059b-4349-8291-68d9ab3af43e) # The resulting visualization of the density and the limiting parameter then looks like this. @@ -227,5 +218,5 @@ plot(sol) # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116) -# ## Target bounds checking +# ## Bounds checking # TODO From daf4c5adf5a26fbcdd85f7e35f66f99eaecb842e Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Mar 2024 13:18:23 +0100 Subject: [PATCH 08/17] Add explanation about limiter, vol integral and solver --- .../src/files/subcell_shock_capturing.jl | 24 +++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 946dad6cfc5..524b19308c5 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -39,7 +39,7 @@ # - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) # - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds) -# # `SubcellLimiterIDP` +# # [`SubcellLimiterIDP`](@id SubcellLimiterIDP) # The IDP limiter supports several options of limiting which are passed very flexible as parameters to # the limiter individually. @@ -85,7 +85,7 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # positivity_variables_nonlinear = [pressure] # ```` -# ### Local bounds (Shock capturing) +# ### Local bounds # Second, Trixi.jl supports the limiting with local bounds for conservative variables. They # allow to avoid spurious oscillations within the global bounds and to improve the # shock-capturing capabilities of the method. The corresponding numerical admissibility @@ -131,16 +131,32 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation end initial_condition = initial_condition_blast_wave; -# TODO: Some explanation +# Since the surface integral is equal for both the DG and the subcell FV method, only the volume +# integral divides between the two methods. +#- +# Note, that the DG method is based on the flux differencing formulation. Hence, you have to use a +# two-point flux, such as [`flux_ranocha`](@ref), [`flux_shima_etal`](@ref), [`flux_chandrashekar`](@ref) +# or [`flux_kennedy_gruber`](@ref), for the DG volume flux. surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) + +# The actuall limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the +# parameters `equations` and `basis`. With additional parameters (described [above](@ref SubcellLimiterIDP) +# or listed in the docstring) you can specify and enabled the wanted limiting options. +# Here, the simulation should contain local limiting for the density using lower and upper bounds. limiter_idp = SubcellLimiterIDP(equations, basis; local_minmax_variables_cons = ["rho"]) + +# The initialized limiter is passed to `VolumeIntegralSubcellLimiting` in addition to the volume +# volume fluxes of the low-order and high-order scheme. volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) -solver = DGSEM(basis, surface_flux, volume_integral); + +# Then, the volume integral is passed to `solver` as it is done for the standard flux-differencing +# DG scheme or the element-wise limiting. +solver = DGSEM(basis, surface_flux, volume_integral) #- coordinates_min = (-2.0, -2.0) From 2a216058b1f536191c088188026da427eee84260 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Mar 2024 14:39:54 +0100 Subject: [PATCH 09/17] Add bounds checking section --- .../src/files/subcell_shock_capturing.jl | 49 +++++++++++++++++-- 1 file changed, 46 insertions(+), 3 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 524b19308c5..c1418894f64 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -205,8 +205,8 @@ summary_callback() # print the timer summary # ## Visualization -# As for a standard simulation in Trixi.jl, it is possible to visualize the solution using Plots -# `plot` routine. +# As for a standard simulation in Trixi.jl, it is possible to visualize the solution using the +# `plot` routine from Plots.jl. using Plots plot(sol) @@ -235,4 +235,47 @@ plot(sol) # ## Bounds checking -# TODO +# Subcell limiting is based on the fulfillment of target bounds - either global or local. +# Although the implementation is correct there are some cases where these bounds are not met. +# On the one hand, the deviations could be in machine precision which is not problematic. +# On the other hand, the target bound can be exceeded by larger deviations. Mostly, this is due +# to too large time step sizes and can be easily fixed by reducing the CFL number. +# Rarely, the larger deviations come from the use of specifc boundary condition of source terms. +#- +# In all this cases the exceeded bounds are not too problematic. Though, it is reasonable to +# monitor them. Because of that Trixi.jl supports a bounds checking routine implemented using +# the stage callback [`BoundsCheckCallback`](@ref). It checks all target bounds for fulfillment +# in every RK stage. If added to the tuple of stage callbacks like +# ````julia +# stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) +# ```` +# and passed to the time integration method, a summary is added to the final console output. +# For the given example this summary shows that all bounds are met at all times. +# ```` +# ──────────────────────────────────────────────────────────────────────────────────────────────────── +# Maximum deviation from bounds: +# ──────────────────────────────────────────────────────────────────────────────────────────────────── +# rho: +# - lower bound: 0.0 +# - upper bound: 0.0 +# ──────────────────────────────────────────────────────────────────────────────────────────────────── +# ```` + +# Moreover, it is also possible to monitor the deviations due to the simulations. To do that use +# the parameter `save_errors = true` and the current deviations are written to `deviations.txt` +# in `output_directory` every `interval` time steps. +# ````julia +# BoundsCheckCallback(save_errors = true, output_directory = "out", interval = 100) +# ```` +# Then, for the given example the deviations file contains all daviations for the current +# timestep and simulation time. +# ```` +# iter, simu_time, rho_min, rho_max +# 1, 0.0, 0.0, 0.0 +# 101, 0.29394033217556337, 0.0, 0.0 +# 201, 0.6012597465597065, 0.0, 0.0 +# 301, 0.9559096690030839, 0.0, 0.0 +# 401, 1.3674274981949077, 0.0, 0.0 +# 501, 1.8395301696603052, 0.0, 0.0 +# 532, 1.9974179806990118, 0.0, 0.0 +# ```` From 67a2d7b6facfff35780fac6629c893c8512eb7bd Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 22 Mar 2024 15:48:31 +0100 Subject: [PATCH 10/17] fix typos --- docs/literate/src/files/subcell_shock_capturing.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index c1418894f64..512e278f66b 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -141,7 +141,7 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) -# The actuall limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the +# The actual limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the # parameters `equations` and `basis`. With additional parameters (described [above](@ref SubcellLimiterIDP) # or listed in the docstring) you can specify and enabled the wanted limiting options. # Here, the simulation should contain local limiting for the density using lower and upper bounds. @@ -240,7 +240,7 @@ plot(sol) # On the one hand, the deviations could be in machine precision which is not problematic. # On the other hand, the target bound can be exceeded by larger deviations. Mostly, this is due # to too large time step sizes and can be easily fixed by reducing the CFL number. -# Rarely, the larger deviations come from the use of specifc boundary condition of source terms. +# Rarely, the larger deviations come from the use of specific boundary condition of source terms. #- # In all this cases the exceeded bounds are not too problematic. Though, it is reasonable to # monitor them. Because of that Trixi.jl supports a bounds checking routine implemented using From 130b7bfdf6698f1a53217a03c48ae938b81766ed Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 5 Apr 2024 15:02:19 +0200 Subject: [PATCH 11/17] Add description in introduction --- docs/literate/src/files/index.jl | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index 9bd75b8aace..dd63d788a24 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -53,7 +53,13 @@ # ### [6 Subcell limiting with the IDP Limiter](@ref subcell_shock_capturing) #- -# TODO +# Trixi.jl also includes subcell-wise limiting using an invariant domain-preserving (IDP) limiting +# approach. It calculates the blending factor between a high-order DG method and a low-order subcell +# finite volume (FV) method for each node within each element independently. Therefore, the limiting +# acts more local and requires less limiting than the element-wise limiting. In addition to the support +# of local bounds (mostly used for shock-capturing), the implementation also supports the application +# of global bounds, which result in the minimal necessary amount of limiting required for physical +# admissibility conditions (e.g. non-negativity of variables). # ### [7 Non-periodic boundary conditions](@ref non_periodic_boundaries) #- From 4f43146baa76ef01c73c47124593b62332524e5d Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 8 Apr 2024 11:00:32 +0200 Subject: [PATCH 12/17] Clear out before hohq mesh tutorial --- docs/literate/src/files/hohqmesh_tutorial.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/literate/src/files/hohqmesh_tutorial.jl b/docs/literate/src/files/hohqmesh_tutorial.jl index b19d363c4bf..dd81f47951e 100644 --- a/docs/literate/src/files/hohqmesh_tutorial.jl +++ b/docs/literate/src/files/hohqmesh_tutorial.jl @@ -35,6 +35,7 @@ # There is a default example for this mesh type that can be executed by using Trixi +rm("out", force = true, recursive = true) #hide #md redirect_stdio(stdout=devnull, stderr=devnull) do # code that prints annoying stuff we don't want to see here #hide #md trixi_include(default_example_unstructured()) end #hide #md From c84fc08101a511e8f3f057875729101475583f89 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 8 Apr 2024 12:24:18 +0200 Subject: [PATCH 13/17] Adapt tutorial --- .../src/files/subcell_shock_capturing.jl | 30 +++++++------------ 1 file changed, 11 insertions(+), 19 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 512e278f66b..78158d3c51e 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -39,9 +39,11 @@ # - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) # - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds) -# # [`SubcellLimiterIDP`](@id SubcellLimiterIDP) -# The IDP limiter supports several options of limiting which are passed very flexible as parameters to -# the limiter individually. +# # [IDP Limiting](@id IDPLimiter) +# The implementation of the invariant domain preserving (IDP) limiting approach ([`SubcellLimiterIDP`](@ref)) +# is based on [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.101/j.compfluid.2022.105627) +# and [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876). +# It supports several types of limiting which are enabled by passing parameters individually. # ### Global bounds # First, there is the use of global bounds. If enabled, they enforce physical admissibility @@ -57,9 +59,7 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) # The quantity name of the density is `rho` shich is how we enable its limiting. -# ````julia -# positivity_variables_cons = ["rho"] -# ```` +positivity_variables_cons = ["rho"] # The quantity names are passed as a vector to allow several quantities. # This is for instance used if you want to limit the density of two different components using @@ -68,22 +68,16 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), gas_constants = (0.287, 1.578)) # Then, we just pass both quantity names. -# ````julia -# positivity_variables_cons = ["rho1", "rho2"] -# ```` +positivity_variables_cons = ["rho1", "rho2"] # Alternatively, it is possible to all limit all density variables with a general command using -# ````julia -# positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)] -# ```` +positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)] # #### Non-linear variables # To allow limitation for all possible non-linear variables including on-the-fly defined ones, # you directly pass function here. # For instance, if you want to enforce non-negativity for the pressure, do as follows. -# ````julia -# positivity_variables_nonlinear = [pressure] -# ```` +positivity_variables_nonlinear = [pressure] # ### Local bounds # Second, Trixi.jl supports the limiting with local bounds for conservative variables. They @@ -96,9 +90,7 @@ equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass # the following. -# ````julia -# local_minmax_variables_cons = ["rho"] -# ```` +local_minmax_variables_cons = ["rho"] # ## Exemplary simulation # How to set up a simulation using the IDP limiting becomes clearer when lokking at a exemplary @@ -142,7 +134,7 @@ volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) # The actual limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the -# parameters `equations` and `basis`. With additional parameters (described [above](@ref SubcellLimiterIDP) +# parameters `equations` and `basis`. With additional parameters (described [above](@ref IDPLimiter) # or listed in the docstring) you can specify and enabled the wanted limiting options. # Here, the simulation should contain local limiting for the density using lower and upper bounds. limiter_idp = SubcellLimiterIDP(equations, basis; From 7cb739b936b192d48419a987be809642318eec41 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Fri, 26 Apr 2024 13:53:42 +0200 Subject: [PATCH 14/17] Implement suggestions. --- docs/literate/src/files/index.jl | 15 +-- .../src/files/subcell_shock_capturing.jl | 99 +++++++++---------- 2 files changed, 57 insertions(+), 57 deletions(-) diff --git a/docs/literate/src/files/index.jl b/docs/literate/src/files/index.jl index e46e2c7d8fd..fd1d884d47d 100644 --- a/docs/literate/src/files/index.jl +++ b/docs/literate/src/files/index.jl @@ -59,13 +59,14 @@ #src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing # ### [{index} Subcell limiting with the IDP Limiter](@ref subcell_shock_capturing) #- -# Trixi.jl also includes subcell-wise limiting using an invariant domain-preserving (IDP) limiting -# approach. It calculates the blending factor between a high-order DG method and a low-order subcell -# finite volume (FV) method for each node within each element independently. Therefore, the limiting -# acts more local and requires less limiting than the element-wise limiting. In addition to the support -# of local bounds (mostly used for shock-capturing), the implementation also supports the application -# of global bounds, which result in the minimal necessary amount of limiting required for physical -# admissibility conditions (e.g. non-negativity of variables). +# Trixi.jl features a subcell-wise limiting strategy utilizing an Invariant Domain-Preserving (IDP) +# approach. This IDP approach computes a blending factor that balances the high-order +# discontinuous Galerkin (DG) method with a low-order subcell finite volume (FV) method for each +# node within an element. This localized approach minimizes the application of dissipation, +# resulting in less limiting compared to the element-wise strategy. Additionally, the framework +# supports both local bounds, which are primarily used for shock capturing, and global bounds. +# The application of global bounds ensures the minimal necessary limiting to meet physical +# admissibility conditions, such as ensuring the non-negativity of variables. #src Note to developers: Use "{ index }" (but without spaces, see next line) to enable automatic indexing # ### [{index} Non-periodic boundary conditions](@ref non_periodic_boundaries) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 78158d3c51e..1251c82b6e5 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -5,26 +5,25 @@ # introduction to the idea and implementation of subcell shock capturing approaches in Trixi.jl, # which is also based on the DGSEM scheme in flux differencing formulation. # Trixi.jl contains the a-posteriori invariant domain-preserving (IDP) limiter which was -# introduced by [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.1016/j.compfluid.2022.105627) -# and [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876). It is a flux-corrected -# transport-type (FCT) limiter and is implemented using [`SubcellLimiterIDP`](@ref) and -# [`VolumeIntegralSubcellLimiting`](@ref). +# introduced by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) and +# [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.1016/j.compfluid.2022.105627). +# It is a flux-corrected transport-type (FCT) limiter and is implemented using [`SubcellLimiterIDP`](@ref) +# and [`VolumeIntegralSubcellLimiting`](@ref). # Since it is an a-posteriori limiter you have to apply a correction stage after each Runge-Kutta # stage. This is done by passing the stage callback [`SubcellLimiterIDPCorrection`](@ref) to the # time integration method. # ## Time integration method # As mentioned before, the IDP limiting is an a-posteriori limiter. Its limiting process -# guarantees the target bounds for a simple Euler evolution. To still achieve a high-order -# approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta method -# which can be written as convex combination of these forward Euler steps. +# guarantees the target bounds for an explicit (forward) Euler time step. To still achieve a +# high-order approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta +# methods, which can be written as convex combinations of forward Euler steps. #- -# Due to this functionality of the limiting procedure the correcting stage and therefore the stage -# callbacks has to be applied to the solution after the forward Euler step and before further -# computation. Unfortunately, the `solve(...)` routines of -# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), which is normally used in -# Trixi.jl for the time integration, does not support calculations via callback at this point -# in the simulation. +# Since IDP/FCT limiting procedure operates on independent forward Euler steps, its +# a-posteriori correction stage is implemented as a stage callback that is triggered after each +# forward Euler step in an SSP Runge-Kutta method. Unfortunately, the `solve(...)` routines in +# [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl), typically employed for time +# integration in Trixi.jl, do not support this type of stage callback. #- # Therefore, subcell limiting with the IDP limiter requires the use of a Trixi-intern # time integration SSPRK method called with @@ -41,19 +40,19 @@ # # [IDP Limiting](@id IDPLimiter) # The implementation of the invariant domain preserving (IDP) limiting approach ([`SubcellLimiterIDP`](@ref)) -# is based on [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.101/j.compfluid.2022.105627) -# and [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876). +# is based on [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) and +# [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.101/j.compfluid.2022.105627). # It supports several types of limiting which are enabled by passing parameters individually. # ### Global bounds -# First, there is the use of global bounds. If enabled, they enforce physical admissibility -# conditions, such as non-negativity of variables. -# This can be done for conservative variables, where the limiter is of a one-sided Zalesak-type, and -# general non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds. +# If enabled, the global bounds enforce physical admissibility conditions, such as non-negativity +# of variables. This can be done for conservative variables, where the limiter is of a one-sided +# Zalesak-type, and general non-linear variables, where a Newton-bisection algorithm is used to +# enforce the bounds. # #### Conservative variables # The procedure to enforce global bounds for a conservative variables is as follows: -# If you want to guarantee non-negativity for the density of compressible Euler equations, +# If you want to guarantee non-negativity for the density of the compressible Euler equations, # you pass the specific quantity name of the conservative variable. using Trixi equations = CompressibleEulerEquations2D(1.4) @@ -62,7 +61,7 @@ equations = CompressibleEulerEquations2D(1.4) positivity_variables_cons = ["rho"] # The quantity names are passed as a vector to allow several quantities. -# This is for instance used if you want to limit the density of two different components using +# This is used, for instance, if you want to limit the density of two different components using # the multicomponent compressible Euler equations. equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648), gas_constants = (0.287, 1.578)) @@ -74,18 +73,19 @@ positivity_variables_cons = ["rho1", "rho2"] positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations)] # #### Non-linear variables -# To allow limitation for all possible non-linear variables including on-the-fly defined ones, -# you directly pass function here. -# For instance, if you want to enforce non-negativity for the pressure, do as follows. +# To allow limitation for all possible non-linear variables, including variables defined +# on-the-fly, you can directly pass the function that computes the quantity for which you want +# to enforce positivity. For instance, if you want to enforce non-negativity for the pressure, +# do as follows. positivity_variables_nonlinear = [pressure] # ### Local bounds # Second, Trixi.jl supports the limiting with local bounds for conservative variables. They -# allow to avoid spurious oscillations within the global bounds and to improve the +# allow to avoid spurious oscillations within the global bounds and to improve the # shock-capturing capabilities of the method. The corresponding numerical admissibility # conditions are frequently formulated as local maximum or minimum principles. -# There are different option to choose local bounds. We calculate them using the low-order FV -# solution. +# The local bounds are computed using the maximum and minimum values of all local neighboring +# nodes. Within this calculation we use the low-order FV solution values for each node. # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass @@ -123,8 +123,8 @@ function initial_condition_blast_wave(x, t, equations::CompressibleEulerEquation end initial_condition = initial_condition_blast_wave; -# Since the surface integral is equal for both the DG and the subcell FV method, only the volume -# integral divides between the two methods. +# Since the surface integral is equal for both the DG and the subcell FV method, the limiting is +# applied only in the volume integral. #- # Note, that the DG method is based on the flux differencing formulation. Hence, you have to use a # two-point flux, such as [`flux_ranocha`](@ref), [`flux_shima_etal`](@ref), [`flux_chandrashekar`](@ref) @@ -133,9 +133,9 @@ surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha basis = LobattoLegendreBasis(3) -# The actual limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the +# The limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the # parameters `equations` and `basis`. With additional parameters (described [above](@ref IDPLimiter) -# or listed in the docstring) you can specify and enabled the wanted limiting options. +# or listed in the docstring) you can specify and enable additional limiting options. # Here, the simulation should contain local limiting for the density using lower and upper bounds. limiter_idp = SubcellLimiterIDP(equations, basis; local_minmax_variables_cons = ["rho"]) @@ -206,8 +206,8 @@ plot(sol) # approach using the [`SaveSolutionCallback`](@ref), [`Trixi2Vtk`](https://github.com/trixi-framework/Trixi2Vtk.jl) # and [ParaView](https://www.paraview.org/download/). More details about this procedure # can be found in the [visualization documentation](@ref visualization). -# Unfortunately, the support for subcell limiting data is not yet merge into the main branch -# of Trixi2Vtk but lies in the branch `bennibolm/node-variables`. +# Unfortunately, the support for subcell limiting data is not yet merged into the main branch +# of Trixi2Vtk but lies in the branch [`bennibolm/node-variables`](https://github.com/bennibolm/Trixi2Vtk.jl/tree/node-variables). #- # With that implementation and the standard procedure used for Trixi2Vtk you get the following # dropdown menu in ParaView. @@ -217,32 +217,31 @@ plot(sol) # The resulting visualization of the density and the limiting parameter then looks like this. # ![blast_wave_paraview](https://github.com/trixi-framework/Trixi.jl/assets/74359358/e5808bed-c8ab-43bf-af7a-050fe43dd630) -# You can see that the limiting coefficient does not lie in the interval [0,1], what actually was -# expected due to its calculation. +# You can see that the limiting coefficient does not lie in the interval [0,1] because Trixi2Vtk +# interpolates all quantities to regular nodes by default. # TODO: Did I write something about this calculation? -# This is due to the reconstruction functionality which is defaultly enabled in Trixi2Vtk. -# You can disabled it with `reinterpolate=false` within the call of `trixi2vtk(...)` and get the -# following visualization. +# You can disable this functionality with `reinterpolate=false` within the call of `trixi2vtk(...)` +# and get the following visualization. # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116) # ## Bounds checking # Subcell limiting is based on the fulfillment of target bounds - either global or local. -# Although the implementation is correct there are some cases where these bounds are not met. -# On the one hand, the deviations could be in machine precision which is not problematic. -# On the other hand, the target bound can be exceeded by larger deviations. Mostly, this is due -# to too large time step sizes and can be easily fixed by reducing the CFL number. -# Rarely, the larger deviations come from the use of specific boundary condition of source terms. +# Although the implementation works and has been thoroughly tested, there are some cases where +# these bounds are not met. +# For instance, the deviations could be in machine precision, which is not problematic. +# Larger deviations can be cause by too large time-step sizes (which can be easily fixed by +# reducing the CFL number), specific boundary conditions or source terms. #- -# In all this cases the exceeded bounds are not too problematic. Though, it is reasonable to -# monitor them. Because of that Trixi.jl supports a bounds checking routine implemented using -# the stage callback [`BoundsCheckCallback`](@ref). It checks all target bounds for fulfillment +# In many cases, it is reasonable to monitor the bounds deviations. +# Because of that, Trixi.jl supports a bounds checking routine implemented using the stage +# callback [`BoundsCheckCallback`](@ref). It checks all target bounds for fulfillment # in every RK stage. If added to the tuple of stage callbacks like # ````julia # stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback()) # ```` # and passed to the time integration method, a summary is added to the final console output. -# For the given example this summary shows that all bounds are met at all times. +# For the given example, this summary shows that all bounds are met at all times. # ```` # ──────────────────────────────────────────────────────────────────────────────────────────────────── # Maximum deviation from bounds: @@ -253,9 +252,9 @@ plot(sol) # ──────────────────────────────────────────────────────────────────────────────────────────────────── # ```` -# Moreover, it is also possible to monitor the deviations due to the simulations. To do that use -# the parameter `save_errors = true` and the current deviations are written to `deviations.txt` -# in `output_directory` every `interval` time steps. +# Moreover, it is also possible to monitor the bounds deviations incurred during the simulations. +# To do that use the parameter `save_errors = true`, such that the instant deviations are written +# to `deviations.txt` in `output_directory` every `interval` time steps. # ````julia # BoundsCheckCallback(save_errors = true, output_directory = "out", interval = 100) # ```` From 387d7871e2476466440e173b7f1d3eb292dd76a8 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 29 Apr 2024 12:17:13 +0200 Subject: [PATCH 15/17] Add section about newton method and correction factor --- .../src/files/subcell_shock_capturing.jl | 26 ++++++++++++++----- 1 file changed, 19 insertions(+), 7 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 1251c82b6e5..1e4769c32e4 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -34,22 +34,33 @@ # Right now, only the canonical three-stage, third-order SSPRK method (Shu-Osher) # [`Trixi.SimpleSSPRK33`](@ref) is implemented. -# TODO: Some comments about -# - parameters of Newton method (max_iterations_newton = 10, newton_tolerances = (1.0e-12, 1.0e-14), gamma_constant_newton = 2 * ndims(equations))) -# - positivity_correction_factor (Maybe show calculation of bounds, also of local bounds) - # # [IDP Limiting](@id IDPLimiter) # The implementation of the invariant domain preserving (IDP) limiting approach ([`SubcellLimiterIDP`](@ref)) # is based on [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) and # [Rueda-Ramírez, Pazner, Gassner (2022)](https://doi.org/10.101/j.compfluid.2022.105627). # It supports several types of limiting which are enabled by passing parameters individually. -# ### Global bounds +# ### [Global bounds](@id global_bounds) # If enabled, the global bounds enforce physical admissibility conditions, such as non-negativity # of variables. This can be done for conservative variables, where the limiter is of a one-sided # Zalesak-type, and general non-linear variables, where a Newton-bisection algorithm is used to # enforce the bounds. +# The Newton-bisection algorithm is an iterative method and requires some parameters. +# It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and +# relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are +# sufficient in most cases and therefore used as default. Additionally, there is the parameter +# `gamma_constant_newton` which approximately scales the flux into one +# to all directions. Here, one can use basically always the default of `2 * ndims(equations)`. + +# Very small non-negative values can be an issue as well. That's why we use an additional +# correction factor in the calculation of the global bounds, +# ```math +# u^{new} \geq \beta * u^{FV}. +# ``` +# By default, $\beta$ (named `positivity_correction_factor`) is set to `0.1` which works properly +# in most of the tested setups. + # #### Conservative variables # The procedure to enforce global bounds for a conservative variables is as follows: # If you want to guarantee non-negativity for the density of the compressible Euler equations, @@ -219,7 +230,6 @@ plot(sol) # You can see that the limiting coefficient does not lie in the interval [0,1] because Trixi2Vtk # interpolates all quantities to regular nodes by default. -# TODO: Did I write something about this calculation? # You can disable this functionality with `reinterpolate=false` within the call of `trixi2vtk(...)` # and get the following visualization. # ![blast_wave_paraview_reinterpolate=false](https://github.com/trixi-framework/Trixi.jl/assets/74359358/39274f18-0064-469c-b4da-bac4b843e116) @@ -231,7 +241,9 @@ plot(sol) # these bounds are not met. # For instance, the deviations could be in machine precision, which is not problematic. # Larger deviations can be cause by too large time-step sizes (which can be easily fixed by -# reducing the CFL number), specific boundary conditions or source terms. +# reducing the CFL number), specific boundary conditions or source terms. Insufficient parameters +# for the Newton-bisection algorithm can also be a reason when limiting non-linear variables. +# There are described [above](@ref global_bounds). #- # In many cases, it is reasonable to monitor the bounds deviations. # Because of that, Trixi.jl supports a bounds checking routine implemented using the stage From 868edcc9647df33a1933ee79722739b922150a9a Mon Sep 17 00:00:00 2001 From: bennibolm Date: Mon, 6 May 2024 11:09:04 +0200 Subject: [PATCH 16/17] Implement suggestion; Fix typos --- .../src/files/subcell_shock_capturing.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 1e4769c32e4..1a49fc73cc4 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -50,8 +50,10 @@ # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and # relative/absolute tolerances (`newton_tolerances = (1.0e-12, 1.0e-14)`). The given values are # sufficient in most cases and therefore used as default. Additionally, there is the parameter -# `gamma_constant_newton` which approximately scales the flux into one -# to all directions. Here, one can use basically always the default of `2 * ndims(equations)`. +# `gamma_constant_newton`, which can be used to scale the antidiffusive flux for the computation +# of the blending coefficients of nonlinear variables. The default value is `2 * ndims(equations)`, +# as it was shown by [Pazner (2020)](https://doi.org/10.1016/j.cma.2021.113876) [Section 4.2.2.] +# that this value guarantees the fulfillment of bounds for a forward-Euler increment. # Very small non-negative values can be an issue as well. That's why we use an additional # correction factor in the calculation of the global bounds, @@ -68,7 +70,7 @@ using Trixi equations = CompressibleEulerEquations2D(1.4) -# The quantity name of the density is `rho` shich is how we enable its limiting. +# The quantity name of the density is `rho` which is how we enable its limiting. positivity_variables_cons = ["rho"] # The quantity names are passed as a vector to allow several quantities. @@ -142,17 +144,17 @@ initial_condition = initial_condition_blast_wave; # or [`flux_kennedy_gruber`](@ref), for the DG volume flux. surface_flux = flux_lax_friedrichs volume_flux = flux_ranocha -basis = LobattoLegendreBasis(3) # The limiter is implemented within [`SubcellLimiterIDP`](@ref). It always requires the # parameters `equations` and `basis`. With additional parameters (described [above](@ref IDPLimiter) # or listed in the docstring) you can specify and enable additional limiting options. # Here, the simulation should contain local limiting for the density using lower and upper bounds. +basis = LobattoLegendreBasis(3) limiter_idp = SubcellLimiterIDP(equations, basis; local_minmax_variables_cons = ["rho"]) # The initialized limiter is passed to `VolumeIntegralSubcellLimiting` in addition to the volume -# volume fluxes of the low-order and high-order scheme. +# fluxes of the low-order and high-order scheme. volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; volume_flux_dg = volume_flux, volume_flux_fv = surface_flux) @@ -195,12 +197,11 @@ callbacks = CallbackSet(summary_callback, # As explained above, the IDP limiter works a-posteriori and requires the additional use of a # correction stage implemented with the stage callback [`SubcellLimiterIDPCorrection`](@ref). # This callback is passed within a tuple to the time integration method. -#- +stage_callbacks = (SubcellLimiterIDPCorrection(),) + # Moreover, as mentioned before as well, simulations with subcell limiting require a Trixi-intern # SSPRK time integration methods with passed stage callbacks and a Trixi-intern `Trixi.solve(...)` # routine. -stage_callbacks = (SubcellLimiterIDPCorrection(),) - sol = Trixi.solve(ode, Trixi.SimpleSSPRK33(stage_callbacks = stage_callbacks); dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback callback = callbacks); From cf49c377c88cb64175335e32815c8f1fab0b4389 Mon Sep 17 00:00:00 2001 From: bennibolm Date: Tue, 7 May 2024 10:35:43 +0200 Subject: [PATCH 17/17] Implement suggestions --- .../src/files/subcell_shock_capturing.jl | 24 ++++++++++--------- 1 file changed, 13 insertions(+), 11 deletions(-) diff --git a/docs/literate/src/files/subcell_shock_capturing.jl b/docs/literate/src/files/subcell_shock_capturing.jl index 1a49fc73cc4..8b8e49a28a8 100644 --- a/docs/literate/src/files/subcell_shock_capturing.jl +++ b/docs/literate/src/files/subcell_shock_capturing.jl @@ -18,6 +18,8 @@ # guarantees the target bounds for an explicit (forward) Euler time step. To still achieve a # high-order approximation, the implementation uses strong-stability preserving (SSP) Runge-Kutta # methods, which can be written as convex combinations of forward Euler steps. +# As such, they preserve the convexity of convex functions and functionals, such as the TVD +# semi-norm and the maximum principle in 1D, for instance. #- # Since IDP/FCT limiting procedure operates on independent forward Euler steps, its # a-posteriori correction stage is implemented as a stage callback that is triggered after each @@ -43,8 +45,8 @@ # ### [Global bounds](@id global_bounds) # If enabled, the global bounds enforce physical admissibility conditions, such as non-negativity # of variables. This can be done for conservative variables, where the limiter is of a one-sided -# Zalesak-type, and general non-linear variables, where a Newton-bisection algorithm is used to -# enforce the bounds. +# Zalesak-type ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2)), and general +# non-linear variables, where a Newton-bisection algorithm is used to enforce the bounds. # The Newton-bisection algorithm is an iterative method and requires some parameters. # It uses a fixed maximum number of iteration steps (`max_iterations_newton = 10`) and @@ -93,12 +95,13 @@ positivity_variables_cons = ["rho" * string(i) for i in eachcomponent(equations) positivity_variables_nonlinear = [pressure] # ### Local bounds -# Second, Trixi.jl supports the limiting with local bounds for conservative variables. They -# allow to avoid spurious oscillations within the global bounds and to improve the -# shock-capturing capabilities of the method. The corresponding numerical admissibility -# conditions are frequently formulated as local maximum or minimum principles. -# The local bounds are computed using the maximum and minimum values of all local neighboring -# nodes. Within this calculation we use the low-order FV solution values for each node. +# Second, Trixi.jl supports the limiting with local bounds for conservative variables using a +# two-sided Zalesak-type limiter ([Zalesak, 1979](https://doi.org/10.1016/0021-9991(79)90051-2)). +# They allow to avoid spurious oscillations within the global bounds and to improve the +# shock-capturing capabilities of the method. The corresponding numerical admissibility conditions +# are frequently formulated as local maximum or minimum principles. The local bounds are computed +# using the maximum and minimum values of all local neighboring nodes. Within this calculation we +# use the low-order FV solution values for each node. # As for the limiting with global bounds you are passing the quantity names of the conservative # variables you want to limit. So, to limit the density with lower and upper local bounds pass @@ -106,10 +109,10 @@ positivity_variables_nonlinear = [pressure] local_minmax_variables_cons = ["rho"] # ## Exemplary simulation -# How to set up a simulation using the IDP limiting becomes clearer when lokking at a exemplary +# How to set up a simulation using the IDP limiting becomes clearer when looking at an exemplary # setup. This will be a simplified version of `tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl`. # Since the setup is mostly very similar to a pure DGSEM setup as in -# `tree_2d_dgsem/elixir_euler_blast_wave.jl`, the equivalent parts are without any explanation +# `tree_2d_dgsem/elixir_euler_blast_wave.jl`, the equivalent parts are used without any explanation # here. using OrdinaryDiffEq using Trixi @@ -162,7 +165,6 @@ volume_integral = VolumeIntegralSubcellLimiting(limiter_idp; # Then, the volume integral is passed to `solver` as it is done for the standard flux-differencing # DG scheme or the element-wise limiting. solver = DGSEM(basis, surface_flux, volume_integral) - #- coordinates_min = (-2.0, -2.0) coordinates_max = (2.0, 2.0)