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
mutable 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
        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
    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 = 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 - 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 = 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 - 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 = 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 - 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 = 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 - 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_ll = get_component(k, u_ll,
879
                                                                              equations)
880
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
881
                                                                              equations)
882

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

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

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

933
            # Fill the fluxes for the mass and momentum equations
934
            f1 = rho_mean * v1_avg
935
            f2 = f1 * v1_avg + p_mean
936
            f3 = f1 * v2_avg
937
            f4 = f1 * v3_avg
938

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

957
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
958
        end
959
    elseif orientation == 2
960
        psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr)
961

962
        # Magnetic field components from f^MHD
963
        f6 = v2_plus_avg * B1_avg - v1_plus_avg * B2_avg
964
        f7 = equations.c_h * psi_avg
965
        f8 = v2_plus_avg * B3_avg - v3_plus_avg * B2_avg
966
        f9 = equations.c_h * B2_avg
967

968
        # Start building the flux
969
        f[1] = f6
970
        f[2] = f7
971
        f[3] = f8
972
        f[end] = f9
973

974
        # Iterate over all components
975
        for k in eachcomponent(equations)
976
            # Unpack left and right states
977
            rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = get_component(k, u_ll,
978
                                                                              equations)
979
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
980
                                                                              equations)
981

982
            rho_inv_ll = 1 / rho_ll
983
            v1_ll = rho_v1_ll * rho_inv_ll
984
            v2_ll = rho_v2_ll * rho_inv_ll
985
            v3_ll = rho_v3_ll * rho_inv_ll
986
            rho_inv_rr = 1 / rho_rr
987
            v1_rr = rho_v1_rr * rho_inv_rr
988
            v2_rr = rho_v2_rr * rho_inv_rr
989
            v3_rr = rho_v3_rr * rho_inv_rr
990
            vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
991
            vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
992

993
            p_ll = (gammas[k] - 1) *
994
                   (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll -
995
                    0.5f0 * psi_ll^2)
996
            p_rr = (gammas[k] - 1) *
997
                   (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr -
998
                    0.5f0 * psi_rr^2)
999
            beta_ll = 0.5f0 * rho_ll / p_ll
1000
            beta_rr = 0.5f0 * rho_rr / p_rr
1001
            # for convenience store vk_plus⋅B
1002
            vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1003
                             vk3_plus_ll[k] * B3_ll
1004
            vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1005
                             vk3_plus_rr[k] * B3_rr
1006

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

1032
            # Fill the fluxes for the mass and momentum equations
1033
            f1 = rho_mean * v2_avg
1034
            f2 = f1 * v1_avg
1035
            f3 = f1 * v2_avg + p_mean
1036
            f4 = f1 * v3_avg
1037

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

1056
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
1057
        end
1058
    else #if orientation == 3
1059
        psi_B3_avg = 0.5f0 * (B3_ll * psi_ll + B3_rr * psi_rr)
1060

1061
        # Magnetic field components from f^MHD
1062
        f6 = v3_plus_avg * B1_avg - v1_plus_avg * B3_avg
1063
        f7 = v3_plus_avg * B2_avg - v2_plus_avg * B3_avg
1064
        f8 = equations.c_h * psi_avg
1065
        f9 = equations.c_h * B3_avg
1066

1067
        # Start building the flux
1068
        f[1] = f6
1069
        f[2] = f7
1070
        f[3] = f8
1071
        f[end] = f9
1072

1073
        # Iterate over all components
1074
        for k in eachcomponent(equations)
1075
            # Unpack left and right states
1076
            rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = get_component(k, u_ll,
1077
                                                                              equations)
1078
            rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
1079
                                                                              equations)
1080

1081
            rho_inv_ll = 1 / rho_ll
1082
            v1_ll = rho_v1_ll * rho_inv_ll
1083
            v2_ll = rho_v2_ll * rho_inv_ll
1084
            v3_ll = rho_v3_ll * rho_inv_ll
1085
            rho_inv_rr = 1 / rho_rr
1086
            v1_rr = rho_v1_rr * rho_inv_rr
1087
            v2_rr = rho_v2_rr * rho_inv_rr
1088
            v3_rr = rho_v3_rr * rho_inv_rr
1089
            vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
1090
            vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
1091

1092
            p_ll = (gammas[k] - 1) *
