First PR for 3D subcell limiting

examples/p4est_3d_dgsem/elixir_euler_sedov_sc_subcell.jl

@@ -0,0 +1,99 @@
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations3D(1.4)
+
+"""
+    initial_condition_medium_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
+                  ode_default_options()..., 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,55 @@
+# 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
+    @unpack antidiffusive_flux1_L, antidiffusive_flux1_R, antidiffusive_flux2_L, antidiffusive_flux2_R, antidiffusive_flux3_L, antidiffusive_flux3_R = cache.antidiffusive_fluxes
+    @unpack alpha1, alpha2, alpha3 = 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: antidiffusive_flux1[v, i, xi, eta, element] = antidiffusive_flux2[v, xi, i, eta, element] = antidiffusive_flux3[v, xi, eta, i, element] = 0 for all i in 1:nnodes and xi, eta in {1, nnodes+1}
+            alpha_flux1 = (1 - alpha1[i, j, k, element]) *
+                          get_node_vars(antidiffusive_flux1_R, equations, dg,
+                                        i, j, k, element)
+            alpha_flux1_ip1 = (1 - alpha1[i + 1, j, k, element]) *
+                              get_node_vars(antidiffusive_flux1_L, equations, dg,
+                                            i + 1, j, k, element)
+            alpha_flux2 = (1 - alpha2[i, j, k, element]) *
+                          get_node_vars(antidiffusive_flux2_R, equations, dg,
+                                        i, j, k, element)
+            alpha_flux2_jp1 = (1 - alpha2[i, j + 1, k, element]) *
+                              get_node_vars(antidiffusive_flux2_L, equations, dg,
+                                            i, j + 1, k, element)
+            alpha_flux3 = (1 - alpha3[i, j, k, element]) *
+                          get_node_vars(antidiffusive_flux3_R, equations, dg,
+                                        i, j, k, element)
+            alpha_flux3_jp1 = (1 - alpha3[i, j, k + 1, element]) *
+                              get_node_vars(antidiffusive_flux3_L, equations, dg,
+                                            i, j, k + 1, element)
+
+            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_jp1[v] - alpha_flux3[v]))
+            end
+        end
+    end
+
+    return nothing
+end
+end # @muladd

src/solvers/dgsem_structured/dg.jl

@@ -156,6 +156,7 @@ include("indicators_3d.jl")
 
 include("subcell_limiters_2d.jl")
 include("dg_2d_subcell_limiters.jl")
+include("dg_3d_subcell_limiters.jl")
 
 # Specialized implementations used to improve performance
 include("dg_2d_compressible_euler.jl")

src/solvers/dgsem_structured/dg_3d_subcell_limiters.jl

@@ -0,0 +1,156 @@
+# 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
+
+# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element
+# (**without non-conservative terms**).
+#
+# See also `flux_differencing_kernel!`.
+@inline function calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, fhat3_L, fhat3_R, u,
+                                mesh::P4estMesh{3},
+                                nonconservative_terms::False, equations,
+                                volume_flux, dg::DGSEM, element, cache)
+    (; contravariant_vectors) = cache.elements
+    (; weights, derivative_split) = dg.basis
+    (; flux_temp_threaded) = cache
+
+    flux_temp = flux_temp_threaded[Threads.threadid()]
+
+    # The FV-form fluxes are calculated in a recursive manner, i.e.:
+    # fhat_(0,1)   = w_0 * FVol_0,
+    # fhat_(j,j+1) = fhat_(j-1,j) + w_j * FVol_j,   for j=1,...,N-1,
+    # with the split form volume fluxes FVol_j = -2 * sum_i=0^N D_ji f*_(j,i).
+
+    # To use the symmetry of the `volume_flux`, the split form volume flux is precalculated
+    # like in `calc_volume_integral!` for the `VolumeIntegralFluxDifferencing`
+    # and saved in in `flux_temp`.
+
+    # Split form volume flux in orientation 1: x direction
+    flux_temp .= zero(eltype(flux_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, k, element) # x direction
+
+        # All diagonal entries of `derivative_split` are zero. Thus, we can skip
+        # the computation of the diagonal terms. In addition, we use the symmetry
+        # of the `volume_flux` to save half of the possible two-point flux
+        # computations.
+
+        # x direction
+        for ii in (i + 1):nnodes(dg)
+            u_node_ii = get_node_vars(u, equations, dg, ii, j, k, element)
+            # pull the contravariant vectors and compute the average
+            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j, k,
+                                                   element)
+            Ja1_avg = 0.5f0 * (Ja1_node + Ja1_node_ii)
+
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde1 = volume_flux(u_node, u_node_ii, Ja1_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[i, ii], fluxtilde1,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,
+                                       equations, dg, ii, j, k)
+        end
+    end
+
+    # FV-form flux `fhat` in x direction
+    fhat1_L[:, 1, :, :] .= zero(eltype(fhat1_L))
+    fhat1_L[:, nnodes(dg) + 1, :, :] .= zero(eltype(fhat1_L))
+    fhat1_R[:, 1, :, :] .= zero(eltype(fhat1_R))
+    fhat1_R[:, nnodes(dg) + 1, :, :] .= zero(eltype(fhat1_R))
+
+    for k in eachnode(dg), j in eachnode(dg), i in 1:(nnodes(dg) - 1),
+        v in eachvariable(equations)
+
+        fhat1_L[v, i + 1, j, k] = fhat1_L[v, i, j, k] +
+                                  weights[i] * flux_temp[v, i, j, k]
+        fhat1_R[v, i + 1, j, k] = fhat1_L[v, i + 1, j, k]
+    end
+
+    # Split form volume flux in orientation 2: y direction
+    flux_temp .= zero(eltype(flux_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, k, element)
+
+        # y direction
+        for jj in (j + 1):nnodes(dg)
+            u_node_jj = get_node_vars(u, equations, dg, i, jj, k, element)
+            # pull the contravariant vectors and compute the average
+            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj, k,
+                                                   element)
+            Ja2_avg = 0.5f0 * (Ja2_node + Ja2_node_jj)
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde2 = volume_flux(u_node, u_node_jj, Ja2_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[j, jj], fluxtilde2,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,
+                                       equations, dg, i, jj, k)
+        end
+    end
+
+    # FV-form flux `fhat` in y direction
+    fhat2_L[:, :, 1, :] .= zero(eltype(fhat2_L))
+    fhat2_L[:, :, nnodes(dg) + 1, :] .= zero(eltype(fhat2_L))
+    fhat2_R[:, :, 1, :] .= zero(eltype(fhat2_R))
+    fhat2_R[:, :, nnodes(dg) + 1, :] .= zero(eltype(fhat2_R))
+
+    for k in eachnode(dg), j in 1:(nnodes(dg) - 1), i in eachnode(dg),
+        v in eachvariable(equations)
+
+        fhat2_L[v, i, j + 1, k] = fhat2_L[v, i, j, k] +
+                                  weights[j] * flux_temp[v, i, j, k]
+        fhat2_R[v, i, j + 1, k] = fhat2_L[v, i, j + 1, k]
+    end
+
+    # Split form volume flux in orientation 3: z direction
+    flux_temp .= zero(eltype(flux_temp))
+
+    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, k, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja3_node = get_contravariant_vector(3, contravariant_vectors, i, j, k, element)
+
+        # y direction
+        for kk in (k + 1):nnodes(dg)
+            u_node_kk = get_node_vars(u, equations, dg, i, j, kk, element)
+            # pull the contravariant vectors and compute the average
+            Ja3_node_kk = get_contravariant_vector(3, contravariant_vectors, i, j, kk,
+                                                   element)
+            Ja3_avg = 0.5f0 * (Ja3_node + Ja3_node_kk)
+            # compute the contravariant sharp flux in the direction of the averaged contravariant vector
+            fluxtilde3 = volume_flux(u_node, u_node_kk, Ja3_avg, equations)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[k, kk], fluxtilde3,
+                                       equations, dg, i, j, k)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[kk, k], fluxtilde3,
+                                       equations, dg, i, j, kk)
+        end
+    end
+
+    # FV-form flux `fhat` in y direction
+    fhat3_L[:, :, :, 1] .= zero(eltype(fhat3_L))
+    fhat3_L[:, :, :, nnodes(dg) + 1] .= zero(eltype(fhat3_L))
+    fhat3_R[:, :, :, 1] .= zero(eltype(fhat3_R))
+    fhat3_R[:, :, :, nnodes(dg) + 1] .= zero(eltype(fhat3_R))
+
+    for k in 1:(nnodes(dg) - 1), j in eachnode(dg), i in eachnode(dg),
+        v in eachvariable(equations)
+
+        fhat3_L[v, i, j, k + 1] = fhat3_L[v, i, j, k] +
+                                  weights[k] * flux_temp[v, i, j, k]
+        fhat3_R[v, i, j, k + 1] = fhat3_L[v, i, j, k + 1]
+    end
+
+    return nothing
+end
+end # @muladd