Add two-sided Zalesak-type IDP subcell limiting

NEWS.md

@@ -12,7 +12,8 @@ for human readability.
 - Non-uniform `TreeMesh` available for hyperbolic-parabolic equations.
 - Capability to set truly discontinuous initial conditions in 1D.
 - Wetting and drying feature and examples for 1D and 2D 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`
 
 #### Changed
 

examples/tree_2d_dgsem/elixir_euler_blast_wave_sc_subcell.jl

@@ -0,0 +1,96 @@
+
+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=[1])
+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(),)
+
+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

examples/tree_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -0,0 +1,92 @@
+
+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=[1])
+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(),)
+
+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

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_minmax.jl

@@ -0,0 +1,139 @@
+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
+  (; 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=[(i+3 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(),)
+
+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

examples/tree_2d_dgsem/elixir_eulermulti_shock_bubble_shockcapturing_subcell_positivity.jl

@@ -20,7 +20,7 @@ function initial_condition_shock_bubble(x, t, equations::CompressibleEulerMultic
   # 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
+  (; gas_constants) = equations
 
   # Positivity Preserving Parameter, can be set to zero if scheme is positivity preserving
   delta   = 0.03
@@ -137,4 +137,4 @@ 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
                   save_everystep=false, callback=callbacks);
-summary_callback() # print the timer summary
+summary_callback() # print the timer summary

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -6,9 +6,9 @@
 #! format: noindent
 
 function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
-    @unpack inverse_weights = dg.basis
-    @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes
-    @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
+    (; antidiffusive_flux1, antidiffusive_flux2) = cache.antidiffusive_fluxes
+    (; inverse_weights) = dg.basis
+    (; alpha1, alpha2) = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
 
     @threaded for element in eachelement(dg, cache)
         # Sign switch as in apply_jacobian!

src/solvers/dg.jl

@@ -181,6 +181,12 @@ end
 A subcell limiting volume integral type for DG methods based on subcell blending approaches
 with a low-order FV method. Used with limiter [`SubcellLimiterIDP`](@ref).
 
+!!! note
+    Subcell limiting methods are not fully functional on non-conforming meshes. This is
+    mainly because the implementation assumes that low- and high-order schemes have the same
+    surface terms, which is not guaranteed for non-conforming meshes. The low-order scheme
+    with a high-order mortar is not invariant domain preserving.
+
 !!! warning "Experimental implementation"
     This is an experimental feature and may change in future releases.
 """

src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl

@@ -190,4 +190,16 @@ end
 
     return nothing
 end
+
+@inline function get_boundary_outer_state(u_inner, cache, t,
+                                          boundary_condition::BoundaryConditionDirichlet,
+                                          orientation_or_normal, direction, equations,
+                                          dg, indices...)
+    (; node_coordinates) = cache.elements
+
+    x = get_node_coords(node_coordinates, equations, dg, indices...)
+    u_outer = boundary_condition.boundary_value_function(x, t, equations)
+
+    return u_outer
+end
 end # @muladd