1093
                   (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll -
1094
                    0.5f0 * psi_ll^2)
1095
            p_rr = (gammas[k] - 1) *
1096
                   (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr -
1097
                    0.5f0 * psi_rr^2)
1098
            beta_ll = 0.5f0 * rho_ll / p_ll
1099
            beta_rr = 0.5f0 * rho_rr / p_rr
1100
            # for convenience store vk_plus⋅B
1101
            vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1102
                             vk3_plus_ll[k] * B3_ll
1103
            vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1104
                             vk3_plus_rr[k] * B3_rr
1105

1106
            # Compute the necessary mean values needed for either direction
1107
            rho_avg = 0.5f0 * (rho_ll + rho_rr)
1108
            rho_mean = ln_mean(rho_ll, rho_rr)
1109
            beta_mean = ln_mean(beta_ll, beta_rr)
1110
            beta_avg = 0.5f0 * (beta_ll + beta_rr)
1111
            p_mean = 0.5f0 * rho_avg / beta_avg
1112
            v1_avg = 0.5f0 * (v1_ll + v1_rr)
1113
            v2_avg = 0.5f0 * (v2_ll + v2_rr)
1114
            v3_avg = 0.5f0 * (v3_ll + v3_rr)
1115
            vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
1116
            vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)
1117
            vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k])
1118
            vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k])
1119
            vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k])
1120
            # v_minus
1121
            vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k]
1122
            vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k]
1123
            vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k]
1124
            vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k]
1125
            vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k]
1126
            vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k]
1127
            vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr)
1128
            vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr)
1129
            vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr)
1130

1131
            # Fill the fluxes for the mass and momentum equations
1132
            f1 = rho_mean * v3_avg
1133
            f2 = f1 * v1_avg
1134
            f3 = f1 * v2_avg
1135
            f4 = f1 * v3_avg + p_mean
1136

1137
            # total energy flux is complicated and involves the previous eight components
1138
            v3_plus_mag_avg = 0.5f0 * (vk3_plus_ll[k] * mag_norm_ll +
1139
                               vk3_plus_rr[k] * mag_norm_rr)
1140
            # Euler part
1141
            f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
1142
                 f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
1143
            # MHD part
1144
            f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v3_plus_mag_avg +
1145
                   B3_avg * vel_dot_mag_avg                                               # Same terms as in Derigs (but with v_plus)
1146
                   + f9 * psi_avg - equations.c_h * psi_B3_avg # GLM term
1147
                   +
1148
                   0.5f0 * vk3_plus_avg * mag_norm_avg -
1149
                   vk1_plus_avg * B3_avg * B1_avg - vk2_plus_avg * B3_avg * B2_avg -
1150
                   vk3_plus_avg * B3_avg * B3_avg   # Additional terms related to the Lorentz non-conservative term (momentum eqs)
1151
                   -
1152
                   B1_avg * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) -
1153
                   B2_avg * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg))             # Terms related to the multi-ion non-conservative term (induction equation!)
1154

1155
            set_component!(f, k, f1, f2, f3, f4, f5, equations)
1156
        end
1157
    end
1158

1159
    return SVector(f)
1160
end
1161

1162
function flux_ruedaramirez_etal(u_ll, u_rr, normal_direction::AbstractVector,
1163
                                equations::IdealGlmMhdMultiIonEquations3D)
1164
    @unpack gammas = equations
1165
    # Unpack left and right states to get the magnetic field
1166
    B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
1167
    B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations)
1168
    psi_ll = divergence_cleaning_field(u_ll, equations)
1169
    psi_rr = divergence_cleaning_field(u_rr, equations)
1170
    B_dot_n_ll = B1_ll * normal_direction[1] +
1171
                 B2_ll * normal_direction[2] +
1172
                 B3_ll * normal_direction[3]
1173
    B_dot_n_rr = B1_rr * normal_direction[1] +
1174
                 B2_rr * normal_direction[2] +
1175
                 B3_rr * normal_direction[3]
1176
    B_dot_n_avg = 0.5f0 * (B_dot_n_ll + B_dot_n_rr)
1177

1178
    v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll,
1179
                                                                                                           equations)
1180
    v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr,
1181
                                                                                                           equations)
1182

1183
    f = zero(MVector{nvariables(equations), eltype(u_ll)})
1184

