WIP: Implement subcell limiting for non-conservative systems

examples/tree_2d_dgsem/elixir_mhd_shockcapturing_subcell.jl

@@ -0,0 +1,108 @@
+
+using OrdinaryDiffEq
+using Trixi
+
+###############################################################################
+# semidiscretization of the compressible ideal GLM-MHD equations
+
+equations = IdealGlmMhdEquations2D(1.4)
+
+"""
+    initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
+
+An MHD blast wave modified from:
+- Dominik Derigs, Gregor J. Gassner, Stefanie Walch & Andrew R. Winters (2018)
+  Entropy Stable Finite Volume Approximations for Ideal Magnetohydrodynamics
+  [doi: 10.1365/s13291-018-0178-9](https://doi.org/10.1365/s13291-018-0178-9)
+This setup needs a positivity limiter for the density.
+"""
+function initial_condition_blast_wave(x, t, equations::IdealGlmMhdEquations2D)
+  # setup taken from Derigs et al. DMV article (2018)
+  # domain must be [-0.5, 0.5] x [-0.5, 0.5], γ = 1.4
+  r = sqrt(x[1]^2 + x[2]^2)
+
+  pmax = 10.0
+  pmin = 1.0
+  rhomax = 1.0
+  rhomin = 0.01
+  if r <= 0.09
+    p = pmax
+    rho = rhomax
+  elseif r >= 0.1
+    p = pmin
+    rho = rhomin
+  else
+    p = pmin + (0.1 - r) * (pmax - pmin) / 0.01
+    rho = rhomin + (0.1 - r) * (rhomax - rhomin) / 0.01
+  end
+  v1 = 0.0
+  v2 = 0.0
+  v3 = 0.0
+  B1 = 1.0/sqrt(4.0*pi)
+  B2 = 0.0
+  B3 = 0.0
+  psi = 0.0
+  return prim2cons(SVector(rho, v1, v2, v3, p, B1, B2, B3, psi), equations)
+end
+initial_condition = initial_condition_blast_wave
+
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_powell2)
+volume_flux  = (flux_derigs_etal, flux_nonconservative_powell2)
+basis = LobattoLegendreBasis(3)
+
+limiter_idp = SubcellLimiterIDP(equations, basis;
+                                positivity_variables_cons=[1],
+                                positivity_correction_factor=0.5)
+volume_integral = VolumeIntegralSubcellLimiting(limiter_idp;
+                                                volume_flux_dg=volume_flux,
+                                                volume_flux_fv=surface_flux)
+solver = DGSEM(basis, surface_flux, volume_integral)
+
+coordinates_min = (-0.5, -0.5)
+coordinates_max = ( 0.5,  0.5)
+mesh = TreeMesh(coordinates_min, coordinates_max,
+                initial_refinement_level=4,
+                n_cells_max=10_000)
+
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)
+
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 0.1)
+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)
+
+cfl = 0.5
+stepsize_callback = StepsizeCallback(cfl=cfl)
+
+glm_speed_callback = GlmSpeedCallback(glm_scale=0.5, cfl=cfl)
+
+callbacks = CallbackSet(summary_callback,
+                        analysis_callback,
+                        alive_callback,
+                        save_solution,
+                        stepsize_callback,
+                        glm_speed_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
+                  save_everystep=false, callback=callbacks);
+summary_callback() # print the timer summary

src/Trixi.jl

@@ -60,7 +60,7 @@ using RecipesBase: RecipesBase
 using Requires: @require
 using Static: Static, One, True, False
 @reexport using StaticArrays: SVector
-using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix
+using StaticArrays: StaticArrays, MVector, MArray, SMatrix, @SMatrix, MMatrix
 using StrideArrays: PtrArray, StrideArray, StaticInt
 @reexport using StructArrays: StructArrays, StructArray
 using TimerOutputs: TimerOutputs, @notimeit, TimerOutput, print_timer, reset_timer!
@@ -161,7 +161,7 @@ export GradientVariablesPrimitive, GradientVariablesEntropy
 export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle,
        flux_godunov,
        flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner,
-       flux_nonconservative_powell,
+       flux_nonconservative_powell, flux_nonconservative_powell2,
        flux_kennedy_gruber, flux_shima_etal, flux_ec,
        flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_es_fjordholm_etal,
        flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal,

