Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add GLM divergence cleaning to the multi-ion MHD equation #14

Merged
merged 23 commits into from
Nov 21, 2024
Merged
Changes from 1 commit
Commits
Show all changes
23 commits
Select commit Hold shift + click to select a range
1ba2092
First attempt to enable save solution with time intervals for SimpleI…
amrueda Oct 18, 2023
e05363d
Importing `add_tstop!`
amrueda Oct 19, 2023
64e9784
First working version of SimpleIntegratorSSP with SaveSolutionCallbac…
amrueda Oct 19, 2023
05fcaaa
Improved formatting and updated reference solution for test
amrueda Oct 19, 2023
71f5c7a
Modified initialization of tstops to ensure a stop at the end of the …
amrueda Oct 19, 2023
df34188
Added missing docstring
amrueda Oct 23, 2023
d04fd98
Removed OrdinaryDiffEq from Trixi's dependencies
amrueda Oct 23, 2023
c3b4d5d
Merge branch 'main' into simplessprk_dt
amrueda Oct 23, 2023
7aa515b
Empty tstops BinaryHeap during the call to `terminate!(integrator::Si…
amrueda Oct 23, 2023
4322f5b
Fixed bug and added explanatory comments
amrueda Oct 23, 2023
772a388
Updated Project.toml
amrueda Oct 23, 2023
7748aba
Merge branch 'main' into simplessprk_dt
amrueda Oct 25, 2023
683e82e
Merge branch 'main' into simplessprk_dt
amrueda Dec 12, 2023
a3fb5e5
format
amrueda Dec 12, 2023
b8d1a2a
Merge branch 'main' into simplessprk_dt
amrueda Dec 13, 2023
a908fd7
Merge remote-tracking branch 'amrueda/simplessprk_dt' into multi_ion_mhd
amrueda Dec 13, 2023
31d27ad
Implemented GLM divergence cleaning for multi-ion MHD
amrueda Dec 19, 2023
a33ad14
Added new ES numerical fluxes
amrueda Jan 24, 2024
a6b831c
Added missing terms to multi-ion ES flux
amrueda Feb 15, 2024
9fbd286
Export magnetic_field and divergence_clenaing_field
amrueda Feb 15, 2024
f0f55a9
Removed 'lumped' ES flux as it is not really useful
amrueda Feb 21, 2024
79ca7c0
Fixed bug in the discretization of the electron pressure!
amrueda Feb 21, 2024
9e1dbff
format
amrueda Feb 21, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Removed 'lumped' ES flux as it is not really useful
  • Loading branch information
amrueda committed Feb 21, 2024
commit f0f55a99549e43a28811a3552e120f20b05cac7e
103 changes: 0 additions & 103 deletions src/equations/ideal_mhd_multiion_2d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1044,109 +1044,6 @@ Computes the sum of the densities times the sum of the pressures
return rho_total * p_total
end

"""
DissipationEntropyStableLump(max_abs_speed=max_abs_speed_naive)

Create a local Lax-Friedrichs dissipation operator where the maximum absolute wave speed
is estimated as
`max_abs_speed(u_ll, u_rr, orientation_or_normal_direction, equations)`,
defaulting to [`max_abs_speed_naive`](@ref).
"""
struct DissipationEntropyStableLump{MaxAbsSpeed}
max_abs_speed::MaxAbsSpeed
end

DissipationEntropyStableLump() = DissipationEntropyStableLump(max_abs_speed_naive)

@inline function (dissipation::DissipationEntropyStableLump)(u_ll, u_rr,
orientation_or_normal_direction,
equations::IdealMhdMultiIonEquations2D)
@unpack gammas = equations
λ = dissipation.max_abs_speed(u_ll, u_rr, orientation_or_normal_direction,
equations)

w_ll = cons2entropy(u_ll, equations)
w_rr = cons2entropy(u_rr, equations)
prim_ll = cons2prim(u_ll, equations)
prim_rr = cons2prim(u_rr, equations)
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)

# Some global averages
B1_avg = 0.5 * (B1_ll + B1_rr)
B2_avg = 0.5 * (B2_ll + B2_rr)
B3_avg = 0.5 * (B3_ll + B3_rr)
psi_avg = 0.5 * (psi_ll + psi_rr)

dissipation = zero(MVector{nvariables(equations), eltype(u_ll)})

beta_plus_ll = 0
beta_plus_rr = 0
# Get the lumped dissipation for all components
for k in eachcomponent(equations)
rho_ll, v1_ll, v2_ll, v3_ll, p_ll = get_component(k, prim_ll, equations)
rho_rr, v1_rr, v2_rr, v3_rr, p_rr = get_component(k, prim_rr, equations)

w1_ll, w2_ll, w3_ll, w4_ll, w5_ll = get_component(k, w_ll, equations)
w1_rr, w2_rr, w3_rr, w4_rr, w5_rr = get_component(k, w_rr, equations)

# Auxiliary variables
beta_ll = 0.5 * rho_ll / p_ll
beta_rr = 0.5 * rho_rr / p_rr
vel_norm_ll = v1_ll^2 + v2_ll^2 + v3_ll^2
vel_norm_rr = v1_rr^2 + v2_rr^2 + v3_rr^2

# Mean variables
rho_ln = ln_mean(rho_ll, rho_rr)
beta_ln = ln_mean(beta_ll, beta_rr)
rho_avg = 0.5 * (rho_ll + rho_rr)
v1_avg = 0.5 * (v1_ll + v1_rr)
v2_avg = 0.5 * (v2_ll + v2_rr)
v3_avg = 0.5 * (v3_ll + v3_rr)
beta_avg = 0.5 * (beta_ll + beta_rr)
tau = 1 / (beta_ll + beta_rr)
p_mean = 0.5 * rho_avg / beta_avg
p_star = 0.5 * rho_ln / beta_ln
vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr)
vel_avg_norm = v1_avg^2 + v2_avg^2 + v3_avg^2
E_bar = p_star / (gammas[k] - 1) + 0.5 * rho_ln * (2 * vel_avg_norm - vel_norm_avg)

h11 = rho_ln
h22 = rho_ln * v1_avg^2 + p_mean
h33 = rho_ln * v2_avg^2 + p_mean
h44 = rho_ln * v3_avg^2 + p_mean
h55 = ((p_star^2 / (gammas[k] - 1) + E_bar * E_bar) / rho_ln
+ vel_avg_norm * p_mean
+ tau * ( B1_avg^2 + B2_avg^2 + B3_avg^2 + psi_avg^2 ))

d1 = -0.5 * λ * h11 * (w1_rr - w1_ll)
d2 = -0.5 * λ * h22 * (w2_rr - w2_ll)
d3 = -0.5 * λ * h33 * (w3_rr - w3_ll)
d4 = -0.5 * λ * h44 * (w4_rr - w4_ll)
d5 = -0.5 * λ * h55 * (w5_rr - w5_ll)

beta_plus_ll += beta_ll
beta_plus_rr += beta_rr

set_component!(dissipation, k, d1, d2, d3, d4, d5, equations)
end

# Set the magnetic field and psi terms
h_B_psi = 1 / (beta_plus_ll + beta_plus_rr)
dissipation[1] = -0.5 * λ * h_B_psi * (w_rr[1] - w_ll[1])
dissipation[2] = -0.5 * λ * h_B_psi * (w_rr[2] - w_ll[2])
dissipation[3] = -0.5 * λ * h_B_psi * (w_rr[3] - w_ll[3])
dissipation[end] = -0.5 * λ * h_B_psi * (w_rr[end] - w_ll[end])

return dissipation
end

function Base.show(io::IO, d::DissipationEntropyStableLump)
print(io, "DissipationEntropyStableLump(", d.max_abs_speed, ")")
end

"""
DissipationEntropyStable(max_abs_speed=max_abs_speed_naive)

Expand Down