CoCalc -- ideal_glm_mhd_multiion

GitHub Repository: trixi-framework/Trixi.jl
Path: blob/main/src/equations/ideal_glm_mhd_multiion_3d.jl
⁵⁵⁸⁶ views
1
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
2
# Since these FMAs can increase the performance of many numerical algorithms,
3
# we need to opt-in explicitly.
4
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
5
@muladd begin
6
#! format: noindent
7

8
@doc raw"""
9
    IdealGlmMhdMultiIonEquations3D(; gammas, charge_to_mass,
10
                                   electron_pressure = electron_pressure_zero,
11
                                   initial_c_h = NaN)
12

13
The ideal compressible multi-ion MHD equations in three space dimensions augmented with a
14
generalized Langange multipliers (GLM) divergence-cleaning technique. This is a
15
multi-species variant of the ideal GLM-MHD equations for calorically perfect plasmas
16
with independent momentum and energy equations for each ion species. This implementation
17
assumes that the equations are non-dimensionalized, such that the vacuum permeability is ``\mu_0 = 1``.
18

19
In case of more than one ion species, the specific heat capacity ratios `gammas` and the charge-to-mass
20
ratios `charge_to_mass` should be passed as tuples, e.g., `gammas=(1.4, 1.667)`.
21

22
The argument `electron_pressure` can be used to pass a function that computes the electron
23
pressure as a function of the state `u` with the signature `electron_pressure(u, equations)`.
24
By default, the electron pressure is zero.
25

26
The argument `initial_c_h` can be used to set the GLM divergence-cleaning speed. Note that
27
`initial_c_h = 0` deactivates the divergence cleaning. The callback [`GlmSpeedCallback`](@ref)
28
can be used to adjust the GLM divergence-cleaning speed according to the time-step size.
29

30
References:
31
- G. Toth, A. Glocer, Y. Ma, D. Najib, Multi-Ion Magnetohydrodynamics 429 (2010). Numerical
32
  Modeling of Space Plasma Flows, 213–218.
33
- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization
34
  of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics.
35
  [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655).
36

37
!!! info "The multi-ion GLM-MHD equations require source terms"
38
    In case of more than one ion species, the multi-ion GLM-MHD equations should ALWAYS be used
39
    with [`source_terms_lorentz`](@ref).
40
"""
41
struct IdealGlmMhdMultiIonEquations3D{NVARS, NCOMP, RealT <: Real,
42
                                      ElectronPressure} <:
43
       AbstractIdealGlmMhdMultiIonEquations{3, NVARS, NCOMP}
44
    gammas::SVector{NCOMP, RealT} # Heat capacity ratios
45
    charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios
46
    electron_pressure::ElectronPressure # Function to compute the electron pressure
47
    c_h::RealT # GLM cleaning speed
48
    function IdealGlmMhdMultiIonEquations3D{NVARS, NCOMP, RealT,
49
                                            ElectronPressure}(gammas
50
                                                              ::SVector{NCOMP, RealT},
51
                                                              charge_to_mass
52
                                                              ::SVector{NCOMP, RealT},
53
                                                              electron_pressure
54
                                                              ::ElectronPressure,
55
                                                              c_h::RealT) where
56
             {NVARS, NCOMP, RealT <: Real, ElectronPressure}
57
        NCOMP >= 1 ||
58
            throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value"))
59

60
        return new(gammas, charge_to_mass, electron_pressure, c_h)
61
    end
62
end
63

64
function IdealGlmMhdMultiIonEquations3D(; gammas, charge_to_mass,
65
                                        electron_pressure = electron_pressure_zero,
66
                                        initial_c_h = convert(eltype(gammas), NaN))
67
    _gammas = promote(gammas...)
68
    _charge_to_mass = promote(charge_to_mass...)
69
    RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass))
70
    __gammas = SVector(map(RealT, _gammas))
71
    __charge_to_mass = SVector(map(RealT, _charge_to_mass))
72

73
    NVARS = length(_gammas) * 5 + 4
74
    NCOMP = length(_gammas)
75

76
    return IdealGlmMhdMultiIonEquations3D{NVARS, NCOMP, RealT,
77
                                          typeof(electron_pressure)}(__gammas,
78
                                                                     __charge_to_mass,
79
                                                                     electron_pressure,
80
                                                                     initial_c_h)
81
end
82

83
# Outer constructor for `@reset` works correctly
84
function IdealGlmMhdMultiIonEquations3D(gammas, charge_to_mass, electron_pressure, c_h)
85
    return IdealGlmMhdMultiIonEquations3D(gammas = gammas,
86
                                          charge_to_mass = charge_to_mass,
87
                                          electron_pressure = electron_pressure,
88
                                          initial_c_h = c_h)
89
end
90

91
@inline function Base.real(::IdealGlmMhdMultiIonEquations3D{NVARS, NCOMP, RealT}) where {
92
                                                                                         NVARS,
93
                                                                                         NCOMP,
94
                                                                                         RealT
95
                                                                                         }
96
    return RealT
97
end
98

99
"""
100
    initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdMultiIonEquations3D)
101

102
A weak blast wave (adapted to multi-ion MHD) from
103
- Hennemann, S., Rueda-Ramírez, A. M., Hindenlang, F. J., & Gassner, G. J. (2021). A provably entropy
104
  stable subcell shock capturing approach for high order split form DG for the compressible Euler equations.
105
  Journal of Computational Physics, 426, 109935. [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044).
106
  [DOI: 10.1016/j.jcp.2020.109935](https://doi.org/10.1016/j.jcp.2020.109935)
107
"""
108
function initial_condition_weak_blast_wave(x, t,
109
                                           equations::IdealGlmMhdMultiIonEquations3D)
110
    # Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
111
    # Same discontinuity in the velocities but with magnetic fields
112
    RealT = eltype(x)
113
    # Set up polar coordinates
114
    inicenter = (0, 0, 0)
115
    x_norm = x[1] - inicenter[1]
116
    y_norm = x[2] - inicenter[2]
117
    z_norm = x[3] - inicenter[3]
118
    r = sqrt(x_norm^2 + y_norm^2 + z_norm^2)
119
    phi = atan(y_norm, x_norm)
120
    theta = iszero(r) ? zero(RealT) : acos(z_norm / r)
121

122
    # Calculate primitive variables
123
    rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691)
124
    v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) * sin(theta)
125
    v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) * sin(theta)
126
    v3 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(theta)
127
    p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245)
128

129
    prim = zero(MVector{nvariables(equations), real(equations)})
130
    prim[1] = 1
131
    prim[2] = 1
132
    prim[3] = 1
133

134
    for k in eachcomponent(equations)
135
        # We initialize each species with a fraction of the total density `rho`, such
136
        # that the sum of the densities is `rho := density(prim, equations)`. The density of
137
        # a species is double the density of the next species.
138
        fraction = 2^(k - 1) * (1 - 2) / (1 - 2^ncomponents(equations))
139
        set_component!(prim, k,
140
                       fraction * rho, v1,
141
                       v2, v3, p, equations)
142
    end
143

144
    return prim2cons(SVector(prim), equations)
145
end
146

147
# 3D flux of the multi-ion GLM-MHD equations in the direction `orientation`
148
@inline function flux(u, orientation::Integer,
149
                      equations::IdealGlmMhdMultiIonEquations3D)
150
    B1, B2, B3 = magnetic_field(u, equations)
151
    psi = divergence_cleaning_field(u, equations)
152

153
    v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u,
154
                                                                                         equations)
155

156
    mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2)
157
    div_clean_energy = 0.5f0 * psi^2
158

159
    f = zero(MVector{nvariables(equations), eltype(u)})
160

161
    if orientation == 1
162
        f[1] = equations.c_h * psi
163
        f[2] = v1_plus * B2 - v2_plus * B1
164
        f[3] = v1_plus * B3 - v3_plus * B1
165

166
        for k in eachcomponent(equations)
167
            rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, u, equations)
168
            rho_inv = 1 / rho
169
            v1 = rho_v1 * rho_inv
170
            v2 = rho_v2 * rho_inv
171
            v3 = rho_v3 * rho_inv
172
            kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
173

174
            gamma = equations.gammas[k]
175
            p = (gamma - 1) * (rho_e_total - kin_en - mag_en - div_clean_energy)
176

177
            f1 = rho_v1
178
            f2 = rho_v1 * v1 + p
179
            f3 = rho_v1 * v2
180
            f4 = rho_v1 * v3
181
            f5 = (kin_en + gamma * p / (gamma - 1)) * v1 + 2 * mag_en * vk1_plus[k] -
182
                 B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) +
183
                 equations.c_h * psi * B1
184

185
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
186
        end
187
        f[end] = equations.c_h * B1
188
    elseif orientation == 2
189
        f[1] = v2_plus * B1 - v1_plus * B2
190
        f[2] = equations.c_h * psi