1185
    # Compute averages for global variables
1186
    v1_plus_avg = 0.5f0 * (v1_plus_ll + v1_plus_rr)
1187
    v2_plus_avg = 0.5f0 * (v2_plus_ll + v2_plus_rr)
1188
    v3_plus_avg = 0.5f0 * (v3_plus_ll + v3_plus_rr)
1189
    B1_avg = 0.5f0 * (B1_ll + B1_rr)
1190
    B2_avg = 0.5f0 * (B2_ll + B2_rr)
1191
    B3_avg = 0.5f0 * (B3_ll + B3_rr)
1192
    mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2
1193
    mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2
1194
    mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
1195
    psi_avg = 0.5f0 * (psi_ll + psi_rr)
1196
    # Averages that depend on normal direction
1197
    psi_B_dot_n_avg = 0.5f0 *
1198
                      ((B1_ll * psi_ll + B1_rr * psi_rr) * normal_direction[1] +
1199
                       (B2_ll * psi_ll + B2_rr * psi_rr) * normal_direction[2] +
1200
                       (B3_ll * psi_ll + B3_rr * psi_rr) * normal_direction[3])
1201

1202
    # Magnetic field components from f^MHD (we store them in f6..f9 for consistency with single-fluid MHD)
1203
    f6 = ((equations.c_h * psi_avg) * normal_direction[1] +
1204
          (v2_plus_avg * B1_avg - v1_plus_avg * B2_avg) * normal_direction[2] +
1205
          (v3_plus_avg * B1_avg - v1_plus_avg * B3_avg) * normal_direction[3])
1206
    f7 = ((v1_plus_avg * B2_avg - v2_plus_avg * B1_avg) * normal_direction[1] +
1207
          (equations.c_h * psi_avg) * normal_direction[2] +
1208
          (v3_plus_avg * B2_avg - v2_plus_avg * B3_avg) * normal_direction[3])
1209
    f8 = ((v1_plus_avg * B3_avg - v3_plus_avg * B1_avg) * normal_direction[1] +
1210
          (v2_plus_avg * B3_avg - v3_plus_avg * B2_avg) * normal_direction[2] +
1211
          (equations.c_h * psi_avg) * normal_direction[3])
1212
    f9 = equations.c_h * B_dot_n_avg
1213

1214
    # Start building the flux
1215
    f[1] = f6
1216
    f[2] = f7
1217
    f[3] = f8
1218
    f[end] = f9
1219

1220
    # Iterate over all components
1221
    for k in eachcomponent(equations)
1222
        # Unpack left and right states
1223
        rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = get_component(k, u_ll,
1224
                                                                          equations)
1225
        rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr,
1226
                                                                          equations)
1227

1228
        rho_inv_ll = 1 / rho_ll
1229
        v1_ll = rho_v1_ll * rho_inv_ll
1230
        v2_ll = rho_v2_ll * rho_inv_ll
1231
        v3_ll = rho_v3_ll * rho_inv_ll
1232
        rho_inv_rr = 1 / rho_rr
1233
        v1_rr = rho_v1_rr * rho_inv_rr
1234
        v2_rr = rho_v2_rr * rho_inv_rr
1235
        v3_rr = rho_v3_rr * rho_inv_rr
1236
        vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
1237
        vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2
1238

1239
        p_ll = (gammas[k] - 1) *
1240
               (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll -
1241
                0.5f0 * psi_ll^2)
1242
        p_rr = (gammas[k] - 1) *
1243
               (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr -
1244
                0.5f0 * psi_rr^2)
1245
        beta_ll = 0.5f0 * rho_ll / p_ll
1246
        beta_rr = 0.5f0 * rho_rr / p_rr
1247

1248
        # for convenience store vk_plus⋅B
1249
        vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll +
1250
                         vk3_plus_ll[k] * B3_ll
1251
        vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr +
1252
                         vk3_plus_rr[k] * B3_rr
1253

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

1296
        # Fill the fluxes for the mass and momentum equations
1297
        f1 = rho_mean * 0.5f0 * (v_dot_n_ll + v_dot_n_rr)
1298
        f2 = f1 * v1_avg + p_mean * normal_direction[1]
1299
        f3 = f1 * v2_avg + p_mean * normal_direction[2]
1300
        f4 = f1 * v3_avg + p_mean * normal_direction[3]
