Add a central non-conservative flux for 3D multi-ion MHD and a convergence test

examples/p4est_3d_dgsem/elixir_mhdmultiion_convergence.jl

@@ -0,0 +1,207 @@
+
+using OrdinaryDiffEqLowStorageRK
+using Trixi
+
+# 3D version of multi-ion MHD convergence test from
+#   - Rueda-Ramírez, A. M., Sikstel, A., & Gassner, G. J. (2025). An entropy-stable discontinuous Galerkin 
+#     discretization of the ideal multi-ion magnetohydrodynamics system. Journal of Computational Physics, 523, 113655.
+
+###############################################################################
+"""
+  electron_pressure_alpha(u, equations::IdealGlmMhdMultiIonEquations3D)
+Returns a fraction (alpha) of the total ion pressure for the electron pressure.
+"""
+function electron_pressure_alpha(u, equations::IdealGlmMhdMultiIonEquations3D)
+    alpha = 0.2
+    prim = cons2prim(u, equations)
+    p_e = zero(u[1])
+    for k in eachcomponent(equations)
+        _, _, _, _, p_k = Trixi.get_component(k, prim, equations)
+        p_e += p_k
+    end
+    return alpha * p_e
+end
+# semidiscretization of the ideal multi-ion MHD equations
+equations = IdealGlmMhdMultiIonEquations3D(gammas = (2.0, 4.0),
+                                           charge_to_mass = (2.0, 1.0),
+                                           electron_pressure = electron_pressure_alpha)
+
+"""
+Initial (and exact) solution for the the manufactured solution test. Runs with 
+* gammas = (2.0, 4.0),
+* charge_to_mass = (2.0, 1.0)
+* Domain size: [-1,1]²
+"""
+function initial_condition_manufactured_solution(x, t,
+                                                 equations::IdealGlmMhdMultiIonEquations3D)
+    am = 0.1
+    om = π
+    h = am * sin(om * (x[1] + x[2] + x[3] - t)) + 2
+    hh1 = am * 0.4 * sin(om * (x[1] + x[2] + x[3] - t)) + 1
+    hh2 = h - hh1
+
+    rho_1 = hh1
+    rhou_1 = hh1
+    rhov_1 = hh1
+    rhow_1 = 0.1 * hh1
+    rhoe_1 = 2 * hh1^2 + hh1
+    rho_2 = hh2
+    rhou_2 = hh2
+    rhov_2 = hh2
+    rhow_2 = 0.1 * hh2
+    rhoe_2 = 2 * hh2^2 + hh2
+    B1 = 0.5 * h
+    B2 = -0.25 * h
+    B3 = -0.25 * h
+
+    return SVector{nvariables(equations), real(equations)}([
+                                                               B1,
+                                                               B2,
+                                                               B3,
+                                                               rho_1,
+                                                               rhou_1,
+                                                               rhov_1,
+                                                               rhow_1,
+                                                               rhoe_1,
+                                                               rho_2,
+                                                               rhou_2,
+                                                               rhov_2,
+                                                               rhow_2,
+                                                               rhoe_2,
+                                                               zero(eltype(x))
+                                                           ])
+end
+
+"""
+Source term that corresponds to the manufactured solution test. Runs with 
+* gammas = (2.0, 4.0),
+* charge_to_mass = (2.0, 1.0)
+* Domain size: [-1,1]²
+"""
+function source_terms_manufactured_solution_pe(u, x, t,
+                                               equations::IdealGlmMhdMultiIonEquations3D)
+    am = 0.1
+    om = pi
+    h1 = am * sin(om * (x[1] + x[2] + x[3] - t))
+    hx = am * om * cos(om * (x[1] + x[2] + x[3] - t))
+
+    s1 = (11 * hx) / 25
+    s2 = (30615 * hx * h1^2 + 156461 * hx * h1 + 191990 * hx) / (35000 * h1 + 75000)
+    s3 = (30615 * hx * h1^2 + 156461 * hx * h1 + 191990 * hx) / (35000 * h1 + 75000)
+    s4 = (30615 * hx * h1^2 + 142601 * hx * h1 + 162290 * hx) / (35000 * h1 + 75000)
+    s5 = (4735167957644739545 * hx * h1^2 + 22683915240114795103 * hx * h1 +
+          26562869799135852870 * hx) / (1863579030125050000 * h1 + 3993383635982250000)
+    s6 = (33 * hx) / 50
+    s7 = (63915 * hx * h1^2 + 245476 * hx * h1 + 233885 * hx) / (17500 * h1 + 37500)
+    s8 = (63915 * hx * h1^2 + 245476 * hx * h1 + 233885 * hx) / (17500 * h1 + 37500)
+    s9 = (63915 * hx * h1^2 + 235081 * hx * h1 + 211610 * hx) / (17500 * h1 + 37500)
+    s10 = (1619415 * hx * h1^2 + 6083946 * hx * h1 + 5629335 * hx) / (175000 * h1 + 375000)
+    s11 = (11 * hx) / 20
+    s12 = -((11 * hx) / 40)
+    s13 = -((11 * hx) / 40)
+
+    s = SVector{nvariables(equations), real(equations)}(s11,
+                                                        s12,
+                                                        s13,
+                                                        s1,
+                                                        s2,
+                                                        s3,
+                                                        s4,
+                                                        s5,
+                                                        s6,
+                                                        s7,
+                                                        s8,
+                                                        s9,
+                                                        s10,
+                                                        zero(eltype(u)))
+    S_std = source_terms_lorentz(u, x, t, equations::IdealGlmMhdMultiIonEquations3D)
+
+    return S_std + s
+end
+
+initial_condition = initial_condition_manufactured_solution
+source_terms = source_terms_manufactured_solution_pe
+
+volume_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal)
+surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)
+
+polydeg = 3
+solver = DGSEM(polydeg = polydeg, surface_flux = surface_flux,
+               volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
+
+coordinates_min = (-1.0, -1.0, -1.0)
+coordinates_max = (1.0, 1.0, 1.0)
+# Mapping as described in https://arxiv.org/abs/2012.12040
+function mapping(xi_, eta_, zeta_)
+    # Transform input variables between -1 and 1 onto [0,3]
+    xi = 1.5 * xi_ + 1.5
+    eta = 1.5 * eta_ + 1.5
+    zeta = 1.5 * zeta_ + 1.5
+
+    y = eta +
+        0.1 * (cos(1.5 * pi * (2 * xi - 3) / 3) *
+         cos(0.5 * pi * (2 * eta - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    x = xi +
+        0.1 * (cos(0.5 * pi * (2 * xi - 3) / 3) *
+         cos(2 * pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    z = zeta +
+        0.1 * (cos(0.5 * pi * (2 * x - 3) / 3) *
+         cos(pi * (2 * y - 3) / 3) *
+         cos(0.5 * pi * (2 * zeta - 3) / 3))
+
+    # Go back to [-1,1]^3
+    x = x * 2 / 3 - 1
+    y = y * 2 / 3 - 1
+    z = z * 2 / 3 - 1
+
+    return SVector(x, y, z)
+end
+
+cells_per_dimension = (8, 8, 8)
+mesh = P4estMesh(cells_per_dimension,
+                 polydeg = polydeg, #initial_refinement_level = 0,
+                 mapping = mapping,
+                 #coordinates_min = coordinates_min, coordinates_max = coordinates_max,
+                 periodicity = true)
+
+semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
+                                    source_terms = source_terms)
+
+###############################################################################
+# ODE solvers, callbacks etc.
+
+tspan = (0.0, 1.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,
+                                     output_directory = joinpath(@__DIR__, "out"))
+
+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
+sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false),
+            dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
+            save_everystep = false, callback = callbacks);

src/equations/ideal_glm_mhd_multiion_3d.jl

@@ -816,6 +816,133 @@ The term is composed of four individual non-conservative terms:
     return SVector(f)
 end
 
+@inline function flux_nonconservative_central(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)
+    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)
+    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)
+    psi_avg = 0.5f0 * (psi_ll + psi_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)
+
+    # 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)
+    v_plus_dot_n_ll = (v1_plus_ll * normal_direction[1] +
+                       v2_plus_ll * normal_direction[2] +
+                       v3_plus_ll * normal_direction[3])
+    f = zero(MVector{nvariables(equations), eltype(u_ll)})
+
+    # 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 * v1_plus_ll * B_dot_n_avg
+    f[2] = 2 * v2_plus_ll * B_dot_n_avg
+    f[3] = 2 * v3_plus_ll * B_dot_n_avg
+
+    for k in eachcomponent(equations)
+        # Compute terms for each species
+        # (we multiply by 2 because the non-conservative flux is multiplied by 0.5 whenever it's used in the Trixi code)
+
+        # Compute term Lorentz term
+        f2_ll = ((0.5f0 * mag_norm_ll + pe_ll) * normal_direction[1] -
+                 B_dot_n_ll * B1_ll)
+        f2_rr = ((0.5f0 * mag_norm_rr + pe_rr) * normal_direction[1] -
+                 B_dot_n_rr * B1_rr)
+        f2 = charge_ratio_ll[k] * (f2_ll + f2_rr)
+
+        f3_ll = ((0.5f0 * mag_norm_ll + pe_ll) * normal_direction[2] -
+                 B_dot_n_ll * B2_ll)
+        f3_rr = ((0.5f0 * mag_norm_rr + pe_rr) * normal_direction[2] -
+                 B_dot_n_rr * B2_rr)
+        f3 = charge_ratio_ll[k] * (f3_ll + f3_rr)
+
+        f4_ll = ((0.5f0 * mag_norm_ll + pe_ll) * normal_direction[3] -
+                 B_dot_n_ll * B3_ll)
+        f4_rr = ((0.5f0 * mag_norm_rr + pe_rr) * normal_direction[3] -
+                 B_dot_n_rr * B3_rr)
+        f4 = charge_ratio_ll[k] * (f4_ll + f4_rr)
+
+        f5 = (vk1_plus_ll[k] * normal_direction[1] +
+              vk2_plus_ll[k] * normal_direction[2] +
+              vk3_plus_ll[k] * normal_direction[3]) * pe_mean * 2
+
+        # Compute multi-ion term (vanishes for NCOMP==1)
+        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]
+        f5 += ((B2_ll * ((vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) +
+                 (vk1_minus_rr * B2_rr - vk2_minus_rr * B1_rr)) +
+                B3_ll * ((vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll) +
+                 (vk1_minus_rr * B3_rr - vk3_minus_rr * B1_rr))) *
+               normal_direction[1] +
+               (B1_ll * ((vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) +
+                 (vk2_minus_rr * B1_rr - vk1_minus_rr * B2_rr)) +
+                B3_ll * ((vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll) +
+                 (vk2_minus_rr * B3_rr - vk3_minus_rr * B2_rr))) *
+               normal_direction[2] +
+               (B1_ll * ((vk3_minus_ll * B1_ll - vk1_minus_ll * B3_ll) +
+                 (vk3_minus_rr * B1_rr - vk1_minus_rr * B3_rr)) +
+                B2_ll * ((vk3_minus_ll * B2_ll - vk2_minus_ll * B3_ll) +
+                 (vk3_minus_rr * B2_rr - vk2_minus_rr * B3_rr))) *
+               normal_direction[3])
+
+        # Compute Godunov-Powell term
+        f2 += charge_ratio_ll[k] * B1_ll * B_dot_n_avg * 2
+        f3 += charge_ratio_ll[k] * B2_ll * B_dot_n_avg * 2
+        f4 += charge_ratio_ll[k] * B3_ll * B_dot_n_avg * 2
+        f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
+              B_dot_n_avg * 2
+
+        # Compute GLM term for the energy
+        f5 += v_plus_dot_n_ll * psi_ll * psi_avg * 2
+
+        # Add to the flux vector
+        set_component!(f, k, 0, f2, f3, f4, 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 * v_plus_dot_n_ll * psi_avg
+
+    return SVector(f)
+end
+
 """
     flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations3D)
 

test/test_p4est_3d_mhdmultiion.jl

@@ -109,5 +109,49 @@ end
         @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
     end
 end
+
+@trixi_testset "elixir_mhdmultiion_convergence.jl" begin
+    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_convergence.jl"),
+                        l2=[
+                            2.7007694451840977e-5,
+                            2.252632596997783e-5,
+                            1.830892822103072e-5,
+                            1.7457386132678724e-5,
+                            3.965825276181703e-5,
+                            6.886878771068099e-5,
+                            3.216774733720572e-5,
+                            0.00013796601797391608,
+                            2.762642533644496e-5,
+                            7.877500410069398e-5,
+                            0.00012184040930856932,
+                            8.918795955887214e-5,
+                            0.0002122739932637704,
+                            1.0532691581216071e-6
+                        ],
+                        linf=[
+                            0.0005846835977684206,
+                            0.00031591380039502903,
+                            0.0002529555339790268,
+                            0.0003873459403432866,
+                            0.0007355557980894822,
+                            0.0012929706727252688,
+                            0.0002558003707378437,
+                            0.0028085112041740246,
+                            0.0006114366794293113,
+                            0.001257825301983151,
+                            0.0018924211424776738,
+                            0.0007347447431757664,
+                            0.004148291057411768,
+                            1.8948511576480304e-5
+                        ], tspan=(0.0, 0.05))
+    # Ensure that we do not have excessive memory allocations
+    # (e.g., from type instabilities)
+    let
+        t = sol.t[end]
+        u_ode = sol.u[end]
+        du_ode = similar(u_ode)
+        @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000
+    end
+end
 end
 end # module