Skip to content

Commit

Permalink
Added 'GLM' to the name of the 2D multi-ion MHD equations
Browse files Browse the repository at this point in the history
  • Loading branch information
amrueda committed Dec 6, 2024
1 parent 0732746 commit ac54ebd
Show file tree
Hide file tree
Showing 6 changed files with 66 additions and 62 deletions.
4 changes: 2 additions & 2 deletions examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@ using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
equations = IdealMhdMultiIonEquations2D(gammas = (2.0, 2.0),
charge_to_mass = (1.0, 1.0))
equations = IdealGlmMhdMultiIonEquations2D(gammas = (2.0, 2.0),
charge_to_mass = (1.0, 1.0))

initial_condition = initial_condition_weak_blast_wave

Expand Down
4 changes: 2 additions & 2 deletions examples/tree_2d_dgsem/elixir_mhdmultiion_es.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,8 +4,8 @@ using Trixi

###############################################################################
# semidiscretization of the ideal MHD equations
equations = IdealMhdMultiIonEquations2D(gammas = (2.0, 2.0),
charge_to_mass = (1.0, 1.0))
equations = IdealGlmMhdMultiIonEquations2D(gammas = (2.0, 2.0),
charge_to_mass = (1.0, 1.0))

initial_condition = initial_condition_weak_blast_wave

Expand Down
2 changes: 1 addition & 1 deletion src/Trixi.jl
Original file line number Diff line number Diff line change
Expand Up @@ -156,7 +156,7 @@ export AcousticPerturbationEquations2D,
CompressibleEulerEquationsQuasi1D,
IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D,
IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D,
IdealMhdMultiIonEquations1D, IdealMhdMultiIonEquations2D,
IdealMhdMultiIonEquations1D, IdealGlmMhdMultiIonEquations2D,
HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D,
HyperbolicDiffusionEquations3D,
LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D,
Expand Down
4 changes: 2 additions & 2 deletions src/callbacks_step/glm_speed_dg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@
function calc_dt_for_cleaning_speed(cfl::Real, mesh,
equations::Union{AbstractIdealGlmMhdEquations,
AbstractIdealGlmMhdMulticomponentEquations,
IdealMhdMultiIonEquations2D},
IdealGlmMhdMultiIonEquations2D},
dg::DG, cache)
# compute time step for GLM linear advection equation with c_h=1 for the DG discretization on
# Cartesian meshes
Expand All @@ -30,7 +30,7 @@ end
function calc_dt_for_cleaning_speed(cfl::Real, mesh,
equations::Union{AbstractIdealGlmMhdEquations,
AbstractIdealGlmMhdMulticomponentEquations,
IdealMhdMultiIonEquations2D},
IdealGlmMhdMultiIonEquations2D},
dg::DGMulti, cache)
rd = dg.basis
md = mesh.md
Expand Down
2 changes: 1 addition & 1 deletion src/equations/equations.jl
Original file line number Diff line number Diff line change
Expand Up @@ -495,7 +495,7 @@ include("ideal_glm_mhd_multicomponent_2d.jl")
abstract type AbstractIdealMhdMultiIonEquations{NDIMS, NVARS, NCOMP} <:
AbstractEquations{NDIMS, NVARS} end
include("ideal_mhd_multiion_1d.jl")
include("ideal_mhd_multiion_2d.jl")
include("ideal_glm_mhd_multiion_2d.jl")

# Retrieve number of components from equation instance for the multicomponent case
@inline function ncomponents(::AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS,
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -6,25 +6,25 @@
#! format: noindent

@doc raw"""
IdealMhdMultiIonEquations2D
IdealGlmMhdMultiIonEquations2D
The ideal compressible multi-ion MHD equations in two space dimensions.
"""
mutable struct IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real,
ElectronPressure} <:
mutable struct IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real,
ElectronPressure} <:
AbstractIdealMhdMultiIonEquations{2, NVARS, NCOMP}
gammas::SVector{NCOMP, RealT} # Heat capacity ratios
charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios
electron_pressure::ElectronPressure # Function to compute the electron pressure
c_h::RealT # GLM cleaning speed
function IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT,
ElectronPressure}(gammas
::SVector{NCOMP, RealT},
charge_to_mass
::SVector{NCOMP, RealT},
electron_pressure
::ElectronPressure,
c_h::RealT) where
function IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT,
ElectronPressure}(gammas
::SVector{NCOMP, RealT},
charge_to_mass
::SVector{NCOMP, RealT},
electron_pressure
::ElectronPressure,
c_h::RealT) where
{NVARS, NCOMP, RealT <: Real, ElectronPressure}
NCOMP >= 1 ||
throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value"))
Expand All @@ -33,9 +33,9 @@ mutable struct IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real,
end
end

