Add subcell limiting support for StructuredMesh

examples/structured_2d_dgsem/elixir_euler_free_stream_sc_subcell.jl

@@ -0,0 +1,85 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_constant
+
+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)])
+
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Mapping as described in https://arxiv.org/abs/2012.12040 but reduced to 2D.
+# This particular mesh is unstructured in the yz-plane, but extruded in x-direction.
+# Apply the warping mapping in the yz-plane to get a curved 2D mesh that is extruded
+# in x-direction to ensure free stream preservation on a non-conforming mesh.
+# See https://doi.org/10.1007/s10915-018-00897-9, Section 6.
+
+# Mapping as described in https://arxiv.org/abs/2012.12040, but reduced to 2D
+function mapping(xi_, eta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+
+    y = eta + 3 / 8 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+                       cos(0.5 * pi * (2 * eta - 3) / 3))
+
+    x = xi + 3 / 8 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+                      cos(2 * pi * (2 * y - 3) / 3))
+
+    return SVector(x, y)
+end
+
+cells_per_dimension = (16, 16)
+mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.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,
+                        stepsize_callback,
+                        save_solution)
+
+###############################################################################
+# 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

examples/structured_2d_dgsem/elixir_euler_source_terms_sc_subcell.jl

@@ -0,0 +1,79 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible Euler equations
+
+equations = CompressibleEulerEquations2D(1.4)
+
+initial_condition = initial_condition_convergence_test
+source_terms = source_terms_convergence_test
+
+boundary_condition = BoundaryConditionDirichlet(initial_condition)
+boundary_conditions = (x_neg = boundary_condition,
+                       x_pos = boundary_condition,
+                       y_neg = boundary_condition,
+                       y_pos = boundary_condition)
+
+# Get the DG approximation space
+surface_flux = flux_lax_friedrichs
+volume_flux = flux_ranocha
+polydeg = 3
+basis = LobattoLegendreBasis(polydeg)
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons = ["rho"],
+                                positivity_variables_nonlinear = [pressure])
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg = volume_flux,
+                                                volume_flux_fv = surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+# Waving flag
+f1(s) = SVector(-1.0, s - 1.0)
+f2(s) = SVector(1.0, s + 1.0)
+f3(s) = SVector(s, -1.0 + sin(0.5 * pi * s))
+f4(s) = SVector(s, 1.0 + sin(0.5 * pi * s))
+mapping = Trixi.transfinite_mapping((f1, f2, f3, f4))
+
+cells_per_dimension = (16, 16)
+mesh = StructuredMesh(cells_per_dimension, mapping, periodicity = false)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    boundary_conditions = boundary_conditions,
+                                    source_terms = source_terms)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 2.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.8)
+
+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_limiter_idp_correction_2d.jl

@@ -6,39 +6,60 @@
 #! format: noindent
 
 function perform_idp_correction!(u, dt, mesh::TreeMesh2D, 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
-
     @threaded for element in eachelement(dg, cache)
         # Sign switch as in apply_jacobian!
         inverse_jacobian = -cache.elements.inverse_jacobian[element]
 
         for j in eachnode(dg), i in eachnode(dg)
-            # 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
+            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::StructuredMesh{2}, equations, dg, cache)
+    @threaded for element in eachelement(dg, cache)
+        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)
+        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_structured/dg.jl

@@ -85,6 +85,9 @@ include("indicators_1d.jl")
 include("indicators_2d.jl")
 include("indicators_3d.jl")
 
+include("subcell_limiters_2d.jl")
+include("dg_2d_subcell_limiters.jl")
+
 # Specialized implementations used to improve performance
 include("dg_2d_compressible_euler.jl")
 include("dg_3d_compressible_euler.jl")

src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl

@@ -0,0 +1,111 @@
+# 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, u,
+                                mesh::StructuredMesh{2},
+                                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 j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, 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, element)
+            # pull the contravariant vectors and compute the average
+            Ja1_node_ii = get_contravariant_vector(1, contravariant_vectors, ii, j,
+                                                   element)
+            Ja1_avg = 0.5 * (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)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,
+                                       equations, dg, ii, j)
+        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 j in eachnode(dg), i in 1:(nnodes(dg) - 1), v in eachvariable(equations)
+        fhat1_L[v, i + 1, j] = fhat1_L[v, i, j] + weights[i] * flux_temp[v, i, j]
+        fhat1_R[v, i + 1, j] = fhat1_L[v, i + 1, j]
+    end
+
+    # Split form volume flux in orientation 2: y direction
+    flux_temp .= zero(eltype(flux_temp))
+
+    for j in eachnode(dg), i in eachnode(dg)
+        u_node = get_node_vars(u, equations, dg, i, j, element)
+
+        # pull the contravariant vectors in each coordinate direction
+        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+
+        # y direction
+        for jj in (j + 1):nnodes(dg)
+            u_node_jj = get_node_vars(u, equations, dg, i, jj, element)
+            # pull the contravariant vectors and compute the average
+            Ja2_node_jj = get_contravariant_vector(2, contravariant_vectors, i, jj,
+                                                   element)
+            Ja2_avg = 0.5 * (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)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,
+                                       equations, dg, i, jj)
+        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 j in 1:(nnodes(dg) - 1), i in eachnode(dg), v in eachvariable(equations)
+        fhat2_L[v, i, j + 1] = fhat2_L[v, i, j] + weights[j] * flux_temp[v, i, j]
+        fhat2_R[v, i, j + 1] = fhat2_L[v, i, j + 1]
+    end
+
+    return nothing
+end
+end # @muladd