src/callbacks_stage/subcell_limiter_idp_correction_2d.jl

@@ -7,7 +7,7 @@
 
 function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
     @unpack inverse_weights = dg.basis
-    @unpack antidiffusive_flux1, antidiffusive_flux2 = cache.antidiffusive_fluxes
+    @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)
@@ -17,16 +17,16 @@ function perform_idp_correction!(u, dt, mesh::TreeMesh2D, equations, dg, cache)
         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, equations, dg, i, j,
+                          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, equations, dg, i + 1,
+                              get_node_vars(antidiffusive_flux1_L, equations, dg, i + 1,
                                             j, element)
             alpha_flux2 = (1 - alpha2[i, j, element]) *
-                          get_node_vars(antidiffusive_flux2, equations, dg, i, j,
+                          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, equations, dg, i,
+                              get_node_vars(antidiffusive_flux2_L, equations, dg, i,
                                             j + 1, element)
 
             for v in eachvariable(equations)

src/equations/equations.jl

@@ -207,7 +207,22 @@ where `x` specifies the coordinates, `t` is the current time, and `equation` is
 struct BoundaryConditionNeumann{B}
     boundary_normal_flux_function::B
 end
+"""
+    NonConservativeLocal()
+
+Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". 
+When the argument `nonconservative_type` is of type `NonConservativeLocal`, 
+the function returns the local part of the non-conservative term.
+"""
+struct NonConservativeLocal end
+"""
+    NonConservativeSymmetric()
 
+Struct used for multiple dispatch on non-conservative flux functions in the format of "local * symmetric". 
+When the argument `nonconservative_type` is of type `NonConservativeSymmetric`, 
+the function returns the symmetric part of the non-conservative term.
+"""
+struct NonConservativeSymmetric end
 # set sensible default values that may be overwritten by specific equations
 """
     have_nonconservative_terms(equations)
@@ -220,6 +235,12 @@ example of equations with nonconservative terms.
 The return value will be `True()` or `False()` to allow dispatching on the return type.
 """
 have_nonconservative_terms(::AbstractEquations) = False()
+"""
+    nnoncons(equations)
+Number of nonconservative terms for a particular equation. The default is 0 and 
+it must be defined for each nonconservative equation independently.
+"""
+nnoncons(::AbstractEquations) = 0
 have_constant_speed(::AbstractEquations) = False()
 
 default_analysis_errors(::AbstractEquations) = (:l2_error, :linf_error)

src/equations/ideal_glm_mhd_2d.jl

@@ -29,6 +29,8 @@ function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN)
 end
 
 have_nonconservative_terms(::IdealGlmMhdEquations2D) = True()
+nnoncons(::IdealGlmMhdEquations2D) = 2
+
 function varnames(::typeof(cons2cons), ::IdealGlmMhdEquations2D)
     ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3", "psi")
 end
@@ -279,6 +281,185 @@ end
     return f
 end
 
