Add subcell limiting support for P4estMesh (trixi-framework#1954)

* First part of P4estMesh support

* Adapt `get_boundary_outer_state` to allow supersoniv cylinder eilixir

* Add test

* Rename routine

* Rename other routine

* Adapt elixir and complete testing

* Adapt tests

* Adapt cfl in elixir

* Fix test

* Add sedov elixir to test periodic boundaries

* fmt

* Activate test for coverage run

* Remove extra lines of code

* Implement suggestions

* Adapt parameters to reduce bounds checking errors

* Implement first suggestions

* Remove `mesh` from `get_boundary_outer_state`

* Use `foreach_enumerate` to remove allocations

* Move `get_boundary_outer_state` to elixir

---------

Co-authored-by: Michael Schlottke-Lakemper <[email protected]>

examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -0,0 +1,102 @@
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+"""
+    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)
+    E = 1.0
+    p0_inner = 3 * (equations.gamma - 1) * E / (3 * pi * r0^2)
+    p0_outer = 1.0e-5 # = true Sedov 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
+
+# Get the DG approximation space
+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 parameters are not sufficient to fulfill bounds properly.
+                                newton_tolerances = (1.0e-14, 1.0e-15))
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+###############################################################################
+
+coordinates_min = (-1.0, -1.0)
+coordinates_max = (1.0, 1.0)
+
+trees_per_dimension = (4, 4)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = polydeg, initial_refinement_level = 2,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+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 = 300
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 300,
+                                     save_initial_solution = true,
+                                     save_final_solution = true)
+
+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())
+
+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
+                  callback = callbacks);
+summary_callback() # print the timer summary

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl

@@ -0,0 +1,162 @@
+# Channel flow around a cylinder at Mach 3
+#
+# Boundary conditions are supersonic Mach 3 inflow at the left portion of the domain
+# and supersonic outflow at the right portion of the domain. The top and bottom of the
+# channel as well as the cylinder are treated as Euler slip wall boundaries.
+# This flow results in strong shock reflections / interactions as well as Kelvin-Helmholtz
+# instabilities at later times as two Mach stems form above and below the cylinder.
+#
+# For complete details on the problem setup see Section 5.7 of the paper:
+# - Jean-Luc Guermond, Murtazo Nazarov, Bojan Popov, and Ignacio Tomas (2018)
+#   Second-Order Invariant Domain Preserving Approximation of the Euler Equations using Convex Limiting.
+#   [DOI: 10.1137/17M1149961](https://doi.org/10.1137/17M1149961)
+#
+# Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+@inline function initial_condition_mach3_flow(x, t, equations::CompressibleEulerEquations2D)
+    # set the freestream flow parameters
+    rho_freestream = 1.4
+    v1 = 3.0
+    v2 = 0.0
+    p_freestream = 1.0
+
+    prim = SVector(rho_freestream, v1, v2, p_freestream)
+    return prim2cons(prim, equations)
+end
+
+initial_condition = initial_condition_mach3_flow
+
+# Supersonic inflow boundary condition.
+# Calculate the boundary flux entirely from the external solution state, i.e., set
+# external solution state values for everything entering the domain.
+@inline function boundary_condition_supersonic_inflow(u_inner,
+                                                      normal_direction::AbstractVector,
+                                                      x, t, surface_flux_function,
+                                                      equations::CompressibleEulerEquations2D)
+    u_boundary = initial_condition_mach3_flow(x, t, equations)
+    flux = Trixi.flux(u_boundary, normal_direction, equations)
+
+    return flux
+end
+
+# For subcell limiting, the calculation of local bounds for non-periodic domains requires the
+# boundary outer state. Those functions return the boundary value for a specific boundary condition
+# at time `t`, for the node with spatial indices `indices` and the given `normal_direction`.
+# only for P4estMesh{2}
+@inline function Trixi.get_boundary_outer_state(u_inner, t,
+                                                boundary_condition::typeof(boundary_condition_supersonic_inflow),
+                                                normal_direction::AbstractVector,
+                                                equations, dg, cache,
+                                                indices...)
+    x = Trixi.get_node_coords(cache.elements.node_coordinates, equations, dg, indices...)
+
+    return initial_condition_mach3_flow(x, t, equations)
+end
+
+# Supersonic outflow boundary condition.
+# Calculate the boundary flux entirely from the internal solution state. Analogous to supersonic inflow
+# except all the solution state values are set from the internal solution as everything leaves the domain
+@inline function boundary_condition_outflow(u_inner, normal_direction::AbstractVector, x, t,
+                                            surface_flux_function,
+                                            equations::CompressibleEulerEquations2D)
+    flux = Trixi.flux(u_inner, normal_direction, equations)
+
+    return flux
+end
+
+# only for P4estMesh{2}
+@inline function Trixi.get_boundary_outer_state(u_inner, t,
+                                                boundary_condition::typeof(boundary_condition_outflow),
+                                                normal_direction::AbstractVector,
+                                                equations, dg, cache,
+                                                indices...)
+    return u_inner
+end
+
+# only for P4estMesh{2}
+@inline function Trixi.get_boundary_outer_state(u_inner, t,
+                                                boundary_condition::typeof(boundary_condition_slip_wall),
+                                                normal_direction::AbstractVector,
+                                                equations::CompressibleEulerEquations2D,
+                                                dg, cache, indices...)
+    factor = (normal_direction[1] * u_inner[2] + normal_direction[2] * u_inner[3])
+    u_normal = (factor / sum(normal_direction .^ 2)) * normal_direction
+
+    return SVector(u_inner[1],
+                   u_inner[2] - 2 * u_normal[1],
+                   u_inner[3] - 2 * u_normal[2],
+                   u_inner[4])
+end
+
+boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,
+                           :Circle => boundary_condition_slip_wall,
+                           :Top => boundary_condition_slip_wall,
+                           :Right => boundary_condition_outflow,
+                           :Left => boundary_condition_supersonic_inflow)
+
+volume_flux = flux_ranocha_turbo
+surface_flux = flux_lax_friedrichs
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                local_twosided_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure],
+                                local_onesided_variables_nonlinear = [(Trixi.entropy_guermond_etal,
+                                                                       min)],
+                                max_iterations_newton = 50) # Default value of 10 iterations is too low to fulfill bounds.
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Get the unstructured quad mesh from a file (downloads the file if not available locally)
+mesh_file = Trixi.download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp",
+                           joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp"))
+
+mesh = P4estMesh{2}(mesh_file, initial_refinement_level = 0)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions)
+
+###############################################################################
+# ODE solvers
+
+tspan = (0.0, 2.0)
+ode = semidiscretize(semi, tspan)
+
+# Callbacks
+
+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.8)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback, alive_callback,
+                        save_solution,
+                        stepsize_callback)
+
+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
+                  callback = callbacks);
+summary_callback() # print the timer summary

src/callbacks_stage/subcell_bounds_check.jl

@@ -171,5 +171,47 @@ end
     return nothing
 end
 
+@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
+
+    # Print to output file
+    open(joinpath(output_directory, "deviations.txt"), "a") do f
+        print(f, iter, ", ", time)
+        if local_twosided
+            for v in limiter.local_twosided_variables_cons
+                v_string = string(v)
+                print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")],
+                      ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
+            end
+        end
+        if local_onesided
+            for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
+                key = Symbol(string(variable), "_", string(min_or_max))
+                print(f, ", ", idp_bounds_delta_local[key])
+            end
+        end
+        if positivity
+            for v in limiter.positivity_variables_cons
+                if v in limiter.local_twosided_variables_cons
+                    continue
+                end
+                print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])
+            end
+            for variable in limiter.positivity_variables_nonlinear
+                print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")])
+            end
+        end
+        println(f)
+    end
+    # Reset local maximum deviations
+    for (key, _) in idp_bounds_delta_local
+        idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key]))
+    end
+
+    return nothing
+end
+
 include("subcell_bounds_check_2d.jl")
 end # @muladd

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -5,8 +5,8 @@
 @muladd begin
 #! format: noindent
 
-@inline function check_bounds(u, equations, solver, cache,
-                              limiter::SubcellLimiterIDP)
+@inline function check_bounds(u, equations::AbstractEquations{2}, # only works for 2D
+                              solver, cache, limiter::SubcellLimiterIDP)
     (; local_twosided, positivity, local_onesided) = solver.volume_integral.limiter
     (; variable_bounds) = limiter.cache.subcell_limiter_coefficients
     (; idp_bounds_delta_local, idp_bounds_delta_global) = limiter.cache
@@ -102,47 +102,4 @@
 
     return nothing
 end
-
-@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
-
-    # Print to output file
-    open(joinpath(output_directory, "deviations.txt"), "a") do f
-        print(f, iter, ", ", time)
-        if local_twosided
-            for v in limiter.local_twosided_variables_cons
-                v_string = string(v)
-                print(f, ", ", idp_bounds_delta_local[Symbol(v_string, "_min")],
-                      ", ", idp_bounds_delta_local[Symbol(v_string, "_max")])
-            end
-        end
-        if local_onesided
-            for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
-                key = Symbol(string(variable), "_", string(min_or_max))
-                print(f, ", ", idp_bounds_delta_local[key])
-            end
-        end
-        if positivity
-            for v in limiter.positivity_variables_cons
-                if v in limiter.local_twosided_variables_cons
-                    continue
-                end
-                print(f, ", ", idp_bounds_delta_local[Symbol(string(v), "_min")])
-            end
-            for variable in limiter.positivity_variables_nonlinear
-                print(f, ", ", idp_bounds_delta_local[Symbol(string(variable), "_min")])
-            end
-        end
-        println(f)
-    end
-
-    # Reset local maximum deviations
-    for (key, _) in idp_bounds_delta_local
-        idp_bounds_delta_local[key] = zero(eltype(idp_bounds_delta_local[key]))
-    end
-
-    return nothing
-end
 end # @muladd

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -6,7 +6,8 @@
 #! format: noindent
 
 function perform_idp_correction!(u, dt,
-                                 mesh::Union{TreeMesh{2}, StructuredMesh{2}},
+                                 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

src/solvers/dgsem_p4est/dg.jl

@@ -53,4 +53,6 @@ include("dg_2d_parabolic.jl")
 include("dg_3d.jl")
 include("dg_3d_parabolic.jl")
 include("dg_parallel.jl")
+
+include("subcell_limiters_2d.jl")
 end # @muladd