function IdealMhdMultiIonEquations2D(; gammas, charge_to_mass,
electron_pressure = electron_pressure_zero,
initial_c_h = convert(eltype(gammas), NaN))
function IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass,
electron_pressure = electron_pressure_zero,
initial_c_h = convert(eltype(gammas), NaN))
_gammas = promote(gammas...)
_charge_to_mass = promote(charge_to_mass...)
RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass))
Expand All @@ -46,23 +46,24 @@ function IdealMhdMultiIonEquations2D(; gammas, charge_to_mass,
__gammas = SVector(map(RealT, _gammas))
__charge_to_mass = SVector(map(RealT, _charge_to_mass))

return IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT, typeof(electron_pressure)}(__gammas,
__charge_to_mass,
electron_pressure,
initial_c_h)
return IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT,
typeof(electron_pressure)}(__gammas,
__charge_to_mass,
electron_pressure,
initial_c_h)
end

@inline function Base.real(::IdealMhdMultiIonEquations2D{NVARS, NCOMP, RealT}) where {
NVARS,
NCOMP,
RealT
}
@inline function Base.real(::IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT}) where {
NVARS,
NCOMP,
RealT
}
RealT
end

have_nonconservative_terms(::IdealMhdMultiIonEquations2D) = True()
have_nonconservative_terms(::IdealGlmMhdMultiIonEquations2D) = True()

