Merge main: Remove dispatching for mesh types

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_sedov_blast_wave_sc_subcell.jl

@@ -0,0 +1,115 @@
+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
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (x_neg = boundary_condition,
+                       x_pos = boundary_condition,
+                       y_neg = boundary_condition,
+                       y_pos = boundary_condition)
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                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 = [],
+                                bar_states = false)
+# 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)
+
+# 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)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 3.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.7)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+###############################################################################
+# run the simulation
+
+stage_callbacks = (SubcellLimiterIDPCorrection(), BoundsCheckCallback())
+
+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

src/callbacks_stage/subcell_bounds_check.jl

@@ -38,7 +38,7 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
     (; t, iter, alg) = integrator
     u = wrap_array(u_ode, mesh, equations, solver, cache)
 
-    @trixi_timeit timer() "check_bounds" check_bounds(u, mesh, equations, solver, cache,
+    @trixi_timeit timer() "check_bounds" check_bounds(u, equations, solver, cache,
                                                       solver.volume_integral)
 
     save_errors = callback.save_errors && (callback.interval > 0) &&
@@ -48,20 +48,20 @@ function (callback::BoundsCheckCallback)(u_ode, integrator, stage)
                    (iter + 1) >= integrator.opts.maxiters)  # Maximum iterations reached
     if save_errors
         @trixi_timeit timer() "save_errors" save_bounds_check_errors(callback.output_directory,
-                                                                     u, t, iter + 1,
+                                                                     t, iter + 1,
                                                                      equations,
                                                                      solver.volume_integral)
     end
 end
 