191
        f[3] = v2_plus * B3 - v3_plus * B2
192

193
        for k in eachcomponent(equations)
194
            rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, u, equations)
195
            rho_inv = 1 / rho
196
            v1 = rho_v1 * rho_inv
197
            v2 = rho_v2 * rho_inv
198
            v3 = rho_v3 * rho_inv
199
            kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
200

201
            gamma = equations.gammas[k]
202
            p = (gamma - 1) * (rho_e_total - kin_en - mag_en - div_clean_energy)
203

204
            f1 = rho_v2
205
            f2 = rho_v2 * v1
206
            f3 = rho_v2 * v2 + p
207
            f4 = rho_v2 * v3
208
            f5 = (kin_en + gamma * p / (gamma - 1)) * v2 + 2 * mag_en * vk2_plus[k] -
209
                 B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) +
210
                 equations.c_h * psi * B2
211

212
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
213
        end
214
        f[end] = equations.c_h * B2
215
    else #if orientation == 3
216
        f[1] = v3_plus * B1 - v1_plus * B3
217
        f[2] = v3_plus * B2 - v2_plus * B3
218
        f[3] = equations.c_h * psi
219

220
        for k in eachcomponent(equations)
221
            rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, u, equations)
222
            rho_inv = 1 / rho
223
            v1 = rho_v1 * rho_inv
224
            v2 = rho_v2 * rho_inv
225
            v3 = rho_v3 * rho_inv
226
            kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
227

228
            gamma = equations.gammas[k]
229
            p = (gamma - 1) * (rho_e_total - kin_en - mag_en - div_clean_energy)
230

231
            f1 = rho_v3
232
            f2 = rho_v3 * v1
233
            f3 = rho_v3 * v2
234
            f4 = rho_v3 * v3 + p
235
            f5 = (kin_en + gamma * p / (gamma - 1)) * v3 + 2 * mag_en * vk3_plus[k] -
236
                 B3 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) +
237
                 equations.c_h * psi * B3
238

239
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
240
        end
241
        f[end] = equations.c_h * B3
242
    end
243

244
    return SVector(f)
245
end
246

247
# Calculate 1D flux for a single point in the normal direction
248
# Note, this directional vector is not normalized
249
@inline function flux(u, normal_direction::AbstractVector,
250
                      equations::IdealGlmMhdMultiIonEquations3D)
251
    B1, B2, B3 = magnetic_field(u, equations)
252
    psi = divergence_cleaning_field(u, equations)
253

254
    v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u,
255
                                                                                         equations)
256

257
    mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2)
258
    div_clean_energy = 0.5f0 * psi^2
259
    B_normal = B1 * normal_direction[1] + B2 * normal_direction[2] +
260
               B3 * normal_direction[3]
261

262
    f = zero(MVector{nvariables(equations), eltype(u)})
263

264
    f[1] = (equations.c_h * psi * normal_direction[1] +
265
            (v2_plus * B1 - v1_plus * B2) * normal_direction[2] +
266
            (v3_plus * B1 - v1_plus * B3) * normal_direction[3])
267
    f[2] = ((v1_plus * B2 - v2_plus * B1) * normal_direction[1] +
268
            equations.c_h * psi * normal_direction[2] +
269
            (v3_plus * B2 - v2_plus * B3) * normal_direction[3])
270
    f[3] = ((v1_plus * B3 - v3_plus * B1) * normal_direction[1] +
271
            (v2_plus * B3 - v3_plus * B2) * normal_direction[2] +
272
            equations.c_h * psi * normal_direction[3])
273

274
    for k in eachcomponent(equations)
275
        rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, u, equations)
276
        rho_inv = 1 / rho
277
        v1 = rho_v1 * rho_inv
278
        v2 = rho_v2 * rho_inv
279
        v3 = rho_v3 * rho_inv
280
        kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
281
        v_normal = v1 * normal_direction[1] + v2 * normal_direction[2] +
282
                   v3 * normal_direction[3]
283
        rho_v_normal = rho * v_normal
284
        vk_plus_normal = vk1_plus[k] * normal_direction[1] +
285
                         vk2_plus[k] * normal_direction[2] +
286
                         vk3_plus[k] * normal_direction[3]
287

288
        gamma = equations.gammas[k]
289
        p = (gamma - 1) * (rho_e_total - kin_en - mag_en - div_clean_energy)
290

291
        f1 = rho_v_normal
292
        f2 = rho_v_normal * v1 + p * normal_direction[1]
293
        f3 = rho_v_normal * v2 + p * normal_direction[2]
294
        f4 = rho_v_normal * v3 + p * normal_direction[3]
295
        f5 = (kin_en + gamma * p / (gamma - 1)) * v_normal +
296
             2 * mag_en * vk_plus_normal -
297
             B_normal * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) +
298
             equations.c_h * psi * B_normal
299

300
        set_component!(f, k, f1, f2, f3, f4, f5, equations)
301
    end
302
    f[end] = equations.c_h * B_normal
303

304
    return SVector(f)
305
end
306

307
"""
308
    flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
309
                                           orientation_or_normal_direction,
310
                                           equations::IdealGlmMhdMultiIonEquations3D)
311

312
Entropy-conserving non-conservative two-point "flux" as described in
313
- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization
314
  of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics.
315
  [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655).
316

317
!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi.jl"
318
    The non-conservative fluxes derived in the reference above are written as the product
319
    of local and symmetric parts and are meant to be used in the same way as the conservative
320
    fluxes (i.e., flux + flux_noncons in both volume and surface integrals). In this routine,
321
    the fluxes are multiplied by 2 because the non-conservative fluxes are always multiplied
322
    by 0.5 whenever they are used in the Trixi.jl code.
323

324
The term is composed of four individual non-conservative terms:
325
1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and
326
   is evaluated as a function of the charge-averaged velocity.
327
2. The Lorentz-force term, which becomes a conservative term in the limit of one ion species for vanishing
328
   electron pressure gradients.
329
3. The "multi-ion" term, which vanishes in the limit of one ion species.
330
4. The GLM term, which is needed for Galilean invariance.
331
"""
332
@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
333
                                                        orientation::Integer,
334
                                                        equations::IdealGlmMhdMultiIonEquations3D)
335
    @unpack charge_to_mass = equations
336
    # Unpack left and right states to get the magnetic field
337
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
338
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
339
    psi_ll = divergence_cleaning_field(u_ll, equations)
340
    psi_rr = divergence_cleaning_field(u_rr, equations)
341

342
    # Compute important averages
343
    B1_avg = 0.5f0 * (B1_ll + B1_rr)
344
    B2_avg = 0.5f0 * (B2_ll + B2_rr)
345
    B3_avg = 0.5f0 * (B3_ll + B3_rr)
346
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
347
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
348
    mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
349
    psi_avg = 0.5f0 * (psi_ll + psi_rr)
350

351
    # Mean electron pressure
352
    pe_ll = equations.electron_pressure(u_ll, equations)
353
    pe_rr = equations.electron_pressure(u_rr, equations)
354
    pe_mean = 0.5f0 * (pe_ll + pe_rr)
355

356
    # Compute charge ratio of u_ll
357
    charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)})
358
    total_electron_charge = zero(eltype(u_ll))
359
    for k in eachcomponent(equations)
360
        rho_k = u_ll[3 + (k - 1) * 5 + 1]
361
        charge_ratio_ll[k] = rho_k * charge_to_mass[k]
362
        total_electron_charge += charge_ratio_ll[k]
363
    end
364
    charge_ratio_ll ./= total_electron_charge
365

366
    # Compute auxiliary variables
367
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
368
                                                                                                           equations)
369
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
370
                                                                                                           equations)
371

372
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
373

374
    if orientation == 1
375
        # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is
376
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
377
        f[1] = 2 * v1_plus_ll * B1_avg
378
        f[2] = 2 * v2_plus_ll * B1_avg
379
        f[3] = 2 * v3_plus_ll * B1_avg
380

381
        for k in eachcomponent(equations)
382
            # Compute term Lorentz term
383
            f2 = charge_ratio_ll[k] * (0.5f0 * mag_norm_avg - B1_avg * B1_avg + pe_mean)
384
            f3 = charge_ratio_ll[k] * (-B1_avg * B2_avg)
385
            f4 = charge_ratio_ll[k] * (-B1_avg * B3_avg)
386
            f5 = vk1_plus_ll[k] * pe_mean
387

388
            # Compute multi-ion term (vanishes for NCOMP==1)
389
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
390
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
391
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
392
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
393
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
394
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
395
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
396
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
397
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
398
            f5 += (B2_ll * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) +
399
                   B3_ll * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg))
400

401
            # Compute Godunov-Powell term
402
            f2 += charge_ratio_ll[k] * B1_ll * B1_avg
403
            f3 += charge_ratio_ll[k] * B2_ll * B1_avg
404
            f4 += charge_ratio_ll[k] * B3_ll * B1_avg