function varnames(::typeof(cons2cons), equations::IdealMhdMultiIonEquations2D)
function varnames(::typeof(cons2cons), equations::IdealGlmMhdMultiIonEquations2D)
cons = ("B1", "B2", "B3")
for i in eachcomponent(equations)
cons = (cons...,
Expand All @@ -74,7 +75,7 @@ function varnames(::typeof(cons2cons), equations::IdealMhdMultiIonEquations2D)
return cons
end

function varnames(::typeof(cons2prim), equations::IdealMhdMultiIonEquations2D)
function varnames(::typeof(cons2prim), equations::IdealGlmMhdMultiIonEquations2D)
prim = ("B1", "B2", "B3")
for i in eachcomponent(equations)
prim = (prim...,
Expand All @@ -86,16 +87,16 @@ function varnames(::typeof(cons2prim), equations::IdealMhdMultiIonEquations2D)
return prim
end

function default_analysis_integrals(::IdealMhdMultiIonEquations2D)
function default_analysis_integrals(::IdealGlmMhdMultiIonEquations2D)
(entropy_timederivative, Val(:l2_divb), Val(:linf_divb))
end

# """
# initial_condition_convergence_test(x, t, equations::IdealMhdMultiIonEquations2D)
# initial_condition_convergence_test(x, t, equations::IdealGlmMhdMultiIonEquations2D)

# An Alfvén wave as smooth initial condition used for convergence tests.
# """
# function initial_condition_convergence_test(x, t, equations::IdealMhdMultiIonEquations2D)
# function initial_condition_convergence_test(x, t, equations::IdealGlmMhdMultiIonEquations2D)
# # smooth Alfvén wave test from Derigs et al. FLASH (2016)
# # domain must be set to [0, 1], γ = 5/3

Expand All @@ -114,14 +115,15 @@ end
# end

"""
initial_condition_weak_blast_wave(x, t, equations::IdealMhdMultiIonEquations2D)
initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdMultiIonEquations2D)
A weak blast wave adapted from
- Sebastian Hennemann, Gregor J. Gassner (2020)
A provably entropy stable subcell shock capturing approach for high order split form DG
[arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
function initial_condition_weak_blast_wave(x, t, equations::IdealMhdMultiIonEquations2D)
function initial_condition_weak_blast_wave(x, t,
equations::IdealGlmMhdMultiIonEquations2D)
# Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
# Same discontinuity in the velocities but with magnetic fields
# Set up polar coordinates
Expand Down Expand Up @@ -158,7 +160,8 @@ end
# TODO: Add initial condition equilibrium

# Calculate 1D flux in for a single point
@inline function flux(u, orientation::Integer, equations::IdealMhdMultiIonEquations2D)
@inline function flux(u, orientation::Integer,
equations::IdealGlmMhdMultiIonEquations2D)
B1, B2, B3 = magnetic_field(u, equations)
psi = divergence_cleaning_field(u, equations)

Expand Down Expand Up @@ -230,7 +233,7 @@ end
"""
Standard source terms of the multi-ion MHD equations
"""
function source_terms_standard(u, x, t, equations::IdealMhdMultiIonEquations2D)
function source_terms_standard(u, x, t, equations::IdealGlmMhdMultiIonEquations2D)
@unpack charge_to_mass = equations
B1, B2, B3 = magnetic_field(u, equations)
v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u,
Expand Down Expand Up @@ -259,10 +262,10 @@ function source_terms_standard(u, x, t, equations::IdealMhdMultiIonEquations2D)
end

"""
electron_pressure_zero(u, equations::IdealMhdMultiIonEquations2D)
electron_pressure_zero(u, equations::IdealGlmMhdMultiIonEquations2D)
Returns the value of zero for the electron pressure. Consistent with the single-fluid MHD equations.
"""
function electron_pressure_zero(u, equations::IdealMhdMultiIonEquations2D)
function electron_pressure_zero(u, equations::IdealGlmMhdMultiIonEquations2D)
return zero(u[1])
end

Expand All @@ -276,7 +279,7 @@ The term is composed of three parts
"""
@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr,
orientation::Integer,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
@unpack charge_to_mass = equations
# Unpack left and right states to get the magnetic field
B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
Expand Down Expand Up @@ -435,7 +438,7 @@ The term is composed of three parts
* The "term 3": Implemented
"""
@inline function flux_nonconservative_central(u_ll, u_rr, orientation::Integer,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
@unpack charge_to_mass = equations
# Unpack left and right states to get the magnetic field
B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
Expand Down Expand Up @@ -545,7 +548,7 @@ The term is composed of three parts
end

"""
flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealMhdMultiIonEquations2D)
flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations2D)
Entropy conserving two-point flux adapted by:
- Rueda-Ramírez et al. (2023)
Expand All @@ -556,7 +559,7 @@ This flux (together with the MHD non-conservative term) is consistent in the cas
[DOI: 10.1016/j.jcp.2018.03.002](https://doi.org/10.1016/j.jcp.2018.03.002)
"""
function flux_ruedaramirez_etal(u_ll, u_rr, orientation::Integer,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
@unpack gammas = equations
# Unpack left and right states to get the magnetic field
B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations)
Expand Down Expand Up @@ -787,7 +790,7 @@ end
!!!ATTENTION: This routine is provisional. TODO: Update with the right max_abs_speed
"""
@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
# Calculate fast magnetoacoustic wave speeds
# left
cf_ll = calc_fast_wavespeed(u_ll, orientation, equations)
Expand Down Expand Up @@ -816,7 +819,7 @@ end
λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr)
end

