Example Zhang-Shu Positivity-Preserving Limiter for MHD (#2637)

* Example Zhang-Shu Positivity-Preserving Limiter for MHD

* coveragge

* mv

* rm

Author: DanielDoehring

examples/tree_1d_dgsem/elixir_mhd_briowu_shock_tube_mod_positivity.jl

@@ -0,0 +1,107 @@
+using OrdinaryDiffEqSSPRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the ideal MHD equations
+gamma = 2
+equations = IdealGlmMhdEquations1D(gamma)
+
+"""
+    initial_condition_briowu_shock_tube_mod(x, t, equations::IdealGlmMhdEquations1D)
+
+Compound shock tube test case for one dimensional ideal MHD equations.
+It is basically an MHD extension of the Sod shock tube.
+Taken from Section V of the article with a modified lower right pressure (1e-3 instead of 0.1)
+- Brio and Wu (1988)
+  An Upwind Differencing Scheme for the Equations of Ideal Magnetohydrodynamics
+  [DOI: 10.1016/0021-9991(88)90120-9](https://doi.org/10.1016/0021-9991(88)90120-9)
+"""
+function initial_condition_briowu_shock_tube_mod(x, t, equations::IdealGlmMhdEquations1D)
+    # domain must be set to [0, 1], γ = 2, final time = 0.12
+    RealT = eltype(x)
+    rho = x[1] < 0.5f0 ? 1.0f0 : 0.125f0
+    v1 = 0
+    v2 = 0
+    v3 = 0
+    p = x[1] < 0.5f0 ? one(RealT) : RealT(1e-3) # 1e-3 instead of 0.1
+    B1 = 0.75f0
+    B2 = x[1] < 0.5f0 ? 1 : -1
+    B3 = 0
+    return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3), equations)
+end
+initial_condition = initial_condition_briowu_shock_tube_mod
+
+boundary_conditions = BoundaryConditionDirichlet(initial_condition)
+
+surface_flux = flux_hlle
+volume_flux = flux_derigs_etal
+basis = LobattoLegendreBasis(4)
+
+indicator_sc = IndicatorHennemannGassner(equations, basis,
+                                         alpha_max = 0.5,
+                                         alpha_min = 0.001,
+                                         alpha_smooth = true,
+                                         variable = density_pressure)
+volume_integral = VolumeIntegralShockCapturingHG(indicator_sc;
+                                                 volume_flux_dg = volume_flux,
+                                                 volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = 0.0
+coordinates_max = 1.0
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level = 4,
+                n_cells_max = 10_000,
+                periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.12)
+ode = semidiscretize(semi, tspan)
+
+summary_callback = SummaryCallback()
+
+analysis_interval = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval,
+                                     extra_analysis_errors = (:l2_error_primitive,
+                                                              :linf_error_primitive))
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_restart = SaveRestartCallback(interval = 100,
+                                   save_final_restart = true)
+
+amr_indicator = IndicatorHennemannGassner(semi,
+                                          alpha_max = 0.5,
+                                          alpha_min = 0.001,
+                                          alpha_smooth = true,
+                                          variable = density_pressure)
+amr_controller = ControllerThreeLevel(semi, amr_indicator,
+                                      base_level = 4,
+                                      max_level = 6, max_threshold = 0.01)
+amr_callback = AMRCallback(semi, amr_controller,
+                           interval = 5,
+                           adapt_initial_condition = true,
+                           adapt_initial_condition_only_refine = true)
+
+stepsize_callback = StepsizeCallback(cfl = 1.1)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_restart,
+                        amr_callback, stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+# crashes without positivity limiter
+stage_limiter! = PositivityPreservingLimiterZhangShu(thresholds = (5.0e-6,),
+                                                     variables = (pressure,))
+
+sol = solve(ode, SSPRK54(stage_limiter! = stage_limiter!);
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            ode_default_options()..., callback = callbacks);

test/test_tree_1d_mhd.jl

@@ -167,6 +167,34 @@ end
     @test_allocations(Trixi.rhs!, semi, sol, 1000)
 end
 
+@trixi_testset "elixir_mhd_briowu_shock_tube_mod_positivity.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR,
+                                 "elixir_mhd_briowu_shock_tube_mod_positivity.jl"),
+                        l2=[
+                            0.17958200783518533,
+                            0.21403217755031365,
+                            0.3528008349880033,
+                            0.0,
+                            0.40156743941762096,
+                            7.603990369710917e-14,
+                            0.34299587882334953,
+                            0.0
+                        ],
+                        linf=[
+                            0.5347779567211823,
+                            0.4651105133406562,
+                            0.9873266889381644,
+                            0.0,
+                            1.000540495846343,
+                            8.615330671091215e-14,
+                            1.4425998855617321,
+                            0.0
+                        ])
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    @test_allocations(Trixi.rhs!, semi, sol, 1000)
+end
+
 @trixi_testset "elixir_mhd_torrilhon_shock_tube.jl" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_torrilhon_shock_tube.jl"),
                         l2=[