Merge branch 'main' into subcell-limiting

examples/p4est_2d_dgsem/elixir_euler_sedov_blast_wave_sc_subcell.jl

@@ -45,7 +45,8 @@ 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.
+                                max_iterations_newton = 40, # Default parameters are not sufficient to fulfill bounds properly.
+                                newton_tolerances = (1.0e-14, 1.0e-15),
                                 bar_states = false)
 
 volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_MCL.jl

@@ -13,7 +13,6 @@
 #
 # Keywords: supersonic flow, shock capturing, AMR, unstructured curved mesh, positivity preservation, compressible Euler, 2D
 
-using Downloads: download
 using OrdinaryDiffEq
 using Trixi
 
@@ -48,6 +47,9 @@ initial_condition = initial_condition_mach3_flow
     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`.
 @inline function Trixi.get_boundary_outer_state(u_inner, t,
                                                 boundary_condition::typeof(boundary_condition_supersonic_inflow),
                                                 normal_direction::AbstractVector,
@@ -77,6 +79,21 @@ end
     return u_inner
 end
 
+# In `compressible_euler_2d.jl`
+# @inline function Trixi.get_boundary_outer_state(u_inner, t,
+#                                                 boundary_condition::typeof(boundary_condition_slip_wall),
+#                                                 normal_direction::AbstractVector,
+#                                                 mesh::P4estMesh{2}, 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_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)
 
 # boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,
@@ -110,11 +127,8 @@ volume_integral = VolumeIntegralSubcellLimiting(limiter_mcl;
 solver = DGSEM(basis, surface_flux, volume_integral)
 
 # Get the unstructured quad mesh from a file (downloads the file if not available locally)
-default_mesh_file = joinpath(@__DIR__, "abaqus_cylinder_in_channel.inp")
-isfile(default_mesh_file) ||
-    download("https://gist.githubusercontent.com/andrewwinters5000/a08f78f6b185b63c3baeff911a63f628/raw/addac716ea0541f588b9d2bd3f92f643eb27b88f/abaqus_cylinder_in_channel.inp",
-             default_mesh_file)
-mesh_file = default_mesh_file
+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)
 

examples/p4est_2d_dgsem/elixir_euler_supersonic_cylinder_sc_subcell.jl

@@ -47,6 +47,9 @@ initial_condition = initial_condition_mach3_flow
     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`.
 @inline function Trixi.get_boundary_outer_state(u_inner, t,
                                                 boundary_condition::typeof(boundary_condition_supersonic_inflow),
                                                 normal_direction::AbstractVector,
@@ -76,6 +79,21 @@ end
     return u_inner
 end
 
+# In `compressible_euler_2d.jl`
+# @inline function Trixi.get_boundary_outer_state(u_inner, t,
+#                                                 boundary_condition::typeof(boundary_condition_slip_wall),
+#                                                 normal_direction::AbstractVector,
+#                                                 mesh::P4estMesh{2}, 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_condition_inflow_outflow = BoundaryConditionCharacteristic(initial_condition)
 
 # boundary_conditions = Dict(:Bottom => boundary_condition_slip_wall,

src/callbacks_stage/subcell_bounds_check_2d.jl

@@ -5,9 +5,8 @@
 @muladd begin
 #! format: noindent
 
-@inline function check_bounds(u::AbstractArray{<:Any, 4},
-                              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
@@ -104,9 +103,8 @@
     return nothing
 end
 
-@inline function check_bounds(u::AbstractArray{<:Any, 4},
-                              equations, solver, cache,
-                              limiter::SubcellLimiterMCL)
+@inline function check_bounds(u, equations::AbstractEquations{2}, # only works for 2D
+                              solver, cache, limiter::SubcellLimiterMCL)
     (; var_min, var_max) = limiter.cache.subcell_limiter_coefficients
     (; bar_states1, bar_states2, lambda1, lambda2) = limiter.cache.container_bar_states
     (; mcl_bounds_delta_local, mcl_bounds_delta_global) = limiter.cache

src/solvers/dgsem_p4est/subcell_limiters_2d.jl

@@ -93,7 +93,7 @@ end
     (; boundaries) = cache
     index_range = eachnode(dg)
 
-    foreach(enumerate(boundary_condition_types)) do (i, boundary_condition)
+    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
         for boundary in boundary_indices[i]
             element = boundaries.neighbor_ids[boundary]
             node_indices = boundaries.node_indices[boundary]
@@ -216,7 +216,7 @@ end
     (; boundaries) = cache
     index_range = eachnode(dg)
 
-    foreach(enumerate(boundary_condition_types)) do (i, boundary_condition)
+    foreach_enumerate(boundary_condition_types) do (i, boundary_condition)
         for boundary in boundary_indices[i]
             element = boundaries.neighbor_ids[boundary]
             node_indices = boundaries.node_indices[boundary]

src/solvers/dgsem_tree/subcell_limiters.jl

@@ -231,40 +231,32 @@ function Base.show(io::IO, ::MIME"text/plain", limiter::SubcellLimiterIDP)
         else
             setup = ["Limiter" => ""]
             if local_twosided
-                setup = [
-                    setup...,
-                    "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)",
-                ]
+                push!(setup,
+                      "" => "Local two-sided limiting for conservative variables $(limiter.local_twosided_variables_cons)")
             end
             if positivity
                 if !isempty(limiter.positivity_variables_cons)
                     string = "conservative variables $(limiter.positivity_variables_cons)"
-                    setup = [setup..., "" => "Positivity limiting for " * string]
+                    push!(setup, "" => "Positivity limiting for " * string)
                 end
                 if !isempty(limiter.positivity_variables_nonlinear)
                     string = "$(limiter.positivity_variables_nonlinear)"
-                    setup = [setup..., "" => "Positivity limiting for " * string]
+                    push!(setup, "" => "Positivity limiting for " * string)
                 end
-                setup = [
-                    setup...,
-                    "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)",
-                ]
+                push!(setup,
+                      "" => "- with positivity correction factor = $(limiter.positivity_correction_factor)")
             end
             if local_onesided
                 for (variable, min_or_max) in limiter.local_onesided_variables_nonlinear
-                    setup = [setup..., "" => "Local $min_or_max limiting for $variable"]
+                    push!(setup, "" => "Local $min_or_max limiting for $variable")
                 end
             end
-            setup = [
-                setup...,
-                "Local bounds with" => (limiter.bar_states ? "Bar States" :
-                                        "FV solution"),
-            ]
+            push!(setup,
+                  "Local bounds with" => (limiter.bar_states ? "Bar States" :
+                                          "FV solution"))
             if limiter.smoothness_indicator
-                setup = [
-                    setup...,
-                    "Smoothness indicator" => "$(limiter.IndicatorHG) using threshold $(limiter.threshold_smoothness_indicator)",
-                ]
+                push!(setup,
+                      "Smoothness indicator" => "$(limiter.IndicatorHG) using threshold $(limiter.threshold_smoothness_indicator)")
             end
         end
         summary_box(io, "SubcellLimiterIDP", setup)

test/test_p4est_2d.jl

@@ -330,16 +330,16 @@ end
     @test_trixi_include(joinpath(EXAMPLES_DIR,
                                  "elixir_euler_sedov_blast_wave_sc_subcell.jl"),
                         l2=[
-                            0.45737877846538905,
-                            0.2852097276261684,
-                            0.28527281809798694,
-                            1.2881460122856072,
+                            0.4573787784168518,
+                            0.28520972760728397,
+                            0.28527281808006966,
+                            1.2881460122982442,
                         ],
                         linf=[
-                            1.6444110425837906,
-                            1.6743368122678752,
-                            1.6760847980983236,
-                            6.268843623083507,
+                            1.644411040701827,
+                            1.6743368119653912,
+                            1.6760847977977988,
+                            6.268843623142863,
                         ],
                         tspan=(0.0, 0.3))
     # Ensure that we do not have excessive memory allocations