405
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
406
                  B1_avg
407

408
            # Compute GLM term for the energy
409
            f5 += v1_plus_ll * psi_ll * psi_avg
410

411
            # Add to the flux vector (multiply by 2 because the non-conservative flux is
412
            # multiplied by 0.5 whenever it's used in the Trixi code)
413
            set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
414
                           equations)
415
        end
416
        # Compute GLM term for psi (multiply by 2 because the non-conservative flux is
417
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
418
        f[end] = 2 * v1_plus_ll * psi_avg
419

420
    elseif orientation == 2
421
        # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is
422
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
423
        f[1] = 2 * v1_plus_ll * B2_avg
424
        f[2] = 2 * v2_plus_ll * B2_avg
425
        f[3] = 2 * v3_plus_ll * B2_avg
426

427
        for k in eachcomponent(equations)
428
            # Compute term Lorentz term
429
            f2 = charge_ratio_ll[k] * (-B2_avg * B1_avg)
430
            f3 = charge_ratio_ll[k] *
431
                 (-B2_avg * B2_avg + 0.5f0 * mag_norm_avg + pe_mean)
432
            f4 = charge_ratio_ll[k] * (-B2_avg * B3_avg)
433
            f5 = vk2_plus_ll[k] * pe_mean
434

435
            # Compute multi-ion term (vanishes for NCOMP==1)
436
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
437
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
438
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
439
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
440
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
441
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
442
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
443
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
444
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
445
            f5 += (B1_ll * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) +
446
                   B3_ll * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg))
447

448
            # Compute Godunov-Powell term
449
            f2 += charge_ratio_ll[k] * B1_ll * B2_avg
450
            f3 += charge_ratio_ll[k] * B2_ll * B2_avg
451
            f4 += charge_ratio_ll[k] * B3_ll * B2_avg
452
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
453
                  B2_avg
454

455
            # Compute GLM term for the energy
456
            f5 += v2_plus_ll * psi_ll * psi_avg
457

458
            # Add to the flux vector (multiply by 2 because the non-conservative flux is
459
            # multiplied by 0.5 whenever it's used in the Trixi.jl code)
460
            set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
461
                           equations)
462
        end
463
        # Compute GLM term for psi (multiply by 2 because the non-conservative flux is
464
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
465
        f[end] = 2 * v2_plus_ll * psi_avg
466
    else #if orientation == 3
467
        # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is
468
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
469
        f[1] = 2 * v1_plus_ll * B3_avg
470
        f[2] = 2 * v2_plus_ll * B3_avg
471
        f[3] = 2 * v3_plus_ll * B3_avg
472

473
        for k in eachcomponent(equations)
474
            # Compute term Lorentz term
475
            f2 = charge_ratio_ll[k] * (-B3_avg * B1_avg)
476
            f3 = charge_ratio_ll[k] * (-B3_avg * B2_avg)
477
            f4 = charge_ratio_ll[k] *
478
                 (-B3_avg * B3_avg + 0.5f0 * mag_norm_avg + pe_mean)
479
            f5 = vk3_plus_ll[k] * pe_mean
480

481
            # Compute multi-ion term (vanishes for NCOMP==1)
482
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
483
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
484
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
485
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
486
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
487
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
488
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
489
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
490
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
491
            f5 += (B1_ll * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) +
492
                   B2_ll * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg))
493

494
            # Compute Godunov-Powell term
495
            f2 += charge_ratio_ll[k] * B1_ll * B3_avg
496
            f3 += charge_ratio_ll[k] * B2_ll * B3_avg
497
            f4 += charge_ratio_ll[k] * B3_ll * B3_avg
498
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
499
                  B3_avg
500

501
            # Compute GLM term for the energy
502
            f5 += v3_plus_ll * psi_ll * psi_avg
503

504
            # Add to the flux vector (multiply by 2 because the non-conservative flux is
505
            # multiplied by 0.5 whenever it's used in the Trixi.jl code)
506
            set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
507
                           equations)
508
        end
509
        # Compute GLM term for psi (multiply by 2 because the non-conservative flux is
510
        # multiplied by 0.5 whenever it's used in the Trixi.jl code)
511
        f[end] = 2 * v3_plus_ll * psi_avg
512
    end
513

514
    return SVector(f)
515
end
516

517
@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
518
                                                        normal_direction::AbstractVector,
519
                                                        equations::IdealGlmMhdMultiIonEquations3D)
520
    @unpack charge_to_mass = equations
521
    # Unpack left and right states to get the magnetic field
522
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
523
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
524
    psi_ll = divergence_cleaning_field(u_ll, equations)
525
    psi_rr = divergence_cleaning_field(u_rr, equations)
526
    B_dot_n_ll = B1_ll * normal_direction[1] +
527
                 B2_ll * normal_direction[2] +
528
                 B3_ll * normal_direction[3]
529
    B_dot_n_rr = B1_rr * normal_direction[1] +
530
                 B2_rr * normal_direction[2] +
531
                 B3_rr * normal_direction[3]
532
    B_dot_n_avg = 0.5f0 * (B_dot_n_ll + B_dot_n_rr)
533

534
    # Compute important averages
535
    B1_avg = 0.5f0 * (B1_ll + B1_rr)
536
    B2_avg = 0.5f0 * (B2_ll + B2_rr)
537
    B3_avg = 0.5f0 * (B3_ll + B3_rr)
538
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
539
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
540
    mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
541
    psi_avg = 0.5f0 * (psi_ll + psi_rr)
542

543
    # Mean electron pressure
544
    pe_ll = equations.electron_pressure(u_ll, equations)
545
    pe_rr = equations.electron_pressure(u_rr, equations)
546
    pe_mean = 0.5f0 * (pe_ll + pe_rr)
547

548
    # Compute charge ratio of u_ll
549
    charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)})
550
    total_electron_charge = zero(eltype(u_ll))
551
    for k in eachcomponent(equations)
552
        rho_k = u_ll[3 + (k - 1) * 5 + 1] # Extract densities from conserved variable vector
553
        charge_ratio_ll[k] = rho_k * charge_to_mass[k]
554
        total_electron_charge += charge_ratio_ll[k]
555
    end
556
    charge_ratio_ll ./= total_electron_charge
557

558
    # Compute auxiliary variables
559
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
560
                                                                                                           equations)
561
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
562
                                                                                                           equations)
563
    v_plus_dot_n_ll = (v1_plus_ll * normal_direction[1] +
564
                       v2_plus_ll * normal_direction[2] +
565
                       v3_plus_ll * normal_direction[3])
566
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
567

568
    # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is
569
    # multiplied by 0.5 whenever it's used in the Trixi code)
570
    f[1] = 2 * v1_plus_ll * B_dot_n_avg
571
    f[2] = 2 * v2_plus_ll * B_dot_n_avg
572
    f[3] = 2 * v3_plus_ll * B_dot_n_avg
573

574
    for k in eachcomponent(equations)
575
        # Compute term Lorentz term
576
        f2 = charge_ratio_ll[k] *
577
             ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[1] -
578
              B_dot_n_avg * B1_avg)
579
        f3 = charge_ratio_ll[k] *
580
             ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[2] -
581
              B_dot_n_avg * B2_avg)
582
        f4 = charge_ratio_ll[k] *
583
             ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[3] -
584
              B_dot_n_avg * B3_avg)
585
        f5 = (vk1_plus_ll[k] * normal_direction[1] +
586
              vk2_plus_ll[k] * normal_direction[2] +
587
              vk3_plus_ll[k] * normal_direction[3]) * pe_mean
588

589
        # Compute multi-ion term (vanishes for NCOMP==1)
590
        vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
591
        vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
592
        vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
593
        vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
594
        vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
595
        vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
596
        vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
597
        vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
598
        vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
599
        f5 += ((B2_ll * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) +
600
                B3_ll * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) *
601
               normal_direction[1] +
602
               (B1_ll * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) +
603
                B3_ll * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg)) *
604
               normal_direction[2] +
605
               (B1_ll * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) +
606
                B2_ll * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg)) *
607
               normal_direction[3])
608

609
        # Compute Godunov-Powell term
610
        f2 += charge_ratio_ll[k] * B1_ll * B_dot_n_avg
611
        f3 += charge_ratio_ll[k] * B2_ll * B_dot_n_avg
612
        f4 += charge_ratio_ll[k] * B3_ll * B_dot_n_avg
613
        f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
614
              B_dot_n_avg
615

616
        # Compute GLM term for the energy
617
        f5 += v_plus_dot_n_ll * psi_ll * psi_avg
618

619
        # Add to the flux vector (multiply by 2 because the non-conservative flux is
620
        # multiplied by 0.5 whenever it's used in the Trixi code)
621
        set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5,
622
                       equations)
623
    end
624
    # Compute GLM term for psi (multiply by 2 because the non-conservative flux is
625
    # multiplied by 0.5 whenever it's used in the Trixi code)
