WIP: Add subcell-limiting-ready 3D fluxes for multi-ion MHD (P4estMesh)

src/equations/ideal_glm_mhd_multiion.jl

@@ -45,6 +45,31 @@ function default_analysis_integrals(::AbstractIdealGlmMhdMultiIonEquations)
     (entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
 end
 
+# For `VolumeIntegralSubcellLimiting` the nonconservative flux is created as a callable struct to 
+# enable dispatch on the type of the nonconservative term (symmetric / jump).
+struct FluxNonConservativeRuedaRamirezEtAl <:
+       FluxNonConservative{NonConservativeSymmetric()}
+end
+
+# We specify 6 non-conservative terms for FluxNonConservativeRuedaRamirezEtAl. This is the number of
+# non-conservative terms in 3D. In 2D, only 5 terms are needed. TODO: Should we create a different struct
+# for the 2D non-conservative term?
+n_nonconservative_terms(::FluxNonConservativeRuedaRamirezEtAl) = 6
+
+const flux_nonconservative_ruedaramirez_etal = FluxNonConservativeRuedaRamirezEtAl()
+
+# State validation for Newton-bisection method of subcell IDP limiting
+@inline function Base.isvalid(u, equations::AbstractIdealGlmMhdMultiIonEquations)
+    p = pressure(u, equations)
+    for k in eachcomponent(equations)
+        u_k = get_component(k, u, equations)
+        if u_k[1] <= 0 || p[k] <= 0
+            return false
+        end
+    end
+    return true
+end
+
 """
     source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations)
 

src/equations/ideal_glm_mhd_multiion_2d.jl

@@ -343,9 +343,9 @@ The term is composed of four individual non-conservative terms:
 3. The "multi-ion" term, which vanishes in the limit of one ion species.
 4. The GLM term, which is needed for Galilean invariance.
 """
-@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::IdealGlmMhdMultiIonEquations2D)
+@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr,
+                                                                     orientation::Integer,
+                                                                     equations::IdealGlmMhdMultiIonEquations2D)
     @unpack charge_to_mass = equations
     # Unpack left and right states to get the magnetic field
     B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)

src/equations/ideal_glm_mhd_multiion_3d.jl

