Add 2D multi-ion GLM-MHD equations for the TreeMesh solver

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

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,9 +229,10 @@ 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
+export get_component
 
 export TreeMesh, StructuredMesh, StructuredMeshView, UnstructuredMesh2D, P4estMesh,
        T8codeMesh

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`

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

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