Implement subcell limiting for non-conservative systems

src/equations/equations.jl

@@ -239,11 +239,13 @@ The return value will be `True()` or `False()` to allow dispatching on the retur
 """
 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.
+    n_nonconservative_terms(equations)
+
+Number of nonconservative terms in the form local * symmetric for a particular equation.
+This function needs to be specialized only if equations with nonconservative terms are
+combined with certain solvers (e.g., subcell limiting).
 """
-nnoncons(::AbstractEquations) = 0
+function n_nonconservative_terms(::AbstractEquations) end
 have_constant_speed(::AbstractEquations) = False()
 
 default_analysis_errors(::AbstractEquations) = (:l2_error, :linf_error)

src/equations/ideal_glm_mhd_2d.jl

@@ -29,7 +29,7 @@ function IdealGlmMhdEquations2D(gamma; initial_c_h = convert(typeof(gamma), NaN)
 end
 
 have_nonconservative_terms(::IdealGlmMhdEquations2D) = True()
-nnoncons(::IdealGlmMhdEquations2D) = 2
+n_nonconservative_terms(::IdealGlmMhdEquations2D) = 2
 
 function varnames(::typeof(cons2cons), ::IdealGlmMhdEquations2D)
     ("rho", "rho_v1", "rho_v2", "rho_v3", "rho_e", "B1", "B2", "B3", "psi")
@@ -295,7 +295,7 @@ of local and symmetric parts. It is equivalent to the non-conservative flux of B
 et al. (`flux_nonconservative_powell`) for conforming meshes but it yields different
 results on non-conforming meshes(!).
 
-The two functions below, which share the same name, return yield either the local 
+The two other flux functions with the same name return either the local 
 or symmetric portion of the non-conservative flux based on the type of the 
 nonconservative_type argument, employing multiple dispatch. They are used to
 compute the subcell fluxes in dg_2d_subcell_limiters.jl.
@@ -317,9 +317,9 @@ compute the subcell fluxes in dg_2d_subcell_limiters.jl.
 
     # 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
+    psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
     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
+        B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
         f = SVector(0,
                     B1_ll * B1_avg,
                     B2_ll * B1_avg,
@@ -330,7 +330,7 @@ compute the subcell fluxes in dg_2d_subcell_limiters.jl.
                     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
+        B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
         f = SVector(0,
                     B1_ll * B2_avg,
                     B2_ll * B2_avg,
@@ -430,7 +430,7 @@ This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl
     if nonconservative_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
+            B1_avg = (B1_ll + B1_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
             f = SVector(0,
                         B1_avg,
                         B1_avg,
@@ -441,7 +441,7 @@ This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl
                         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
+            B2_avg = (B2_ll + B2_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
             f = SVector(0,
                         B2_avg,
                         B2_avg,
@@ -454,7 +454,7 @@ This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl
         end
     else #nonconservative_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
+        psi_avg = (psi_ll + psi_rr)#* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
         f = SVector(0,
                     0,
                     0,

src/solvers/dgsem_tree/dg_2d_subcell_limiters.jl

@@ -26,30 +26,33 @@ function create_cache(mesh::TreeMesh{2}, equations,
                                        nnodes(dg) + 1) for _ in 1:Threads.nthreads()]
     flux_temp_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg), nnodes(dg))
                              for _ in 1:Threads.nthreads()]
-    flux_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations),
-                                                 nnoncons(equations),
-                                                 nnodes(dg), nnodes(dg))
-                                             for _ in 1:Threads.nthreads()]
     fhat_temp_threaded = A3d[A3d(undef, nvariables(equations), nnodes(dg),
                                  nnodes(dg))
                              for _ in 1:Threads.nthreads()]
-    fhat_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations),
-                                                 nnoncons(equations),
-                                                 nnodes(dg), nnodes(dg))
-                                             for _ in 1:Threads.nthreads()]
-
-    phi_threaded = A4d[A4d(undef, nvariables(equations),
-                           nnoncons(equations),
-                           nnodes(dg), nnodes(dg))
-                       for _ in 1:Threads.nthreads()]
-
     antidiffusive_fluxes = Trixi.ContainerAntidiffusiveFlux2D{uEltype}(0,
                                                                        nvariables(equations),
                                                                        nnodes(dg))