@@ -329,9 +329,9 @@ The term is composed of four individual non-conservative terms:
 3. The "multi-ion" term, which vanishes in the limit of one ion species.
 4. The GLM term, which is needed for Galilean invariance.
 """
-@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
-                                                        orientation::Integer,
-                                                        equations::IdealGlmMhdMultiIonEquations3D)
+@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr,
+                                                                     orientation::Integer,
+                                                                     equations::IdealGlmMhdMultiIonEquations3D)
     @unpack charge_to_mass = equations
     # Unpack left and right states to get the magnetic field
     B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
@@ -514,9 +514,9 @@ The term is composed of four individual non-conservative terms:
     return SVector(f)
 end
 
-@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
-                                                        normal_direction::AbstractVector,
-                                                        equations::IdealGlmMhdMultiIonEquations3D)
+@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr,
+                                                                     normal_direction::AbstractVector,
+                                                                     equations::IdealGlmMhdMultiIonEquations3D)
     @unpack charge_to_mass = equations
     # Unpack left and right states to get the magnetic field
     B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
@@ -628,6 +628,291 @@ end
     return SVector(f)
 end
 
+"""
+    flux_nonconservative_ruedaramirez_etal(u_ll, normal_direction_ll::AbstractVector,
+                                           equations::IdealGlmMhdMultiIonEquations3D,
+                                           nonconservative_type::NonConservativeLocal,
+                                           nonconservative_term::Integer)
+
+Non-symmetric local part of the non-conservative terms for multi-ion MHD. 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.
+This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl.
+
+On curvilinear meshes this formulation applies the local normal direction compared to the averaged one used 
+in [`flux_nonconservative_ruedaramirez_etal`](@ref) when used for volume fluxes. This is done to reduce the 
+number of operations to obtain the subcell limiting fluxes. However, this decision causes this flux to be 
+"slightly" different when used for subcell limiting or for a flux-differencing DG method.
+"""
+@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll,
+                                                                     normal_direction::AbstractVector,
+                                                                     equations::IdealGlmMhdMultiIonEquations3D,
+                                                                     nonconservative_type::NonConservativeLocal,
+                                                                     nonconservative_term::Integer)
+    @unpack charge_to_mass = equations
+    # Unpack left and right states to get the magnetic field
+    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
+    psi_ll = divergence_cleaning_field(u_ll, equations)
+
+    # Compute charge ratio of u_ll
+    charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)})
+    total_electron_charge = zero(eltype(u_ll))
+    for k in eachcomponent(equations)
+        rho_k = u_ll[3 + (k - 1) * 5 + 1] # Extract densities from conserved variable vector
+        charge_ratio_ll[k] = rho_k * charge_to_mass[k]
+        total_electron_charge += charge_ratio_ll[k]
+    end
+    charge_ratio_ll ./= total_electron_charge
+
+    # Compute auxiliary variables
+    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
+                                                                                                           equations)
+
+    f = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    if nonconservative_term == 1
+        f[1] = v1_plus_ll
+        f[2] = v2_plus_ll
+        f[3] = v3_plus_ll
+
+        for k in eachcomponent(equations)
+            # Compute Godunov-Powell term
+            f2 = charge_ratio_ll[k] * B1_ll
+            f3 = charge_ratio_ll[k] * B2_ll
+            f4 = charge_ratio_ll[k] * B3_ll
+            f5 = (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll)
+
+            set_component!(f, k, 0, f2, f3, f4, f5, equations)
+        end
+    elseif nonconservative_term == 2
+        # Compute Lorentz term
+
+        for k in eachcomponent(equations)
+            f2 = charge_ratio_ll[k]
+            f3 = charge_ratio_ll[k]
+            f4 = charge_ratio_ll[k]
+            f5 = (vk1_plus_ll[k] * normal_direction[1] +
+                  vk2_plus_ll[k] * normal_direction[2] +
+                  vk3_plus_ll[k] * normal_direction[3])
+
+            set_component!(f, k, 0, f2, f3, f4, f5, equations)
+        end
+    elseif nonconservative_term == 3
+        # Compute GLM term
+        v_plus_dot_n_ll = (v1_plus_ll * normal_direction[1] +
+                           v2_plus_ll * normal_direction[2] +
+                           v3_plus_ll * normal_direction[3])
+        for k in eachcomponent(equations)
+            f5 = v_plus_dot_n_ll * psi_ll
+            set_component!(f, k, 0, 0, 0, 0, f5, equations)
+        end
+        # Compute GLM term for psi
+        f[end] = v_plus_dot_n_ll
+    elseif nonconservative_term == 4
+        # Multi-ion term (vanishes for NCOMP==1) for B1
+        for k in eachcomponent(equations)
+            f5 = B1_ll
+            set_component!(f, k, 0, 0, 0, 0, f5, equations)
+        end
+    elseif nonconservative_term == 5
+        # Multi-ion term (vanishes for NCOMP==1) for B2
+        for k in eachcomponent(equations)
+            f5 = B2_ll
+            set_component!(f, k, 0, 0, 0, 0, f5, equations)
+        end
+    elseif nonconservative_term == 6
+        # Multi-ion term (vanishes for NCOMP==1) for B3
+        for k in eachcomponent(equations)
+            f5 = B3_ll
+            set_component!(f, k, 0, 0, 0, 0, f5, equations)
+        end
+    end
+
+    return SVector(f)
+end
+
+"""
+    flux_nonconservative_ruedaramirez_etal(u_ll, normal_direction_avg::AbstractVector,
+                                           equations::IdealGlmMhdMultiIonEquations3D,
+                                           nonconservative_type::NonConservativeSymmetric,
+                                           nonconservative_term::Integer)
+
+Symmetric part of the multi-ion 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.
+This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl.
+"""
+@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr,
+                                                                     normal_direction::AbstractVector,
+                                                                     equations::IdealGlmMhdMultiIonEquations3D,
+                                                                     nonconservative_type::NonConservativeSymmetric,
+                                                                     nonconservative_term::Integer)
+    @unpack charge_to_mass = equations
+    # Unpack left and right states to get the magnetic field
+    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
+    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
+    psi_ll = divergence_cleaning_field(u_ll, equations)
+    psi_rr = divergence_cleaning_field(u_rr, equations)
+    B_dot_n_ll = B1_ll * normal_direction[1] +
+                 B2_ll * normal_direction[2] +
+                 B3_ll * normal_direction[3]
+    B_dot_n_rr = B1_rr * normal_direction[1] +
+                 B2_rr * normal_direction[2] +
+                 B3_rr * normal_direction[3]
+    B_dot_n_avg = 0.5f0 * (B_dot_n_ll + B_dot_n_rr)
+
+    # Compute important averages
+    B1_avg = 0.5f0 * (B1_ll + B1_rr)
+    B2_avg = 0.5f0 * (B2_ll + B2_rr)
+    B3_avg = 0.5f0 * (B3_ll + B3_rr)
+
+    # Compute charge ratio of u_ll
+    charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)})
+    total_electron_charge = zero(eltype(u_ll))
+    for k in eachcomponent(equations)
+        rho_k = u_ll[3 + (k - 1) * 5 + 1] # Extract densities from conserved variable vector
+        charge_ratio_ll[k] = rho_k * charge_to_mass[k]
+        total_electron_charge += charge_ratio_ll[k]
+    end
+    charge_ratio_ll ./= total_electron_charge
+
+    # Compute auxiliary variables
+    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
+                                                                                                           equations)
+    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
+                                                                                                           equations)
+
+    f = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    if nonconservative_term == 1
+        # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is 
+        # multiplied by 0.5 whenever it's used in the Trixi code)
+        f[1] = 2 * B_dot_n_avg
+        f[2] = 2 * B_dot_n_avg
+        f[3] = 2 * B_dot_n_avg
+
+        for k in eachcomponent(equations)
+            # Compute Godunov-Powell term
+            f2 = B_dot_n_avg
+            f3 = B_dot_n_avg
+            f4 = B_dot_n_avg
+            f5 = B_dot_n_avg
+            # Add to the flux vector (multiply by 2 because the non-conservative flux is 
+            # multiplied by 0.5 whenever it's used in the Trixi code)
+            set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
+                           equations)
+        end
+    elseif nonconservative_term == 2
+        # Compute Lorentz term
+
+        # Important averages
+        mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
+        mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
+        mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
+
+        # Mean electron pressure
+        pe_ll = equations.electron_pressure(u_ll, equations)
+        pe_rr = equations.electron_pressure(u_rr, equations)
+        pe_mean = 0.5f0 * (pe_ll + pe_rr)
+
+        for k in eachcomponent(equations)
+            f2 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[1] -
+                  B_dot_n_avg * B1_avg)
+            f3 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[2] -
+                  B_dot_n_avg * B2_avg)
+            f4 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[3] -
+                  B_dot_n_avg * B3_avg)
+            f5 = pe_mean
+
+            set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
+                           equations)
+        end
+    elseif nonconservative_term == 3
+        # Compute GLM term
+        psi_avg = 0.5f0 * (psi_ll + psi_rr)
+        for k in eachcomponent(equations)
+            f5 = psi_avg
+            # (multiply by 2 because the non-conservative flux is 
+            # multiplied by 0.5 whenever it's used in the Trixi code)
+            set_component!(f, k, 0, 0, 0, 0, 2 * f5, equations)
+        end
+        # Compute GLM term for psi (multiply by 2 because the non-conservative flux is 
+        # multiplied by 0.5 whenever it's used in the Trixi code)
+        f[end] = 2 * psi_avg
+    elseif nonconservative_term == 4
+        # Multi-ion term (vanishes for NCOMP==1) for B1
+        for k in eachcomponent(equations)
+            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
+            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
+            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
+            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
+            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
+            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
+            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
+            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
+            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
+
+            f5 = ((vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) *
+                  normal_direction[2] +
+                  (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) *
+                  normal_direction[3])
+
+            # Add to the flux vector (multiply by 2 because the non-conservative flux is 
+            # multiplied by 0.5 whenever it's used in the Trixi code)
+            set_component!(f, k, 0, 0, 0, 0, 2 * f5,
+                           equations)
+        end
+    elseif nonconservative_term == 5
+        # Multi-ion term (vanishes for NCOMP==1) for B2
+        for k in eachcomponent(equations)
+            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
+            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
+            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
+            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
+            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
+            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
+            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
+            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
+            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
+
+            f5 = ((vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) *
+                  normal_direction[1] +
+                  (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg) *
+                  normal_direction[3])
+            # Add to the flux vector (multiply by 2 because the non-conservative flux is 
+            # multiplied by 0.5 whenever it's used in the Trixi code)
+            set_component!(f, k, 0, 0, 0, 0, 2 * f5,
+                           equations)
+        end
+    elseif nonconservative_term == 6
+        # Multi-ion term (vanishes for NCOMP==1) for B3
+        for k in eachcomponent(equations)
+            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
+            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
+            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
+            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
+            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
+            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
+            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
+            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
+            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
+
+            f5 = ((vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg) *
+                  normal_direction[1] +
+                  (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg) *
+                  normal_direction[2])
+            # Add to the flux vector (multiply by 2 because the non-conservative flux is 
+            # multiplied by 0.5 whenever it's used in the Trixi code)
+            set_component!(f, k, 0, 0, 0, 0, 2 * f5,
+                           equations)
+        end
+    end
+
+    return SVector(f)
+end
+
 """
     flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
                                  equations::IdealGlmMhdMultiIonEquations3D)