626
    f[end] = 2 * v_plus_dot_n_ll * psi_avg
627

628
    return SVector(f)
629
end
630

631
"""
632
    flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
633
                                 equations::IdealGlmMhdMultiIonEquations3D)
634

635
Central non-conservative two-point "flux", where the symmetric parts are computed with standard averages.
636
The use of this term together with [`flux_central`](@ref)
637
with [`VolumeIntegralFluxDifferencing`](@ref) yields a "standard"
638
(weak-form) DGSEM discretization of the multi-ion GLM-MHD system. This flux can also be used to construct a
639
standard local Lax-Friedrichs flux using `surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)`.
640

641
!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi.jl"
642
    The central non-conservative fluxes implemented in this function are written as the product
643
    of local and symmetric parts, where the symmetric part is a standard average. These fluxes
644
    are meant to be used in the same way as the conservative fluxes (i.e., flux + flux_noncons
645
    in both volume and surface integrals). In this routine, the fluxes are multiplied by 2 because
646
    the non-conservative fluxes are always multiplied by 0.5 whenever they are used in the Trixi.jl code.
647

648
The term is composed of four individual non-conservative terms:
649
1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and
650
   is evaluated as a function of the charge-averaged velocity.
651
2. The Lorentz-force term, which becomes a conservative term in the limit of one ion species for vanishing
652
   electron pressure gradients.
653
3. The "multi-ion" term, which vanishes in the limit of one ion species.
654
4. The GLM term, which is needed for Galilean invariance.
655
"""
656
@inline function flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
657
                                              equations::IdealGlmMhdMultiIonEquations3D)
658
    @unpack charge_to_mass = equations
659
    # Unpack left and right states to get the magnetic field
660
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
661
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
662
    psi_ll = divergence_cleaning_field(u_ll, equations)
663
    psi_rr = divergence_cleaning_field(u_rr, equations)
664

665
    # Compute important averages
666
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
667
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
668

669
    # Electron pressure
670
    pe_ll = equations.electron_pressure(u_ll, equations)
671
    pe_rr = equations.electron_pressure(u_rr, equations)
672

673
    # Compute charge ratio of u_ll
674
    charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)})
675
    total_electron_charge = zero(real(equations))
676
    for k in eachcomponent(equations)
677
        rho_k = u_ll[3 + (k - 1) * 5 + 1]
678
        charge_ratio_ll[k] = rho_k * charge_to_mass[k]
679
        total_electron_charge += charge_ratio_ll[k]
680
    end
681
    charge_ratio_ll ./= total_electron_charge
682

683
    # Compute auxiliary variables
684
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
685
                                                                                                           equations)
686
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
687
                                                                                                           equations)
688

689
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
690

691
    if orientation == 1
692
        # Entries of Godunov-Powell term for induction equation
693
        f[1] = v1_plus_ll * (B1_ll + B1_rr)
694
        f[2] = v2_plus_ll * (B1_ll + B1_rr)
695
        f[3] = v3_plus_ll * (B1_ll + B1_rr)
696
        for k in eachcomponent(equations)
697
            # Compute Lorentz term
698
            f2 = charge_ratio_ll[k] * ((0.5f0 * mag_norm_ll - B1_ll * B1_ll + pe_ll) +
699
                  (0.5f0 * mag_norm_rr - B1_rr * B1_rr + pe_rr))
700
            f3 = charge_ratio_ll[k] * ((-B1_ll * B2_ll) + (-B1_rr * B2_rr))
701
            f4 = charge_ratio_ll[k] * ((-B1_ll * B3_ll) + (-B1_rr * B3_rr))
702
            f5 = vk1_plus_ll[k] * (pe_ll + pe_rr)
703

704
            # Compute multi-ion term, which vanishes for NCOMP==1
705
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
706
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
707
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
708
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
709
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
710
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
711
            f5 += (B2_ll * ((vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) +
712
                    (vk1_minus_rr * B2_rr - vk2_minus_rr * B1_rr)) +
713
                   B3_ll * ((vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll) +
714
                    (vk1_minus_rr * B3_rr - vk3_minus_rr * B1_rr)))
715

716
            # Compute Godunov-Powell term
717
            f2 += charge_ratio_ll[k] * B1_ll * (B1_ll + B1_rr)
718
            f3 += charge_ratio_ll[k] * B2_ll * (B1_ll + B1_rr)
719
            f4 += charge_ratio_ll[k] * B3_ll * (B1_ll + B1_rr)
720
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
721
                  (B1_ll + B1_rr)
722

723
            # Compute GLM term for the energy
724
            f5 += v1_plus_ll * psi_ll * (psi_ll + psi_rr)
725

726
            # Append to the flux vector
727
            set_component!(f, k, 0, f2, f3, f4, f5, equations)
728
        end
729
        # Compute GLM term for psi
730
        f[end] = v1_plus_ll * (psi_ll + psi_rr)
731

732
    elseif orientation == 2
733
        # Entries of Godunov-Powell term for induction equation
734
        f[1] = v1_plus_ll * (B2_ll + B2_rr)
735
        f[2] = v2_plus_ll * (B2_ll + B2_rr)
736
        f[3] = v3_plus_ll * (B2_ll + B2_rr)
737

738
        for k in eachcomponent(equations)
739
            # Compute Lorentz term
740
            f2 = charge_ratio_ll[k] * ((-B2_ll * B1_ll) + (-B2_rr * B1_rr))
741
            f3 = charge_ratio_ll[k] * ((-B2_ll * B2_ll + 0.5f0 * mag_norm_ll + pe_ll) +
742
                  (-B2_rr * B2_rr + 0.5f0 * mag_norm_rr + pe_rr))
743
            f4 = charge_ratio_ll[k] * ((-B2_ll * B3_ll) + (-B2_rr * B3_rr))
744
            f5 = vk2_plus_ll[k] * (pe_ll + pe_rr)
745

746
            # Compute multi-ion term (vanishes for NCOMP==1)
747
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
748
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
749
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
750
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
751
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
752
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
753
            f5 += (B1_ll * ((vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) +
754
                    (vk2_minus_rr * B1_rr - vk1_minus_rr * B2_rr)) +
755
                   B3_ll * ((vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll) +
756
                    (vk2_minus_rr * B3_rr - vk3_minus_rr * B2_rr)))
757

758
            # Compute Godunov-Powell term
759
            f2 += charge_ratio_ll[k] * B1_ll * (B2_ll + B2_rr)
760
            f3 += charge_ratio_ll[k] * B2_ll * (B2_ll + B2_rr)
761
            f4 += charge_ratio_ll[k] * B3_ll * (B2_ll + B2_rr)
762
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
763
                  (B2_ll + B2_rr)
764

765
            # Compute GLM term for the energy
766
            f5 += v2_plus_ll * psi_ll * (psi_ll + psi_rr)
767

768
            # Append to the flux vector
769
            set_component!(f, k, 0, f2, f3, f4, f5, equations)
770
        end
771
        # Compute GLM term for psi
772
        f[end] = v2_plus_ll * (psi_ll + psi_rr)
773
    else #if orientation == 3
774
        # Entries of Godunov-Powell term for induction equation
775
        f[1] = v1_plus_ll * (B3_ll + B3_rr)
776
        f[2] = v2_plus_ll * (B3_ll + B3_rr)
777
        f[3] = v3_plus_ll * (B3_ll + B3_rr)
778

779
        for k in eachcomponent(equations)
780
            # Compute Lorentz term
781
            f2 = charge_ratio_ll[k] * ((-B3_ll * B1_ll) + (-B3_rr * B1_rr))
782
            f3 = charge_ratio_ll[k] * ((-B3_ll * B2_ll) + (-B3_rr * B2_rr))
783
            f4 = charge_ratio_ll[k] * ((-B3_ll * B3_ll + 0.5f0 * mag_norm_ll + pe_ll) +
784
                  (-B3_rr * B3_rr + 0.5f0 * mag_norm_rr + pe_rr))
785
            f5 = vk3_plus_ll[k] * (pe_ll + pe_rr)
786

787
            # Compute multi-ion term (vanishes for NCOMP==1)
788
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
789
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
790
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
791
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
792
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
793
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
794
            f5 += (B1_ll * ((vk3_minus_ll * B1_ll - vk1_minus_ll * B3_ll) +
795
                    (vk3_minus_rr * B1_rr - vk1_minus_rr * B3_rr)) +
796
                   B2_ll * ((vk3_minus_ll * B2_ll - vk2_minus_ll * B3_ll) +
797
                    (vk3_minus_rr * B2_rr - vk2_minus_rr * B3_rr)))
798

799
            # Compute Godunov-Powell term
800
            f2 += charge_ratio_ll[k] * B1_ll * (B3_ll + B3_rr)
801
            f3 += charge_ratio_ll[k] * B2_ll * (B3_ll + B3_rr)