+
+    if typeof(have_nonconservative_terms(equations)) == True
+        flux_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations),
+                                                     n_nonconservative_terms(equations),
+                                                     nnodes(dg), nnodes(dg))
+                                                 for _ in 1:Threads.nthreads()]
+        fhat_nonconservative_temp_threaded = A4d[A4d(undef, nvariables(equations),
+                                                     n_nonconservative_terms(equations),
+                                                     nnodes(dg), nnodes(dg))
+                                                 for _ in 1:Threads.nthreads()]
+        phi_threaded = A4d[A4d(undef, nvariables(equations),
+                               n_nonconservative_terms(equations),
+                               nnodes(dg), nnodes(dg))
+                           for _ in 1:Threads.nthreads()]
+        cache = (; cache..., flux_nonconservative_temp_threaded,
+                 fhat_nonconservative_temp_threaded, phi_threaded)
+    end
+
     return (; cache..., antidiffusive_fluxes,
             fhat1_L_threaded, fhat2_L_threaded, fhat1_R_threaded, fhat2_R_threaded,
-            flux_temp_threaded, flux_nonconservative_temp_threaded, fhat_temp_threaded,
-            fhat_nonconservative_temp_threaded, phi_threaded)
+            flux_temp_threaded, fhat_temp_threaded)
 end
 
 function calc_volume_integral!(du, u,
@@ -253,7 +256,7 @@ end
                                        equations, dg, i, j)
             multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], flux1,
                                        equations, dg, ii, j)
-            for noncons in 1:nnoncons(equations)
+            for noncons in 1:n_nonconservative_terms(equations)
                 # We multiply by 0.5 because that is done in other parts of Trixi
                 flux1_noncons = volume_flux_noncons(u_node, u_node_ii, 1, equations,
                                                     NonConservativeSymmetric(), noncons)
@@ -281,7 +284,7 @@ end
     # Compute local contribution to non-conservative flux
     for j in eachnode(dg), i in eachnode(dg)
         u_local = get_node_vars(u, equations, dg, i, j, element)
-        for noncons in 1:nnoncons(equations)
+        for noncons in 1:n_nonconservative_terms(equations)
             set_node_vars!(phi,
                            volume_flux_noncons(u_local, 1, equations,
                                                NonConservativeLocal(), noncons),
@@ -298,7 +301,9 @@ end
             fhat1_R[v, i + 1, j] = value
         end
         # Nonconservative part
-        for noncons in 1:nnoncons(equations), v in eachvariable(equations)
+        for noncons in 1:n_nonconservative_terms(equations),
+            v in eachvariable(equations)
+
             value = fhat_noncons_temp[v, noncons, i, j] +
                     weights[i] * flux_noncons_temp[v, noncons, i, j]
             fhat_noncons_temp[v, noncons, i + 1, j] = value
@@ -322,7 +327,7 @@ end
                                        equations, dg, i, j)
             multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], flux2,
                                        equations, dg, i, jj)
-            for noncons in 1:nnoncons(equations)
+            for noncons in 1:n_nonconservative_terms(equations)
                 # We multiply by 0.5 because that is done in other parts of Trixi
                 flux2_noncons = volume_flux_noncons(u_node, u_node_jj, 2, equations,
                                                     NonConservativeSymmetric(), noncons)
@@ -350,7 +355,7 @@ end
     # Compute local contribution to non-conservative flux
     for j in eachnode(dg), i in eachnode(dg)
         u_local = get_node_vars(u, equations, dg, i, j, element)
-        for noncons in 1:nnoncons(equations)
+        for noncons in 1:n_nonconservative_terms(equations)
             set_node_vars!(phi,
                            volume_flux_noncons(u_local, 2, equations,
                                                NonConservativeLocal(), noncons),
@@ -367,7 +372,9 @@ end
             fhat2_R[v, i, j + 1] = value
         end
         # Nonconservative part
-        for noncons in 1:nnoncons(equations), v in eachvariable(equations)
+        for noncons in 1:n_nonconservative_terms(equations),
+            v in eachvariable(equations)
+
             value = fhat_noncons_temp[v, noncons, i, j] +
                     weights[j] * flux_noncons_temp[v, noncons, i, j]
             fhat_noncons_temp[v, noncons, i, j + 1] = value