diff --git a/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl new file mode 100644 index 00000000000..d9b66496b19 --- /dev/null +++ b/examples/tree_2d_dgsem/elixir_mhdmultiion_ec.jl @@ -0,0 +1,61 @@ + +using OrdinaryDiffEq +using Trixi + +############################################################################### +# semidiscretization of the ideal multi-ion MHD equations +equations = IdealGlmMhdMultiIonEquations2D(gammas = (1.4, 1.667), + charge_to_mass = (1.0, 2.0)) + +initial_condition = initial_condition_weak_blast_wave + +volume_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal) +surface_flux = (flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal) +solver = DGSEM(polydeg = 3, surface_flux = surface_flux, + volume_integral = VolumeIntegralFluxDifferencing(volume_flux)) + +coordinates_min = (-2.0, -2.0) +coordinates_max = (2.0, 2.0) +mesh = TreeMesh(coordinates_min, coordinates_max, + initial_refinement_level = 4, + n_cells_max = 10_000) + +semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, + source_terms = source_terms_lorentz) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 0.4) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 10 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval) +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(dt = 0.1, # interval=100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +cfl = 0.5 + +stepsize_callback = StepsizeCallback(cfl = cfl) + +glm_speed_callback = GlmSpeedCallback(glm_scale = 0.5, cfl = cfl) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_callback, + glm_speed_callback) + +############################################################################### +# run the simulation + +sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false), + dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback + save_everystep = false, callback = callbacks); +summary_callback() # print the timer summary diff --git a/src/Trixi.jl b/src/Trixi.jl index 7d557ddde38..cb50ca2e7df 100644 --- a/src/Trixi.jl +++ b/src/Trixi.jl @@ -156,6 +156,7 @@ export AcousticPerturbationEquations2D, CompressibleEulerEquationsQuasi1D, IdealGlmMhdEquations1D, IdealGlmMhdEquations2D, IdealGlmMhdEquations3D, IdealGlmMhdMulticomponentEquations1D, IdealGlmMhdMulticomponentEquations2D, + IdealGlmMhdMultiIonEquations2D, HyperbolicDiffusionEquations1D, HyperbolicDiffusionEquations2D, HyperbolicDiffusionEquations3D, LinearScalarAdvectionEquation1D, LinearScalarAdvectionEquation2D, @@ -179,6 +180,8 @@ export flux, flux_central, flux_lax_friedrichs, flux_hll, flux_hllc, flux_hlle, flux_godunov, flux_chandrashekar, flux_ranocha, flux_derigs_etal, flux_hindenlang_gassner, flux_nonconservative_powell, flux_nonconservative_powell_local_symmetric, + flux_ruedaramirez_etal, flux_nonconservative_ruedaramirez_etal, + flux_nonconservative_central, flux_kennedy_gruber, flux_shima_etal, flux_ec, flux_fjordholm_etal, flux_nonconservative_fjordholm_etal, flux_wintermeyer_etal, flux_nonconservative_wintermeyer_etal, @@ -212,7 +215,8 @@ export boundary_condition_do_nothing, BoundaryConditionNavierStokesWall, NoSlip, Adiabatic, Isothermal, BoundaryConditionCoupled -export initial_condition_convergence_test, source_terms_convergence_test +export initial_condition_convergence_test, source_terms_convergence_test, + source_terms_lorentz export source_terms_harmonic export initial_condition_poisson_nonperiodic, source_terms_poisson_nonperiodic, boundary_condition_poisson_nonperiodic @@ -225,7 +229,7 @@ export density, pressure, density_pressure, velocity, global_mean_vars, equilibrium_distribution, waterheight_pressure export entropy, energy_total, energy_kinetic, energy_internal, energy_magnetic, cross_helicity, - enstrophy + enstrophy, magnetic_field, divergence_cleaning_field export lake_at_rest_error export ncomponents, eachcomponent diff --git a/src/callbacks_step/glm_speed.jl b/src/callbacks_step/glm_speed.jl index 0bd87c189bc..075e084949a 100644 --- a/src/callbacks_step/glm_speed.jl +++ b/src/callbacks_step/glm_speed.jl @@ -9,7 +9,8 @@ GlmSpeedCallback(; glm_scale=0.5, cfl, semi_indices=()) Update the divergence cleaning wave speed `c_h` according to the time step -computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations. +computed in [`StepsizeCallback`](@ref) for the ideal GLM-MHD equations, the multi-component +GLM-MHD equations, and the multi-ion GLM-MHD equations. The `cfl` number should be set to the same value as for the time step size calculation. The `glm_scale` ensures that the GLM wave speed is lower than the fastest physical waves in the MHD solution and should thus be set to a value within the interval [0,1]. Note that `glm_scale = 0` diff --git a/src/callbacks_step/glm_speed_dg.jl b/src/callbacks_step/glm_speed_dg.jl index 4d9e4eee23f..dc2579a3359 100644 --- a/src/callbacks_step/glm_speed_dg.jl +++ b/src/callbacks_step/glm_speed_dg.jl @@ -7,7 +7,8 @@ function calc_dt_for_cleaning_speed(cfl::Real, mesh, equations::Union{AbstractIdealGlmMhdEquations, - AbstractIdealGlmMhdMulticomponentEquations}, + AbstractIdealGlmMhdMulticomponentEquations, + AbstractIdealGlmMhdMultiIonEquations}, dg::DG, cache) # compute time step for GLM linear advection equation with c_h=1 for the DG discretization on # Cartesian meshes @@ -28,7 +29,8 @@ end function calc_dt_for_cleaning_speed(cfl::Real, mesh, equations::Union{AbstractIdealGlmMhdEquations, - AbstractIdealGlmMhdMulticomponentEquations}, + AbstractIdealGlmMhdMulticomponentEquations, + AbstractIdealGlmMhdMultiIonEquations}, dg::DGMulti, cache) rd = dg.basis md = mesh.md diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 0304ea16cfe..0799a040bb6 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -491,6 +491,12 @@ abstract type AbstractIdealGlmMhdMulticomponentEquations{NDIMS, NVARS, NCOMP} <: include("ideal_glm_mhd_multicomponent_1d.jl") include("ideal_glm_mhd_multicomponent_2d.jl") +# IdealMhdMultiIonEquations +abstract type AbstractIdealGlmMhdMultiIonEquations{NDIMS, NVARS, NCOMP} <: + AbstractEquations{NDIMS, NVARS} end +include("ideal_glm_mhd_multiion.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, NCOMP}) where { @@ -511,6 +517,27 @@ In particular, not the components themselves are returned. Base.OneTo(ncomponents(equations)) end +# Retrieve number of components from equation instance for the multi-ion case +@inline function ncomponents(::AbstractIdealGlmMhdMultiIonEquations{NDIMS, NVARS, + NCOMP}) where { + NDIMS, + NVARS, + NCOMP + } + NCOMP +end + +""" + eachcomponent(equations::AbstractIdealGlmMhdMultiIonEquations) + +Return an iterator over the indices that specify the location in relevant data structures +for the components in `AbstractIdealGlmMhdMultiIonEquations`. +In particular, not the components themselves are returned. +""" +@inline function eachcomponent(equations::AbstractIdealGlmMhdMultiIonEquations) + Base.OneTo(ncomponents(equations)) +end + # Diffusion equation: first order hyperbolic system abstract type AbstractHyperbolicDiffusionEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end diff --git a/src/equations/ideal_glm_mhd_multiion.jl b/src/equations/ideal_glm_mhd_multiion.jl new file mode 100644 index 00000000000..c76351a3750 --- /dev/null +++ b/src/equations/ideal_glm_mhd_multiion.jl @@ -0,0 +1,266 @@ +# This file includes functions that are common for all AbstractIdealGlmMhdMultiIonEquations + +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +have_nonconservative_terms(::AbstractIdealGlmMhdMultiIonEquations) = True() + +function varnames(::typeof(cons2cons), equations::AbstractIdealGlmMhdMultiIonEquations) + cons = ("B1", "B2", "B3") + for i in eachcomponent(equations) + cons = (cons..., + tuple("rho_" * string(i), "rho_v1_" * string(i), "rho_v2_" * string(i), + "rho_v3_" * string(i), "rho_e_" * string(i))...) + end + cons = (cons..., "psi") + + return cons +end + +function varnames(::typeof(cons2prim), equations::AbstractIdealGlmMhdMultiIonEquations) + prim = ("B1", "B2", "B3") + for i in eachcomponent(equations) + prim = (prim..., + tuple("rho_" * string(i), "v1_" * string(i), "v2_" * string(i), + "v3_" * string(i), "p_" * string(i))...) + end + prim = (prim..., "psi") + + return prim +end + +function default_analysis_integrals(::AbstractIdealGlmMhdMultiIonEquations) + (entropy_timederivative, Val(:l2_divb), Val(:linf_divb)) +end + +""" + source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations) + +Source terms due to the Lorentz' force for plasmas with more than one ion species. These source +terms are a fundamental, inseparable part of the multi-ion GLM-MHD equations, and vanish for +a single-species plasma. In particular, they have to be used for every +simulation of [`IdealGlmMhdMultiIonEquations2D`](@ref). +""" +function source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations) + @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, + equations) + + s = zero(MVector{nvariables(equations), eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations) + rho_inv = 1 / rho + v1 = rho_v1 * rho_inv + v2 = rho_v2 * rho_inv + v3 = rho_v3 * rho_inv + v1_diff = v1_plus - v1 + v2_diff = v2_plus - v2 + v3_diff = v3_plus - v3 + r_rho = charge_to_mass[k] * rho + s2 = r_rho * (v2_diff * B3 - v3_diff * B2) + s3 = r_rho * (v3_diff * B1 - v1_diff * B3) + s4 = r_rho * (v1_diff * B2 - v2_diff * B1) + s5 = v1 * s2 + v2 * s3 + v3 * s4 + + set_component!(s, k, 0, s2, s3, s4, s5, equations) + end + + return SVector(s) +end + +""" + electron_pressure_zero(u, equations::AbstractIdealGlmMhdMultiIonEquations) + +Returns the value of zero for the electron pressure. Needed for consistency with the +single-fluid MHD equations in the limit of one ion species. +""" +function electron_pressure_zero(u, equations::AbstractIdealGlmMhdMultiIonEquations) + return zero(u[1]) +end + +""" + v1, v2, v3, vk1, vk2, vk3 = charge_averaged_velocities(u, + equations::AbstractIdealGlmMhdMultiIonEquations) + + +Compute the charge-averaged velocities (`v1`, `v2`, and `v3`) and each ion species' contribution +to the charge-averaged velocities (`vk1`, `vk2`, and `vk3`). The output variables `vk1`, `vk2`, and `vk3` +are `SVectors` of size `ncomponents(equations)`. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +@inline function charge_averaged_velocities(u, + equations::AbstractIdealGlmMhdMultiIonEquations) + total_electron_charge = zero(real(equations)) + + vk1_plus = zero(MVector{ncomponents(equations), eltype(u)}) + vk2_plus = zero(MVector{ncomponents(equations), eltype(u)}) + vk3_plus = zero(MVector{ncomponents(equations), eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, _ = get_component(k, u, equations) + + total_electron_charge += rho * equations.charge_to_mass[k] + vk1_plus[k] = rho_v1 * equations.charge_to_mass[k] + vk2_plus[k] = rho_v2 * equations.charge_to_mass[k] + vk3_plus[k] = rho_v3 * equations.charge_to_mass[k] + end + vk1_plus ./= total_electron_charge + vk2_plus ./= total_electron_charge + vk3_plus ./= total_electron_charge + v1_plus = sum(vk1_plus) + v2_plus = sum(vk2_plus) + v3_plus = sum(vk3_plus) + + return v1_plus, v2_plus, v3_plus, SVector(vk1_plus), SVector(vk2_plus), + SVector(vk3_plus) +end + +""" + get_component(k, u, equations::AbstractIdealGlmMhdMultiIonEquations) + +Get the hydrodynamic variables of component (ion species) `k`. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +@inline function get_component(k, u, equations::AbstractIdealGlmMhdMultiIonEquations) + return SVector(u[3 + (k - 1) * 5 + 1], + u[3 + (k - 1) * 5 + 2], + u[3 + (k - 1) * 5 + 3], + u[3 + (k - 1) * 5 + 4], + u[3 + (k - 1) * 5 + 5]) +end + +""" + set_component!(u, k, u1, u2, u3, u4, u5, + equations::AbstractIdealGlmMhdMultiIonEquations) + +Set the hydrodynamic variables (`u1` to `u5`) of component (ion species) `k`. + +!!! warning "Experimental implementation" + This is an experimental feature and may change in future releases. +""" +@inline function set_component!(u, k, u1, u2, u3, u4, u5, + equations::AbstractIdealGlmMhdMultiIonEquations) + u[3 + (k - 1) * 5 + 1] = u1 + u[3 + (k - 1) * 5 + 2] = u2 + u[3 + (k - 1) * 5 + 3] = u3 + u[3 + (k - 1) * 5 + 4] = u4 + u[3 + (k - 1) * 5 + 5] = u5 + + return u +end + +# Extract magnetic field from solution vector +magnetic_field(u, equations::AbstractIdealGlmMhdMultiIonEquations) = SVector(u[1], u[2], + u[3]) + +# Extract GLM divergence-cleaning field from solution vector +divergence_cleaning_field(u, equations::AbstractIdealGlmMhdMultiIonEquations) = u[end] + +# Get total density as the sum of the individual densities of the ion species +@inline function density(u, equations::AbstractIdealGlmMhdMultiIonEquations) + rho = zero(real(equations)) + for k in eachcomponent(equations) + rho += u[3 + (k - 1) * 5 + 1] + end + return rho +end + +#Convert conservative variables to primitive +function cons2prim(u, equations::AbstractIdealGlmMhdMultiIonEquations) + @unpack gammas = equations + B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(u, equations) + + prim = zero(MVector{nvariables(equations), eltype(u)}) + prim[1] = B1 + prim[2] = B2 + prim[3] = B3 + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations) + rho_inv = 1 / rho + v1 = rho_inv * rho_v1 + v2 = rho_inv * rho_v2 + v3 = rho_inv * rho_v3 + + p = (gammas[k] - 1) * (rho_e - + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3 + + B1 * B1 + B2 * B2 + B3 * B3 + + psi * psi)) + + set_component!(prim, k, rho, v1, v2, v3, p, equations) + end + prim[end] = psi + + return SVector(prim) +end + +#Convert conservative variables to entropy variables +@inline function cons2entropy(u, equations::AbstractIdealGlmMhdMultiIonEquations) + @unpack gammas = equations + B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(u, equations) + + prim = cons2prim(u, equations) + entropy = zero(MVector{nvariables(equations), eltype(u)}) + rho_p_plus = zero(real(equations)) + for k in eachcomponent(equations) + rho, v1, v2, v3, p = get_component(k, prim, equations) + s = log(p) - gammas[k] * log(rho) + rho_p = rho / p + w1 = (gammas[k] - s) / (gammas[k] - 1) - 0.5f0 * rho_p * (v1^2 + v2^2 + v3^2) + w2 = rho_p * v1 + w3 = rho_p * v2 + w4 = rho_p * v3 + w5 = -rho_p + rho_p_plus += rho_p + + set_component!(entropy, k, w1, w2, w3, w4, w5, equations) + end + + # Additional non-conservative variables + entropy[1] = rho_p_plus * B1 + entropy[2] = rho_p_plus * B2 + entropy[3] = rho_p_plus * B3 + entropy[end] = rho_p_plus * psi + + return SVector(entropy) +end + +# Convert primitive to conservative variables +@inline function prim2cons(prim, equations::AbstractIdealGlmMhdMultiIonEquations) + @unpack gammas = equations + B1, B2, B3 = magnetic_field(prim, equations) + psi = divergence_cleaning_field(prim, equations) + + cons = zero(MVector{nvariables(equations), eltype(prim)}) + cons[1] = B1 + cons[2] = B2 + cons[3] = B3 + for k in eachcomponent(equations) + rho, v1, v2, v3, p = get_component(k, prim, equations) + rho_v1 = rho * v1 + rho_v2 = rho * v2 + rho_v3 = rho * v3 + + rho_e = p / (gammas[k] - 1) + + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) + + 0.5f0 * (B1^2 + B2^2 + B3^2) + + 0.5f0 * psi^2 + + set_component!(cons, k, rho, rho_v1, rho_v2, rho_v3, rho_e, equations) + end + cons[end] = psi + + return SVector(cons) +end +end diff --git a/src/equations/ideal_glm_mhd_multiion_2d.jl b/src/equations/ideal_glm_mhd_multiion_2d.jl new file mode 100644 index 00000000000..3bebdefc0e3 --- /dev/null +++ b/src/equations/ideal_glm_mhd_multiion_2d.jl @@ -0,0 +1,859 @@ +# By default, Julia/LLVM does not use fused multiply-add operations (FMAs). +# Since these FMAs can increase the performance of many numerical algorithms, +# we need to opt-in explicitly. +# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details. +@muladd begin +#! format: noindent + +@doc raw""" + IdealGlmMhdMultiIonEquations2D(; gammas, charge_to_mass, + electron_pressure = electron_pressure_zero, + initial_c_h = NaN) + +The ideal compressible multi-ion MHD equations in two space dimensions augmented with a +generalized Langange multipliers (GLM) divergence-cleaning technique. This is a +multi-species variant of the ideal GLM-MHD equations for calorically perfect plasmas +with independent momentum and energy equations for each ion species. This implementation +assumes that the equations are non-dimensionalized such, that the vacuum permeability is ``\mu_0 = 1``. + +In case of more than one ion species, the specific heat capacity ratios `gammas` and the charge-to-mass +ratios `charge_to_mass` should be passed as tuples, e.g., `gammas=(1.4, 1.667)`. + +The argument `electron_pressure` can be used to pass a function that computes the electron +pressure as a function of the state `u` with the signature `electron_pressure(u, equations::IdealGlmMhdMultiIonEquations2D)`. +By default, the electron pressure is zero. + +The argument `initial_c_h` can be used to set the GLM divergence-cleaning speed. Note that +`initial_c_h = 0` deactivates the divergence cleaning. The callback [`GlmSpeedCallback`](@ref) +can be used to adjust the GLM divergence-cleaning speed according to the time-step size. + +References: +- G. Toth, A. Glocer, Y. Ma, D. Najib, Multi-Ion Magnetohydrodynamics 429 (2010). Numerical + Modeling of Space Plasma Flows, 213–218. +- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization + of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. + [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). + +!!! info "The multi-ion GLM-MHD equations require source terms" + In case of more than one ion species, the multi-ion GLM-MHD equations should ALWAYS be used + with [`source_terms_lorentz`](@ref). +""" +mutable struct IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT <: Real, + ElectronPressure} <: + AbstractIdealGlmMhdMultiIonEquations{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 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")) + + new(gammas, charge_to_mass, electron_pressure, c_h) + end +end + +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)) + __gammas = SVector(map(RealT, _gammas)) + __charge_to_mass = SVector(map(RealT, _charge_to_mass)) + + NVARS = length(_gammas) * 5 + 4 + NCOMP = length(_gammas) + + return IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT, + typeof(electron_pressure)}(__gammas, + __charge_to_mass, + electron_pressure, + initial_c_h) +end + +# Outer constructor for `@reset` works correctly +function IdealGlmMhdMultiIonEquations2D(gammas, charge_to_mass, electron_pressure, c_h) + return IdealGlmMhdMultiIonEquations2D(gammas = gammas, + charge_to_mass = charge_to_mass, + electron_pressure = electron_pressure, + initial_c_h = c_h) +end + +@inline function Base.real(::IdealGlmMhdMultiIonEquations2D{NVARS, NCOMP, RealT}) where { + NVARS, + NCOMP, + RealT + } + RealT +end + +""" + initial_condition_weak_blast_wave(x, t, equations::IdealGlmMhdMultiIonEquations2D) + +A weak blast wave (adapted to multi-ion MHD) from +- Hennemann, S., Rueda-Ramírez, A. M., Hindenlang, F. J., & Gassner, G. J. (2021). A provably entropy + stable subcell shock capturing approach for high order split form DG for the compressible Euler equations. + Journal of Computational Physics, 426, 109935. [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044). + [DOI: 10.1016/j.jcp.2020.109935](https://doi.org/10.1016/j.jcp.2020.109935) +""" +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 + RealT = real(equations) + inicenter = (0, 0) + x_norm = x[1] - inicenter[1] + y_norm = x[2] - inicenter[2] + r = sqrt(x_norm^2 + y_norm^2) + phi = atan(y_norm, x_norm) + + # Calculate primitive variables + rho = r > 0.5f0 ? one(RealT) : convert(RealT, 1.1691) + v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi) + v2 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * sin(phi) + p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245) + + prim = zero(MVector{nvariables(equations), RealT}) + prim[1] = 1 + prim[2] = 1 + prim[3] = 1 + for k in eachcomponent(equations) + set_component!(prim, k, + 2^(k - 1) * (1 - 2) / (1 - 2^ncomponents(equations)) * rho, v1, + v2, 0, p, equations) + end + + return prim2cons(SVector(prim), equations) +end + +# 2D flux of the multi-ion GLM-MHD equations in the direction `orientation` +@inline function flux(u, orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + B1, B2, B3 = magnetic_field(u, equations) + psi = divergence_cleaning_field(u, equations) + + v1_plus, v2_plus, v3_plus, vk1_plus, vk2_plus, vk3_plus = charge_averaged_velocities(u, + equations) + + mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2) + div_clean_energy = 0.5f0 * psi^2 + + f = zero(MVector{nvariables(equations), eltype(u)}) + + if orientation == 1 + f[1] = equations.c_h * psi + f[2] = v1_plus * B2 - v2_plus * B1 + f[3] = v1_plus * B3 - v3_plus * B1 + + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations) + rho_inv = 1 / rho + v1 = rho_v1 * rho_inv + v2 = rho_v2 * rho_inv + v3 = rho_v3 * rho_inv + kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2) + + gamma = equations.gammas[k] + p = (gamma - 1) * (rho_e - kin_en - mag_en - div_clean_energy) + + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = rho_v1 * v2 + f4 = rho_v1 * v3 + f5 = (kin_en + gamma * p / (gamma - 1)) * v1 + 2 * mag_en * vk1_plus[k] - + B1 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + + equations.c_h * psi * B1 + + set_component!(f, k, f1, f2, f3, f4, f5, equations) + end + f[end] = equations.c_h * B1 + else #if orientation == 2 + f[1] = v2_plus * B1 - v1_plus * B2 + f[2] = equations.c_h * psi + f[3] = v2_plus * B3 - v3_plus * B2 + + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, u, equations) + rho_inv = 1 / rho + v1 = rho_v1 * rho_inv + v2 = rho_v2 * rho_inv + v3 = rho_v3 * rho_inv + kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2) + + gamma = equations.gammas[k] + p = (gamma - 1) * (rho_e - kin_en - mag_en - div_clean_energy) + + f1 = rho_v2 + f2 = rho_v2 * v1 + f3 = rho_v2 * v2 + p + f4 = rho_v2 * v3 + f5 = (kin_en + gamma * p / (gamma - 1)) * v2 + 2 * mag_en * vk2_plus[k] - + B2 * (vk1_plus[k] * B1 + vk2_plus[k] * B2 + vk3_plus[k] * B3) + + equations.c_h * psi * B2 + + set_component!(f, k, f1, f2, f3, f4, f5, equations) + end + f[end] = equations.c_h * B2 + end + + return SVector(f) +end + +""" + flux_nonconservative_ruedaramirez_etal(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + +Entropy-conserving non-conservative two-point "flux"" as described in +- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization + of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. + [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). + +!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi.jl" + The non-conservative fluxes derived in the reference above are written as the product + of local and symmetric parts and are meant to be used in the same way as the conservative + fluxes (i.e., flux + flux_noncons in both volume and surface integrals). In this routine, + the fluxes are multiplied by 2 because the non-conservative fluxes are always multiplied + by 0.5 whenever they are used in the Trixi code. + +The term is composed of four individual non-conservative terms: +1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and + is evaluated as a function of the charge-averaged velocity. +2. The Lorentz-force term, which becomes a conservative term in the limit of one ion species for vanishing + electron pressure gradients. +3. The "multi-ion" term, which vanishes in the limit of one ion species. +4. The GLM term, which is needed for Galilean invariance. +""" +@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + @unpack charge_to_mass = equations + # Unpack left and right states to get the magnetic field + B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + # Compute important averages + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 + mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + + # Mean electron pressure + pe_ll = equations.electron_pressure(u_ll, equations) + pe_rr = equations.electron_pressure(u_rr, equations) + pe_mean = 0.5f0 * (pe_ll + pe_rr) + + # Compute charge ratio of u_ll + charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)}) + total_electron_charge = zero(eltype(u_ll)) + for k in eachcomponent(equations) + rho_k = u_ll[3 + (k - 1) * 5 + 1] + charge_ratio_ll[k] = rho_k * charge_to_mass[k] + total_electron_charge += charge_ratio_ll[k] + end + charge_ratio_ll ./= total_electron_charge + + # Compute auxiliary variables + v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll, + equations) + v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr, + equations) + + f = zero(MVector{nvariables(equations), eltype(u_ll)}) + + if orientation == 1 + # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + f[1] = 2 * v1_plus_ll * B1_avg + f[2] = 2 * v2_plus_ll * B1_avg + f[3] = 2 * v3_plus_ll * B1_avg + + for k in eachcomponent(equations) + # Compute term Lorentz term + f2 = charge_ratio_ll[k] * (0.5f0 * mag_norm_avg - B1_avg * B1_avg + pe_mean) + f3 = charge_ratio_ll[k] * (-B1_avg * B2_avg) + f4 = charge_ratio_ll[k] * (-B1_avg * B3_avg) + f5 = vk1_plus_ll[k] * pe_mean + + # Compute multi-ion term (vanishes for NCOMP==1) + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr) + vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr) + vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr) + f5 += (B2_ll * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) + + B3_ll * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) + + # Compute Godunov-Powell term + f2 += charge_ratio_ll[k] * B1_ll * B1_avg + f3 += charge_ratio_ll[k] * B2_ll * B1_avg + f4 += charge_ratio_ll[k] * B3_ll * B1_avg + f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * + B1_avg + + # Compute GLM term for the energy + f5 += v1_plus_ll * psi_ll * psi_avg + + # Add to the flux vector (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5, + equations) + end + # Compute GLM term for psi (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + f[end] = 2 * v1_plus_ll * psi_avg + + else #if orientation == 2 + # Entries of Godunov-Powell term for induction equation (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + f[1] = 2 * v1_plus_ll * B2_avg + f[2] = 2 * v2_plus_ll * B2_avg + f[3] = 2 * v3_plus_ll * B2_avg + + for k in eachcomponent(equations) + # Compute term Lorentz term + f2 = charge_ratio_ll[k] * (-B2_avg * B1_avg) + f3 = charge_ratio_ll[k] * + (-B2_avg * B2_avg + 0.5f0 * mag_norm_avg + pe_mean) + f4 = charge_ratio_ll[k] * (-B2_avg * B3_avg) + f5 = vk2_plus_ll[k] * pe_mean + + # Compute multi-ion term (vanishes for NCOMP==1) + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr) + vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr) + vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr) + f5 += (B1_ll * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) + + B3_ll * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg)) + + # Compute Godunov-Powell term + f2 += charge_ratio_ll[k] * B1_ll * B2_avg + f3 += charge_ratio_ll[k] * B2_ll * B2_avg + f4 += charge_ratio_ll[k] * B3_ll * B2_avg + f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * + B2_avg + + # Compute GLM term for the energy + f5 += v2_plus_ll * psi_ll * psi_avg + + # Add to the flux vector (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5, + equations) + end + # Compute GLM term for psi (multiply by 2 because the non-conservative flux is + # multiplied by 0.5 whenever it's used in the Trixi code) + f[end] = 2 * v2_plus_ll * psi_avg + end + + return SVector(f) +end + +""" + flux_nonconservative_central(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + +Central non-conservative two-point "flux", where the symmetric parts are computed with standard averages. +The use of this term together with [`flux_central`](@ref) +with [`VolumeIntegralFluxDifferencing`](@ref) yields a "standard" +(weak-form) DGSEM discretization of the multi-ion GLM-MHD system. This flux can also be used to construct a +standard local Lax-Friedrichs flux using `surface_flux = (flux_lax_friedrichs, flux_nonconservative_central)`. + +!!! info "Usage and Scaling of Non-Conservative Fluxes in Trixi" + The central non-conservative fluxes implemented in this function are written as the product + of local and symmetric parts, where the symmetric part is a standard average. These fluxes + are meant to be used in the same way as the conservative fluxes (i.e., flux + flux_noncons + in both volume and surface integrals). In this routine, the fluxes are multiplied by 2 because + the non-conservative fluxes are always multiplied by 0.5 whenever they are used in the Trixi code. + +The term is composed of four individual non-conservative terms: +1. The Godunov-Powell term, which arises for plasmas with non-vanishing magnetic field divergence, and + is evaluated as a function of the charge-averaged velocity. +2. The Lorentz-force term, which becomes a conservative term in the limit of one ion species for vanishing + electron pressure gradients. +3. The "multi-ion" term, which vanishes in the limit of one ion species. +4. The GLM term, which is needed for Galilean invariance. +""" +@inline function flux_nonconservative_central(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + @unpack charge_to_mass = equations + # Unpack left and right states to get the magnetic field + B1_ll, B2_ll, B3_ll = magnetic_field(u_ll, equations) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + # Compute important averages + mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 + mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 + + # Electron pressure + pe_ll = equations.electron_pressure(u_ll, equations) + pe_rr = equations.electron_pressure(u_rr, equations) + + # Compute charge ratio of u_ll + charge_ratio_ll = zero(MVector{ncomponents(equations), eltype(u_ll)}) + total_electron_charge = zero(real(equations)) + for k in eachcomponent(equations) + rho_k = u_ll[3 + (k - 1) * 5 + 1] + charge_ratio_ll[k] = rho_k * charge_to_mass[k] + total_electron_charge += charge_ratio_ll[k] + end + charge_ratio_ll ./= total_electron_charge + + # Compute auxiliary variables + v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll, + equations) + v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr, + equations) + + f = zero(MVector{nvariables(equations), eltype(u_ll)}) + + if orientation == 1 + # Entries of Godunov-Powell term for induction equation + f[1] = v1_plus_ll * (B1_ll + B1_rr) + f[2] = v2_plus_ll * (B1_ll + B1_rr) + f[3] = v3_plus_ll * (B1_ll + B1_rr) + for k in eachcomponent(equations) + # Compute Lorentz term + f2 = charge_ratio_ll[k] * ((0.5f0 * mag_norm_ll - B1_ll * B1_ll + pe_ll) + + (0.5f0 * mag_norm_rr - B1_rr * B1_rr + pe_rr)) + f3 = charge_ratio_ll[k] * ((-B1_ll * B2_ll) + (-B1_rr * B2_rr)) + f4 = charge_ratio_ll[k] * ((-B1_ll * B3_ll) + (-B1_rr * B3_rr)) + f5 = vk1_plus_ll[k] * (pe_ll + pe_rr) + + # Compute multi-ion term, which vanishes for NCOMP==1 + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + f5 += (B2_ll * ((vk1_minus_ll * B2_ll - vk2_minus_ll * B1_ll) + + (vk1_minus_rr * B2_rr - vk2_minus_rr * B1_rr)) + + B3_ll * ((vk1_minus_ll * B3_ll - vk3_minus_ll * B1_ll) + + (vk1_minus_rr * B3_rr - vk3_minus_rr * B1_rr))) + + # Compute Godunov-Powell term + f2 += charge_ratio_ll[k] * B1_ll * (B1_ll + B1_rr) + f3 += charge_ratio_ll[k] * B2_ll * (B1_ll + B1_rr) + f4 += charge_ratio_ll[k] * B3_ll * (B1_ll + B1_rr) + f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * + (B1_ll + B1_rr) + + # Compute GLM term for the energy + f5 += v1_plus_ll * psi_ll * (psi_ll + psi_rr) + + # Append to the flux vector + set_component!(f, k, 0, f2, f3, f4, f5, equations) + end + # Compute GLM term for psi + f[end] = v1_plus_ll * (psi_ll + psi_rr) + + else #if orientation == 2 + # Entries of Godunov-Powell term for induction equation + f[1] = v1_plus_ll * (B2_ll + B2_rr) + f[2] = v2_plus_ll * (B2_ll + B2_rr) + f[3] = v3_plus_ll * (B2_ll + B2_rr) + + for k in eachcomponent(equations) + # Compute Lorentz term + f2 = charge_ratio_ll[k] * ((-B2_ll * B1_ll) + (-B2_rr * B1_rr)) + f3 = charge_ratio_ll[k] * ((-B2_ll * B2_ll + 0.5f0 * mag_norm_ll + pe_ll) + + (-B2_rr * B2_rr + 0.5f0 * mag_norm_rr + pe_rr)) + f4 = charge_ratio_ll[k] * ((-B2_ll * B3_ll) + (-B2_rr * B3_rr)) + f5 = vk2_plus_ll[k] * (pe_ll + pe_rr) + + # Compute multi-ion term (vanishes for NCOMP==1) + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + f5 += (B1_ll * ((vk2_minus_ll * B1_ll - vk1_minus_ll * B2_ll) + + (vk2_minus_rr * B1_rr - vk1_minus_rr * B2_rr)) + + B3_ll * ((vk2_minus_ll * B3_ll - vk3_minus_ll * B2_ll) + + (vk2_minus_rr * B3_rr - vk3_minus_rr * B2_rr))) + + # Compute Godunov-Powell term + f2 += charge_ratio_ll[k] * B1_ll * (B2_ll + B2_rr) + f3 += charge_ratio_ll[k] * B2_ll * (B2_ll + B2_rr) + f4 += charge_ratio_ll[k] * B3_ll * (B2_ll + B2_rr) + f5 += (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) * + (B2_ll + B2_rr) + + # Compute GLM term for the energy + f5 += v2_plus_ll * psi_ll * (psi_ll + psi_rr) + + # Append to the flux vector + set_component!(f, k, 0, f2, f3, f4, f5, equations) + end + # Compute GLM term for psi + f[end] = v2_plus_ll * (psi_ll + psi_rr) + end + + return SVector(f) +end + +""" + flux_ruedaramirez_etal(u_ll, u_rr, orientation, equations::IdealGlmMhdMultiIonEquations2D) + +Entropy conserving two-point flux for the multi-ion GLM-MHD equations from +- A. Rueda-Ramírez, A. Sikstel, G. Gassner, An Entropy-Stable Discontinuous Galerkin Discretization + of the Ideal Multi-Ion Magnetohydrodynamics System (2024). Journal of Computational Physics. + [DOI: 10.1016/j.jcp.2024.113655](https://doi.org/10.1016/j.jcp.2024.113655). + +This flux (together with the MHD non-conservative term) is consistent in the case of one ion species with the flux of: +- Derigs et al. (2018). Ideal GLM-MHD: About the entropy consistent nine-wave magnetic field + divergence diminishing ideal magnetohydrodynamics equations for multi-ion + [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::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) + B1_rr, B2_rr, B3_rr = magnetic_field(u_rr, equations) + psi_ll = divergence_cleaning_field(u_ll, equations) + psi_rr = divergence_cleaning_field(u_rr, equations) + + v1_plus_ll, v2_plus_ll, v3_plus_ll, vk1_plus_ll, vk2_plus_ll, vk3_plus_ll = charge_averaged_velocities(u_ll, + equations) + v1_plus_rr, v2_plus_rr, v3_plus_rr, vk1_plus_rr, vk2_plus_rr, vk3_plus_rr = charge_averaged_velocities(u_rr, + equations) + + f = zero(MVector{nvariables(equations), eltype(u_ll)}) + + # Compute averages for global variables + v1_plus_avg = 0.5f0 * (v1_plus_ll + v1_plus_rr) + v2_plus_avg = 0.5f0 * (v2_plus_ll + v2_plus_rr) + v3_plus_avg = 0.5f0 * (v3_plus_ll + v3_plus_rr) + B1_avg = 0.5f0 * (B1_ll + B1_rr) + B2_avg = 0.5f0 * (B2_ll + B2_rr) + B3_avg = 0.5f0 * (B3_ll + B3_rr) + mag_norm_ll = B1_ll^2 + B2_ll^2 + B3_ll^2 + mag_norm_rr = B1_rr^2 + B2_rr^2 + B3_rr^2 + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr) + psi_avg = 0.5f0 * (psi_ll + psi_rr) + + if orientation == 1 + psi_B1_avg = 0.5f0 * (B1_ll * psi_ll + B1_rr * psi_rr) + + # Magnetic field components from f^MHD + f6 = equations.c_h * psi_avg + f7 = v1_plus_avg * B2_avg - v2_plus_avg * B1_avg + f8 = v1_plus_avg * B3_avg - v3_plus_avg * B1_avg + f9 = equations.c_h * B1_avg + + # Start building the flux + f[1] = f6 + f[2] = f7 + f[3] = f8 + f[end] = f9 + + # Iterate over all components + for k in eachcomponent(equations) + # Unpack left and right states + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = get_component(k, u_ll, + equations) + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr, + equations) + rho_inv_ll = 1 / rho_ll + v1_ll = rho_v1_ll * rho_inv_ll + v2_ll = rho_v2_ll * rho_inv_ll + v3_ll = rho_v3_ll * rho_inv_ll + rho_inv_rr = 1 / rho_rr + v1_rr = rho_v1_rr * rho_inv_rr + v2_rr = rho_v2_rr * rho_inv_rr + v3_rr = rho_v3_rr * rho_inv_rr + vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 + vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 + + p_ll = (gammas[k] - 1) * + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll - + 0.5f0 * psi_ll^2) + p_rr = (gammas[k] - 1) * + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr - + 0.5f0 * psi_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + # for convenience store vk_plus⋅B + vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll + + vk3_plus_ll[k] * B3_ll + vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr + + vk3_plus_rr[k] * B3_rr + + # Compute the necessary mean values needed for either direction + rho_avg = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) + vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k]) + vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k]) + vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k]) + # v_minus + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr) + vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr) + vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr) + + # Fill the fluxes for the mass and momentum equations + f1 = rho_mean * v1_avg + f2 = f1 * v1_avg + p_mean + f3 = f1 * v2_avg + f4 = f1 * v3_avg + + # total energy flux is complicated and involves the previous eight components + v1_plus_mag_avg = 0.5f0 * (vk1_plus_ll[k] * mag_norm_ll + + vk1_plus_rr[k] * mag_norm_rr) + # Euler part + f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + # MHD part + f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v1_plus_mag_avg + + B1_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + + f9 * psi_avg - equations.c_h * psi_B1_avg # GLM term + + + 0.5f0 * vk1_plus_avg * mag_norm_avg - + vk1_plus_avg * B1_avg * B1_avg - vk2_plus_avg * B1_avg * B2_avg - + vk3_plus_avg * B1_avg * B3_avg # Additional terms related to the Lorentz non-conservative term (momentum eqs) + - + B2_avg * (vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) - + B3_avg * (vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg)) # Terms related to the multi-ion non-conservative term (induction equation!) + + set_component!(f, k, f1, f2, f3, f4, f5, equations) + end + else #if orientation == 2 + psi_B2_avg = 0.5f0 * (B2_ll * psi_ll + B2_rr * psi_rr) + + # Magnetic field components from f^MHD + f6 = v2_plus_avg * B1_avg - v1_plus_avg * B2_avg + f7 = equations.c_h * psi_avg + f8 = v2_plus_avg * B3_avg - v3_plus_avg * B2_avg + f9 = equations.c_h * B2_avg + + # Start building the flux + f[1] = f6 + f[2] = f7 + f[3] = f8 + f[end] = f9 + + # Iterate over all components + for k in eachcomponent(equations) + # Unpack left and right states + rho_ll, rho_v1_ll, rho_v2_ll, rho_v3_ll, rho_e_ll = get_component(k, u_ll, + equations) + rho_rr, rho_v1_rr, rho_v2_rr, rho_v3_rr, rho_e_rr = get_component(k, u_rr, + equations) + + rho_inv_ll = 1 / rho_ll + v1_ll = rho_v1_ll * rho_inv_ll + v2_ll = rho_v2_ll * rho_inv_ll + v3_ll = rho_v3_ll * rho_inv_ll + rho_inv_rr = 1 / rho_rr + v1_rr = rho_v1_rr * rho_inv_rr + v2_rr = rho_v2_rr * rho_inv_rr + v3_rr = rho_v3_rr * rho_inv_rr + vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2 + vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2 + + p_ll = (gammas[k] - 1) * + (rho_e_ll - 0.5f0 * rho_ll * vel_norm_ll - 0.5f0 * mag_norm_ll - + 0.5f0 * psi_ll^2) + p_rr = (gammas[k] - 1) * + (rho_e_rr - 0.5f0 * rho_rr * vel_norm_rr - 0.5f0 * mag_norm_rr - + 0.5f0 * psi_rr^2) + beta_ll = 0.5f0 * rho_ll / p_ll + beta_rr = 0.5f0 * rho_rr / p_rr + # for convenience store vk_plus⋅B + vel_dot_mag_ll = vk1_plus_ll[k] * B1_ll + vk2_plus_ll[k] * B2_ll + + vk3_plus_ll[k] * B3_ll + vel_dot_mag_rr = vk1_plus_rr[k] * B1_rr + vk2_plus_rr[k] * B2_rr + + vk3_plus_rr[k] * B3_rr + + # Compute the necessary mean values needed for either direction + rho_avg = 0.5f0 * (rho_ll + rho_rr) + rho_mean = ln_mean(rho_ll, rho_rr) + beta_mean = ln_mean(beta_ll, beta_rr) + beta_avg = 0.5f0 * (beta_ll + beta_rr) + p_mean = 0.5f0 * rho_avg / beta_avg + v1_avg = 0.5f0 * (v1_ll + v1_rr) + v2_avg = 0.5f0 * (v2_ll + v2_rr) + v3_avg = 0.5f0 * (v3_ll + v3_rr) + vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr) + vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr) + vk1_plus_avg = 0.5f0 * (vk1_plus_ll[k] + vk1_plus_rr[k]) + vk2_plus_avg = 0.5f0 * (vk2_plus_ll[k] + vk2_plus_rr[k]) + vk3_plus_avg = 0.5f0 * (vk3_plus_ll[k] + vk3_plus_rr[k]) + # v_minus + vk1_minus_ll = v1_plus_ll - vk1_plus_ll[k] + vk2_minus_ll = v2_plus_ll - vk2_plus_ll[k] + vk3_minus_ll = v3_plus_ll - vk3_plus_ll[k] + vk1_minus_rr = v1_plus_rr - vk1_plus_rr[k] + vk2_minus_rr = v2_plus_rr - vk2_plus_rr[k] + vk3_minus_rr = v3_plus_rr - vk3_plus_rr[k] + vk1_minus_avg = 0.5f0 * (vk1_minus_ll + vk1_minus_rr) + vk2_minus_avg = 0.5f0 * (vk2_minus_ll + vk2_minus_rr) + vk3_minus_avg = 0.5f0 * (vk3_minus_ll + vk3_minus_rr) + + # Fill the fluxes for the mass and momentum equations + f1 = rho_mean * v2_avg + f2 = f1 * v1_avg + f3 = f1 * v2_avg + p_mean + f4 = f1 * v3_avg + + # total energy flux is complicated and involves the previous eight components + v2_plus_mag_avg = 0.5f0 * (vk2_plus_ll[k] * mag_norm_ll + + vk2_plus_rr[k] * mag_norm_rr) + # Euler part + f5 = f1 * 0.5f0 * (1 / (gammas[k] - 1) / beta_mean - vel_norm_avg) + + f2 * v1_avg + f3 * v2_avg + f4 * v3_avg + # MHD part + f5 += (f6 * B1_avg + f7 * B2_avg + f8 * B3_avg - 0.5f0 * v2_plus_mag_avg + + B2_avg * vel_dot_mag_avg # Same terms as in Derigs (but with v_plus) + + f9 * psi_avg - equations.c_h * psi_B2_avg # GLM term + + + 0.5f0 * vk2_plus_avg * mag_norm_avg - + vk1_plus_avg * B2_avg * B1_avg - vk2_plus_avg * B2_avg * B2_avg - + vk3_plus_avg * B2_avg * B3_avg # Additional terms related to the Lorentz non-conservative term (momentum eqs) + - + B1_avg * (vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) - + B3_avg * (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg)) # Terms related to the multi-ion non-conservative term (induction equation!) + + set_component!(f, k, f1, f2, f3, f4, f5, equations) + end + end + + return SVector(f) +end + +# Calculate maximum wave speed for local Lax-Friedrichs-type dissipation +# This routine approximates the maximum wave speed as sum of the maximum ion velocity +# for all species and the maximum magnetosonic speed. +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + # Calculate fast magnetoacoustic wave speeds + # left + cf_ll = calc_fast_wavespeed(u_ll, orientation, equations) + # right + cf_rr = calc_fast_wavespeed(u_rr, orientation, equations) + + # Calculate velocities + v_ll = zero(eltype(u_ll)) + v_rr = zero(eltype(u_rr)) + if orientation == 1 + for k in eachcomponent(equations) + rho, rho_v1, _ = get_component(k, u_ll, equations) + v_ll = max(v_ll, abs(rho_v1 / rho)) + rho, rho_v1, _ = get_component(k, u_rr, equations) + v_rr = max(v_rr, abs(rho_v1 / rho)) + end + else #if orientation == 2 + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, _ = get_component(k, u_ll, equations) + v_ll = max(v_ll, abs(rho_v2 / rho)) + rho, rho_v1, rho_v2, _ = get_component(k, u_rr, equations) + v_rr = max(v_rr, abs(rho_v2 / rho)) + end + end + + λ_max = max(abs(v_ll), abs(v_rr)) + max(cf_ll, cf_rr) +end + +@inline function max_abs_speeds(u, equations::IdealGlmMhdMultiIonEquations2D) + v1 = zero(real(equations)) + v2 = zero(real(equations)) + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, _ = get_component(k, u, equations) + v1 = max(v1, abs(rho_v1 / rho)) + v2 = max(v2, abs(rho_v2 / rho)) + end + + cf_x_direction = calc_fast_wavespeed(u, 1, equations) + cf_y_direction = calc_fast_wavespeed(u, 2, equations) + + return (abs(v1) + cf_x_direction, abs(v2) + cf_y_direction) +end + +# Compute the fastest wave speed for ideal multi-ion GLM-MHD equations: c_f, the fast +# magnetoacoustic eigenvalue. This routine computes the fast magnetosonic speed for each ion +# species using the single-fluid MHD expressions and approximates the multi-ion c_f as +# the maximum of these individual magnetosonic speeds. +@inline function calc_fast_wavespeed(cons, orientation::Integer, + equations::IdealGlmMhdMultiIonEquations2D) + B1, B2, B3 = magnetic_field(cons, equations) + psi = divergence_cleaning_field(cons, equations) + + c_f = zero(real(equations)) + for k in eachcomponent(equations) + rho, rho_v1, rho_v2, rho_v3, rho_e = get_component(k, cons, equations) + + rho_inv = 1 / rho + v1 = rho_v1 * rho_inv + v2 = rho_v2 * rho_inv + v3 = rho_v3 * rho_inv + v_mag = sqrt(v1^2 + v2^2 + v3^2) + gamma = equations.gammas[k] + p = (gamma - 1) * + (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2) - + 0.5f0 * psi^2) + a_square = gamma * p * rho_inv + inv_sqrt_rho = 1 / sqrt(rho) + + b1 = B1 * inv_sqrt_rho + b2 = B2 * inv_sqrt_rho + b3 = B3 * inv_sqrt_rho + b_square = b1^2 + b2^2 + b3^2 + + if orientation == 1 + c_f = max(c_f, + sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * + sqrt((a_square + b_square)^2 - 4 * a_square * b1^2))) + else #if orientation == 2 + c_f = max(c_f, + sqrt(0.5f0 * (a_square + b_square) + + 0.5f0 * + sqrt((a_square + b_square)^2 - 4 * a_square * b2^2))) + end + end + + return c_f +end +end # @muladd diff --git a/test/test_tree_2d_mhdmultiion.jl b/test/test_tree_2d_mhdmultiion.jl new file mode 100644 index 00000000000..2ef8b56e9e8 --- /dev/null +++ b/test/test_tree_2d_mhdmultiion.jl @@ -0,0 +1,104 @@ +module TestExamples2DIdealGlmMhdMultiIon + +using Test +using Trixi + +include("test_trixi.jl") + +# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory +EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem") + +@testset "MHD Multi-ion" begin +#! format: noindent + +@trixi_testset "elixir_mhdmultiion_ec.jl" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"), + l2=[ + 0.018116158127836963, + 0.018242057606900185, + 0.02842532060540296, + 0.015805662470907644, + 0.018125738668209504, + 0.01841652763250858, + 0.005540867360495384, + 0.2091437538608106, + 0.021619156799420326, + 0.030561821442897923, + 0.029778012392807182, + 0.018651133866567658, + 0.1204153285953923, + 0.00017580766470956546 + ], + linf=[ + 0.08706810449620284, + 0.085845022468381, + 0.15301576281555795, + 0.07924698349506726, + 0.11185098353068998, + 0.10902813348518392, + 0.03512713073214368, + 1.018607282408099, + 0.11897673517549667, + 0.18757934071639615, + 0.19089795413679128, + 0.10989483685170078, + 0.6078381519649727, + 0.00673110606965085 + ]) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end + +@trixi_testset "elixir_mhdmultiion_ec.jl with local Lax-Friedrichs at the surface" begin + @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmultiion_ec.jl"), + l2=[ + 0.017668026737187294, + 0.017797845988889206, + 0.02784171335711194, + 0.015603472482130114, + 0.017849085712788024, + 0.018142011483140937, + 0.005478230098832851, + 0.20585547149049335, + 0.021301307197244185, + 0.030185081369414384, + 0.029385190432285654, + 0.01837280802521888, + 0.11810313415151609, + 0.0002962121353575025 + ], + linf=[ + 0.06594404679380572, + 0.0658741409805832, + 0.09451390288403083, + 0.06787235653416557, + 0.0891017651119973, + 0.08828137594061974, + 0.02364750099643367, + 0.8059321345611803, + 0.12243877249999213, + 0.15930897967901766, + 0.1538327998380914, + 0.08694877996306204, + 0.49493751138636366, + 0.003287414714660175 + ], + surface_flux=(flux_lax_friedrichs, flux_nonconservative_central)) + # Ensure that we do not have excessive memory allocations + # (e.g., from type instabilities) + let + t = sol.t[end] + u_ode = sol.u[end] + du_ode = similar(u_ode) + @test (@allocated Trixi.rhs!(du_ode, u_ode, semi, t)) < 1000 + end +end +end + +end # module diff --git a/test/test_tree_2d_part3.jl b/test/test_tree_2d_part3.jl index 0eff564132c..42633ceae17 100644 --- a/test/test_tree_2d_part3.jl +++ b/test/test_tree_2d_part3.jl @@ -20,6 +20,9 @@ isdir(outdir) && rm(outdir, recursive = true) # MHD Multicomponent include("test_tree_2d_mhdmulti.jl") + # MHD Multi-ion + include("test_tree_2d_mhdmultiion.jl") + # Lattice-Boltzmann include("test_tree_2d_lbm.jl") diff --git a/test/test_type.jl b/test/test_type.jl index 131f5f3e1a1..e26b34cdfc6 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1503,6 +1503,75 @@ isdir(outdir) && rm(outdir, recursive = true) end end + @timed_testset "Ideal Glm Mhd MultiIon 2D" begin + for RealT in (Float32, Float64) + gammas = (RealT(2), RealT(2)) + charge_to_mass = (RealT(2), RealT(2)) + equations = @inferred IdealGlmMhdMultiIonEquations2D(gammas = gammas, + charge_to_mass = charge_to_mass) + + x = SVector(zero(RealT), zero(RealT)) + t = zero(RealT) + u = u_ll = u_rr = cons = SVector(one(RealT), one(RealT), one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT), + one(RealT)) + orientations = [1, 2] + + @test eltype(@inferred initial_condition_weak_blast_wave(x, t, equations)) == + RealT + + @test eltype(@inferred source_terms_lorentz(u, x, t, equations)) == + RealT + + for orientation in orientations + @test eltype(@inferred flux(u, orientation, equations)) == RealT + @test eltype(@inferred flux_nonconservative_ruedaramirez_etal(u_ll, u_rr, + orientation, + equations)) == + RealT + @test eltype(@inferred flux_nonconservative_central(u_ll, u_rr, orientation, + equations)) == + RealT + @test eltype(@inferred flux_ruedaramirez_etal(u_ll, u_rr, orientation, + equations)) == RealT + + @test typeof(@inferred max_abs_speed_naive(u_ll, u_rr, orientation, + equations)) == + RealT + end + + @test eltype(@inferred Trixi.max_abs_speeds(u, equations)) == RealT + @test eltype(@inferred cons2prim(u, equations)) == RealT + @test eltype(@inferred prim2cons(u, equations)) == RealT + @test eltype(@inferred cons2entropy(u, equations)) == RealT + @test typeof(@inferred density(u, equations)) == RealT + @test eltype(@inferred magnetic_field(u, equations)) == RealT + @test typeof(@inferred divergence_cleaning_field(u, equations)) == RealT + @test typeof(@inferred Trixi.electron_pressure_zero(u, equations)) == RealT + + @test typeof(@inferred Trixi.charge_averaged_velocities(u, equations)) == + Tuple{RealT, RealT, RealT, SVector{2, RealT}, SVector{2, RealT}, + SVector{2, RealT}} + + for k in 1:2 + @test eltype(@inferred Trixi.get_component(k, u, equations)) == RealT + end + + for direction in orientations + @test typeof(Trixi.calc_fast_wavespeed(cons, direction, equations)) == RealT + end + end + end + @timed_testset "Inviscid Burgers 1D" begin for RealT in (Float32, Float64) equations = @inferred InviscidBurgersEquation1D()