+"""
+    flux_nonconservative_powell2(u_ll, u_rr, orientation::Integer,
+                                equations::IdealGlmMhdEquations2D)
+
+Non-symmetric two-point flux discretizing the nonconservative (source) term of
+Powell and the Galilean nonconservative term associated with the GLM multiplier
+of the [`IdealGlmMhdEquations2D`](@ref).
+
+This implementation uses a non-conservative term that can be written as the product
+of local and symmetric parts. It is equivalent to the non-conservative flux of Bohm 
+et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different 
+results on non-conforming meshes(!).
+## References
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.
+"""
+@inline function flux_nonconservative_powell2(u_ll, u_rr, orientation::Integer,
+                                              equations::IdealGlmMhdEquations2D)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
+
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+    # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+    # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+    psi_avg = (psi_ll + psi_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+    if orientation == 1
+        B1_avg = (B1_ll + B1_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+        f = SVector(0,
+                    B1_ll * B1_avg,
+                    B2_ll * B1_avg,
+                    B3_ll * B1_avg,
+                    v_dot_B_ll * B1_avg + v1_ll * psi_ll * psi_avg,
+                    v1_ll * B1_avg,
+                    v2_ll * B1_avg,
+                    v3_ll * B1_avg,
+                    v1_ll * psi_avg)
+    else # orientation == 2
+        B2_avg = (B2_ll + B2_rr) #* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+        f = SVector(0,
+                    B1_ll * B2_avg,
+                    B2_ll * B2_avg,
+                    B3_ll * B2_avg,
+                    v_dot_B_ll * B2_avg + v2_ll * psi_ll * psi_avg,
+                    v1_ll * B2_avg,
+                    v2_ll * B2_avg,
+                    v3_ll * B2_avg,
+                    v2_ll * psi_avg)
+    end
+
+    return f
+end
+"""
+    flux_nonconservative_powell2(u_ll, orientation::Integer,
+                                equations::IdealGlmMhdEquations2D, 
+                                nonconservative_type::NonConservativeLocal,
+                                noncons_term::Integer)
+
+Local part of the Powell and GLM non-conservative terms. Needed for the calculation of 
+the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.
+"""
+@inline function flux_nonconservative_powell2(u_ll, orientation::Integer,
+                                              equations::IdealGlmMhdEquations2D,
+                                              nonconservative_type::NonConservativeLocal,
+                                              noncons_term::Integer)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+    if noncons_term == 1
+        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+        f = SVector(0,
+                    B1_ll,
+                    B2_ll,
+                    B3_ll,
+                    v_dot_B_ll,
+                    v1_ll,
+                    v2_ll,
+                    v3_ll,
+                    0)
+    else #noncons_term ==2
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        if orientation == 1
+            f = SVector(0,
+                        0,
+                        0,
+                        0,
+                        v1_ll * psi_ll,
+                        0,
+                        0,
+                        0,
+                        v1_ll)
+        else #orientation == 2
+            f = SVector(0,
+                        0,
+                        0,
+                        0,
+                        v2_ll * psi_ll,
+                        0,
+                        0,
+                        0,
+                        v2_ll)
+        end
+    end
+    return f
+end
+"""
+    flux_nonconservative_powell2(u_ll, orientation::Integer,
+                                equations::IdealGlmMhdEquations2D, 
+                                nonconservative_type::NonConservativeSymmetric,
+                                noncons_term::Integer)
+
+Symmetric part of the Powell and GLM non-conservative terms. Needed for the calculation of 
+the non-conservative staggered "fluxes" for subcell limiting. See, e.g.,
+- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts
+  Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf.
+"""
+@inline function flux_nonconservative_powell2(u_ll, u_rr, orientation::Integer,
+                                              equations::IdealGlmMhdEquations2D,
+                                              nonconservative_type::NonConservativeSymmetric,
+                                              noncons_term::Integer)
+    rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll, B1_ll, B2_ll, B3_ll, psi_ll = u_ll
+    rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr, B1_rr, B2_rr, B3_rr, psi_rr = u_rr
+
+    v1_ll = rho_v1_ll / rho_ll
+    v2_ll = rho_v2_ll / rho_ll
+    v3_ll = rho_v3_ll / rho_ll
+    v_dot_B_ll = v1_ll * B1_ll + v2_ll * B2_ll + v3_ll * B3_ll
+
+    if noncons_term == 1
+        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+        if orientation == 1
+            B1_avg = (B1_ll + B1_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+            f = SVector(0,
+                        B1_avg,
+                        B1_avg,
+                        B1_avg,
+                        B1_avg,
+                        B1_avg,
+                        B1_avg,
+                        B1_avg,
+                        0)
+        else # orientation == 2
+            B2_avg = (B2_ll + B2_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+            f = SVector(0,
+                        B2_avg,
+                        B2_avg,
+                        B2_avg,
+                        B2_avg,
+                        B2_avg,
+                        B2_avg,
+                        B2_avg,
+                        0)
+        end
+    else #noncons_term == 2
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        psi_avg = (psi_ll + psi_rr)#* 0.5 # We remove the 0.5 because the flux is always multiplied by 0.5
+        f = SVector(0,
+                    0,
+                    0,
+                    0,
+                    psi_avg,
+                    0,
+                    0,
+                    0,
+                    psi_avg)
+    end
+
+    return f
+end
+
 """
     flux_derigs_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdEquations2D)