Add subcell limiting support for StructuredMesh

NEWS.md

@@ -14,6 +14,7 @@ for human readability.
 - Subcell local one-sided limiting support for nonlinear variables in 2D for `TreeMesh` ([#1792]).
 - New time integrator `PairedExplicitRK2`, implementing the second-order paired explicit Runge-Kutta
   method with [Convex.jl](https://github.com/jump-dev/Convex.jl) and [ECOS.jl](https://github.com/jump-dev/ECOS.jl) ([#1908])
+- Add subcell limiting support for `StructuredMesh` ([#1946]).
 
 ## Changes when updating to v0.7 from v0.6.x
 

examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl → examples/structured_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -1,57 +1,87 @@
-
 using OrdinaryDiffEq
 using Trixi
 
 ###############################################################################
 # semidiscretization of the compressible Euler equations
-
-equations = CompressibleEulerEquations2D(1.4)
-
-initial_condition = initial_condition_convergence_test
-source_terms = source_terms_convergence_test
+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
 
 boundary_condition = BoundaryConditionDirichlet(initial_condition)
 boundary_conditions = (x_neg = boundary_condition,
                        x_pos = boundary_condition,
                        y_neg = boundary_condition,
                        y_pos = boundary_condition)
 
-# Get the DG approximation space
 surface_flux = flux_lax_friedrichs
 volume_flux = flux_ranocha
 polydeg = 3
 basis = LobattoLegendreBasis(polydeg)
 limiter_idp = SubcellLimiterIDP(equations, basis;
-                                positivity_variables_cons = ["rho"],
-                                positivity_variables_nonlinear = [pressure],
-                                local_twosided_variables_cons = [],
-                                local_onesided_variables_nonlinear = [])
-# Variables for local limiting (`local_twosided_variables_cons` and
-# `local_onesided_variables_nonlinear`) are overwritten and used in the tests.
+                                local_twosided_variables_cons = ["rho"],
+                                local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
+                                                                       min)],
+                                max_iterations_newton = 40, # Default value of 10 iterations is too low to fulfill bounds.
+                                positivity_variables_cons = [],
+                                positivity_variables_nonlinear = [])
+# Variables for global limiting (`positivity_variables_cons` and
+# `positivity_variables_nonlinear`) are overwritten and used in the tests.
 
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
                                                 volume_flux_dg = volume_flux,
                                                 volume_flux_fv = surface_flux)
 solver = DGSEM(basis, surface_flux, volume_integral)
 
-# Waving flag
-f1(s) = SVector(-1.0, s - 1.0)
-f2(s) = SVector(1.0, s + 1.0)
-f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
-f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s))
-mapping = Trixi.transfinite_mapping((f1, f2, f3, f4))
+# Get the curved quad mesh from a mapping function
+# Mapping as described in https://arxiv.org/abs/2012.12040
+function mapping(xi, eta)
+    y = eta + 0.125 * (cos(1.5 * pi * xi) * cos(0.5 * pi * eta))
+
+    x = xi + 0.125 * (cos(0.5 * pi * xi) * cos(2 * pi * y))
+
+    return SVector(x, y)
+end
 
 cells_per_dimension = (16, 16)
 mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false)
 
 semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
-                                    boundary_conditions = boundary_conditions,
-                                    source_terms = source_terms)
+                                    boundary_conditions = boundary_conditions)
 
 ###############################################################################
 # ODE solvers, callbacks etc.
 
-tspan = (0.0, 2.0)
+tspan = (0.0, 3.0)
 ode = semidiscretize(semi, tspan)
 
 summary_callback = SummaryCallback()
@@ -66,7 +96,7 @@ save_solution = SaveSolutionCallback(interval = 100,
                                      save_final_solution = true,
                                      solution_variables = cons2prim)
 
