First PR for 3D subcell limiting

NEWS.md

@@ -6,6 +6,14 @@ used in the Julia ecosystem. Notable changes will be documented in this file
 for human readability.
 
 
+## Changes in the v0.13 lifecycle
+
+#### Added
+- Initial 3D support for subcell limiting with `P4estMesh` was added ([#2582]).
+  In the new version, positivity limiting for conservative variables (using `positivity_variables_cons`)
+  is supported. `BoundsCheckCallback` is not supported in 3D yet.
+
+
 ## Changes when updating to v0.13 from v0.12.x
 
 #### Changed

examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl

@@ -0,0 +1,98 @@
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+"""
+    initial_condition_sedov_blast_wave(x, t, equations::CompressibleEulerEquations3D)
+
+The Sedov blast wave setup based on Flash
+- https://flash.rochester.edu/site/flashcode/user_support/flash_ug_devel/node187.html#SECTION010114000000000000000
+with smaller strength of the initial discontinuity.
+"""
+function initial_condition_sedov_blast_wave(x, t,
+                                            equations::CompressibleEulerEquations3D)
+    # Set up polar coordinates
+    inicenter = SVector(0.0, 0.0, 0.0)
+    x_norm = x[1] - inicenter[1]
+    y_norm = x[2] - inicenter[2]
+    z_norm = x[3] - inicenter[3]
+    r = sqrt(x_norm^2 + y_norm^2 + z_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 / (4 * pi * r0^2)
+    p0_outer = 1.0e-1 # "simpler" setup since positivity limiter for pressure is not yet supported in 3D
+
+    # Calculate primitive variables
+    rho = 1.0
+    v1 = 0.0
+    v2 = 0.0
+    v3 = 0.0
+    p = r > r0 ? p0_outer : p0_inner
+
+    return prim2cons(SVector(rho, v1, v2, v3, p), equations)
+end
+
+initial_condition = initial_condition_sedov_blast_wave
+
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho"])
+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, -1.0)
+coordinates_max = (1.0, 1.0, 1.0)
+
+trees_per_dimension = (8, 8, 8)
+mesh = P4estMesh(trees_per_dimension,
+                 polydeg = 1, initial_refinement_level = 0,
+                 coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+# create the semi discretization object
+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 = 100
+analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
+
+alive_callback = AliveCallback(analysis_interval = analysis_interval)
+
+save_solution = SaveSolutionCallback(interval = 10,
+                                     save_initial_solution = true,
+                                     save_final_solution = true,
+                                     extra_node_variables = (:limiting_coefficient,))
+
+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
+                  callback = callbacks);

src/callbacks_stage/subcell_limiter_idp_correction.jl

@@ -64,4 +64,5 @@ init_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
 finalize_callback(limiter!::SubcellLimiterIDPCorrection, semi) = nothing
 
 include("subcell_limiter_idp_correction_2d.jl")
+include("subcell_limiter_idp_correction_3d.jl")
 end # @muladd

src/callbacks_stage/subcell_limiter_idp_correction_3d.jl

@@ -0,0 +1,87 @@
+# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
+# Since these FMAs can increase the performance of many numerical algorithms,
+# we need to opt-in explicitly.
+# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
+@muladd begin
+#! format: noindent
+
+function perform_idp_correction!(u, dt,
+                                 mesh::P4estMesh{3},
+                                 equations, dg, cache)
+    @unpack inverse_weights = dg.basis # Plays role of DG subcell sizes
+    @unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes
+    @unpack alpha = dg.volume_integral.limiter.cache.subcell_limiter_coefficients
+
+    @threaded for element in eachelement(dg, cache)
+        for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+            # Sign switch as in apply_jacobian!
+            inverse_jacobian = -get_inverse_jacobian(cache.elements.inverse_jacobian,
+                                                     mesh, i, j, k, element)
+
+            # Note: For LGL nodes, the high-order and low-order fluxes at element interfaces are equal.
+            # Therefore, the antidiffusive fluxes are zero.
+            # To avoid accessing zero entries, we directly use zero vectors instead.
+            if i > 1 # Not at "left" boundary node
+                alpha1 = max(alpha[i - 1, j, k, element], alpha[i, j, k, element])
+                alpha_flux1 = (1 - alpha1) *
+                              get_node_vars(antidiffusive_flux1_R, equations, dg,
+                                            i, j, k, element)
+            else # At "left" boundary node
+                alpha_flux1 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+            if i < nnodes(dg) # Not at "right" boundary node
+                alpha1_ip1 = max(alpha[i, j, k, element], alpha[i + 1, j, k, element])
+                alpha_flux1_ip1 = (1 - alpha1_ip1) *
+                                  get_node_vars(antidiffusive_flux1_L, equations, dg,
+                                                i + 1, j, k, element)
+            else # At "right" boundary node
+                alpha_flux1_ip1 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+            if j > 1 # Not at "bottom" boundary node
+                alpha2 = max(alpha[i, j - 1, k, element], alpha[i, j, k, element])
+                alpha_flux2 = (1 - alpha2) *
+                              get_node_vars(antidiffusive_flux2_R, equations, dg,
+                                            i, j, k, element)
+            else # At "bottom" boundary node
+                alpha_flux2 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+            if j < nnodes(dg) # Not at "top" boundary node
+                alpha2_jp1 = max(alpha[i, j, k, element], alpha[i, j + 1, k, element])
+                alpha_flux2_jp1 = (1 - alpha2_jp1) *
+                                  get_node_vars(antidiffusive_flux2_L, equations, dg,
+                                                i, j + 1, k, element)
+            else # At "top" boundary node
+                alpha_flux2_jp1 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+            if k > 1 # Not at "front" boundary node
+                alpha3 = max(alpha[i, j, k - 1, element], alpha[i, j, k, element])
+                alpha_flux3 = (1 - alpha3) *
+                              get_node_vars(antidiffusive_flux3_R, equations, dg,
+                                            i, j, k, element)
+            else # At "front" boundary node
+                alpha_flux3 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+            if k < nnodes(dg) # Not at "back" boundary node
+                alpha3_kp1 = max(alpha[i, j, k, element], alpha[i, j, k + 1, element])
+                alpha_flux3_kp1 = (1 - alpha3_kp1) *
+                                  get_node_vars(antidiffusive_flux3_L, equations, dg,
+                                                i, j, k + 1, element)
+            else # At "back" boundary node
+                alpha_flux3_kp1 = zero(SVector{nvariables(equations), eltype(u)})
+            end
+
+            for v in eachvariable(equations)
+                u[v, i, j, k, element] += dt * inverse_jacobian *
+                                          (inverse_weights[i] *
+                                           (alpha_flux1_ip1[v] - alpha_flux1[v]) +
+                                           inverse_weights[j] *
+                                           (alpha_flux2_jp1[v] - alpha_flux2[v]) +
+                                           inverse_weights[k] *
+                                           (alpha_flux3_kp1[v] - alpha_flux3[v]))
+            end
+        end
+    end
+
+    return nothing
+end
+end # @muladd

src/solvers/dgsem_p4est/dg.jl

@@ -91,4 +91,6 @@ include("dg_3d_parabolic.jl")
 include("dg_parallel.jl")
 
 include("subcell_limiters_2d.jl")
+include("subcell_limiters_3d.jl")
+include("dg_3d_subcell_limiters.jl")
 end # @muladd