802
            f4 += charge_ratio_ll[k] * B3_ll * (B3_ll + B3_rr)
803
            f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) *
804
                  (B3_ll + B3_rr)
805

806
            # Compute GLM term for the energy
807
            f5 += v3_plus_ll * psi_ll * (psi_ll + psi_rr)
808

809
            # Append to the flux vector
810
            set_component!(f, k, 0, f2, f3, f4, f5, equations)
811
        end
812
        # Compute GLM term for psi
813
        f[end] = v3_plus_ll * (psi_ll + psi_rr)
814
    end
815

816
    return SVector(f)
817
end
818

819
"""
820
    flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations3D)
821

822
Entropy conserving two-point flux for the multi-ion GLM-MHD equations from
823
- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization
824
  of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics.
825
  [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655).
826

827
This flux (together with the MHD non-conservative term) is consistent in the case of one ion species with the flux of:
828
- Derigs et al. (2018). Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field
829
  divergence diminishing ideal magnetohydrodynamics equations for multi-ion
830
  [DOI: 10.1016/j.jcp.2018.03.002](https://doi.org/10.1016/j.jcp.2018.03.002)
831
"""
832
function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer,
833
                                equations::IdealGlmMhdMultiIonEquations3D)
834
    @unpack gammas = equations
835
    # Unpack left and right states to get the magnetic field
836
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
837
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
838
    psi_ll = divergence_cleaning_field(u_ll, equations)
839
    psi_rr = divergence_cleaning_field(u_rr, equations)
840

841
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
842
                                                                                                           equations)
843
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
844
                                                                                                           equations)
845

846
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
847

848
    # Compute averages for global variables
849
    v1_plus_avg = 0.5f0 * (v1_plus_ll + v1_plus_rr)
850
    v2_plus_avg = 0.5f0 * (v2_plus_ll + v2_plus_rr)
851
    v3_plus_avg = 0.5f0 * (v3_plus_ll + v3_plus_rr)
852
    B1_avg = 0.5f0 * (B1_ll + B1_rr)
853
    B2_avg = 0.5f0 * (B2_ll + B2_rr)
854
    B3_avg = 0.5f0 * (B3_ll + B3_rr)
855
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
856
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
857
    mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
858
    psi_avg = 0.5f0 * (psi_ll + psi_rr)
859

860
    if orientation == 1
861
        psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr)
862

863
        # Magnetic field components from f^MHD
864
        f6 = equations.c_h * psi_avg
865
        f7 = v1_plus_avg * B2_avg - v2_plus_avg * B1_avg
866
        f8 = v1_plus_avg * B3_avg - v3_plus_avg * B1_avg
867
        f9 = equations.c_h * B1_avg
868

869
        # Start building the flux
870
        f[1] = f6
871
        f[2] = f7
872
        f[3] = f8
873
        f[end] = f9
874

875
        # Iterate over all components
876
        for k in eachcomponent(equations)
877
            # Unpack left and right states
878
            rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll = get_component(k,
879
                                                                                    u_ll,
880
                                                                                    equations)
881
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr = get_component(k,
882
                                                                                    u_rr,
883
                                                                                    equations)
884

885
            rho_inv_ll = 1 / rho_ll
886
            v1_ll = rho_v1_ll * rho_inv_ll
887
            v2_ll = rho_v2_ll * rho_inv_ll
888
            v3_ll = rho_v3_ll * rho_inv_ll
889
            rho_inv_rr = 1 / rho_rr
890
            v1_rr = rho_v1_rr * rho_inv_rr
891
            v2_rr = rho_v2_rr * rho_inv_rr
892
            v3_rr = rho_v3_rr * rho_inv_rr
893
            vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
894
            vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
895

896
            p_ll = (gammas[k] - 1) *
897
                   (rho_e_total_ll - 0.5f0 * rho_ll * vel_norm_ll -
898
                    0.5f0 * mag_norm_ll -
899
                    0.5f0 * psi_ll^2)
900
            p_rr = (gammas[k] - 1) *
901
                   (rho_e_total_rr - 0.5f0 * rho_rr * vel_norm_rr -
902
                    0.5f0 * mag_norm_rr -
903
                    0.5f0 * psi_rr^2)
904
            beta_ll = 0.5f0 * rho_ll / p_ll
905
            beta_rr = 0.5f0 * rho_rr / p_rr
906
            # for convenience store vk_plus⋅B
907
            vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
908
                             vk3_plus_ll[k] * B3_ll
909
            vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
910
                             vk3_plus_rr[k] * B3_rr
911

912
            # Compute the necessary mean values needed for either direction
913
            rho_avg = 0.5f0 * (rho_ll + rho_rr)
914
            rho_mean = ln_mean(rho_ll, rho_rr)
915
            beta_mean = ln_mean(beta_ll, beta_rr)
916
            beta_avg = 0.5f0 * (beta_ll + beta_rr)
917
            p_mean = 0.5f0 * rho_avg / beta_avg
918
            v1_avg = 0.5f0 * (v1_ll + v1_rr)
919
            v2_avg = 0.5f0 * (v2_ll + v2_rr)
920
            v3_avg = 0.5f0 * (v3_ll + v3_rr)
921
            vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
922
            vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)
923
            vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k])
924
            vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k])
925
            vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k])
926
            # v_minus
927
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
928
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
929
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
930
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
931
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
932
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
933
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
934
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
935
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
936

937
            # Fill the fluxes for the mass and momentum equations
938
            f1 = rho_mean * v1_avg
939
            f2 = f1 * v1_avg + p_mean
940
            f3 = f1 * v2_avg
941
            f4 = f1 * v3_avg
942

943
            # total energy flux is complicated and involves the previous eight components
944
            v1_plus_mag_avg = 0.5f0 * (vk1_plus_ll[k] * mag_norm_ll +
945
                               vk1_plus_rr[k] * mag_norm_rr)
946
            # Euler part
947
            f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
948
                 f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
949
            # MHD part
950
            f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v1_plus_mag_avg +
951
                   B1_avg * vel_dot_mag_avg                                               # Same terms as in Derigs (but with v_plus)
952
                   + f9 * psi_avg - equations.c_h * psi_B1_avg # GLM term
953
                   +
954
                   0.5f0 * vk1_plus_avg * mag_norm_avg -
955
                   vk1_plus_avg * B1_avg * B1_avg - vk2_plus_avg * B1_avg * B2_avg -
956
                   vk3_plus_avg * B1_avg * B3_avg   # Additional terms related to the Lorentz non-conservative term (momentum eqs)
957
                   -
958
                   B2_avg * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) -
959
                   B3_avg * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg))             # Terms related to the multi-ion non-conservative term (induction equation!)
960

961
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
962
        end
963
    elseif orientation == 2
964
        psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr)
965

966
        # Magnetic field components from f^MHD
967
        f6 = v2_plus_avg * B1_avg - v1_plus_avg * B2_avg
968
        f7 = equations.c_h * psi_avg
969
        f8 = v2_plus_avg * B3_avg - v3_plus_avg * B2_avg
970
        f9 = equations.c_h * B2_avg
971

972
        # Start building the flux
973
        f[1] = f6
974
        f[2] = f7
975
        f[3] = f8
976
        f[end] = f9
977

978
        # Iterate over all components
979
        for k in eachcomponent(equations)
980
            # Unpack left and right states
981
            rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll = get_component(k,
982
                                                                                    u_ll,
983
                                                                                    equations)
984
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr = get_component(k,
985
                                                                                    u_rr,
986
                                                                                    equations)
987

988
            rho_inv_ll = 1 / rho_ll
989
            v1_ll = rho_v1_ll * rho_inv_ll
990
            v2_ll = rho_v2_ll * rho_inv_ll
991
            v3_ll = rho_v3_ll * rho_inv_ll
992
            rho_inv_rr = 1 / rho_rr
993
            v1_rr = rho_v1_rr * rho_inv_rr
994
            v2_rr = rho_v2_rr * rho_inv_rr
995
            v3_rr = rho_v3_rr * rho_inv_rr
996
            vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
997
            vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
998

999
            p_ll = (gammas[k] - 1) *
1000
                   (rho_e_total_ll - 0.5f0 * rho_ll * vel_norm_ll -
1001
                    0.5f0 * mag_norm_ll -
1002
                    0.5f0 * psi_ll^2)
1003
            p_rr = (gammas[k] - 1) *
1004
                   (rho_e_total_rr - 0.5f0 * rho_rr * vel_norm_rr -
1005
                    0.5f0 * mag_norm_rr -
1006
                    0.5f0 * psi_rr^2)
1007
            beta_ll = 0.5f0 * rho_ll / p_ll
1008
            beta_rr = 0.5f0 * rho_rr / p_rr
1009
            # for convenience store vk_plus⋅B
1010
            vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1011
                             vk3_plus_ll[k] * B3_ll
1012
            vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1013
                             vk3_plus_rr[k] * B3_rr