-@inline function check_bounds(u, mesh, equations, solver, cache,
+@inline function check_bounds(u, equations, solver, cache,
                               volume_integral::VolumeIntegralSubcellLimiting)
-    check_bounds(u, mesh, equations, solver, cache, volume_integral.limiter)
+    check_bounds(u, equations, solver, cache, volume_integral.limiter)
 end
 
-@inline function save_bounds_check_errors(output_directory, u, t, iter, equations,
+@inline function save_bounds_check_errors(output_directory, t, iter, equations,
                                           volume_integral::VolumeIntegralSubcellLimiting)
-    save_bounds_check_errors(output_directory, u, t, iter, equations,
+    save_bounds_check_errors(output_directory, t, iter, equations,
                              volume_integral.limiter)
 end
 

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -5,7 +5,8 @@
 @muladd begin
 #! format: noindent
 
-@inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
+@inline function check_bounds(u::AbstractArray{<:Any, 4},
+                              equations, solver, cache,
                               limiter::SubcellLimiterIDP)
     (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
@@ -103,7 +104,7 @@
     return nothing
 end
 
-@inline function save_bounds_check_errors(output_directory, u, time, iter, equations,
+@inline function save_bounds_check_errors(output_directory, time, iter, equations,
                                           limiter::SubcellLimiterIDP)
     (; local_twosided, positivity, local_onesided) = limiter
     (; idp_bounds_delta_local) = limiter.cache
@@ -145,7 +146,8 @@ end
     return nothing
 end
 
-@inline function check_bounds(u, mesh::AbstractMesh{2}, equations, solver, cache,
+@inline function check_bounds(u::AbstractArray{<:Any, 4},
+                              equations, solver, cache,
                               limiter::SubcellLimiterMCL)
     (; var_min, var_max) = limiter.cache.subcell_limiter_coefficients
     (; bar_states1, bar_states2, lambda1, lambda2) = limiter.cache.container_bar_states
@@ -625,14 +627,14 @@ end
     return nothing
 end
 
-@inline function save_bounds_check_errors(output_directory, u, time, iter, equations,
+@inline function save_bounds_check_errors(output_directory, time, iter, equations,
                                           limiter::SubcellLimiterMCL)
     (; mcl_bounds_delta_local) = limiter.cache
 
     n_vars = nvariables(equations)
 
     # Print errors to output file
-    open("$output_directory/deviations.txt", "a") do f
+    open(joinpath(output_directory, "deviations.txt"), "a") do f
         print(f, iter, ", ", time)
         for v in eachvariable(equations)
             print(f, ", ", mcl_bounds_delta_local[1, v], ", ",

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -5,29 +5,14 @@
 @muladd begin
 #! format: noindent
 
-function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
-    if dg.volume_integral.limiter.smoothness_indicator
-        elements = cache.element_ids_dgfv
-    else
-        elements = eachelement(dg, cache)
-    end
-
-    # Loop over blended DG-FV elements
-    @threaded for element in elements
-        # Sign switch as in apply_jacobian!
-        inverse_jacobian = -cache.elements.inverse_jacobian[element]
-
-        for j in eachnode(dg), i in eachnode(dg)
-            perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache,
-                                          i, j, element)
-        end
-    end
-
-    return nothing
-end
-
-function perform_idp_correction!(u, dt, mesh::Union{StructuredMesh{2}, P4estMesh{2}},
+function perform_idp_correction!(u, dt,
+                                 mesh::Union{TreeMesh{2}, StructuredMesh{2},
+                                             P4estMesh{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
+
     if dg.volume_integral.limiter.smoothness_indicator
         elements = cache.element_ids_dgfv
     else
@@ -37,43 +22,33 @@ function perform_idp_correction!(u, dt, mesh::Union{StructuredMesh{2}, P4estMesh
     @threaded for element in elements
         for j in eachnode(dg), i in eachnode(dg)
             # Sign switch as in apply_jacobian!
-            inverse_jacobian = -cache.elements.inverse_jacobian[i, j, element]
-
-            perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg, cache,
-                                          i, j, element)
+            inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                     mesh, i, j, element)
+
+            # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1}
+            alpha_flux1 = (1 - alpha1[i, j, element]) *
+                          get_node_vars(antidiffusive_flux1_R, equations, dg,
+                                        i, j, element)
+            alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) *
+                              get_node_vars(antidiffusive_flux1_L, equations, dg,
+                                            i + 1, j, element)
+            alpha_flux2 = (1 - alpha2[i, j, element]) *
+                          get_node_vars(antidiffusive_flux2_R, equations, dg,
+                                        i, j, element)
+            alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) *
+                              get_node_vars(antidiffusive_flux2_L, equations, dg,
+                                            i, j + 1, element)
+
+            for v in eachvariable(equations)
+                u[v, i, j, element] += dt * inverse_jacobian *
+                                       (inverse_weights[i] *
+                                        (alpha_flux1_ip1[v] - alpha_flux1[v]) +
+                                        inverse_weights[j] *
+                                        (alpha_flux2_jp1[v] - alpha_flux2[v]))
+            end
         end
     end
 
     return nothing
 end
-
-# Function barrier to dispatch outer function by mesh type
-@inline function perform_idp_correction_inner!(u, dt, inverse_jacobian, equations, dg,
-                                               cache, i, j, element)
-    (; inverse_weights) = dg.basis
-    (; antidiffusive_flux1_L, antidiffusive_flux2_L, antidiffusive_flux1_R, antidiffusive_flux2_R) = cache.antidiffusive_fluxes
-    (; alpha1, alpha2) = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
-
-    # Note: antidiffusive_flux1[v, i, xi, element] = antidiffusive_flux2[v, xi, i, element] = 0 for all i in 1:nnodes and xi in {1, nnodes+1}
-    alpha_flux1 = (1 - alpha1[i, j, element]) *
-                  get_node_vars(antidiffusive_flux1_R, equations, dg, i, j, element)
-    alpha_flux1_ip1 = (1 - alpha1[i + 1, j, element]) *
-                      get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1, j,
-                                    element)
-    alpha_flux2 = (1 - alpha2[i, j, element]) *
-                  get_node_vars(antidiffusive_flux2_R, equations, dg, i, j, element)
-    alpha_flux2_jp1 = (1 - alpha2[i, j + 1, element]) *
-                      get_node_vars(antidiffusive_flux2_L, equations, dg, i, j + 1,
-                                    element)
-
-    for v in eachvariable(equations)
-        u[v, i, j, element] += dt * inverse_jacobian *
-                               (inverse_weights[i] *
-                                (alpha_flux1_ip1[v] - alpha_flux1[v]) +
-                                inverse_weights[j] *
-                                (alpha_flux2_jp1[v] - alpha_flux2[v]))
-    end
-
-    return nothing
-end
 end # @muladd

src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl

@@ -508,7 +508,9 @@ 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::Union{TreeMesh{2}, StructuredMesh{2},
+                                                     P4estMesh{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
@@ -546,7 +548,9 @@ 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::Union{TreeMesh{2}, StructuredMesh{2},
+                                                     P4estMesh{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