-stepsize_callback = StepsizeCallback(cfl = 0.8)
+stepsize_callback = StepsizeCallback(cfl = 0.7)
 
 callbacks = CallbackSet(summary_callback,
                         analysis_callback, alive_callback,

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -5,7 +5,7 @@
 @muladd begin
 #! format: noindent
 
-function perform_idp_correction!(u, dt, mesh, equations, dg, cache)
+function perform_idp_correction!(u, dt, mesh::AbstractMesh{2}, equations, dg, cache)
     @unpack inverse_weights = dg.basis
     @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes
     @unpack alpha1, alpha2 = dg.volume_integral.limiter.cache.subcell_limiter_coefficients

src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl

@@ -5,7 +5,7 @@
 @muladd begin
 #! format: noindent
 
-function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+function create_cache(mesh::AbstractMesh{2},
                       equations, volume_integral::VolumeIntegralSubcellLimiting,
                       dg::DG, uEltype)
     cache = create_cache(mesh, equations,
@@ -57,7 +57,7 @@ function create_cache(mesh::Union{TreeMesh{2}, StructuredMesh{2}},
 end
 
 function calc_volume_integral!(du, u,
-                               mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+                               mesh::Union{TreeMesh, StructuredMesh},
                                nonconservative_terms, equations,
                                volume_integral::VolumeIntegralSubcellLimiting,
                                dg::DGSEM, cache)
@@ -392,7 +392,7 @@ end
 # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems.
 @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R,
                                          fstar1_L, fstar1_R, fstar2_L, fstar2_R,
-                                         u, mesh,
+                                         u, mesh::AbstractMesh{2},
                                          nonconservative_terms::False, equations,
                                          limiter::SubcellLimiterIDP, dg, element, cache)
     @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes
@@ -430,7 +430,7 @@ end
 # Calculate the antidiffusive flux `antidiffusive_flux` as the subtraction between `fhat` and `fstar` for conservative systems.
 @inline function calcflux_antidiffusive!(fhat1_L, fhat1_R, fhat2_L, fhat2_R,
                                          fstar1_L, fstar1_R, fstar2_L, fstar2_R,
-                                         u, mesh,
+                                         u, mesh::AbstractMesh{2},
                                          nonconservative_terms::True, equations,
                                          limiter::SubcellLimiterIDP, dg, element, cache)
     @unpack antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R = cache.antidiffusive_fluxes

test/test_structured_2d.jl

@@ -475,20 +475,20 @@ end
     end
 end
 
-@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (global bounds)" begin
+@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (local bounds)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
-                                 "elixir_euler_source_terms_sc_subcell.jl"),
+                                 "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
                         l2=[
-                            7.939813971944535e-5,
-                            3.6805976512737965e-5,
-                            7.50142174035037e-5,
-                            0.0001576144969663713,
+                            0.6337774834710513,
+                            0.30377119245852724,
+                            0.3111372568571772,
+                            1.2976221893997268,
                         ],
                         linf=[
-                            0.0007421863054295486,
-                            0.00037635704053817776,
-                            0.0007278708058917616,
-                            0.001550366737408826,
+                            2.2064877103138207,
+                            1.541067099687334,
+                            1.5487587769900337,
+                            6.271271639873466,
                         ],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations
@@ -501,26 +501,24 @@ end
     end
 end
 
-@trixi_testset "elixir_euler_source_terms_sc_subcell.jl (local bounds)" begin
+@trixi_testset "elixir_euler_sedov_blast_wave_sc_subcell.jl (global bounds)" begin
     @test_trixi_include(joinpath(EXAMPLES_DIR,
-                                 "elixir_euler_source_terms_sc_subcell.jl"),
-                        positivity_variables_cons=[],
-                        positivity_variables_nonlinear=[],
-                        local_twosided_variables_cons=["rho"],
-                        local_onesided_variables_nonlinear=[(Trixi.entropy_guermond_etal,
-                                                             min)],
-                        cfl=0.5,
+                                 "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
+                        positivity_variables_cons=["rho"],
+                        positivity_variables_nonlinear=[pressure],
+                        local_twosided_variables_cons=[],
+                        local_onesided_variables_nonlinear=[],
                         l2=[
-                            0.007788374804566746,
-                            0.006564110929038488,
-                            0.00841157867996223,
-                            0.023360113888099734,
+                            0.7869912572385168,
+                            0.39170886758882073,
+                            0.39613257454431977,
+                            1.2951760266455101,
                         ],
                         linf=[
-                            0.033816854019882214,
-                            0.03938807953377843,
-                            0.044093017726981376,
-                            0.08006778793724179,
+                            5.156044534854053,
+                            3.6261667239538986,
+                            3.1807681416546085,
+                            6.3028422220287235,
                         ],
                         tspan=(0.0, 0.5))
     # Ensure that we do not have excessive memory allocations