1014

1015
            # Compute the necessary mean values needed for either direction
1016
            rho_avg = 0.5f0 * (rho_ll + rho_rr)
1017
            rho_mean = ln_mean(rho_ll, rho_rr)
1018
            beta_mean = ln_mean(beta_ll, beta_rr)
1019
            beta_avg = 0.5f0 * (beta_ll + beta_rr)
1020
            p_mean = 0.5f0 * rho_avg / beta_avg
1021
            v1_avg = 0.5f0 * (v1_ll + v1_rr)
1022
            v2_avg = 0.5f0 * (v2_ll + v2_rr)
1023
            v3_avg = 0.5f0 * (v3_ll + v3_rr)
1024
            vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
1025
            vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)
1026
            vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k])
1027
            vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k])
1028
            vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k])
1029
            # v_minus
1030
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
1031
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
1032
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
1033
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
1034
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
1035
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
1036
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
1037
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
1038
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
1039

1040
            # Fill the fluxes for the mass and momentum equations
1041
            f1 = rho_mean * v2_avg
1042
            f2 = f1 * v1_avg
1043
            f3 = f1 * v2_avg + p_mean
1044
            f4 = f1 * v3_avg
1045

1046
            # total energy flux is complicated and involves the previous eight components
1047
            v2_plus_mag_avg = 0.5f0 * (vk2_plus_ll[k] * mag_norm_ll +
1048
                               vk2_plus_rr[k] * mag_norm_rr)
1049
            # Euler part
1050
            f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
1051
                 f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
1052
            # MHD part
1053
            f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v2_plus_mag_avg +
1054
                   B2_avg * vel_dot_mag_avg                                               # Same terms as in Derigs (but with v_plus)
1055
                   + f9 * psi_avg - equations.c_h * psi_B2_avg # GLM term
1056
                   +
1057
                   0.5f0 * vk2_plus_avg * mag_norm_avg -
1058
                   vk1_plus_avg * B2_avg * B1_avg - vk2_plus_avg * B2_avg * B2_avg -
1059
                   vk3_plus_avg * B2_avg * B3_avg   # Additional terms related to the Lorentz non-conservative term (momentum eqs)
1060
                   -
1061
                   B1_avg * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) -
1062
                   B3_avg * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg))             # Terms related to the multi-ion non-conservative term (induction equation!)
1063

1064
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
1065
        end
1066
    else #if orientation == 3
1067
        psi_B3_avg = 0.5f0 * (B3_ll * psi_ll + B3_rr * psi_rr)
1068

1069
        # Magnetic field components from f^MHD
1070
        f6 = v3_plus_avg * B1_avg - v1_plus_avg * B3_avg
1071
        f7 = v3_plus_avg * B2_avg - v2_plus_avg * B3_avg
1072
        f8 = equations.c_h * psi_avg
1073
        f9 = equations.c_h * B3_avg
1074

1075
        # Start building the flux
1076
        f[1] = f6
1077
        f[2] = f7
1078
        f[3] = f8
1079
        f[end] = f9
1080

1081
        # Iterate over all components
1082
        for k in eachcomponent(equations)
1083
            # Unpack left and right states
1084
            rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll = get_component(k,
1085
                                                                                    u_ll,
1086
                                                                                    equations)
1087
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr = get_component(k,
1088
                                                                                    u_rr,
1089
                                                                                    equations)
1090

1091
            rho_inv_ll = 1 / rho_ll
1092
            v1_ll = rho_v1_ll * rho_inv_ll
1093
            v2_ll = rho_v2_ll * rho_inv_ll
1094
            v3_ll = rho_v3_ll * rho_inv_ll
1095
            rho_inv_rr = 1 / rho_rr
1096
            v1_rr = rho_v1_rr * rho_inv_rr
1097
            v2_rr = rho_v2_rr * rho_inv_rr
1098
            v3_rr = rho_v3_rr * rho_inv_rr
1099
            vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
1100
            vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
1101

1102
            p_ll = (gammas[k] - 1) *
1103
                   (rho_e_total_ll - 0.5f0 * rho_ll * vel_norm_ll -
1104
                    0.5f0 * mag_norm_ll -
1105
                    0.5f0 * psi_ll^2)
1106
            p_rr = (gammas[k] - 1) *
1107
                   (rho_e_total_rr - 0.5f0 * rho_rr * vel_norm_rr -
1108
                    0.5f0 * mag_norm_rr -
1109
                    0.5f0 * psi_rr^2)
1110
            beta_ll = 0.5f0 * rho_ll / p_ll
1111
            beta_rr = 0.5f0 * rho_rr / p_rr
1112
            # for convenience store vk_plus⋅B
1113
            vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1114
                             vk3_plus_ll[k] * B3_ll
1115
            vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1116
                             vk3_plus_rr[k] * B3_rr
1117

1118
            # Compute the necessary mean values needed for either direction
1119
            rho_avg = 0.5f0 * (rho_ll + rho_rr)
1120
            rho_mean = ln_mean(rho_ll, rho_rr)
1121
            beta_mean = ln_mean(beta_ll, beta_rr)
1122
            beta_avg = 0.5f0 * (beta_ll + beta_rr)
1123
            p_mean = 0.5f0 * rho_avg / beta_avg
1124
            v1_avg = 0.5f0 * (v1_ll + v1_rr)
1125
            v2_avg = 0.5f0 * (v2_ll + v2_rr)
1126
            v3_avg = 0.5f0 * (v3_ll + v3_rr)
1127
            vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
1128
            vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)
1129
            vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k])
1130
            vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k])
1131
            vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k])
1132
            # v_minus
1133
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
1134
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
1135
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
1136
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
1137
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
1138
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
1139
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
1140
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
1141
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
1142

1143
            # Fill the fluxes for the mass and momentum equations
1144
            f1 = rho_mean * v3_avg
1145
            f2 = f1 * v1_avg
1146
            f3 = f1 * v2_avg
1147
            f4 = f1 * v3_avg + p_mean
1148

1149
            # total energy flux is complicated and involves the previous eight components
1150
            v3_plus_mag_avg = 0.5f0 * (vk3_plus_ll[k] * mag_norm_ll +
1151
                               vk3_plus_rr[k] * mag_norm_rr)
1152
            # Euler part
1153
            f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
1154
                 f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
1155
            # MHD part
1156
            f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v3_plus_mag_avg +
1157
                   B3_avg * vel_dot_mag_avg                                               # Same terms as in Derigs (but with v_plus)
1158
                   + f9 * psi_avg - equations.c_h * psi_B3_avg # GLM term
1159
                   +
1160
                   0.5f0 * vk3_plus_avg * mag_norm_avg -
1161
                   vk1_plus_avg * B3_avg * B1_avg - vk2_plus_avg * B3_avg * B2_avg -
1162
                   vk3_plus_avg * B3_avg * B3_avg   # Additional terms related to the Lorentz non-conservative term (momentum eqs)
1163
                   -
1164
                   B1_avg * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) -
1165
                   B2_avg * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg))             # Terms related to the multi-ion non-conservative term (induction equation!)
1166

1167
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
1168
        end
1169
    end
1170

1171
    return SVector(f)
1172
end
1173

1174
function flux_ruedaramirez_etal(u_ll, u_rr, normal_direction::AbstractVector,
1175
                                equations::IdealGlmMhdMultiIonEquations3D)
1176
    @unpack gammas = equations
1177
    # Unpack left and right states to get the magnetic field
1178
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
1179
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
1180
    psi_ll = divergence_cleaning_field(u_ll, equations)
1181
    psi_rr = divergence_cleaning_field(u_rr, equations)
1182
    B_dot_n_ll = B1_ll * normal_direction[1] +
1183
                 B2_ll * normal_direction[2] +
1184
                 B3_ll * normal_direction[3]
1185
    B_dot_n_rr = B1_rr * normal_direction[1] +
1186
                 B2_rr * normal_direction[2] +
1187
                 B3_rr * normal_direction[3]
1188
    B_dot_n_avg = 0.5f0 * (B_dot_n_ll + B_dot_n_rr)
1189

1190
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
1191
                                                                                                           equations)
1192
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
1193
                                                                                                           equations)
1194

1195
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
1196

1197
    # Compute averages for global variables
1198
    v1_plus_avg = 0.5f0 * (v1_plus_ll + v1_plus_rr)
1199
    v2_plus_avg = 0.5f0 * (v2_plus_ll + v2_plus_rr)
1200
    v3_plus_avg = 0.5f0 * (v3_plus_ll + v3_plus_rr)
1201
    B1_avg = 0.5f0 * (B1_ll + B1_rr)
1202
    B2_avg = 0.5f0 * (B2_ll + B2_rr)
1203
    B3_avg = 0.5f0 * (B3_ll + B3_rr)
1204
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
1205
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
1206
    mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
1207
    psi_avg = 0.5f0 * (psi_ll + psi_rr)
