Add two-sided Zalesak-type IDP subcell limiting (#1648)

* Add two-sided limiting for conservative variables * Fix visualization routines * Add bounds calculation for BoundaryConditionDirichlet * Reduce cfl in elixir * Fix test * Add comment about subcell limiting with non-conforming meshes * Remove subcell visualization * Fix last commit * Remove @unpack * Add comment to `News.md` * Fix source for sedov blast setup; Formatting * Reduce allocations * Replace construction of Symbols * Add bounds check for local limiting * Implement suggestions * Fix format * Add subcell allocation tests; Add changes to minmax limiter * Undo changes in elixirs * Implement suggestions * Skip positivity limiting if local limiting is more restrictive * Reduce allocations * Pass variables as strings instead of ints * Add `_nonperiodic` to elixir name * Fix unit test * Implement suggestions * Add missing comma in export of bounds check deviation * Implement suggestions --------- Co-authored-by: Michael Schlottke-Lakemper <[email protected]> Co-authored-by: Andrés Rueda-Ramírez <[email protected]>
trixi-framework · Nov 13, 2023 · 7bb3b46 · 7bb3b46
1 parent d772bf8
commit 7bb3b46
Show file tree

Hide file tree

Showing 17 changed files with 744 additions and 30 deletions.
diff --git a/NEWS.md b/NEWS.md
@@ -38,7 +38,8 @@ for human readability.
 - Wetting and drying feature and examples for 1D and 2D shallow water equations
 - Implementation of the polytropic Euler equations in 2D
 - Implementation of the quasi-1D shallow water equations
-- Subcell positivity limiting support for conservative variables in 2D for `TreeMesh`
+- Subcell (positivity and local min/max) limiting support for conservative variables
+  in 2D for `TreeMesh`
 - AMR for hyperbolic-parabolic equations on 2D/3D `TreeMesh`
 - Added `GradientVariables` type parameter to `AbstractEquationsParabolic`
 

diff --git a/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl b/examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell_nonperiodic.jl
@@ -0,0 +1,93 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+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
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+
+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 = 6,
+                n_cells_max = 10_000,
+                periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_condition)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 100,
+                                     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)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+
+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
+                  save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl
@@ -0,0 +1,91 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+gamma = 1.4
+equations = CompressibleEulerEquations2D(gamma)
+
+"""
+    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+
+The Sedov blast wave setup based on Flash
+- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
+"""
+function initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations2D)
+    # 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)
+
+    # Setup based on https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
+    r0 = 0.21875 # = 3.5 * smallest dx (for domain length=4 and max-ref=6)
+    # r0 = 0.5 # = more reasonable setup
+    E = 1.0
+    p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
+    p0_outer = 1.0e-5 # = true Sedov setup
+    # p0_outer = 1.0e-3 # = more reasonable setup
+
+    # Calculate primitive variables
+    rho = 1.0
+    v1 = 0.0
+    v2 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, p), equations)
+end
+initial_condition = initial_condition_sedov_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_chandrashekar
+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 = 3,
+                n_cells_max = 100_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.0)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 1000
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 1000,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.6)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        stepsize_callback,
+                        save_solution)
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+
+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
+                  save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_euler_shockcapturing_subcell.jl
@@ -39,7 +39,7 @@ surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
 basis = LobattoLegendreBasis(3)
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                positivity_variables_cons = [1],
+                                positivity_variables_cons = ["rho"],
                                 positivity_correction_factor = 0.5)
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,

diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl
@@ -0,0 +1,142 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler multicomponent equations
+
+# 1) Dry Air  2) Helium + 28% Air
+equations = CompressibleEulerMulticomponentEquations2D(gammas = (1.4, 1.648),
+                                                       gas_constants = (0.287, 1.578))
+
+"""
+    initial_condition_shock_bubble(x, t, equations::CompressibleEulerMulticomponentEquations2D{5, 2})
+
+A shock-bubble testcase for multicomponent Euler equations
+- Ayoub Gouasmi, Karthik Duraisamy, Scott Murman
+  Formulation of Entropy-Stable schemes for the multicomponent compressible Euler equations
+  [arXiv: 1904.00972](https://arxiv.org/abs/1904.00972)
+"""
+function initial_condition_shock_bubble(x, t,
+                                        equations::CompressibleEulerMulticomponentEquations2D{
+                                                                                              5,
+                                                                                              2
+                                                                                              })
+    # bubble test case, see Gouasmi et al. https://arxiv.org/pdf/1904.00972
+    # other reference: https://www.researchgate.net/profile/Pep_Mulet/publication/222675930_A_flux-split_algorithm_applied_to_conservative_models_for_multicomponent_compressible_flows/links/568da54508aeaa1481ae7af0.pdf
+    # typical domain is rectangular, we change it to a square, as Trixi can only do squares
+    @unpack gas_constants = equations
+
+    # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
+    delta = 0.03
+
+    # Region I
+    rho1_1 = delta
+    rho2_1 = 1.225 * gas_constants[1] / gas_constants[2] - delta
+    v1_1 = zero(delta)
+    v2_1 = zero(delta)
+    p_1 = 101325
+
+    # Region II
+    rho1_2 = 1.225 - delta
+    rho2_2 = delta
+    v1_2 = zero(delta)
+    v2_2 = zero(delta)
+    p_2 = 101325
+
+    # Region III
+    rho1_3 = 1.6861 - delta
+    rho2_3 = delta
+    v1_3 = -113.5243
+    v2_3 = zero(delta)
+    p_3 = 159060
+
+    # Set up Region I & II:
+    inicenter = SVector(zero(delta), zero(delta))
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    r = sqrt(x_norm^2 + y_norm^2)
+
+    if (x[1] > 0.50)
+        # Set up Region III
+        rho1 = rho1_3
+        rho2 = rho2_3
+        v1 = v1_3
+        v2 = v2_3
+        p = p_3
+    elseif (r < 0.25)
+        # Set up Region I
+        rho1 = rho1_1
+        rho2 = rho2_1
+        v1 = v1_1
+        v2 = v2_1
+        p = p_1
+    else
+        # Set up Region II
+        rho1 = rho1_2
+        rho2 = rho2_2
+        v1 = v1_2
+        v2 = v2_2
+        p = p_2
+    end
+
+    return prim2cons(SVector(v1, v2, p, rho1, rho2), equations)
+end
+initial_condition = initial_condition_shock_bubble
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+basis = LobattoLegendreBasis(3)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                local_minmax_variables_cons = ["rho" * string(i)
+                                                               for i in eachcomponent(equations)])
+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.25, -2.225)
+coordinates_max = (2.20, 2.225)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 3,
+                n_cells_max = 1_000_000)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.01)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 300
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_integrals = (Trixi.density,))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 600,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     solution_variables = cons2prim)
+
+stepsize_callback = StepsizeCallback(cfl = 0.5)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback(save_errors = false))
+
+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
+                  save_everystep = false, callback = callbacks);
+summary_callback() # print the timer summary
diff --git a/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl b/examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl
@@ -88,9 +88,8 @@ volume_flux = flux_ranocha
 basis = LobattoLegendreBasis(3)
 
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                positivity_variables_cons = [
-                                    (i + 3 for i in eachcomponent(equations))...,
-                                ])
+                                positivity_variables_cons = ["rho" * string(i)
+                                                             for i in eachcomponent(equations)])
 
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,

diff --git a/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl b/examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl
@@ -51,7 +51,7 @@ volume_flux = (flux_derigs_etal, flux_nonconservative_powell_local_symmetric)
 basis = LobattoLegendreBasis(3)
 
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                positivity_variables_cons = [1],
+                                positivity_variables_cons = ["rho"],
                                 positivity_correction_factor = 0.5)
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,