From 64d92bf87ed590a0163f2f93029339afa5393c2d Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 13 Jan 2025 22:03:20 +0100 Subject: [PATCH 1/6] add new equation with two different species of ion and plasma --- .../elixir_euler_ion_electron.jl | 55 +++ src/equations/compressible_euler_plasma.jl | 323 ++++++++++++++++++ src/equations/equations.jl | 1 + 3 files changed, 379 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_euler_ion_electron.jl create mode 100644 src/equations/compressible_euler_plasma.jl diff --git a/examples/tree_1d_dgsem/elixir_euler_ion_electron.jl b/examples/tree_1d_dgsem/elixir_euler_ion_electron.jl new file mode 100644 index 0000000000..ec8da32fee --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_ion_electron.jl @@ -0,0 +1,55 @@ +using Trixi +using OrdinaryDiffEq + +############################################################################### +# semidiscretization of the compressible Euler equations +equations = Trixi.CompressibleEulerPlasmaEquations1D(1.4) + +initial_condition = initial_condition_convergence_test + +solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, + volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc)) + +coordinates_min = 0.0 +coordinates_max = 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_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_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 \ No newline at end of file diff --git a/src/equations/compressible_euler_plasma.jl b/src/equations/compressible_euler_plasma.jl new file mode 100644 index 0000000000..69a33af7de --- /dev/null +++ b/src/equations/compressible_euler_plasma.jl @@ -0,0 +1,323 @@ +struct CompressibleEulerPlasmaEquations1D{RealT <: Real} <: + AbstractCompressibleEulerEquations{1, 6} + gamma::RealT # ratio of specific heats + inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + function CompressibleEulerPlasmaEquations1D(gamma) + γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) + new{typeof(γ)}(γ, inv_gamma_minus_one) + end +end + +function varnames(::typeof(cons2cons), ::CompressibleEulerPlasmaEquations1D) + ("rho_ion", "rho_v1_ion", "rho_e_ion", "rho_el", "rho_v1_el", "rho_e_el") +end +varnames(::typeof(cons2prim), ::CompressibleEulerPlasmaEquations1D) = ("rho_ion", "v1_ion", "p_ion", "rho_el", "v1_el", "p_el") + +""" + single_species_flux(u, gamma) + +Calculate the flux vector for a single species (ion or electron). +""" +@inline function single_species_flux(rho, rho_v1, rho_e, gamma) + v1 = rho_v1 / rho + p = (gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + return SVector(rho_v1, rho_v1 * v1 + p, (rho_e + p) * v1) +end + +""" + flux(u, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) + +Calculate the flux for the full plasma system with ions and electrons. +""" +@inline function flux(u, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) + rho_ion, rho_v1_ion, rho_e_ion, rho_el, rho_v1_el, rho_e_el = u + + f_ion = single_species_flux(rho_ion, rho_v1_ion, rho_e_ion, equations.gamma) + + f_el = single_species_flux(rho_el, rho_v1_el, rho_e_el, equations.gamma) + + return SVector(f_ion[1], f_ion[2], f_ion[3], f_el[1], f_el[2], f_el[3]) +end + +""" + initial_condition_constant(x, t, equations::CompressibleEulerPlasmaEquations1D) + +A constant initial condition to test free-stream preservation. +""" +function initial_condition_constant(x, t, equations::CompressibleEulerPlasmaEquations1D) + RealT = eltype(x) + rho_ion = 1 + rho_el = 1 + + # Ion and electron velocities + rho_v1_ion = convert(RealT, 0.01) + rho_v1_el = convert(RealT, 0.1) + + # Ion and electron internal energies + rho_e_ion = 10.0 + rho_e_el = 10.0 + + return SVector(rho_ion, rho_v1_ion, rho_e_ion, rho_el, rho_v1_el, rho_e_el) +end + +""" + initial_condition_convergence_test(x, t, equations::CompressibleEulerPlasmaEquations1D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function initial_condition_convergence_test(x, t, equations::CompressibleEulerPlasmaEquations1D) + RealT = eltype(x) + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f + ini = c + A * sin(ω * (x[1] - t)) + + rho = ini + rho_v1 = ini + rho_e = ini^2 + + return SVector(rho, rho_v1, rho_e, rho, rho_v1, rho_e) +end + + +""" + source_terms_convergence_test(u, x, t, equations::CompressibleEulerPlasmaEquations1D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::CompressibleEulerPlasmaEquations1D) + # Same settings as in `initial_condition` + RealT = eltype(u) + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f + γ = equations.gamma + + x1, = x + + si, co = sincos(ω * (x1 - t)) + rho = c + A * si + rho_x = ω * A * co + + # Note that d/dt rho = -d/dx rho. + # This yields du2 = du3 = d/dx p (derivative of pressure). + # Other terms vanish because of v = 1. + du1 = 0 + du2 = rho_x * (2 * rho - 0.5f0) * (γ - 1) + du3 = du2 + + return SVector(du1, du2, du3, du1, du2, du3) +end + +#TODO: Maybe there is a better way to do this idk +function single_species_flux_hllc(u_ll, u_rr, orientation::Integer, + gamma) + # Calculate primitive variables and speed of sound + rho_ll, rho_v1_ll, rho_e_ll = u_ll + rho_rr, rho_v1_rr, rho_e_rr = u_rr + + v1_ll = rho_v1_ll / rho_ll + e_ll = rho_e_ll / rho_ll + p_ll = (gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * v1_ll^2) + c_ll = sqrt(gamma * p_ll / rho_ll) + + v1_rr = rho_v1_rr / rho_rr + e_rr = rho_e_rr / rho_rr + p_rr = (gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * v1_rr^2) + c_rr = sqrt(gamma * p_rr / rho_rr) + + # Obtain left and right fluxes + f_ll = single_species_flux(rho_ll, rho_v1_ll, rho_e_ll, gamma) + f_rr = single_species_flux(rho_rr, rho_v1_rr, rho_e_rr, gamma) + + # Compute Roe averages + sqrt_rho_ll = sqrt(rho_ll) + sqrt_rho_rr = sqrt(rho_rr) + sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr + vel_L = v1_ll + vel_R = v1_rr + vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho + ekin_roe = 0.5f0 * vel_roe^2 + H_ll = (rho_e_ll + p_ll) / rho_ll + H_rr = (rho_e_rr + p_rr) / rho_rr + H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho + c_roe = sqrt((gamma - 1) * (H_roe - ekin_roe)) + + Ssl = min(vel_L - c_ll, vel_roe - c_roe) + Ssr = max(vel_R + c_rr, vel_roe + c_roe) + sMu_L = Ssl - vel_L + sMu_R = Ssr - vel_R + if Ssl >= 0 + f1 = f_ll[1] + f2 = f_ll[2] + f3 = f_ll[3] + elseif Ssr <= 0 + f1 = f_rr[1] + f2 = f_rr[2] + f3 = f_rr[3] + else + SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / + (rho_ll * sMu_L - rho_rr * sMu_R) + if Ssl <= 0 <= SStar + densStar = rho_ll * sMu_L / (Ssl - SStar) + enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) + UStar1 = densStar + UStar2 = densStar * SStar + UStar3 = densStar * enerStar + + f1 = f_ll[1] + Ssl * (UStar1 - rho_ll) + f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) + f3 = f_ll[3] + Ssl * (UStar3 - rho_e_ll) + else + densStar = rho_rr * sMu_R / (Ssr - SStar) + enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) + UStar1 = densStar + UStar2 = densStar * SStar + UStar3 = densStar * enerStar + + #end + f1 = f_rr[1] + Ssr * (UStar1 - rho_rr) + f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) + f3 = f_rr[3] + Ssr * (UStar3 - rho_e_rr) + end + end + return SVector(f1, f2, f3) +end + +""" + flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerPlasmaEquations1D) + +Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro +[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) +Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140) +""" +function flux_hllc(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) + flux_hllc_ion = single_species_flux_hllc(u_ll[1:3], u_rr[1:3], orientation, equations.gamma) + flux_hllc_el = single_species_flux_hllc(u_ll[4:6], u_rr[4:6], orientation, equations.gamma) + return vcat(flux_hllc_ion, flux_hllc_el) +end + +@inline function max_abs_speeds(u, equations::CompressibleEulerPlasmaEquations1D) + rho_ion, v1_ion, p_ion, rho_el, v1_el, p_el = cons2prim(u, equations) + c_ion = sqrt(equations.gamma * p_ion / rho_ion) + c_el = sqrt(equations.gamma * p_el / rho_el) + + return (abs(v1_ion) + c_ion, abs(v1_el) + c_el) # Is this the right way? I've seen multicomponent where (abs(v1) + c,) is returned +end + +@inline function single_species_cons2prim(u, gamma) + rho, rho_v1, rho_e = u + v1 = rho_v1 / rho + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v1^2) + return SVector(rho, v1, p) +end + +@inline function cons2prim(u, equations::CompressibleEulerPlasmaEquations1D) + prim_ion = single_species_cons2prim(u[1:3], equations.gamma) + + prim_el = single_species_cons2prim(u[4:6], equations.gamma) + + return vcat(prim_ion, prim_el) +end + +@inline function single_species_cons2entropy(u, gamma, inv_gamma_minus_one) + rho, rho_v1, rho_e = u + + v1 = rho_v1 / rho + v_square = v1^2 + p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) + s = log(p) - gamma * log(rho) + rho_p = rho / p + + w1 = (gamma - s) * inv_gamma_minus_one - + 0.5f0 * rho_p * v_square + w2 = rho_p * v1 + w3 = -rho_p + + return SVector(w1, w2, w3) +end + +# Convert conservative variables to entropy +@inline function cons2entropy(u, equations::CompressibleEulerPlasmaEquations1D) + w_ion = single_species_cons2entropy(u[1:3], equations.gamma, equations.inv_gamma_minus_one) + w_el = single_species_cons2entropy(u[4:6], equations.gamma, equations.inv_gamma_minus_one) + return vcat(w_ion, w_el) +end + +@inline function single_species_entropy2cons(w, gamma, inv_gamma_minus_one) + # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD + # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + + # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986) + # instead of `-rho * s / (gamma - 1)` + V1, V2, V5 = w .* (gamma - 1) + + # specific entropy, eq. (53) + s = gamma - V1 + 0.5f0 * (V2^2) / V5 + + # eq. (52) + energy_internal = ((gamma - 1) / (-V5)^gamma)^(inv_gamma_minus_one) * + exp(-s * inv_gamma_minus_one) + + # eq. (51) + rho = -V5 * energy_internal + rho_v1 = V2 * energy_internal + rho_e = (1 - 0.5f0 * (V2^2) / V5) * energy_internal + return SVector(rho, rho_v1, rho_e) +end + +@inline function entropy2cons(w, equations::CompressibleEulerPlasmaEquations1D) + u_ion = single_species_entropy2cons(w[1:3], equations.gamma, equations.inv_gamma_minus_one) + u_el = single_species_entropy2cons(w[4:6], equations.gamma, equations.inv_gamma_minus_one) + return vcat(u_ion, u_el) +end + +# Convert primitive to conservative variables +@inline function single_species_prim2cons(prim, inv_gamma_minus_one) + rho, v1, p = prim + rho_v1 = rho * v1 + rho_e = p * inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1) + return SVector(rho, rho_v1, rho_e) +end + +@inline function prim2cons(prim, equations::CompressibleEulerPlasmaEquations1D) + u_ion = single_species_prim2cons(prim[1:3], equations.inv_gamma_minus_one) + u_el = single_species_prim2cons(prim[4:6], equations.inv_gamma_minus_one) + return vcat(u_ion, u_el) +end + +@inline function density(u, equations::CompressibleEulerPlasmaEquations1D) + rho_ion = u[1] + rho_el = u[4] + return rho_ion, rho_el +end + +@inline function velocity(u, equations::CompressibleEulerPlasmaEquations1D) + rho_ion = u[1] + rho_el = u[4] + v1_ion = u[2] / rho_ion + v1_el = u[5] / rho_el + return v1_ion, v1_el +end + +@inline function single_species_pressure(u) + rho, rho_v1, rho_e = u + p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) + return p +end + +@inline function pressure(u, equations::CompressibleEulerPlasmaEquations1D) + p_ion = single_species_pressure(u[1:3]) + p_el = single_species_pressure(u[4:6]) + return p_ion, p_el +end \ No newline at end of file diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 2cecce6759..4051b56d65 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -442,6 +442,7 @@ include("compressible_euler_1d.jl") include("compressible_euler_2d.jl") include("compressible_euler_3d.jl") include("compressible_euler_quasi_1d.jl") +include("compressible_euler_plasma.jl") # CompressibleEulerMulticomponentEquations abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: From 1220fb663846feed4f78f7300ce60b954a8e5042 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 14 Jan 2025 17:32:37 +0100 Subject: [PATCH 2/6] multiion system implementation in 1D --- ...n_electron.jl => elixir_euler_multiion.jl} | 7 +- src/equations/compressible_euler_multiion.jl | 424 ++++++++++++++++++ src/equations/compressible_euler_plasma.jl | 323 ------------- src/equations/equations.jl | 26 +- 4 files changed, 453 insertions(+), 327 deletions(-) rename examples/tree_1d_dgsem/{elixir_euler_ion_electron.jl => elixir_euler_multiion.jl} (90%) create mode 100644 src/equations/compressible_euler_multiion.jl delete mode 100644 src/equations/compressible_euler_plasma.jl diff --git a/examples/tree_1d_dgsem/elixir_euler_ion_electron.jl b/examples/tree_1d_dgsem/elixir_euler_multiion.jl similarity index 90% rename from examples/tree_1d_dgsem/elixir_euler_ion_electron.jl rename to examples/tree_1d_dgsem/elixir_euler_multiion.jl index ec8da32fee..a1c7e2cd64 100644 --- a/examples/tree_1d_dgsem/elixir_euler_ion_electron.jl +++ b/examples/tree_1d_dgsem/elixir_euler_multiion.jl @@ -3,12 +3,13 @@ using OrdinaryDiffEq ############################################################################### # semidiscretization of the compressible Euler equations -equations = Trixi.CompressibleEulerPlasmaEquations1D(1.4) +gammas = SVector(1.4, 1.76, 1.4) # Heat capacity ratios for two species +equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas) initial_condition = initial_condition_convergence_test -solver = DGSEM(polydeg = 3, surface_flux = flux_hllc, - volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hllc)) +solver = DGSEM(polydeg = 3, surface_flux = flux_hll, + volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hll)) coordinates_min = 0.0 coordinates_max = 2.0 diff --git a/src/equations/compressible_euler_multiion.jl b/src/equations/compressible_euler_multiion.jl new file mode 100644 index 0000000000..f04c41501c --- /dev/null +++ b/src/equations/compressible_euler_multiion.jl @@ -0,0 +1,424 @@ +# 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 + + +mutable struct CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT <: Real} <: + AbstractCompressibleEulerMultiIonEquations{1, NVARS, NCOMP} + gammas::SVector{NCOMP, RealT} # Heat capacity ratios + inv_gammas_minus_one::SVector{NCOMP, RealT} # = inv(gamma - 1); can be used to write slow divisions as fast multiplications + + function CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(gammas) where {NVARS, NCOMP, RealT <: Real} + + NCOMP >= 1 || + throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value")) + + inv_gammas_minus_one = SVector{NCOMP, RealT}(inv.(gammas .- 1)) + new(gammas, inv_gammas_minus_one) + end +end + +function CompressibleEulerMultiIonEquations1D(; gammas) + _gammas = promote(gammas...) + NCOMP = length(_gammas) + NVARS = 3 * NCOMP + RealT = promote_type(eltype(_gammas)) + + return CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(_gammas) +end + +# Outer constructor for `@reset` works correctly +function CompressibleEulerMultiIonEquations1D(gammas) + return CompressibleEulerMultiIonEquations1D(gammas = gammas) +end + +@inline function Base.real(::CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}) where {NVARS, NCOMP, RealT <: Real} + RealT +end + + +function varnames(::typeof(cons2cons), equations::CompressibleEulerMultiIonEquations1D) + cons = () + for i in eachcomponent(equations) + cons = (cons..., + tuple("rho_" * string(i), "rho_v1_" * string(i), "rho_e_" * string(i))...) + end + + return cons +end +function varnames(::typeof(cons2prim), equations::CompressibleEulerMultiIonEquations1D) + prim = () + for i in eachcomponent(equations) + prim = (prim..., + tuple("rho_" * string(i), "v1_" * string(i), "p_" * string(i))...) + end + return prim +end + +""" + get_component(k, u, equations::CompressibleEulerMultiIonEquations1D) + +Get the 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::CompressibleEulerMultiIonEquations1D) + return SVector(u[(k - 1) * 3 + 1], + u[(k - 1) * 3 + 2], + u[(k - 1) * 3 + 3]) +end + +""" + set_component!(u, k, u1, u2, u3, + equations::CompressibleEulerMultiIonEquations1D) + +Set the variables (`u1` to `u3`) 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, + equations::CompressibleEulerMultiIonEquations1D) + u[(k - 1) * 3 + 1] = u1 + u[(k - 1) * 3 + 2] = u2 + u[(k - 1) * 3 + 3] = u3 + + return u +end + +""" + flux(u, orientation::Integer, equations::CompressibleEulerMultiIonEquations1D) + +Calculate the flux for the multiion system. The flux is calculated for each ion species separately. +""" +@inline function flux(u, orientation::Integer, equations::CompressibleEulerMultiIonEquations1D) + f = zero(MVector{nvariables(equations), eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, rho_e = get_component(k, u, equations) + rho_inv = 1 / rho + v1 = rho_v1 * rho_inv + + gamma = equations.gammas[k] + p = (gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) + + f1 = rho_v1 + f2 = rho_v1 * v1 + p + f3 = (rho_e + p) * v1 + + set_component!(f, k, f1, f2, f3, equations) + end + + return SVector(f) +end + +""" + initial_condition_constant(x, t, equations::CompressibleEulerMultiIonEquations1D) + +A constant initial condition to test free-stream preservation. +""" +function initial_condition_constant(x, t, equations::CompressibleEulerMultiIonEquations1D) + cons = zero(MVector{nvariables(equations), eltype(x)}) + + rho = 0.1 + rho_v1 = 1 + rho_e = 10 + + for k in eachcomponent(equations) + set_component!(cons, k, rho, rho_v1, rho_e, equations) + end + + return SVector(cons) +end + +""" + initial_condition_convergence_test(x, t, equations::CompressibleEulerMultiIonEquations1D) + +A smooth initial condition used for convergence tests in combination with +[`source_terms_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +function initial_condition_convergence_test(x, t, equations::CompressibleEulerMultiIonEquations1D) + RealT = eltype(x) + cons = zero(MVector{nvariables(equations), RealT}) + + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f + ini = c + A * sin(ω * (x[1] - t)) + + rho = ini + rho_v1 = ini + rho_e = ini^2 + + for k in eachcomponent(equations) + set_component!(cons, k, rho, rho_v1, rho_e, equations) + end + + return SVector(cons) +end + + +""" + source_terms_convergence_test(u, x, t, equations::CompressibleEulerMultiIonEquations1D) + +Source terms used for convergence tests in combination with +[`initial_condition_convergence_test`](@ref) +(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). +""" +@inline function source_terms_convergence_test(u, x, t, + equations::CompressibleEulerMultiIonEquations1D) + # Same settings as in `initial_condition` + RealT = eltype(u) + cons = zero(MVector{nvariables(equations), RealT}) + + c = 2 + A = convert(RealT, 0.1) + L = 2 + f = 1.0f0 / L + ω = 2 * convert(RealT, pi) * f + + x1, = x + + si, co = sincos(ω * (x1 - t)) + + for k in eachcomponent(equations) + rho, _, _ = get_component(k, u, equations) + gamma = equations.gammas[k] + rho = c + A * si + rho_x = ω * A * co + + # Note that d/dt rho = -d/dx rho. + # This yields du2 = du3 = d/dx p (derivative of pressure). + # Other terms vanish because of v = 1. + du1 = 0 + du2 = rho_x * (2 * rho - 0.5f0) * (gamma - 1) + du3 = du2 + set_component!(cons, k, du1, du2, du3, equations) + end + + return SVector(cons) +end + +# Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the +# maximum velocity magnitude plus the maximum speed of sound +@inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerMultiIonEquations1D) + + # Calculate velocities + v_mag_ll = zero(eltype(u_ll)) + v_mag_rr = zero(eltype(u_rr)) + c_ll = zero(eltype(u_ll)) + c_rr = zero(eltype(u_rr)) + + for k in eachcomponent(equations) + rho_ll, rho_v1_ll, rho_e_ll = get_component(k, u_ll, equations) + rho_rr, rho_v1_rr, rho_e_rr = get_component(k, u_rr, equations) + gamma = equations.gammas[k] + + # Calculate primitive variables and speed of sound + v1_ll = rho_v1_ll / rho_ll + p_ll = (gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * abs(v1_ll)^2) + v_mag_ll = max(v_mag_ll,abs(v1_ll)) + c_ll = max(c_ll, sqrt(gamma * p_ll / rho_ll)) + + v1_rr = rho_v1_rr / rho_rr + p_rr = (gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * abs(v1_rr)^2) + v_mag_rr = max(v_mag_rr,abs(v1_rr)) + c_rr = max(c_rr, sqrt(gamma * p_rr / rho_rr)) + end + + λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr) +end + +@inline function max_abs_speeds(u, + equations::CompressibleEulerMultiIonEquations1D) + prim = cons2prim(u, equations) + v1_max = zero(eltype(u)) + for k in eachcomponent(equations) + rho, v1, p = get_component(k, prim, equations) + gamma = equations.gammas[k] + c = sqrt(gamma * p / rho) + v1_max = max(v1_max, abs(v1) + c) + end + return (v1_max,) +end + + +# Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerMultiIonEquations1D) + prim_ll = cons2prim(u_ll, equations) + prim_rr = cons2prim(u_rr, equations) + + λ_min = oftype(u_ll[1], Inf) + λ_max = oftype(u_ll[1], -Inf) + + for k in eachcomponent(equations) + rho_ll, v1_ll, p_ll = get_component(k, prim_ll, equations) + rho_rr, v1_rr, p_rr = get_component(k, prim_rr, equations) + gamma = equations.gammas[k] + + λ_min = min(λ_min, v1_ll - sqrt(gamma * p_ll / rho_ll)) + λ_max = max(λ_max, v1_rr + sqrt(gamma * p_rr / rho_rr)) + end + + #Assert that λ_min and λ_max are not Inf + @assert !isinf(λ_min) "λ_min is Inf" + @assert !isinf(λ_max) "λ_max is Inf" + + return λ_min, λ_max +end + +# More refined estimates for minimum and maximum wave speeds for HLL-type fluxes +@inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, + equations::CompressibleEulerMultiIonEquations1D) + prim_ll = cons2prim(u_ll, equations) + prim_rr = cons2prim(u_rr, equations) + + λ_min = oftype(u_ll[1], Inf) + λ_max = oftype(u_ll[1], -Inf) + + for k in eachcomponent(equations) + rho_ll, v1_ll, p_ll = get_component(k, prim_ll, equations) + rho_rr, v1_rr, p_rr = get_component(k, prim_rr, equations) + gamma = equations.gammas[k] + + c_ll = sqrt(gamma * p_ll / rho_ll) + c_rr = sqrt(gamma * p_rr / rho_rr) + + λ_min = min(λ_min, min(v1_ll - c_ll, v1_rr - c_rr)) + λ_max = max(λ_max, max(v1_ll + c_ll, v1_rr + c_rr)) + end + + return λ_min, λ_max +end + +@inline function cons2prim(u, equations::CompressibleEulerMultiIonEquations1D) + prim = zero(MVector{nvariables(equations), eltype(u)}) + for k in eachcomponent(equations) + rho, rho_v1, rho_e = get_component(k, u, equations) + v1 = rho_v1 / rho + p = (equations.gammas[k] - 1) * (rho_e - 0.5 * rho * v1^2) + set_component!(prim, k, rho, v1, p, equations) + end + + return SVector(prim) +end + +# Convert conservative variables to entropy +@inline function cons2entropy(u, equations::CompressibleEulerMultiIonEquations1D) + w = zero(MVector{nvariables(equations), eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, rho_e = get_component(k, u, equations) + gamma = equations.gammas[k] + inv_gamma_minus_one = equations.inv_gammas_minus_one[k] + v1 = rho_v1 / rho + v_square = v1^2 + p = (gamma - 1) * (rho_e - 0.5 * rho * v_square) + s = log(p) - gamma * log(rho) + rho_p = rho / p + + w1 = (gamma - s) * inv_gamma_minus_one - + 0.5 * rho_p * v_square + w2 = rho_p * v1 + w3 = -rho_p + set_component!(w, k, w1, w2, w3, equations) + end + + return SVector(w) +end + +@inline function entropy2cons(w, equations::CompressibleEulerMultiIonEquations1D) + cons = zero(MVector{nvariables(equations), eltype(w)}) + # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD + # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) + + # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986) + # instead of `-rho * s / (gamma - 1)` + for k in eachcomponent(equations) + w_k = get_component(k, w, equations) + V1, V2, V5 = w .* (gamma - 1) + + # specific entropy, eq. (53) + s = gamma - V1 + 0.5f0 * (V2^2) / V5 + + # eq. (52) + energy_internal = ((gamma - 1) / (-V5)^gamma)^(inv_gamma_minus_one) * + exp(-s * inv_gamma_minus_one) + + # eq. (51) + rho = -V5 * energy_internal + rho_v1 = V2 * energy_internal + rho_e = (1 - 0.5f0 * (V2^2) / V5) * energy_internal + set_component!(cons, k, rho, rho_v1, rho_e, equations) + end + return SVector(cons) +end + +# Convert primitive to conservative variables +@inline function single_species_prim2cons(prim, inv_gamma_minus_one) + rho, v1, p = prim + rho_v1 = rho * v1 + rho_e = p * inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1) + return SVector(rho, rho_v1, rho_e) +end + +@inline function prim2cons(prim, equations::CompressibleEulerMultiIonEquations1D) + cons = zero(MVector{nvariables(equations), eltype(prim)}) + + for k in eachcomponent(equations) + rho, v1, p = get_component(k, prim, equations) + inv_gamma_minus_one = equations.inv_gammas_minus_one[k] + rho_v1 = rho * v1 + rho_e = p * inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1) + set_component!(cons, k, rho, rho_v1, rho_e, equations) + end + + return SVector(cons) +end + +@inline function density(u, equations::CompressibleEulerMultiIonEquations1D) + num_components = div(nvariables(equations), 3) + rhos = zero(MVector{num_components, eltype(u)}) + + for k in eachcomponent(equations) + rho, _, _ = get_component(k, u, equations) + rhos[k] = rho + end + + return rhos +end + +@inline function velocity(u, equations::CompressibleEulerMultiIonEquations1D) + num_components = div(nvariables(equations), 3) + velocities = zero(MVector{num_components, eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, _ = get_component(k, u, equations) + velocities[k] = rho_v1 / rho + end + + return velocities +end + +@inline function pressure(u, equations::CompressibleEulerMultiIonEquations1D) + num_components = div(nvariables(equations), 3) + pressures = zero(MVector{num_components, eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, rho_e = get_component(k, u, equations) + p = (equations.gammas[k] - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) + pressures[k] = p + end + return pressure +end + +end # @muladd \ No newline at end of file diff --git a/src/equations/compressible_euler_plasma.jl b/src/equations/compressible_euler_plasma.jl deleted file mode 100644 index 69a33af7de..0000000000 --- a/src/equations/compressible_euler_plasma.jl +++ /dev/null @@ -1,323 +0,0 @@ -struct CompressibleEulerPlasmaEquations1D{RealT <: Real} <: - AbstractCompressibleEulerEquations{1, 6} - gamma::RealT # ratio of specific heats - inv_gamma_minus_one::RealT # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - - function CompressibleEulerPlasmaEquations1D(gamma) - γ, inv_gamma_minus_one = promote(gamma, inv(gamma - 1)) - new{typeof(γ)}(γ, inv_gamma_minus_one) - end -end - -function varnames(::typeof(cons2cons), ::CompressibleEulerPlasmaEquations1D) - ("rho_ion", "rho_v1_ion", "rho_e_ion", "rho_el", "rho_v1_el", "rho_e_el") -end -varnames(::typeof(cons2prim), ::CompressibleEulerPlasmaEquations1D) = ("rho_ion", "v1_ion", "p_ion", "rho_el", "v1_el", "p_el") - -""" - single_species_flux(u, gamma) - -Calculate the flux vector for a single species (ion or electron). -""" -@inline function single_species_flux(rho, rho_v1, rho_e, gamma) - v1 = rho_v1 / rho - p = (gamma - 1) * (rho_e - 0.5 * rho_v1 * v1) - return SVector(rho_v1, rho_v1 * v1 + p, (rho_e + p) * v1) -end - -""" - flux(u, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) - -Calculate the flux for the full plasma system with ions and electrons. -""" -@inline function flux(u, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) - rho_ion, rho_v1_ion, rho_e_ion, rho_el, rho_v1_el, rho_e_el = u - - f_ion = single_species_flux(rho_ion, rho_v1_ion, rho_e_ion, equations.gamma) - - f_el = single_species_flux(rho_el, rho_v1_el, rho_e_el, equations.gamma) - - return SVector(f_ion[1], f_ion[2], f_ion[3], f_el[1], f_el[2], f_el[3]) -end - -""" - initial_condition_constant(x, t, equations::CompressibleEulerPlasmaEquations1D) - -A constant initial condition to test free-stream preservation. -""" -function initial_condition_constant(x, t, equations::CompressibleEulerPlasmaEquations1D) - RealT = eltype(x) - rho_ion = 1 - rho_el = 1 - - # Ion and electron velocities - rho_v1_ion = convert(RealT, 0.01) - rho_v1_el = convert(RealT, 0.1) - - # Ion and electron internal energies - rho_e_ion = 10.0 - rho_e_el = 10.0 - - return SVector(rho_ion, rho_v1_ion, rho_e_ion, rho_el, rho_v1_el, rho_e_el) -end - -""" - initial_condition_convergence_test(x, t, equations::CompressibleEulerPlasmaEquations1D) - -A smooth initial condition used for convergence tests in combination with -[`source_terms_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). -""" -function initial_condition_convergence_test(x, t, equations::CompressibleEulerPlasmaEquations1D) - RealT = eltype(x) - c = 2 - A = convert(RealT, 0.1) - L = 2 - f = 1.0f0 / L - ω = 2 * convert(RealT, pi) * f - ini = c + A * sin(ω * (x[1] - t)) - - rho = ini - rho_v1 = ini - rho_e = ini^2 - - return SVector(rho, rho_v1, rho_e, rho, rho_v1, rho_e) -end - - -""" - source_terms_convergence_test(u, x, t, equations::CompressibleEulerPlasmaEquations1D) - -Source terms used for convergence tests in combination with -[`initial_condition_convergence_test`](@ref) -(and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). -""" -@inline function source_terms_convergence_test(u, x, t, - equations::CompressibleEulerPlasmaEquations1D) - # Same settings as in `initial_condition` - RealT = eltype(u) - c = 2 - A = convert(RealT, 0.1) - L = 2 - f = 1.0f0 / L - ω = 2 * convert(RealT, pi) * f - γ = equations.gamma - - x1, = x - - si, co = sincos(ω * (x1 - t)) - rho = c + A * si - rho_x = ω * A * co - - # Note that d/dt rho = -d/dx rho. - # This yields du2 = du3 = d/dx p (derivative of pressure). - # Other terms vanish because of v = 1. - du1 = 0 - du2 = rho_x * (2 * rho - 0.5f0) * (γ - 1) - du3 = du2 - - return SVector(du1, du2, du3, du1, du2, du3) -end - -#TODO: Maybe there is a better way to do this idk -function single_species_flux_hllc(u_ll, u_rr, orientation::Integer, - gamma) - # Calculate primitive variables and speed of sound - rho_ll, rho_v1_ll, rho_e_ll = u_ll - rho_rr, rho_v1_rr, rho_e_rr = u_rr - - v1_ll = rho_v1_ll / rho_ll - e_ll = rho_e_ll / rho_ll - p_ll = (gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * v1_ll^2) - c_ll = sqrt(gamma * p_ll / rho_ll) - - v1_rr = rho_v1_rr / rho_rr - e_rr = rho_e_rr / rho_rr - p_rr = (gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * v1_rr^2) - c_rr = sqrt(gamma * p_rr / rho_rr) - - # Obtain left and right fluxes - f_ll = single_species_flux(rho_ll, rho_v1_ll, rho_e_ll, gamma) - f_rr = single_species_flux(rho_rr, rho_v1_rr, rho_e_rr, gamma) - - # Compute Roe averages - sqrt_rho_ll = sqrt(rho_ll) - sqrt_rho_rr = sqrt(rho_rr) - sum_sqrt_rho = sqrt_rho_ll + sqrt_rho_rr - vel_L = v1_ll - vel_R = v1_rr - vel_roe = (sqrt_rho_ll * vel_L + sqrt_rho_rr * vel_R) / sum_sqrt_rho - ekin_roe = 0.5f0 * vel_roe^2 - H_ll = (rho_e_ll + p_ll) / rho_ll - H_rr = (rho_e_rr + p_rr) / rho_rr - H_roe = (sqrt_rho_ll * H_ll + sqrt_rho_rr * H_rr) / sum_sqrt_rho - c_roe = sqrt((gamma - 1) * (H_roe - ekin_roe)) - - Ssl = min(vel_L - c_ll, vel_roe - c_roe) - Ssr = max(vel_R + c_rr, vel_roe + c_roe) - sMu_L = Ssl - vel_L - sMu_R = Ssr - vel_R - if Ssl >= 0 - f1 = f_ll[1] - f2 = f_ll[2] - f3 = f_ll[3] - elseif Ssr <= 0 - f1 = f_rr[1] - f2 = f_rr[2] - f3 = f_rr[3] - else - SStar = (p_rr - p_ll + rho_ll * vel_L * sMu_L - rho_rr * vel_R * sMu_R) / - (rho_ll * sMu_L - rho_rr * sMu_R) - if Ssl <= 0 <= SStar - densStar = rho_ll * sMu_L / (Ssl - SStar) - enerStar = e_ll + (SStar - vel_L) * (SStar + p_ll / (rho_ll * sMu_L)) - UStar1 = densStar - UStar2 = densStar * SStar - UStar3 = densStar * enerStar - - f1 = f_ll[1] + Ssl * (UStar1 - rho_ll) - f2 = f_ll[2] + Ssl * (UStar2 - rho_v1_ll) - f3 = f_ll[3] + Ssl * (UStar3 - rho_e_ll) - else - densStar = rho_rr * sMu_R / (Ssr - SStar) - enerStar = e_rr + (SStar - vel_R) * (SStar + p_rr / (rho_rr * sMu_R)) - UStar1 = densStar - UStar2 = densStar * SStar - UStar3 = densStar * enerStar - - #end - f1 = f_rr[1] + Ssr * (UStar1 - rho_rr) - f2 = f_rr[2] + Ssr * (UStar2 - rho_v1_rr) - f3 = f_rr[3] + Ssr * (UStar3 - rho_e_rr) - end - end - return SVector(f1, f2, f3) -end - -""" - flux_hllc(u_ll, u_rr, orientation, equations::CompressibleEulerPlasmaEquations1D) - -Computes the HLLC flux (HLL with Contact) for compressible Euler equations developed by E.F. Toro -[Lecture slides](http://www.prague-sum.com/download/2012/Toro_2-HLLC-RiemannSolver.pdf) -Signal speeds: [DOI: 10.1137/S1064827593260140](https://doi.org/10.1137/S1064827593260140) -""" -function flux_hllc(u_ll, u_rr, orientation::Integer, equations::CompressibleEulerPlasmaEquations1D) - flux_hllc_ion = single_species_flux_hllc(u_ll[1:3], u_rr[1:3], orientation, equations.gamma) - flux_hllc_el = single_species_flux_hllc(u_ll[4:6], u_rr[4:6], orientation, equations.gamma) - return vcat(flux_hllc_ion, flux_hllc_el) -end - -@inline function max_abs_speeds(u, equations::CompressibleEulerPlasmaEquations1D) - rho_ion, v1_ion, p_ion, rho_el, v1_el, p_el = cons2prim(u, equations) - c_ion = sqrt(equations.gamma * p_ion / rho_ion) - c_el = sqrt(equations.gamma * p_el / rho_el) - - return (abs(v1_ion) + c_ion, abs(v1_el) + c_el) # Is this the right way? I've seen multicomponent where (abs(v1) + c,) is returned -end - -@inline function single_species_cons2prim(u, gamma) - rho, rho_v1, rho_e = u - v1 = rho_v1 / rho - p = (gamma - 1) * (rho_e - 0.5f0 * rho * v1^2) - return SVector(rho, v1, p) -end - -@inline function cons2prim(u, equations::CompressibleEulerPlasmaEquations1D) - prim_ion = single_species_cons2prim(u[1:3], equations.gamma) - - prim_el = single_species_cons2prim(u[4:6], equations.gamma) - - return vcat(prim_ion, prim_el) -end - -@inline function single_species_cons2entropy(u, gamma, inv_gamma_minus_one) - rho, rho_v1, rho_e = u - - v1 = rho_v1 / rho - v_square = v1^2 - p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square) - s = log(p) - gamma * log(rho) - rho_p = rho / p - - w1 = (gamma - s) * inv_gamma_minus_one - - 0.5f0 * rho_p * v_square - w2 = rho_p * v1 - w3 = -rho_p - - return SVector(w1, w2, w3) -end - -# Convert conservative variables to entropy -@inline function cons2entropy(u, equations::CompressibleEulerPlasmaEquations1D) - w_ion = single_species_cons2entropy(u[1:3], equations.gamma, equations.inv_gamma_minus_one) - w_el = single_species_cons2entropy(u[4:6], equations.gamma, equations.inv_gamma_minus_one) - return vcat(w_ion, w_el) -end - -@inline function single_species_entropy2cons(w, gamma, inv_gamma_minus_one) - # See Hughes, Franca, Mallet (1986) A new finite element formulation for CFD - # [DOI: 10.1016/0045-7825(86)90127-1](https://doi.org/10.1016/0045-7825(86)90127-1) - - # convert to entropy `-rho * s` used by Hughes, France, Mallet (1986) - # instead of `-rho * s / (gamma - 1)` - V1, V2, V5 = w .* (gamma - 1) - - # specific entropy, eq. (53) - s = gamma - V1 + 0.5f0 * (V2^2) / V5 - - # eq. (52) - energy_internal = ((gamma - 1) / (-V5)^gamma)^(inv_gamma_minus_one) * - exp(-s * inv_gamma_minus_one) - - # eq. (51) - rho = -V5 * energy_internal - rho_v1 = V2 * energy_internal - rho_e = (1 - 0.5f0 * (V2^2) / V5) * energy_internal - return SVector(rho, rho_v1, rho_e) -end - -@inline function entropy2cons(w, equations::CompressibleEulerPlasmaEquations1D) - u_ion = single_species_entropy2cons(w[1:3], equations.gamma, equations.inv_gamma_minus_one) - u_el = single_species_entropy2cons(w[4:6], equations.gamma, equations.inv_gamma_minus_one) - return vcat(u_ion, u_el) -end - -# Convert primitive to conservative variables -@inline function single_species_prim2cons(prim, inv_gamma_minus_one) - rho, v1, p = prim - rho_v1 = rho * v1 - rho_e = p * inv_gamma_minus_one + 0.5f0 * (rho_v1 * v1) - return SVector(rho, rho_v1, rho_e) -end - -@inline function prim2cons(prim, equations::CompressibleEulerPlasmaEquations1D) - u_ion = single_species_prim2cons(prim[1:3], equations.inv_gamma_minus_one) - u_el = single_species_prim2cons(prim[4:6], equations.inv_gamma_minus_one) - return vcat(u_ion, u_el) -end - -@inline function density(u, equations::CompressibleEulerPlasmaEquations1D) - rho_ion = u[1] - rho_el = u[4] - return rho_ion, rho_el -end - -@inline function velocity(u, equations::CompressibleEulerPlasmaEquations1D) - rho_ion = u[1] - rho_el = u[4] - v1_ion = u[2] / rho_ion - v1_el = u[5] / rho_el - return v1_ion, v1_el -end - -@inline function single_species_pressure(u) - rho, rho_v1, rho_e = u - p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2) / rho) - return p -end - -@inline function pressure(u, equations::CompressibleEulerPlasmaEquations1D) - p_ion = single_species_pressure(u[1:3]) - p_el = single_species_pressure(u[4:6]) - return p_ion, p_el -end \ No newline at end of file diff --git a/src/equations/equations.jl b/src/equations/equations.jl index 4051b56d65..a375347b82 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -442,7 +442,10 @@ include("compressible_euler_1d.jl") include("compressible_euler_2d.jl") include("compressible_euler_3d.jl") include("compressible_euler_quasi_1d.jl") -include("compressible_euler_plasma.jl") + +abstract type AbstractCompressibleEulerMultiIonEquations{NDIMS, NVARS, NCOMP} <: + AbstractEquations{NDIMS, NVARS} end +include("compressible_euler_multiion.jl") # CompressibleEulerMulticomponentEquations abstract type AbstractCompressibleEulerMulticomponentEquations{NDIMS, NVARS, NCOMP} <: @@ -476,6 +479,27 @@ In particular, not the components themselves are returned. Base.OneTo(ncomponents(equations)) end +# Retrieve number of components from equation instance for the multicomponent case +@inline function ncomponents(::AbstractCompressibleEulerMultiIonEquations{NDIMS, + NVARS, + NCOMP}) where { + NDIMS, + NVARS, + NCOMP + } + NCOMP +end +""" + eachcomponent(equations::AbstractCompressibleEulerMultiIonEquationsEquations) + +Return an iterator over the indices that specify the location in relevant data structures +for the components in `AbstractCompressibleEulerMultiIonEquationsEquations`. +In particular, not the components themselves are returned. +""" +@inline function eachcomponent(equations::AbstractCompressibleEulerMultiIonEquations) + Base.OneTo(ncomponents(equations)) +end + # Ideal MHD abstract type AbstractIdealGlmMhdEquations{NDIMS, NVARS} <: AbstractEquations{NDIMS, NVARS} end From 53705b5e5125a218a2d79d128c44ea0f573e4aa6 Mon Sep 17 00:00:00 2001 From: Warisa Date: Tue, 14 Jan 2025 17:35:54 +0100 Subject: [PATCH 3/6] fmt and a comment --- .../tree_1d_dgsem/elixir_euler_multiion.jl | 4 +- src/equations/compressible_euler_multiion.jl | 62 +++++++++++-------- src/equations/equations.jl | 14 ++--- 3 files changed, 44 insertions(+), 36 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_euler_multiion.jl b/examples/tree_1d_dgsem/elixir_euler_multiion.jl index a1c7e2cd64..643d0c5b08 100644 --- a/examples/tree_1d_dgsem/elixir_euler_multiion.jl +++ b/examples/tree_1d_dgsem/elixir_euler_multiion.jl @@ -4,7 +4,7 @@ using OrdinaryDiffEq ############################################################################### # semidiscretization of the compressible Euler equations gammas = SVector(1.4, 1.76, 1.4) # Heat capacity ratios for two species -equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas) +equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas) # Since I haven't exported this yet initial_condition = initial_condition_convergence_test @@ -53,4 +53,4 @@ callbacks = CallbackSet(summary_callback, 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 \ No newline at end of file +summary_callback() # print the timer summary diff --git a/src/equations/compressible_euler_multiion.jl b/src/equations/compressible_euler_multiion.jl index f04c41501c..25572d722e 100644 --- a/src/equations/compressible_euler_multiion.jl +++ b/src/equations/compressible_euler_multiion.jl @@ -5,17 +5,20 @@ @muladd begin #! format: noindent - mutable struct CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT <: Real} <: - AbstractCompressibleEulerMultiIonEquations{1, NVARS, NCOMP} + AbstractCompressibleEulerMultiIonEquations{1, NVARS, NCOMP} gammas::SVector{NCOMP, RealT} # Heat capacity ratios inv_gammas_minus_one::SVector{NCOMP, RealT} # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - function CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(gammas) where {NVARS, NCOMP, RealT <: Real} - + function CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(gammas) where { + NVARS, + NCOMP, + RealT <: + Real + } NCOMP >= 1 || throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value")) - + inv_gammas_minus_one = SVector{NCOMP, RealT}(inv.(gammas .- 1)) new(gammas, inv_gammas_minus_one) end @@ -35,16 +38,21 @@ function CompressibleEulerMultiIonEquations1D(gammas) return CompressibleEulerMultiIonEquations1D(gammas = gammas) end -@inline function Base.real(::CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}) where {NVARS, NCOMP, RealT <: Real} +@inline function Base.real(::CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}) where { + NVARS, + NCOMP, + RealT <: + Real + } RealT end - function varnames(::typeof(cons2cons), equations::CompressibleEulerMultiIonEquations1D) cons = () for i in eachcomponent(equations) cons = (cons..., - tuple("rho_" * string(i), "rho_v1_" * string(i), "rho_e_" * string(i))...) + tuple("rho_" * string(i), "rho_v1_" * string(i), + "rho_e_" * string(i))...) end return cons @@ -95,7 +103,8 @@ end Calculate the flux for the multiion system. The flux is calculated for each ion species separately. """ -@inline function flux(u, orientation::Integer, equations::CompressibleEulerMultiIonEquations1D) +@inline function flux(u, orientation::Integer, + equations::CompressibleEulerMultiIonEquations1D) f = zero(MVector{nvariables(equations), eltype(u)}) for k in eachcomponent(equations) @@ -121,9 +130,10 @@ end A constant initial condition to test free-stream preservation. """ -function initial_condition_constant(x, t, equations::CompressibleEulerMultiIonEquations1D) +function initial_condition_constant(x, t, + equations::CompressibleEulerMultiIonEquations1D) cons = zero(MVector{nvariables(equations), eltype(x)}) - + rho = 0.1 rho_v1 = 1 rho_e = 10 @@ -142,10 +152,11 @@ A smooth initial condition used for convergence tests in combination with [`source_terms_convergence_test`](@ref) (and [`BoundaryConditionDirichlet(initial_condition_convergence_test)`](@ref) in non-periodic domains). """ -function initial_condition_convergence_test(x, t, equations::CompressibleEulerMultiIonEquations1D) +function initial_condition_convergence_test(x, t, + equations::CompressibleEulerMultiIonEquations1D) RealT = eltype(x) cons = zero(MVector{nvariables(equations), RealT}) - + c = 2 A = convert(RealT, 0.1) L = 2 @@ -164,7 +175,6 @@ function initial_condition_convergence_test(x, t, equations::CompressibleEulerMu return SVector(cons) end - """ source_terms_convergence_test(u, x, t, equations::CompressibleEulerMultiIonEquations1D) @@ -209,8 +219,8 @@ end # Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the # maximum velocity magnitude plus the maximum speed of sound @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMultiIonEquations1D) - + equations::CompressibleEulerMultiIonEquations1D) + # Calculate velocities v_mag_ll = zero(eltype(u_ll)) v_mag_rr = zero(eltype(u_rr)) @@ -225,12 +235,12 @@ end # Calculate primitive variables and speed of sound v1_ll = rho_v1_ll / rho_ll p_ll = (gamma - 1) * (rho_e_ll - 0.5f0 * rho_ll * abs(v1_ll)^2) - v_mag_ll = max(v_mag_ll,abs(v1_ll)) + v_mag_ll = max(v_mag_ll, abs(v1_ll)) c_ll = max(c_ll, sqrt(gamma * p_ll / rho_ll)) v1_rr = rho_v1_rr / rho_rr p_rr = (gamma - 1) * (rho_e_rr - 0.5f0 * rho_rr * abs(v1_rr)^2) - v_mag_rr = max(v_mag_rr,abs(v1_rr)) + v_mag_rr = max(v_mag_rr, abs(v1_rr)) c_rr = max(c_rr, sqrt(gamma * p_rr / rho_rr)) end @@ -238,7 +248,7 @@ end end @inline function max_abs_speeds(u, - equations::CompressibleEulerMultiIonEquations1D) + equations::CompressibleEulerMultiIonEquations1D) prim = cons2prim(u, equations) v1_max = zero(eltype(u)) for k in eachcomponent(equations) @@ -250,10 +260,9 @@ end return (v1_max,) end - # Calculate estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_naive(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMultiIonEquations1D) + equations::CompressibleEulerMultiIonEquations1D) prim_ll = cons2prim(u_ll, equations) prim_rr = cons2prim(u_rr, equations) @@ -268,7 +277,7 @@ end λ_min = min(λ_min, v1_ll - sqrt(gamma * p_ll / rho_ll)) λ_max = max(λ_max, v1_rr + sqrt(gamma * p_rr / rho_rr)) end - + #Assert that λ_min and λ_max are not Inf @assert !isinf(λ_min) "λ_min is Inf" @assert !isinf(λ_max) "λ_max is Inf" @@ -278,7 +287,7 @@ end # More refined estimates for minimum and maximum wave speeds for HLL-type fluxes @inline function min_max_speed_davis(u_ll, u_rr, orientation::Integer, - equations::CompressibleEulerMultiIonEquations1D) + equations::CompressibleEulerMultiIonEquations1D) prim_ll = cons2prim(u_ll, equations) prim_rr = cons2prim(u_rr, equations) @@ -327,7 +336,7 @@ end rho_p = rho / p w1 = (gamma - s) * inv_gamma_minus_one - - 0.5 * rho_p * v_square + 0.5 * rho_p * v_square w2 = rho_p * v1 w3 = -rho_p set_component!(w, k, w1, w2, w3, equations) @@ -352,7 +361,7 @@ end # eq. (52) energy_internal = ((gamma - 1) / (-V5)^gamma)^(inv_gamma_minus_one) * - exp(-s * inv_gamma_minus_one) + exp(-s * inv_gamma_minus_one) # eq. (51) rho = -V5 * energy_internal @@ -420,5 +429,4 @@ end end return pressure end - -end # @muladd \ No newline at end of file +end # @muladd diff --git a/src/equations/equations.jl b/src/equations/equations.jl index a375347b82..07705ab14e 100644 --- a/src/equations/equations.jl +++ b/src/equations/equations.jl @@ -444,7 +444,7 @@ include("compressible_euler_3d.jl") include("compressible_euler_quasi_1d.jl") abstract type AbstractCompressibleEulerMultiIonEquations{NDIMS, NVARS, NCOMP} <: - AbstractEquations{NDIMS, NVARS} end + AbstractEquations{NDIMS, NVARS} end include("compressible_euler_multiion.jl") # CompressibleEulerMulticomponentEquations @@ -481,12 +481,12 @@ end # Retrieve number of components from equation instance for the multicomponent case @inline function ncomponents(::AbstractCompressibleEulerMultiIonEquations{NDIMS, - NVARS, - NCOMP}) where { - NDIMS, - NVARS, - NCOMP - } + NVARS, + NCOMP}) where { + NDIMS, + NVARS, + NCOMP + } NCOMP end """ From 3fcd4f8b5b376cb888431cdeeb8605cec0f99339 Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 20 Jan 2025 09:49:57 +0100 Subject: [PATCH 4/6] add collision source terms --- .../tree_1d_dgsem/elixir_euler_multiion.jl | 99 +++++- src/equations/compressible_euler_multiion.jl | 288 ++++++++++++++++-- 2 files changed, 363 insertions(+), 24 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_euler_multiion.jl b/examples/tree_1d_dgsem/elixir_euler_multiion.jl index 643d0c5b08..78086c69b3 100644 --- a/examples/tree_1d_dgsem/elixir_euler_multiion.jl +++ b/examples/tree_1d_dgsem/elixir_euler_multiion.jl @@ -2,21 +2,110 @@ using Trixi using OrdinaryDiffEq ############################################################################### -# semidiscretization of the compressible Euler equations -gammas = SVector(1.4, 1.76, 1.4) # Heat capacity ratios for two species -equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas) # Since I haven't exported this yet +# This elixir describes the frictional slowing of an ionized carbon fluid (C6+) with respect to another species +# of a background ionized carbon fluid with an initially nonzero relative velocity. It is the second slow-down +# test (fluids with different densities) described in: +# - Ghosh, D., Chapman, T. D., Berger, R. L., Dimits, A., & Banks, J. W. (2019). A +# multispecies, multifluid model for laser–induced counterstreaming plasma simulations. +# Computers & Fluids, 186, 38-57. +# +# This is effectively a zero-dimensional case because the spatial gradients are zero, and we use it to test the +# collision source terms. +# +# To run this physically relevant test, we use the following characteristic quantities to non-dimensionalize +# the equations: +# Characteristic length: L_inf = 1.00E-03 m (domain size) +# Characteristic density: rho_inf = 1.99E+00 kg/m^3 (corresponds to a number density of 1e20 cm^{-3}) +# Characteristic vacuum permeability: mu0_inf = 1.26E-06 N/A^2 (for equations with mu0 = 1) +# Characteristic gas constant: R_inf = 6.92237E+02 J/kg/K (specific gas constant for a Carbon fluid) +# Characteristic velocity: V_inf = 1.00E+06 m/s +# +# The results of the paper can be reproduced using `source_terms = source_terms_collision_ion_ion` (i.e., only +# taking into account ion-ion collisions). However, we include ion-electron collisions assuming a constant +# electron temperature of 1 keV in this elixir to test the function `source_terms_collision_ion_electron` -initial_condition = initial_condition_convergence_test +# Return the electron pressure for a constant electron temperature Te = 1 keV +function electron_pressure_constantTe(u, equations::Trixi.CompressibleEulerMultiIonEquations1D) + @unpack charge_to_mass = equations + Te = 0.008029953773 # [nondimensional] = 1 [keV] + total_electron_charge = zero(eltype(u)) + for k in eachcomponent(equations) + rho_k = u[(k - 1) * 3 + 1] + total_electron_charge += rho_k * charge_to_mass[k] + end + + # Boltzmann constant divided by elementary charge + kB_e = 7.86319034E-02 #[nondimensional] + + return total_electron_charge * kB_e * Te +end + +# Return the constant electron temperature Te = 1 keV +function electron_temperature_constantTe(u, equations::Trixi.CompressibleEulerMultiIonEquations1D) + return 0.008029953773 # [nondimensional] = 1 [keV] +end + +equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas = (5 / 3, 5 / 3), + charge_to_mass = (76.3049060157692000, + 76.3049060157692000), # [nondimensional] + gas_constants = (1.0, 1.0), # [nondimensional] + molar_masses = (1.0, 1.0), # [nondimensional] + ion_ion_collision_constants = [0.0 0.4079382480442680; + 0.4079382480442680 0.0], # [nondimensional] (computed with eq (4.142) of Schunk & Nagy (2009)) + ion_electron_collision_constants = (8.56368379833E-06, + 8.56368379833E-06), # [nondimensional] (computed with eq (9) of Ghosh et al. (2019)) + electron_pressure = electron_pressure_constantTe, + electron_temperature = electron_temperature_constantTe) + +# Frictional slowing of an ionized carbon fluid with respect to another background carbon fluid in motion +function initial_condition_slow_down(x, t, equations::Trixi.CompressibleEulerMultiIonEquations1D) + v11 = 0.65508770000000 + v21 = 0.0 + rho1 = 0.1 + rho2 = 1.0 + + p1 = 0.00040170535986 + p2 = 0.00401705359856 + + return prim2cons(SVector(rho1, v11, p1, rho2, v21, p2), + equations) +end + +# Temperature of ion 1 +function temperature1(u, equations::Trixi.CompressibleEulerMultiIonEquations1D) + rho_1, _ = Trixi.get_component(1, u, equations) + p = pressure(u, equations) + + return p[1] / (rho_1 * equations.gas_constants[1]) +end + +# Temperature of ion 2 +function temperature2(u, equations::Trixi.CompressibleEulerMultiIonEquations1D) + rho_2, _ = Trixi.get_component(2, u, equations) + p = pressure(u, equations) + + return p[2] / (rho_2 * equations.gas_constants[2]) +end + +initial_condition = initial_condition_slow_down +tspan = (0.0, 0.1) # 100 [ps] solver = DGSEM(polydeg = 3, surface_flux = flux_hll, volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hll)) coordinates_min = 0.0 -coordinates_max = 2.0 +coordinates_max = 1.0 mesh = TreeMesh(coordinates_min, coordinates_max, initial_refinement_level = 4, n_cells_max = 10_000) +# Ion-ion and ion-electron collision source terms +# In this particular case, we can omit source_terms_lorentz because the magnetic field is zero! +function source_terms(u, x, t, equations::Trixi.CompressibleEulerMultiIonEquations1D) + source_terms_collision_ion_ion(u, x, t, equations) + + source_terms_collision_ion_electron(u, x, t, equations) +end + semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, source_terms = source_terms_convergence_test) diff --git a/src/equations/compressible_euler_multiion.jl b/src/equations/compressible_euler_multiion.jl index 25572d722e..8684a73cee 100644 --- a/src/equations/compressible_euler_multiion.jl +++ b/src/equations/compressible_euler_multiion.jl @@ -5,37 +5,100 @@ @muladd begin #! format: noindent -mutable struct CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT <: Real} <: +mutable struct CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT <: Real, + ElectronPressure, + ElectronTemperature} <: AbstractCompressibleEulerMultiIonEquations{1, NVARS, NCOMP} gammas::SVector{NCOMP, RealT} # Heat capacity ratios - inv_gammas_minus_one::SVector{NCOMP, RealT} # = inv(gamma - 1); can be used to write slow divisions as fast multiplications - - function CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(gammas) where { - NVARS, - NCOMP, - RealT <: - Real - } + inv_gammas_minus_one::SVector{NCOMP, RealT} # = inv(gamma - 1) + charge_to_mass::SVector{NCOMP, RealT} # Charge to mass ratios + gas_constants::SVector{NCOMP, RealT} # Specific gas constants + molar_masses::SVector{NCOMP, RealT} # Molar masses + ion_ion_collision_constants::Array{RealT, 2} # Collision coefficients + ion_electron_collision_constants::SVector{NCOMP, RealT} # Ion-electron collision constants + electron_pressure::ElectronPressure # Electron pressure function + electron_temperature::ElectronTemperature # Electron temperature function + + # Inner Constructor + function CompressibleEulerMultiIonEquations1D(gammas::SVector{NCOMP, RealT}, + charge_to_mass::SVector{NCOMP, RealT}, + gas_constants::SVector{NCOMP, RealT}, + molar_masses::SVector{NCOMP, RealT}, + ion_ion_collision_constants::Array{RealT, + 2}, + ion_electron_collision_constants::SVector{NCOMP, + RealT}, + electron_pressure::ElectronPressure, + electron_temperature::ElectronTemperature) where { + NCOMP, + RealT <: + Real, + ElectronPressure, + ElectronTemperature + } NCOMP >= 1 || - throw(DimensionMismatch("`gammas` and `charge_to_mass` have to be filled with at least one value")) + throw(DimensionMismatch("`gammas` and `charge_to_mass` must contain at least one value")) + # Precompute inverse gamma - 1 inv_gammas_minus_one = SVector{NCOMP, RealT}(inv.(gammas .- 1)) - new(gammas, inv_gammas_minus_one) + + NVARS = 3 * NCOMP + + return new{NVARS, NCOMP, RealT, ElectronPressure, ElectronTemperature}(gammas, + inv_gammas_minus_one, + charge_to_mass, + gas_constants, + molar_masses, + ion_ion_collision_constants, + ion_electron_collision_constants, + electron_pressure, + electron_temperature) end end -function CompressibleEulerMultiIonEquations1D(; gammas) +# Outer Constructor Delegating to Inner Constructor +function CompressibleEulerMultiIonEquations1D(; gammas, charge_to_mass, + gas_constants = zero(SVector{length(gammas), + eltype(gammas)}), + molar_masses = zero(SVector{length(gammas), + eltype(gammas)}), + ion_ion_collision_constants = zeros(eltype(gammas), + length(gammas), + length(gammas)), + ion_electron_collision_constants = zero(SVector{length(gammas), + eltype(gammas)}), + electron_pressure = electron_pressure_zero, + electron_temperature = electron_pressure_zero) + # Promote input types _gammas = promote(gammas...) + _charge_to_mass = promote(charge_to_mass...) + _gas_constants = promote(gas_constants...) + _molar_masses = promote(molar_masses...) + _ion_electron_collision_constants = promote(ion_electron_collision_constants...) + RealT = promote_type(eltype(_gammas), eltype(_charge_to_mass), + eltype(_gas_constants), eltype(_molar_masses), + eltype(ion_ion_collision_constants), + eltype(_ion_electron_collision_constants)) + __gammas = SVector(map(RealT, _gammas)) + __charge_to_mass = SVector(map(RealT, _charge_to_mass)) + __gas_constants = SVector(map(RealT, _gas_constants)) + __molar_masses = SVector(map(RealT, _molar_masses)) + __ion_ion_collision_constants = map(RealT, ion_ion_collision_constants) + __ion_electron_collision_constants = SVector(map(RealT, + _ion_electron_collision_constants)) NCOMP = length(_gammas) NVARS = 3 * NCOMP - RealT = promote_type(eltype(_gammas)) - - return CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}(_gammas) -end -# Outer constructor for `@reset` works correctly -function CompressibleEulerMultiIonEquations1D(gammas) - return CompressibleEulerMultiIonEquations1D(gammas = gammas) + return CompressibleEulerMultiIonEquations1D( + __gammas, + __charge_to_mass, + __gas_constants, + __molar_masses, + __ion_ion_collision_constants, + __ion_electron_collision_constants, + electron_pressure, + electron_temperature + ) end @inline function Base.real(::CompressibleEulerMultiIonEquations1D{NVARS, NCOMP, RealT}) where { @@ -66,6 +129,35 @@ function varnames(::typeof(cons2prim), equations::CompressibleEulerMultiIonEquat return prim end +""" + v1, vk1 = charge_averaged_velocities(u, equations::CompressibleEulerMultiIonEquations1D) + + +Compute the charge-averaged velocities (`v1`) and each ion species' contribution +to the charge-averaged velocities (`vk1`). The output variables `vk1` +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::CompressibleEulerMultiIonEquations1D) + total_electron_charge = zero(real(equations)) + + vk1_plus = zero(MVector{ncomponents(equations), eltype(u)}) + + for k in eachcomponent(equations) + rho, rho_v1, _ = get_component(k, u, equations) + + total_electron_charge += rho * equations.charge_to_mass[k] + vk1_plus[k] = rho_v1 * equations.charge_to_mass[k] + end + vk1_plus ./= total_electron_charge + v1_plus = sum(vk1_plus) + + return v1_plus, SVector(vk1_plus) +end + """ get_component(k, u, equations::CompressibleEulerMultiIonEquations1D) @@ -216,6 +308,164 @@ Source terms used for convergence tests in combination with return SVector(cons) end +@doc raw""" + source_terms_collision_ion_ion(u, x, t, + equations::CompressibleEulerMultiIonEquations1D) + +Compute the ion-ion collision source terms for the momentum and energy equations of each ion species as +```math +\begin{aligned} + \vec{s}_{\rho_k \vec{v}_k} =& \rho_k\sum_{k'}\bar{\nu}_{kk'}(\vec{v}_{k'} - \vec{v}_k),\\ + s_{E_k} =& + 3 \sum_{k'} \left( + \bar{\nu}_{kk'} \frac{\rho_k M_1}{M_{k'} + M_k} R_1 (T_{k'} - T_k) + \right) + + \sum_{k'} \left( + \bar{\nu}_{kk'} \rho_k \frac{M_{k'}}{M_{k'} + M_k} \|\vec{v}_{k'} - \vec{v}_k\|^2 + \right) + + + \vec{v}_k \cdot \vec{s}_{\rho_k \vec{v}_k}, +\end{aligned} +``` +where ``M_k`` is the molar mass of ion species `k` provided in `equations.molar_masses`, +``R_k`` is the specific gas constant of ion species `k` provided in `equations.gas_constants`, and + ``\bar{\nu}_{kk'}`` is the effective collision frequency of species `k` with species `k'`, which is computed as +```math +\begin{aligned} + \bar{\nu}_{kk'} = \bar{\nu}^1_{kk'} \tilde{B}_{kk'} \frac{\rho_{k'}}{T_{k k'}^{3/2}}, +\end{aligned} +``` +with the so-called reduced temperature ``T_{k k'}`` and the ion-ion collision constants ``\tilde{B}_{kk'}`` provided +in `equations.ion_electron_collision_constants` (see [`CompressibleEulerMultiIonEquations1D`](@ref)). + +The additional coefficient ``\bar{\nu}^1_{kk'}`` is a non-dimensional drift correction factor proposed by Rambo and Denavit. + +References: +- P. Rambo, J. Denavit, Interpenetration and ion separation in colliding plasmas, Physics of Plasmas 1 (1994) 4050–4060. +- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. + Cambridge university press. +""" +@inline function source_terms_collision_ion_ion(u, x, t, + equations::CompressibleEulerMultiIonEquations1D) + s = zero(MVector{nvariables(equations), eltype(u)}) + @unpack gas_constants, molar_masses, ion_ion_collision_constants = equations + + prim = cons2prim(u, equations) + + for k in eachcomponent(equations) + rho_k, v1_k, p_k = get_component(k, prim, equations) + T_k = p_k / (rho_k * gas_constants[k]) + + S_q1 = zero(eltype(u)) + S_E = zero(eltype(u)) + for l in eachcomponent(equations) + # Do not compute collisions of an ion species with itself + k == l && continue + + rho_l, v1_l, p_l = get_component(l, prim, equations) + T_l = p_l / (rho_l * gas_constants[l]) + + # Reduced temperature + T_kl = (molar_masses[l] * T_k + molar_masses[k] * T_l) / + (molar_masses[k] + molar_masses[l]) + + delta_v = (v1_l - v1_k)^2 + + # Compute collision frequency without drifting correction + v_kl = ion_ion_collision_constants[k, l] * rho_l / T_kl^(3 / 2) + + # Correct the collision frequency with the drifting effect + z = delta_v / (p_l / rho_l + p_k / rho_k) + v_kl /= (1 + (2 / (9 * pi))^(1 / 3) * z)^(3 / 2) + + S_q1 += rho_k * v_kl * (v1_l - v1_k) + S_E += (3 * molar_masses[1] * gas_constants[1] * (T_l - T_k) + + + molar_masses[l] * delta_v) * v_kl * rho_k / + (molar_masses[k] + molar_masses[l]) + end + + S_E += v1_k * S_q1 + + set_component!(s, k, 0, S_q1, S_E, equations) + end + return SVector{nvariables(equations), real(equations)}(s) +end + +@doc raw""" + source_terms_collision_ion_electron(u, x, t, + equations::CompressibleEulerMultiIonEquations1D) + +Compute the ion-electron collision source terms for the momentum and energy equations of each ion species. We assume ``v_e = v^+`` +(no effect of currents on the electron velocity). + +The collision sources read as +```math +\begin{aligned} + \vec{s}_{\rho_k \vec{v}_k} =& \rho_k \bar{\nu}_{ke} (\vec{v}_{e} - \vec{v}_k), + \\ + s_{E_k} =& + 3 \left( + \bar{\nu}_{ke} \frac{\rho_k M_{1}}{M_k} R_1 (T_{e} - T_k) + \right) + + + \vec{v}_k \cdot \vec{s}_{\rho_k \vec{v}_k}, +\end{aligned} +``` +where ``T_e`` is the electron temperature computed with the function `equations.electron_temperature`, +``M_k`` is the molar mass of ion species `k` provided in `equations.molar_masses`, +``R_k`` is the specific gas constant of ion species `k` provided in `equations.gas_constants`, and +``\bar{\nu}_{kk'}`` is the collision frequency of species `k` with the electrons, which is computed as +```math +\begin{aligned} + \bar{\nu}_{ke} = \tilde{B}_{ke} \frac{e n_e}{T_e^{3/2}}, +\end{aligned} +``` +with the total electron charge ``e n_e`` (computed assuming quasi-neutrality), and the +ion-electron collision coefficient ``\tilde{B}_{ke}`` provided in `equations.ion_electron_collision_constants`, +which is scaled with the elementary charge (see [`CompressibleEulerMultiIonEquations1D`](@ref)). + +References: +- P. Rambo, J. Denavit, Interpenetration and ion separation in colliding plasmas, Physics of Plasmas 1 (1994) 4050–4060. +- Schunk, R. W., Nagy, A. F. (2000). Ionospheres: Physics, plasma physics, and chemistry. + Cambridge university press. +""" +function source_terms_collision_ion_electron(u, x, t, + equations::CompressibleEulerMultiIonEquations1D) + s = zero(MVector{nvariables(equations), eltype(u)}) + @unpack gas_constants, molar_masses, ion_electron_collision_constants, electron_temperature = equations + + prim = cons2prim(u, equations) + T_e = electron_temperature(u, equations) + T_e32 = T_e^(3 / 2) + v1_plus, vk1_plus = charge_averaged_velocities(u, equations) + + # Compute total electron charge + total_electron_charge = zero(real(equations)) + for k in eachcomponent(equations) + rho, _ = get_component(k, u, equations) + total_electron_charge += rho * equations.charge_to_mass[k] + end + + for k in eachcomponent(equations) + rho_k, v1_k, p_k = get_component(k, prim, equations) + T_k = p_k / (rho_k * gas_constants[k]) + + # Compute effective collision frequency + v_ke = ion_electron_collision_constants[k] * total_electron_charge / T_e32 + + S_q1 = rho_k * v_ke * (v1_plus - v1_k) + + S_E = 3 * molar_masses[1] * gas_constants[1] * (T_e - T_k) * v_ke * rho_k / + molar_masses[k] + + S_E += (v1_k * S_q1) + + set_component!(s, k, 0, S_q1, S_E, equations) + end + return SVector{nvariables(equations), real(equations)}(s) +end + # Calculate estimates for maximum wave speed for local Lax-Friedrichs-type dissipation as the # maximum velocity magnitude plus the maximum speed of sound @inline function max_abs_speed_naive(u_ll, u_rr, orientation::Integer, From 6b9a97ca36ad0512e052ac02f22cbbf008fe7750 Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 20 Jan 2025 15:20:40 +0100 Subject: [PATCH 5/6] Correct the example --- examples/tree_1d_dgsem/elixir_euler_multiion.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/tree_1d_dgsem/elixir_euler_multiion.jl b/examples/tree_1d_dgsem/elixir_euler_multiion.jl index 78086c69b3..585b5a1ed6 100644 --- a/examples/tree_1d_dgsem/elixir_euler_multiion.jl +++ b/examples/tree_1d_dgsem/elixir_euler_multiion.jl @@ -102,17 +102,17 @@ mesh = TreeMesh(coordinates_min, coordinates_max, # Ion-ion and ion-electron collision source terms # In this particular case, we can omit source_terms_lorentz because the magnetic field is zero! function source_terms(u, x, t, equations::Trixi.CompressibleEulerMultiIonEquations1D) - source_terms_collision_ion_ion(u, x, t, equations) + - source_terms_collision_ion_electron(u, x, t, equations) + Trixi.source_terms_collision_ion_ion(u, x, t, equations) + + Trixi.source_terms_collision_ion_electron(u, x, t, equations) end semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver, - source_terms = source_terms_convergence_test) + source_terms = source_terms) ############################################################################### # ODE solvers, callbacks etc. -tspan = (0.0, 2.0) +tspan = (0.0, 0.1) ode = semidiscretize(semi, tspan) summary_callback = SummaryCallback() @@ -129,7 +129,7 @@ save_solution = SaveSolutionCallback(interval = 100, save_final_solution = true, solution_variables = cons2prim) -stepsize_callback = StepsizeCallback(cfl = 0.8) +stepsize_callback = StepsizeCallback(cfl = 0.01) callbacks = CallbackSet(summary_callback, analysis_callback, alive_callback, From a2871f78acb56f470a52459341fde7169f9f1c75 Mon Sep 17 00:00:00 2001 From: Warisa Date: Mon, 20 Jan 2025 18:44:03 +0100 Subject: [PATCH 6/6] add a sanity check example --- .../elixir_euler_multiion_convergence_test.jl | 56 +++++++++++++++++++ 1 file changed, 56 insertions(+) create mode 100644 examples/tree_1d_dgsem/elixir_euler_multiion_convergence_test.jl diff --git a/examples/tree_1d_dgsem/elixir_euler_multiion_convergence_test.jl b/examples/tree_1d_dgsem/elixir_euler_multiion_convergence_test.jl new file mode 100644 index 0000000000..6f2b15556e --- /dev/null +++ b/examples/tree_1d_dgsem/elixir_euler_multiion_convergence_test.jl @@ -0,0 +1,56 @@ +using Trixi +using OrdinaryDiffEq + +############################################################################### +# semidiscretization of the compressible Euler equations +gammas = (1.4, 1.76, 1.4) # Heat capacity ratios for two species +equations = Trixi.CompressibleEulerMultiIonEquations1D(gammas=gammas, charge_to_mass = (0.0, 0.0, 0.0)) # Since I haven't exported this yet + +initial_condition = initial_condition_convergence_test + +solver = DGSEM(polydeg = 3, surface_flux = flux_hll, + volume_integral = VolumeIntegralPureLGLFiniteVolume(flux_hll)) + +coordinates_min = 0.0 +coordinates_max = 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_convergence_test) + +############################################################################### +# ODE solvers, callbacks etc. + +tspan = (0.0, 2.0) +ode = semidiscretize(semi, tspan) + +summary_callback = SummaryCallback() + +analysis_interval = 100 +analysis_callback = AnalysisCallback(semi, interval = analysis_interval, + extra_analysis_errors = (:l2_error_primitive, + :linf_error_primitive)) + +alive_callback = AliveCallback(analysis_interval = analysis_interval) + +save_solution = SaveSolutionCallback(interval = 100, + save_initial_solution = true, + save_final_solution = true, + solution_variables = cons2prim) + +stepsize_callback = StepsizeCallback(cfl = 0.8) + +callbacks = CallbackSet(summary_callback, + analysis_callback, alive_callback, + save_solution, + stepsize_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 \ No newline at end of file