1208
    # Averages that depend on normal direction
1209
    psi_B_dot_n_avg = 0.5f0 *
1210
                      ((B1_ll * psi_ll + B1_rr * psi_rr) * normal_direction[1] +
1211
                       (B2_ll * psi_ll + B2_rr * psi_rr) * normal_direction[2] +
1212
                       (B3_ll * psi_ll + B3_rr * psi_rr) * normal_direction[3])
1213

1214
    # Magnetic field components from f^MHD (we store them in f6..f9 for consistency with single-fluid MHD)
1215
    f6 = ((equations.c_h * psi_avg) * normal_direction[1] +
1216
          (v2_plus_avg * B1_avg - v1_plus_avg * B2_avg) * normal_direction[2] +
1217
          (v3_plus_avg * B1_avg - v1_plus_avg * B3_avg) * normal_direction[3])
1218
    f7 = ((v1_plus_avg * B2_avg - v2_plus_avg * B1_avg) * normal_direction[1] +
1219
          (equations.c_h * psi_avg) * normal_direction[2] +
1220
          (v3_plus_avg * B2_avg - v2_plus_avg * B3_avg) * normal_direction[3])
1221
    f8 = ((v1_plus_avg * B3_avg - v3_plus_avg * B1_avg) * normal_direction[1] +
1222
          (v2_plus_avg * B3_avg - v3_plus_avg * B2_avg) * normal_direction[2] +
1223
          (equations.c_h * psi_avg) * normal_direction[3])
1224
    f9 = equations.c_h * B_dot_n_avg
1225

1226
    # Start building the flux
1227
    f[1] = f6
1228
    f[2] = f7
1229
    f[3] = f8
1230
    f[end] = f9
1231

1232
    # Iterate over all components
1233
    for k in eachcomponent(equations)
1234
        # Unpack left and right states
1235
        rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_total_ll = get_component(k, u_ll,
1236
                                                                                equations)
1237
        rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_total_rr = get_component(k, u_rr,
1238
                                                                                equations)
1239

1240
        rho_inv_ll = 1 / rho_ll
1241
        v1_ll = rho_v1_ll * rho_inv_ll
1242
        v2_ll = rho_v2_ll * rho_inv_ll
1243
        v3_ll = rho_v3_ll * rho_inv_ll
1244
        rho_inv_rr = 1 / rho_rr
1245
        v1_rr = rho_v1_rr * rho_inv_rr
1246
        v2_rr = rho_v2_rr * rho_inv_rr
1247
        v3_rr = rho_v3_rr * rho_inv_rr
1248
        vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
1249
        vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
1250

1251
        p_ll = (gammas[k] - 1) *
1252
               (rho_e_total_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll -
1253
                0.5f0 * psi_ll^2)
1254
        p_rr = (gammas[k] - 1) *
1255
               (rho_e_total_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr -
1256
                0.5f0 * psi_rr^2)
1257
        beta_ll = 0.5f0 * rho_ll / p_ll
1258
        beta_rr = 0.5f0 * rho_rr / p_rr
1259

1260
        # for convenience store vk_plus⋅B
1261
        vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1262
                         vk3_plus_ll[k] * B3_ll
1263
        vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1264
                         vk3_plus_rr[k] * B3_rr
1265

1266
        # Compute the necessary mean values needed for either direction
1267
        rho_avg = 0.5f0 * (rho_ll + rho_rr)
1268
        rho_mean = ln_mean(rho_ll, rho_rr)
1269
        beta_mean = ln_mean(beta_ll, beta_rr)
1270
        beta_avg = 0.5f0 * (beta_ll + beta_rr)
1271
        p_mean = 0.5f0 * rho_avg / beta_avg
1272
        v1_avg = 0.5f0 * (v1_ll + v1_rr)
1273
        v2_avg = 0.5f0 * (v2_ll + v2_rr)
1274
        v3_avg = 0.5f0 * (v3_ll + v3_rr)
1275
        vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
1276
        vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)
1277
        vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k])
1278
        vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k])
1279
        vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k])
1280
        # v_minus
1281
        vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
1282
        vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
1283
        vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
1284
        vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
1285
        vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
1286
        vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
1287
        vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
1288
        vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
1289
        vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
1290
        # Compute terms that depend on normal direction
1291
        v_dot_n_ll = v1_ll * normal_direction[1] + v2_ll * normal_direction[2] +
1292
                     v3_ll * normal_direction[3]
1293
        v_dot_n_rr = v1_rr * normal_direction[1] + v2_rr * normal_direction[2] +
1294
                     v3_rr * normal_direction[3]
1295
        vk_plus_mag_avg_dot_n = 0.5f0 * (normal_direction[1] *
1296
                                 (vk1_plus_ll[k] * mag_norm_ll +
1297
                                  vk1_plus_rr[k] * mag_norm_rr) +
1298
                                 normal_direction[2] *
1299
                                 (vk2_plus_ll[k] * mag_norm_ll +
1300
                                  vk2_plus_rr[k] * mag_norm_rr) +
1301
                                 normal_direction[3] *
1302
                                 (vk3_plus_ll[k] * mag_norm_ll +
1303
                                  vk3_plus_rr[k] * mag_norm_rr))
1304
        vk_plus_dot_n_avg = vk1_plus_avg * normal_direction[1] +
1305
                            vk2_plus_avg * normal_direction[2] +
1306
                            vk3_plus_avg * normal_direction[3]
1307

1308
        # Fill the fluxes for the mass and momentum equations
1309
        f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
1310
        f2 = f1 * v1_avg + p_mean * normal_direction[1]
1311
        f3 = f1 * v2_avg + p_mean * normal_direction[2]
1312
        f4 = f1 * v3_avg + p_mean * normal_direction[3]
1313

1314
        # total energy flux is complicated and involves the previous eight components
1315
        # Euler part
1316
        f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
1317
             f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
1318
        # MHD part
1319
        f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * vk_plus_mag_avg_dot_n +
1320
               B_dot_n_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus)
1321
               + f9 * psi_avg - equations.c_h * psi_B_dot_n_avg # GLM terms
1322
               +
1323
               0.5f0 * vk_plus_dot_n_avg * mag_norm_avg -
1324
               vk1_plus_avg * B_dot_n_avg * B1_avg -
1325
               vk2_plus_avg * B_dot_n_avg * B2_avg -
1326
               vk3_plus_avg * B_dot_n_avg * B3_avg) # Additional terms related to the Lorentz non-conservative term (momentum eqs)
1327

1328
        # Curl terms related to the multi-ion non-conservative term that depend on vk_minus
1329
        # These terms vanish in the limit of one ion species and come from the induction equation!
1330
        f5 -= (normal_direction[1] *
1331
               (B2_avg * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) +
1332
                B3_avg * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) +
1333
               normal_direction[2] *
1334
               (B1_avg * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) +
1335
                B3_avg * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg)) +
1336
               normal_direction[3] *
1337
               (B1_avg * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) +
1338
                B2_avg * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg)))
1339

1340
        set_component!(f, k, f1, f2, f3, f4, f5, equations)
1341
    end
1342

1343
    return SVector(f)
1344
end
1345

1346
# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
1347
# This routine approximates the maximum wave speed as sum of the maximum ion velocity
1348
# for all species and the maximum magnetosonic speed.
1349
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
1350
                                     equations::IdealGlmMhdMultiIonEquations3D)
1351
    # Calculate fast magnetoacoustic wave speeds
1352
    # left
1353
    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
1354
    # right
1355
    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
1356

1357
    # Calculate velocities
1358
    v_ll = zero(eltype(u_ll))
1359
    v_rr = zero(eltype(u_rr))
1360
    if orientation == 1
1361
        for k in eachcomponent(equations)
1362
            rho, rho_v1, _ = get_component(k, u_ll, equations)
1363
            v_ll = max(v_ll, abs(rho_v1 / rho))
1364
            rho, rho_v1, _ = get_component(k, u_rr, equations)
1365
            v_rr = max(v_rr, abs(rho_v1 / rho))
1366
        end
1367
    elseif orientation == 2
1368
        for k in eachcomponent(equations)
1369
            rho, rho_v1, rho_v2, _ = get_component(k, u_ll, equations)
1370
            v_ll = max(v_ll, abs(rho_v2 / rho))
1371
            rho, rho_v1, rho_v2, _ = get_component(k, u_rr, equations)
1372
            v_rr = max(v_rr, abs(rho_v2 / rho))
1373
        end
1374
    else #if orientation == 3
1375
        for k in eachcomponent(equations)
1376
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1377
            v_ll = max(v_ll, abs(rho_v3 / rho))
1378
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1379
            v_rr = max(v_rr, abs(rho_v3 / rho))
1380
        end
1381
    end
1382

1383
    return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
1384
end
1385

1386
@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
1387
                                     equations::IdealGlmMhdMultiIonEquations3D)