@inline function max_abs_speeds(u, equations::IdealMhdMultiIonEquations2D)
@inline function max_abs_speeds(u, equations::IdealGlmMhdMultiIonEquations2D)
v1 = zero(real(equations))
v2 = zero(real(equations))
for k in eachcomponent(equations)
Expand All @@ -834,7 +837,7 @@ end
"""
Convert conservative variables to primitive
"""
function cons2prim(u, equations::IdealMhdMultiIonEquations2D)
function cons2prim(u, equations::IdealGlmMhdMultiIonEquations2D)
@unpack gammas = equations
B1, B2, B3 = magnetic_field(u, equations)
psi = divergence_cleaning_field(u, equations)
Expand Down Expand Up @@ -865,7 +868,7 @@ end
"""
Convert conservative variables to entropy
"""
@inline function cons2entropy(u, equations::IdealMhdMultiIonEquations2D)
@inline function cons2entropy(u, equations::IdealGlmMhdMultiIonEquations2D)
@unpack gammas = equations
B1, B2, B3 = magnetic_field(u, equations)
psi = divergence_cleaning_field(u, equations)
Expand Down Expand Up @@ -899,7 +902,7 @@ end
"""
Convert primitive to conservative variables
"""
@inline function prim2cons(prim, equations::IdealMhdMultiIonEquations2D)
@inline function prim2cons(prim, equations::IdealGlmMhdMultiIonEquations2D)
@unpack gammas = equations
B1, B2, B3 = magnetic_field(prim, equations)
psi = divergence_cleaning_field(prim, equations)
Expand Down Expand Up @@ -931,7 +934,7 @@ Compute the fastest wave speed for ideal MHD equations: c_f, the fast magnetoaco
!!! ATTENTION: This routine is provisional.. Change once the fastest wave speed is known!!
"""
@inline function calc_fast_wavespeed(cons, orientation::Integer,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
B1, B2, B3 = magnetic_field(cons, equations)
psi = divergence_cleaning_field(cons, equations)

Expand Down Expand Up @@ -973,7 +976,8 @@ Routine to compute the Charge-averaged velocities:
* v*_plus: Charge-averaged velocity
* vk*_plus: Contribution of each species to the charge-averaged velocity
"""
@inline function charge_averaged_velocities(u, equations::IdealMhdMultiIonEquations2D)
@inline function charge_averaged_velocities(u,
equations::IdealGlmMhdMultiIonEquations2D)
total_electron_charge = zero(real(equations))

vk1_plus = zero(MVector{ncomponents(equations), eltype(u)})
Expand All @@ -982,7 +986,7 @@ Routine to compute the Charge-averaged velocities:

for k in eachcomponent(equations)
rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)

total_electron_charge += rho * equations.charge_to_mass[k]
vk1_plus[k] = rho_v1 * equations.charge_to_mass[k]
Expand All @@ -1003,7 +1007,7 @@ end
"""
Get the flow variables of component k
"""
@inline function get_component(k, u, equations::IdealMhdMultiIonEquations2D)
@inline function get_component(k, u, equations::IdealGlmMhdMultiIonEquations2D)
# The first 3 entries of u contain the magnetic field. The following entries contain the density, momentum (3 entries), and energy of each component.
return SVector(u[3 + (k - 1) * 5 + 1],
u[3 + (k - 1) * 5 + 2],
Expand All @@ -1016,7 +1020,7 @@ end
Set the flow variables of component k
"""
@inline function set_component!(u, k, u1, u2, u3, u4, u5,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
# The first 3 entries of u contain the magnetic field. The following entries contain the density, momentum (3 entries), and energy of each component.
u[3 + (k - 1) * 5 + 1] = u1
u[3 + (k - 1) * 5 + 2] = u2
Expand All @@ -1027,10 +1031,10 @@ Set the flow variables of component k
return u
end

magnetic_field(u, equations::IdealMhdMultiIonEquations2D) = SVector(u[1], u[2], u[3])
divergence_cleaning_field(u, equations::IdealMhdMultiIonEquations2D) = u[end]
magnetic_field(u, equations::IdealGlmMhdMultiIonEquations2D) = SVector(u[1], u[2], u[3])
divergence_cleaning_field(u, equations::IdealGlmMhdMultiIonEquations2D) = u[end]

@inline function density(u, equations::IdealMhdMultiIonEquations2D)
@inline function density(u, equations::IdealGlmMhdMultiIonEquations2D)
rho = zero(real(equations))
for k in eachcomponent(equations)
rho += u[3 + (k - 1) * 5 + 1]
Expand All @@ -1041,7 +1045,7 @@ end
"""
Computes the sum of the densities times the sum of the pressures
"""
@inline function density_pressure(u, equations::IdealMhdMultiIonEquations2D)
@inline function density_pressure(u, equations::IdealGlmMhdMultiIonEquations2D)
B1, B2, B3 = magnetic_field(u, equations)
psi = divergence_cleaning_field(cons, equations)

Expand Down Expand Up @@ -1083,7 +1087,7 @@ DissipationEntropyStable() = DissipationEntropyStable(max_abs_speed_naive)

@inline function (dissipation::DissipationEntropyStable)(u_ll, u_rr,
orientation_or_normal_direction,
equations::IdealMhdMultiIonEquations2D)
equations::IdealGlmMhdMultiIonEquations2D)
@unpack gammas = equations
λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
equations)
Expand Down

0 comments on commit ac54ebd

Please sign in to comment.