1301

1302
        # total energy flux is complicated and involves the previous eight components
1303
        # Euler part
1304
        f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) +
1305
             f2 * v1_avg + f3 * v2_avg + f4 * v3_avg
1306
        # MHD part
1307
        f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * vk_plus_mag_avg_dot_n +
1308
               B_dot_n_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus)
1309
               + f9 * psi_avg - equations.c_h * psi_B_dot_n_avg # GLM terms
1310
               +
1311
               0.5f0 * vk_plus_dot_n_avg * mag_norm_avg -
1312
               vk1_plus_avg * B_dot_n_avg * B1_avg -
1313
               vk2_plus_avg * B_dot_n_avg * B2_avg -
1314
               vk3_plus_avg * B_dot_n_avg * B3_avg) # Additional terms related to the Lorentz non-conservative term (momentum eqs)
1315

1316
        # Curl terms related to the multi-ion non-conservative term that depend on vk_minus
1317
        # These terms vanish in the limit of one ion species and come from the induction equation!
1318
        f5 -= (normal_direction[1] *
1319
               (B2_avg * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) +
1320
                B3_avg * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) +
1321
               normal_direction[2] *
1322
               (B1_avg * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) +
1323
                B3_avg * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg)) +
1324
               normal_direction[3] *
1325
               (B1_avg * (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) +
1326
                B2_avg * (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg)))
1327

1328
        set_component!(f, k, f1, f2, f3, f4, f5, equations)
1329
    end
1330

1331
    return SVector(f)
1332
end
1333

1334
# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation
1335
# This routine approximates the maximum wave speed as sum of the maximum ion velocity 
1336
# for all species and the maximum magnetosonic speed.
1337
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
1338
                                     equations::IdealGlmMhdMultiIonEquations3D)
1339
    # Calculate fast magnetoacoustic wave speeds
1340
    # left
1341
    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
1342
    # right
1343
    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
1344

1345
    # Calculate velocities
1346
    v_ll = zero(eltype(u_ll))
1347
    v_rr = zero(eltype(u_rr))
1348
    if orientation == 1
1349
        for k in eachcomponent(equations)
1350
            rho, rho_v1, _ = get_component(k, u_ll, equations)
1351
            v_ll = max(v_ll, abs(rho_v1 / rho))
1352
            rho, rho_v1, _ = get_component(k, u_rr, equations)
1353
            v_rr = max(v_rr, abs(rho_v1 / rho))
1354
        end
1355
    elseif orientation == 2
1356
        for k in eachcomponent(equations)
1357
            rho, rho_v1, rho_v2, _ = get_component(k, u_ll, equations)
1358
            v_ll = max(v_ll, abs(rho_v2 / rho))
1359
            rho, rho_v1, rho_v2, _ = get_component(k, u_rr, equations)
1360
            v_rr = max(v_rr, abs(rho_v2 / rho))
1361
        end
1362
    else #if orientation == 3
1363
        for k in eachcomponent(equations)
1364
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1365
            v_ll = max(v_ll, abs(rho_v3 / rho))
1366
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1367
            v_rr = max(v_rr, abs(rho_v3 / rho))
1368
        end
1369
    end
1370

1371
    return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
1372
end
1373

1374
@inline function max_abs_speed_naive(u_ll, u_rr, normal_direction::AbstractVector,
1375
                                     equations::IdealGlmMhdMultiIonEquations3D)
1376
    # Calculate fast magnetoacoustic wave speeds
1377
    # left
1378
    cf_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
1379
    # right
1380
    cf_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
1381

1382
    # Calculate velocities
1383
    v_ll = zero(eltype(u_ll))
1384
    v_rr = zero(eltype(u_rr))
1385
    for k in eachcomponent(equations)
1386
        # Left state
1387
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1388
        inv_rho = 1 / rho
1389
        v_ll = max(v_ll,
1390
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1391
                       rho_v2 * inv_rho * normal_direction[2] +
1392
                       rho_v3 * inv_rho * normal_direction[3]))
1393
        # Right state
1394
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1395
        inv_rho = 1 / rho
1396
        v_rr = max(v_rr,
1397
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1398
                       rho_v2 * inv_rho * normal_direction[2] +
1399
                       rho_v3 * inv_rho * normal_direction[3]))
1400
    end
1401

1402
    return max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