1388
    # Calculate fast magnetoacoustic wave speeds
1389
    # left
1390
    cf_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
1391
    # right
1392
    cf_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
1393

1394
    # Calculate velocities
1395
    v_ll = zero(eltype(u_ll))
1396
    v_rr = zero(eltype(u_rr))
1397
    for k in eachcomponent(equations)
1398
        # Left state
1399
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1400
        inv_rho = 1 / rho
1401
        v_ll = max(v_ll,
1402
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1403
                       rho_v2 * inv_rho * normal_direction[2] +
1404
                       rho_v3 * inv_rho * normal_direction[3]))
1405
        # Right state
1406
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1407
        inv_rho = 1 / rho
1408
        v_rr = max(v_rr,
1409
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1410
                       rho_v2 * inv_rho * normal_direction[2] +
1411
                       rho_v3 * inv_rho * normal_direction[3]))
1412
    end
1413

1414
    return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
1415
end
1416

1417
# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
1418
@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
1419
                               equations::IdealGlmMhdMultiIonEquations3D)
1420
    # Calculate fast magnetoacoustic wave speeds
1421
    # left
1422
    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
1423
    # right
1424
    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
1425

1426
    # Calculate velocities
1427
    v_ll = zero(eltype(u_ll))
1428
    v_rr = zero(eltype(u_rr))
1429
    if orientation == 1
1430
        for k in eachcomponent(equations)
1431
            rho, rho_v1, _ = get_component(k, u_ll, equations)
1432
            v_ll = max(v_ll, abs(rho_v1 / rho))
1433
            rho, rho_v1, _ = get_component(k, u_rr, equations)
1434
            v_rr = max(v_rr, abs(rho_v1 / rho))
1435
        end
1436
    elseif orientation == 2
1437
        for k in eachcomponent(equations)
1438
            rho, rho_v1, rho_v2, _ = get_component(k, u_ll, equations)
1439
            v_ll = max(v_ll, abs(rho_v2 / rho))
1440
            rho, rho_v1, rho_v2, _ = get_component(k, u_rr, equations)
1441
            v_rr = max(v_rr, abs(rho_v2 / rho))
1442
        end
1443
    else #if orientation == 3
1444
        for k in eachcomponent(equations)
1445
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1446
            v_ll = max(v_ll, abs(rho_v3 / rho))
1447
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1448
            v_rr = max(v_rr, abs(rho_v3 / rho))
1449
        end
1450
    end
1451

1452
    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
1453
end
1454

1455
# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
1456
@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
1457
                               equations::IdealGlmMhdMultiIonEquations3D)
1458
    # Calculate fast magnetoacoustic wave speeds
1459
    # left
1460
    cf_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
1461
    # right
1462
    cf_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
1463

1464
    # Calculate velocities
1465
    v_ll = zero(eltype(u_ll))
1466
    v_rr = zero(eltype(u_rr))
1467
    for k in eachcomponent(equations)
1468
        # Left state
1469
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1470
        inv_rho = 1 / rho
1471
        v_ll = max(v_ll,
1472
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1473
                       rho_v2 * inv_rho * normal_direction[2] +
1474
                       rho_v3 * inv_rho * normal_direction[3]))
1475
        # Right state
1476
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1477
        inv_rho = 1 / rho
1478
        v_rr = max(v_rr,
1479
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1480
                       rho_v2 * inv_rho * normal_direction[2] +
1481
                       rho_v3 * inv_rho * normal_direction[3]))
1482
    end
1483

1484
    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
1485
end
1486

1487
@inline function max_abs_speeds(u, equations::IdealGlmMhdMultiIonEquations3D)
1488
    v1 = zero(real(equations))
1489
    v2 = zero(real(equations))
1490
    v3 = zero(real(equations))
1491
    for k in eachcomponent(equations)
1492
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u, equations)
1493
        rho_inv = 1 / rho
1494
        v1 = max(v1, abs(rho_v1 * rho_inv))
1495
        v2 = max(v2, abs(rho_v2 * rho_inv))
1496
        v3 = max(v3, abs(rho_v3 * rho_inv))
1497
    end
1498

1499
    cf_x_direction = calc_fast_wavespeed(u, 1, equations)
1500
    cf_y_direction = calc_fast_wavespeed(u, 2, equations)
1501
    cf_z_direction = calc_fast_wavespeed(u, 3, equations)
1502

1503
    return (abs(v1) + cf_x_direction, abs(v2) + cf_y_direction,
1504
            abs(v3) + cf_z_direction)
1505
end
1506

1507
# Compute the fastest wave speed for ideal multi-ion GLM-MHD equations: c_f, the fast
1508
# magnetoacoustic eigenvalue. This routine computes the fast magnetosonic speed for each ion
1509
# species using the single-fluid MHD expressions and approximates the multi-ion c_f as
1510
# the maximum of these individual magnetosonic speeds.
1511
@inline function calc_fast_wavespeed(cons, orientation::Integer,
1512
                                     equations::IdealGlmMhdMultiIonEquations3D)
1513
    B1, B2, B3 = magnetic_field(cons, equations)
1514
    psi = divergence_cleaning_field(cons, equations)
1515

1516
    c_f = zero(real(equations))
1517
    for k in eachcomponent(equations)
1518
        rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, cons, equations)
1519

1520
        rho_inv = 1 / rho
1521
        v1 = rho_v1 * rho_inv
1522
        v2 = rho_v2 * rho_inv
1523
        v3 = rho_v3 * rho_inv
1524
        gamma = equations.gammas[k]
1525
        p = (gamma - 1) *
1526
            (rho_e_total - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) -
1527
             0.5f0 * (B1^2 + B2^2 + B3^2) -
1528
             0.5f0 * psi^2)
1529
        a_square = gamma * p * rho_inv
1530
        inv_sqrt_rho = 1 / sqrt(rho)
1531

1532
        b1 = B1 * inv_sqrt_rho
1533
        b2 = B2 * inv_sqrt_rho
1534
        b3 = B3 * inv_sqrt_rho
1535
        b_square = b1^2 + b2^2 + b3^2
1536

1537
        if orientation == 1
1538
            c_f = max(c_f,
1539
                      sqrt(0.5f0 * (a_square + b_square) +
1540
                           0.5f0 *
1541
                           sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)))
1542
        elseif orientation == 2
1543
            c_f = max(c_f,
1544
                      sqrt(0.5f0 * (a_square + b_square) +
1545
                           0.5f0 *
1546
                           sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)))
1547
        else #if orientation == 3
1548
            c_f = max(c_f,
1549
                      sqrt(0.5f0 * (a_square + b_square) +
1550
                           0.5f0 *
1551
                           sqrt((a_square + b_square)^2 - 4 * a_square * b3^2)))
1552
        end
1553
    end
1554

1555
    return c_f
1556
end
1557

1558
@inline function calc_fast_wavespeed(cons, normal_direction::AbstractVector,
1559
                                     equations::IdealGlmMhdMultiIonEquations3D)
1560
    B1, B2, B3 = magnetic_field(cons, equations)
1561
    psi = divergence_cleaning_field(cons, equations)
1562

1563
    norm_squared = (normal_direction[1] * normal_direction[1] +
1564
                    normal_direction[2] * normal_direction[2] +
1565
                    normal_direction[3] * normal_direction[3])
1566

1567
    c_f = zero(real(equations))
1568
    for k in eachcomponent(equations)
1569
        rho, rho_v1, rho_v2, rho_v3, rho_e_total = get_component(k, cons, equations)
1570

1571
        rho_inv = 1 / rho
1572
        v1 = rho_v1 * rho_inv
1573
        v2 = rho_v2 * rho_inv
1574
        v3 = rho_v3 * rho_inv
1575
        gamma = equations.gammas[k]
1576
        p = (gamma - 1) *
1577
            (rho_e_total - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) -
1578
             0.5f0 * (B1^2 + B2^2 + B3^2) -
1579
             0.5f0 * psi^2)
1580
        a_square = gamma * p * rho_inv
1581
        inv_sqrt_rho = 1 / sqrt(rho)
1582

1583
        b1 = B1 * inv_sqrt_rho
1584
        b2 = B2 * inv_sqrt_rho
1585
        b3 = B3 * inv_sqrt_rho
1586
        b_square = b1^2 + b2^2 + b3^2
1587
        b_dot_n_squared = (b1 * normal_direction[1] +
1588
                           b2 * normal_direction[2] +
1589
                           b3 * normal_direction[3])^2 / norm_squared
1590

1591
        c_f = max(c_f,
1592
                  sqrt((0.5f0 * (a_square + b_square) +
1593
                        0.5f0 * sqrt((a_square + b_square)^2 -
1594
                             4 * a_square * b_dot_n_squared)) *
1595
                       norm_squared))
1596
    end
1597

1598
    return c_f
1599
end
1600
end # @muladd
1601

1602
Product

Resources

Company
