WIP: Implement subcell limiting for non-conservative systems with the DGSEM structured solver

src/equations/ideal_glm_mhd_2d.jl

@@ -472,6 +472,143 @@ This function is used to compute the subcell fluxes in dg_2d_subcell_limiters.jl
     return f
 end
 
+@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
+                                                             normal_direction_ll::AbstractVector,
+                                                             normal_direction_average::AbstractVector,
+                                                             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
+
+    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
+    # at the same node location) while `B_dot_n_rr` uses the averaged normal
+    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
+    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
+    # the divergence of B.
+    B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+    B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+    psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+    v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
+    B_dot_n_avg = B1_avg * normal_direction_average[1] +
+                  B2_avg * normal_direction_average[2]
+
+    # 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})
+    f = SVector(0,
+                B1_ll * B_dot_n_avg,
+                B2_ll * B_dot_n_avg,
+                B3_ll * B_dot_n_avg,
+                v_dot_B_ll * B_dot_n_avg + v_dot_n_ll * psi_ll * psi_avg,
+                v1_ll * B_dot_n_avg,
+                v2_ll * B_dot_n_avg,
+                v3_ll * B_dot_n_avg,
+                v_dot_n_ll * psi_avg)
+
+    return f
+end
+
+@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
+                                                             normal_direction_ll::AbstractVector,
+                                                             equations::IdealGlmMhdEquations2D,
+                                                             nonconservative_type::NonConservativeLocal,
+                                                             nonconservative_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
+
+    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
+    # at the same node location) while `B_dot_n_rr` uses the averaged normal
+    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
+    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
+    # the divergence of B.
+    if nonconservative_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 #nonconservative_term == 2
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        v_dot_n_ll = v1_ll * normal_direction_ll[1] + v2_ll * normal_direction_ll[2]
+        f = SVector(0,
+                    0,
+                    0,
+                    0,
+                    v_dot_n_ll * psi_ll,
+                    0,
+                    0,
+                    0,
+                    v_dot_n_ll)
+    end
+
+    return f
+end
+
+@inline function flux_nonconservative_powell_local_symmetric(u_ll, u_rr,
+                                                             normal_direction_ll::AbstractVector,
+                                                             normal_direction_average::AbstractVector,
+                                                             equations::IdealGlmMhdEquations2D,
+                                                             nonconservative_type::NonConservativeSymmetric,
+                                                             nonconservative_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
+
+    # Note that `v_dot_n_ll` uses the `normal_direction_ll` (contravariant vector
+    # at the same node location) while `B_dot_n_rr` uses the averaged normal
+    # direction. The reason for this is that `v_dot_n_ll` depends only on the left
+    # state and multiplies some gradient while `B_dot_n_rr` is used to compute
+    # the divergence of B.
+    if nonconservative_term == 1
+        # Powell nonconservative term:   (0, B_1, B_2, B_3, v⋅B, v_1, v_2, v_3, 0)
+        B1_avg = (B1_ll + B1_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+        B2_avg = (B2_ll + B2_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+        B_dot_n_avg = B1_avg * normal_direction_average[1] +
+                      B2_avg * normal_direction_average[2]
+
+        f = SVector(0,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    B_dot_n_avg,
+                    0)
+    else #nonconservative_term == 2
+        psi_avg = (psi_ll + psi_rr) #* 0.5 # The flux is already multiplied by 0.5 wherever it is used in the code
+        # Galilean nonconservative term: (0, 0, 0, 0, ψ v_{1,2}, 0, 0, 0, v_{1,2})
+        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)
 

src/solvers/dgsem_structured/dg_2d_subcell_limiters.jl

@@ -108,4 +108,219 @@
 
     return nothing
 end
+
+# Calculate the DG staggered volume fluxes `fhat` in subcell FV-form inside the element
+# (**with non-conservative terms**).
+#
+# See also `flux_differencing_kernel!`.
+#
+# The calculation of the non-conservative staggered "fluxes" requires non-conservative
+# terms that can be written as a product of local and a symmetric contributions. 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 calcflux_fhat!(fhat1_L, fhat1_R, fhat2_L, fhat2_R, u,
+                                mesh::StructuredMesh{2},
+                                nonconservative_terms::True, equations,
+                                volume_flux, dg::DGSEM, element, cache)
+    (; contravariant_vectors) = cache.elements
+    (; weights, derivative_split) = dg.basis
+    (; flux_temp_threaded, flux_nonconservative_temp_threaded) = cache
+    (; fhat_temp_threaded, fhat_nonconservative_temp_threaded, phi_threaded) = cache
+
+    volume_flux_cons, volume_flux_noncons = volume_flux
+
+    flux_temp = flux_temp_threaded[Threads.threadid()]
+    flux_noncons_temp = flux_nonconservative_temp_threaded[Threads.threadid()]
+
+    fhat_temp = fhat_temp_threaded[Threads.threadid()]
+    fhat_noncons_temp = fhat_nonconservative_temp_threaded[Threads.threadid()]
+    phi = phi_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))
+    flux_noncons_temp .= zero(eltype(flux_noncons_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.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)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[ii, i], fluxtilde1,
+                                       equations, dg, ii, j)
+            for noncons in 1:n_nonconservative_terms(equations)
+                # We multiply by 0.5 because that is done in other parts of Trixi
+                fluxtilde1_noncons = volume_flux_noncons(u_node, u_node_ii, Ja1_node,
+                                                         Ja1_avg, equations,
+                                                         NonConservativeSymmetric(),
+                                                         noncons)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[i, ii],
+                                           fluxtilde1_noncons,
+                                           equations, dg, noncons, i, j)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[ii, i],
+                                           fluxtilde1_noncons,
+                                           equations, dg, noncons, ii, j)
+            end
+        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))
+
+    fhat_temp[:, 1, :] .= zero(eltype(fhat1_L))
+    fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L))
+
+    # 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)
+        Ja1_node = get_contravariant_vector(1, contravariant_vectors, i, j, element) # x direction
+        for noncons in 1:n_nonconservative_terms(equations)
+            set_node_vars!(phi,
+                           volume_flux_noncons(u_local, Ja1_node, equations,
+                                               NonConservativeLocal(), noncons),
+                           equations, dg, noncons, i, j)
+        end
+    end
+
+    for j in eachnode(dg), i in 1:(nnodes(dg) - 1)
+        # Conservative part
+        for v in eachvariable(equations)
+            value = fhat_temp[v, i, j] + weights[i] * flux_temp[v, i, j]
+            fhat_temp[v, i + 1, j] = value
+            fhat1_L[v, i + 1, j] = value
+            fhat1_R[v, i + 1, j] = value
+        end
+        # Nonconservative part
+        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
+
+            fhat1_L[v, i + 1, j] = fhat1_L[v, i + 1, j] + phi[v, noncons, i, j] * value
+            fhat1_R[v, i + 1, j] = fhat1_R[v, i + 1, j] +
+                                   phi[v, noncons, i + 1, j] * value
+        end
+    end
+
+    # Split form volume flux in orientation 2: y direction
+    flux_temp .= zero(eltype(flux_temp))
+    flux_noncons_temp .= zero(eltype(flux_noncons_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.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)
+            multiply_add_to_node_vars!(flux_temp, derivative_split[jj, j], fluxtilde2,
+                                       equations, dg, i, jj)
+            for noncons in 1:n_nonconservative_terms(equations)
+                # We multiply by 0.5 because that is done in other parts of Trixi
+                fluxtilde2_noncons = volume_flux_noncons(u_node, u_node_jj, Ja2_node,
+                                                         Ja2_avg, equations,
+                                                         NonConservativeSymmetric(),
+                                                         noncons)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[j, jj],
+                                           fluxtilde2_noncons,
+                                           equations, dg, noncons, i, j)
+                multiply_add_to_node_vars!(flux_noncons_temp,
+                                           0.5f0 * derivative_split[jj, j],
+                                           fluxtilde2_noncons,
+                                           equations, dg, noncons, i, jj)
+            end
+        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))
+
+    fhat_temp[:, 1, :] .= zero(eltype(fhat1_L))
+    fhat_noncons_temp[:, :, 1, :] .= zero(eltype(fhat1_L))
+
+    # 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)
+        Ja2_node = get_contravariant_vector(2, contravariant_vectors, i, j, element)
+        for noncons in 1:n_nonconservative_terms(equations)
+            set_node_vars!(phi,
+                           volume_flux_noncons(u_local, Ja2_node, equations,
+                                               NonConservativeLocal(), noncons),
+                           equations, dg, noncons, i, j)
+        end
+    end
+
+    for j in 1:(nnodes(dg) - 1), i in eachnode(dg)
+        # Conservative part
+        for v in eachvariable(equations)
+            value = fhat_temp[v, i, j] + weights[j] * flux_temp[v, i, j]
+            fhat_temp[v, i, j + 1] = value
+            fhat1_L[v, i, j + 1] = value
+            fhat1_R[v, i, j + 1] = value
+        end
+        # Nonconservative part
+        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
+
+            fhat1_L[v, i, j + 1] = fhat1_L[v, i, j + 1] + phi[v, noncons, i, j] * value
+            fhat1_R[v, i, j + 1] = fhat1_R[v, i, j + 1] +
+                                   phi[v, noncons, i, j + 1] * value
+        end
+    end
+
+    return nothing
+end
 end # @muladd