1403
end
1404

1405
# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
1406
@inline function max_abs_speed(u_ll, u_rr, orientation::Integer,
1407
                               equations::IdealGlmMhdMultiIonEquations3D)
1408
    # Calculate fast magnetoacoustic wave speeds
1409
    # left
1410
    cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
1411
    # right
1412
    cf_rr = calc_fast_wavespeed(u_rr, orientation, equations)
1413

1414
    # Calculate velocities
1415
    v_ll = zero(eltype(u_ll))
1416
    v_rr = zero(eltype(u_rr))
1417
    if orientation == 1
1418
        for k in eachcomponent(equations)
1419
            rho, rho_v1, _ = get_component(k, u_ll, equations)
1420
            v_ll = max(v_ll, abs(rho_v1 / rho))
1421
            rho, rho_v1, _ = get_component(k, u_rr, equations)
1422
            v_rr = max(v_rr, abs(rho_v1 / rho))
1423
        end
1424
    elseif orientation == 2
1425
        for k in eachcomponent(equations)
1426
            rho, rho_v1, rho_v2, _ = get_component(k, u_ll, equations)
1427
            v_ll = max(v_ll, abs(rho_v2 / rho))
1428
            rho, rho_v1, rho_v2, _ = get_component(k, u_rr, equations)
1429
            v_rr = max(v_rr, abs(rho_v2 / rho))
1430
        end
1431
    else #if orientation == 3
1432
        for k in eachcomponent(equations)
1433
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1434
            v_ll = max(v_ll, abs(rho_v3 / rho))
1435
            rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1436
            v_rr = max(v_rr, abs(rho_v3 / rho))
1437
        end
1438
    end
1439

1440
    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
1441
end
1442

1443
# Less "cautious", i.e., less overestimating `λ_max` compared to `max_abs_speed_naive`
1444
@inline function max_abs_speed(u_ll, u_rr, normal_direction::AbstractVector,
1445
                               equations::IdealGlmMhdMultiIonEquations3D)
1446
    # Calculate fast magnetoacoustic wave speeds
1447
    # left
1448
    cf_ll = calc_fast_wavespeed(u_ll, normal_direction, equations)
1449
    # right
1450
    cf_rr = calc_fast_wavespeed(u_rr, normal_direction, equations)
1451

1452
    # Calculate velocities
1453
    v_ll = zero(eltype(u_ll))
1454
    v_rr = zero(eltype(u_rr))
1455
    for k in eachcomponent(equations)
1456
        # Left state
1457
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_ll, equations)
1458
        inv_rho = 1 / rho
1459
        v_ll = max(v_ll,
1460
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1461
                       rho_v2 * inv_rho * normal_direction[2] +
1462
                       rho_v3 * inv_rho * normal_direction[3]))
1463
        # Right state
1464
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u_rr, equations)
1465
        inv_rho = 1 / rho
1466
        v_rr = max(v_rr,
1467
                   abs(rho_v1 * inv_rho * normal_direction[1] +
1468
                       rho_v2 * inv_rho * normal_direction[2] +
1469
                       rho_v3 * inv_rho * normal_direction[3]))
1470
    end
1471

1472
    return max(abs(v_ll) + cf_ll, abs(v_rr) + cf_rr)
1473
end
1474

1475
@inline function max_abs_speeds(u, equations::IdealGlmMhdMultiIonEquations3D)
1476
    v1 = zero(real(equations))
1477
    v2 = zero(real(equations))
1478
    v3 = zero(real(equations))
1479
    for k in eachcomponent(equations)
1480
        rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u, equations)
1481
        rho_inv = 1 / rho
1482
        v1 = max(v1, abs(rho_v1 * rho_inv))
1483
        v2 = max(v2, abs(rho_v2 * rho_inv))
1484
        v3 = max(v3, abs(rho_v3 * rho_inv))
1485
    end
1486

1487
    cf_x_direction = calc_fast_wavespeed(u, 1, equations)
1488
    cf_y_direction = calc_fast_wavespeed(u, 2, equations)
1489
    cf_z_direction = calc_fast_wavespeed(u, 3, equations)
1490

1491
    return (abs(v1) + cf_x_direction, abs(v2) + cf_y_direction,
1492
            abs(v3) + cf_z_direction)
1493
end
1494

