diff --git a/src/equations/ideal_glm_mhd_multiion.jl b/src/equations/ideal_glm_mhd_multiion.jl index 0d0229e8750..37ba7c52a48 100644 --- a/src/equations/ideal_glm_mhd_multiion.jl +++ b/src/equations/ideal_glm_mhd_multiion.jl @@ -45,6 +45,31 @@ function default_analysis_integrals(::AbstractIdealGlmMhdMultiIonEquations) (entropy_timederivative, Val(:l2_divb), Val(:linf_divb)) end +# For `VolumeIntegralSubcellLimiting` the nonconservative flux is created as a callable struct to +# enable dispatch on the type of the nonconservative term (symmetric / jump). +struct FluxNonConservativeRuedaRamirezEtAl <: + FluxNonConservative{NonConservativeSymmetric()} +end + +# We specify 6 non-conservative terms for FluxNonConservativeRuedaRamirezEtAl. This is the number of +# non-conservative terms in 3D. In 2D, only 5 terms are needed. TODO: Should we create a different struct +# for the 2D non-conservative term? +n_nonconservative_terms(::FluxNonConservativeRuedaRamirezEtAl) = 6 + +const flux_nonconservative_ruedaramirez_etal = FluxNonConservativeRuedaRamirezEtAl() + +# State validation for Newton-bisection method of subcell IDP limiting +@inline function Base.isvalid(u, equations::AbstractIdealGlmMhdMultiIonEquations) + p = pressure(u, equations) + for k in eachcomponent(equations) + u_k = get_component(k, u, equations) + if u_k[1] <= 0 || p[k] <= 0 + return false + end + end + return true +end + """ source_terms_lorentz(u, x, t, equations::AbstractIdealGlmMhdMultiIonEquations) diff --git a/src/equations/ideal_glm_mhd_multiion_2d.jl b/src/equations/ideal_glm_mhd_multiion_2d.jl index 38b3e52fa49..870ceb6bae6 100644 --- a/src/equations/ideal_glm_mhd_multiion_2d.jl +++ b/src/equations/ideal_glm_mhd_multiion_2d.jl @@ -343,9 +343,9 @@ The term is composed of four individual non-conservative terms: 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) +@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(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) diff --git a/src/equations/ideal_glm_mhd_multiion_3d.jl b/src/equations/ideal_glm_mhd_multiion_3d.jl index f5a88ca4192..3c2988c89fc 100644 --- a/src/equations/ideal_glm_mhd_multiion_3d.jl +++ b/src/equations/ideal_glm_mhd_multiion_3d.jl @@ -329,9 +329,9 @@ The term is composed of four individual non-conservative terms: 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::IdealGlmMhdMultiIonEquations3D) +@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr, + orientation::Integer, + equations::IdealGlmMhdMultiIonEquations3D) @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) @@ -514,9 +514,9 @@ The term is composed of four individual non-conservative terms: return SVector(f) end -@inline function flux_nonconservative_ruedaramirez_etal(u_ll, u_rr, - normal_direction::AbstractVector, - equations::IdealGlmMhdMultiIonEquations3D) +@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr, + normal_direction::AbstractVector, + equations::IdealGlmMhdMultiIonEquations3D) @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) @@ -628,6 +628,291 @@ end return SVector(f) end +""" + flux_nonconservative_ruedaramirez_etal(u_ll, normal_direction_ll::AbstractVector, + equations::IdealGlmMhdMultiIonEquations3D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + +Non-symmetric local part of the non-conservative terms for multi-ion MHD. Needed for the calculation of +the non-conservative staggered "fluxes" for subcell limiting. See, e.g., +- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts + Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl. + +On curvilinear meshes this formulation applies the local normal direction compared to the averaged one used +in [`flux_nonconservative_ruedaramirez_etal`](@ref) when used for volume fluxes. This is done to reduce the +number of operations to obtain the subcell limiting fluxes. However, this decision causes this flux to be +"slightly" different when used for subcell limiting or for a flux-differencing DG method. +""" +@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, + normal_direction::AbstractVector, + equations::IdealGlmMhdMultiIonEquations3D, + nonconservative_type::NonConservativeLocal, + nonconservative_term::Integer) + @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) + psi_ll = divergence_cleaning_field(u_ll, equations) + + # 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] # Extract densities from conserved variable vector + 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) + + f = zero(MVector{nvariables(equations), eltype(u_ll)}) + + if nonconservative_term == 1 + f[1] = v1_plus_ll + f[2] = v2_plus_ll + f[3] = v3_plus_ll + + for k in eachcomponent(equations) + # Compute Godunov-Powell term + f2 = charge_ratio_ll[k] * B1_ll + f3 = charge_ratio_ll[k] * B2_ll + f4 = charge_ratio_ll[k] * B3_ll + f5 = (v1_plus_ll * B1_ll + v2_plus_ll * B2_ll + v3_plus_ll * B3_ll) + + set_component!(f, k, 0, f2, f3, f4, f5, equations) + end + elseif nonconservative_term == 2 + # Compute Lorentz term + + for k in eachcomponent(equations) + f2 = charge_ratio_ll[k] + f3 = charge_ratio_ll[k] + f4 = charge_ratio_ll[k] + f5 = (vk1_plus_ll[k] * normal_direction[1] + + vk2_plus_ll[k] * normal_direction[2] + + vk3_plus_ll[k] * normal_direction[3]) + + set_component!(f, k, 0, f2, f3, f4, f5, equations) + end + elseif nonconservative_term == 3 + # Compute GLM term + v_plus_dot_n_ll = (v1_plus_ll * normal_direction[1] + + v2_plus_ll * normal_direction[2] + + v3_plus_ll * normal_direction[3]) + for k in eachcomponent(equations) + f5 = v_plus_dot_n_ll * psi_ll + set_component!(f, k, 0, 0, 0, 0, f5, equations) + end + # Compute GLM term for psi + f[end] = v_plus_dot_n_ll + elseif nonconservative_term == 4 + # Multi-ion term (vanishes for NCOMP==1) for B1 + for k in eachcomponent(equations) + f5 = B1_ll + set_component!(f, k, 0, 0, 0, 0, f5, equations) + end + elseif nonconservative_term == 5 + # Multi-ion term (vanishes for NCOMP==1) for B2 + for k in eachcomponent(equations) + f5 = B2_ll + set_component!(f, k, 0, 0, 0, 0, f5, equations) + end + elseif nonconservative_term == 6 + # Multi-ion term (vanishes for NCOMP==1) for B3 + for k in eachcomponent(equations) + f5 = B3_ll + set_component!(f, k, 0, 0, 0, 0, f5, equations) + end + end + + return SVector(f) +end + +""" + flux_nonconservative_ruedaramirez_etal(u_ll, normal_direction_avg::AbstractVector, + equations::IdealGlmMhdMultiIonEquations3D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + +Symmetric part of the multi-ion non-conservative terms. Needed for the calculation of +the non-conservative staggered "fluxes" for subcell limiting. See, e.g., +- Rueda-Ramírez, Gassner (2023). A Flux-Differencing Formula for Split-Form Summation By Parts + Discretizations of Non-Conservative Systems. https://arxiv.org/pdf/2211.14009.pdf. +This function is used to compute the subcell fluxes in dg_3d_subcell_limiters.jl. +""" +@inline function (noncons_flux::FluxNonConservativeRuedaRamirezEtAl)(u_ll, u_rr, + normal_direction::AbstractVector, + equations::IdealGlmMhdMultiIonEquations3D, + nonconservative_type::NonConservativeSymmetric, + nonconservative_term::Integer) + @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) + B_dot_n_ll = B1_ll * normal_direction[1] + + B2_ll * normal_direction[2] + + B3_ll * normal_direction[3] + B_dot_n_rr = B1_rr * normal_direction[1] + + B2_rr * normal_direction[2] + + B3_rr * normal_direction[3] + B_dot_n_avg = 0.5f0 * (B_dot_n_ll + B_dot_n_rr) + + # 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) + + # 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] # Extract densities from conserved variable vector + 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 nonconservative_term == 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 * B_dot_n_avg + f[2] = 2 * B_dot_n_avg + f[3] = 2 * B_dot_n_avg + + for k in eachcomponent(equations) + # Compute Godunov-Powell term + f2 = B_dot_n_avg + f3 = B_dot_n_avg + f4 = B_dot_n_avg + f5 = B_dot_n_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 + elseif nonconservative_term == 2 + # Compute Lorentz term + + # 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 + mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_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) + + for k in eachcomponent(equations) + f2 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[1] - + B_dot_n_avg * B1_avg) + f3 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[2] - + B_dot_n_avg * B2_avg) + f4 = ((0.5f0 * mag_norm_avg + pe_mean) * normal_direction[3] - + B_dot_n_avg * B3_avg) + f5 = pe_mean + + set_component!(f, k, 0, 2 * f2, 2 * f3, 2 * f4, 2 * f5, + equations) + end + elseif nonconservative_term == 3 + # Compute GLM term + psi_avg = 0.5f0 * (psi_ll + psi_rr) + for k in eachcomponent(equations) + f5 = psi_avg + # (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, 0, 0, 0, 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 * psi_avg + elseif nonconservative_term == 4 + # Multi-ion term (vanishes for NCOMP==1) for B1 + for k in eachcomponent(equations) + 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 = ((vk2_minus_avg * B1_avg - vk1_minus_avg * B2_avg) * + normal_direction[2] + + (vk3_minus_avg * B1_avg - vk1_minus_avg * B3_avg) * + normal_direction[3]) + + # 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, 0, 0, 0, 2 * f5, + equations) + end + elseif nonconservative_term == 5 + # Multi-ion term (vanishes for NCOMP==1) for B2 + for k in eachcomponent(equations) + 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 = ((vk1_minus_avg * B2_avg - vk2_minus_avg * B1_avg) * + normal_direction[1] + + (vk3_minus_avg * B2_avg - vk2_minus_avg * B3_avg) * + normal_direction[3]) + # 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, 0, 0, 0, 2 * f5, + equations) + end + elseif nonconservative_term == 6 + # Multi-ion term (vanishes for NCOMP==1) for B3 + for k in eachcomponent(equations) + 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 = ((vk1_minus_avg * B3_avg - vk3_minus_avg * B1_avg) * + normal_direction[1] + + (vk2_minus_avg * B3_avg - vk3_minus_avg * B2_avg) * + normal_direction[2]) + # 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, 0, 0, 0, 2 * f5, + equations) + end + end + + return SVector(f) +end + """ flux_nonconservative_central(u_ll, u_rr, orientation::Integer, equations::IdealGlmMhdMultiIonEquations3D)