1495
# Compute the fastest wave speed for ideal multi-ion GLM-MHD equations: c_f, the fast 
1496
# magnetoacoustic eigenvalue. This routine computes the fast magnetosonic speed for each ion
1497
# species using the single-fluid MHD expressions and approximates the multi-ion c_f as 
1498
# the maximum of these individual magnetosonic speeds.
1499
@inline function calc_fast_wavespeed(cons, orientation::Integer,
1500
                                     equations::IdealGlmMhdMultiIonEquations3D)
1501
    B1, B2, B3 = magnetic_field(cons, equations)
1502
    psi = divergence_cleaning_field(cons, equations)
1503

1504
    c_f = zero(real(equations))
1505
    for k in eachcomponent(equations)
1506
        rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, cons, equations)
1507

1508
        rho_inv = 1 / rho
1509
        v1 = rho_v1 * rho_inv
1510
        v2 = rho_v2 * rho_inv
1511
        v3 = rho_v3 * rho_inv
1512
        gamma = equations.gammas[k]
1513
        p = (gamma - 1) *
1514
            (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2) -
1515
             0.5f0 * psi^2)
1516
        a_square = gamma * p * rho_inv
1517
        inv_sqrt_rho = 1 / sqrt(rho)
1518

1519
        b1 = B1 * inv_sqrt_rho
1520
        b2 = B2 * inv_sqrt_rho
1521
        b3 = B3 * inv_sqrt_rho
1522
        b_square = b1^2 + b2^2 + b3^2
1523

1524
        if orientation == 1
1525
            c_f = max(c_f,
1526
                      sqrt(0.5f0 * (a_square + b_square) +
1527
                           0.5f0 *
1528
                           sqrt((a_square + b_square)^2 - 4 * a_square * b1^2)))
1529
        elseif orientation == 2
1530
            c_f = max(c_f,
1531
                      sqrt(0.5f0 * (a_square + b_square) +
1532
                           0.5f0 *
1533
                           sqrt((a_square + b_square)^2 - 4 * a_square * b2^2)))
1534
        else #if orientation == 3
1535
            c_f = max(c_f,
1536
                      sqrt(0.5f0 * (a_square + b_square) +
1537
                           0.5f0 *
1538
                           sqrt((a_square + b_square)^2 - 4 * a_square * b3^2)))
1539
        end
1540
    end
1541

1542
    return c_f
1543
end
1544

1545
@inline function calc_fast_wavespeed(cons, normal_direction::AbstractVector,
1546
                                     equations::IdealGlmMhdMultiIonEquations3D)
1547
    B1, B2, B3 = magnetic_field(cons, equations)
1548
    psi = divergence_cleaning_field(cons, equations)
1549

1550
    norm_squared = (normal_direction[1] * normal_direction[1] +
1551
                    normal_direction[2] * normal_direction[2] +
1552
                    normal_direction[3] * normal_direction[3])
1553

1554
    c_f = zero(real(equations))
1555
    for k in eachcomponent(equations)
1556
        rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, cons, equations)
1557

1558
        rho_inv = 1 / rho
1559
        v1 = rho_v1 * rho_inv
1560
        v2 = rho_v2 * rho_inv
1561
        v3 = rho_v3 * rho_inv
1562
        gamma = equations.gammas[k]
1563
        p = (gamma - 1) *
1564
            (rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2) -
1565
             0.5f0 * psi^2)
1566
        a_square = gamma * p * rho_inv
1567
        inv_sqrt_rho = 1 / sqrt(rho)
1568

1569
        b1 = B1 * inv_sqrt_rho
1570
        b2 = B2 * inv_sqrt_rho
1571
        b3 = B3 * inv_sqrt_rho
1572
        b_square = b1^2 + b2^2 + b3^2
1573
        b_dot_n_squared = (b1 * normal_direction[1] +
1574
                           b2 * normal_direction[2] +
1575
                           b3 * normal_direction[3])^2 / norm_squared
1576

1577
        c_f = max(c_f,
1578
                  sqrt((0.5f0 * (a_square + b_square) +
1579
                        0.5f0 * sqrt((a_square + b_square)^2 -
1580
                             4 * a_square * b_dot_n_squared)) *
1581
                       norm_squared))
1582
    end
1583

1584
    return c_f
1585
end
1586
end # @muladd
1587

1588
Product

Resources

Company
