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.

Sign up for GitHub

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 numerical support of other real types (multi) #1988

Merged
merged 11 commits into from
Jul 7, 2024
Merged

Add numerical support of other real types (multi) #1988

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
11 commits
Select commit Hold shift + click to select a range
fdbd012
start
huiyuxie Jun 21, 2024
6c003ad
Merge branch 'main' into main-multi
huiyuxie Jun 30, 2024
772c63f
continue
huiyuxie Jul 1, 2024
7b7b249
Merge branch 'main' into main-multi
huiyuxie Jul 2, 2024
7e148e9
complete equations
huiyuxie Jul 2, 2024
23c118a
Merge branch 'main' into main-multi
huiyuxie Jul 2, 2024
34def58
Merge branch 'main' into main-multi
huiyuxie Jul 3, 2024
24f4464
unit test 1D
huiyuxie Jul 4, 2024
3655e1d
unit test 2D
huiyuxie Jul 4, 2024
7160962
complete unit tests
huiyuxie Jul 5, 2024
4a406ac
Merge branch 'main' into main-multi
ranocha Jul 7, 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
2 changes: 1 addition & 1 deletion src/equations/compressible_euler_multicomponent_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -206,7 +206,7 @@ function initial_condition_weak_blast_wave(x, t,
2^(i - 1) * (1 - 2) /
(1 -
2^ncomponents(equations)) *
1 :
one(RealT) :
ranocha marked this conversation as resolved.
Show resolved Hide resolved
2^(i - 1) * (1 - 2) /
(1 -
2^ncomponents(equations)) *
Expand Down
2 changes: 1 addition & 1 deletion src/equations/compressible_euler_multicomponent_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -224,7 +224,7 @@ function initial_condition_weak_blast_wave(x, t,
2^(i - 1) * (1 - 2) /
(1 -
2^ncomponents(equations)) *
1 :
one(RealT) :
ranocha marked this conversation as resolved.
Show resolved Hide resolved
2^(i - 1) * (1 - 2) /
(1 -
2^ncomponents(equations)) *
Expand Down
173 changes: 92 additions & 81 deletions src/equations/ideal_glm_mhd_multicomponent_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -83,21 +83,24 @@ function initial_condition_convergence_test(x, t,
# smooth Alfvén wave test from Derigs et al. FLASH (2016)
# domain must be set to [0, 1], γ = 5/3

rho = 1.0
RealT = eltype(x)
rho = 1
prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) * (1 - 2) /
(1 -
2^ncomponents(equations)) *
rho
for i in eachcomponent(equations))
v1 = 0.0
si, co = sincos(2 * pi * x[1])
v2 = 0.1 * si
v3 = 0.1 * co
p = 0.1
B1 = 1.0
v1 = 0
# TODO: sincospi
si, co = sincos(2 * convert(RealT, pi) * x[1])
v2 = convert(RealT, 0.1) * si
v3 = convert(RealT, 0.1) * co
p = convert(RealT, 0.1)
B1 = 1
B2 = v2
B3 = v3
prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3)
prim_other = SVector(v1, v2, v3, p, B1, B2, B3)

return prim2cons(vcat(prim_other, prim_rho), equations)
end

Expand All @@ -114,31 +117,32 @@ function initial_condition_weak_blast_wave(x, t,
# Adapted MHD version of the weak blast wave from Hennemann & Gassner JCP paper 2020 (Sec. 6.3)
# Same discontinuity in the velocities but with magnetic fields
# Set up polar coordinates
RealT = eltype(x)
inicenter = (0)
x_norm = x[1] - inicenter[1]
r = sqrt(x_norm^2)
phi = atan(x_norm)

# Calculate primitive variables
if r > 0.5
rho = 1.0
if r > 0.5f0
rho = one(RealT)
prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) *
(1 - 2) / (1 -
2^ncomponents(equations)) *
rho
for i in eachcomponent(equations))
else
rho = 1.1691
rho = convert(RealT, 1.1691)
prim_rho = SVector{ncomponents(equations), real(equations)}(2^(i - 1) *
(1 - 2) / (1 -
2^ncomponents(equations)) *
rho
for i in eachcomponent(equations))
end
v1 = r > 0.5 ? 0.0 : 0.1882 * cos(phi)
p = r > 0.5 ? 1.0 : 1.245
v1 = r > 0.5f0 ? zero(RealT) : convert(RealT, 0.1882) * cos(phi)
p = r > 0.5f0 ? one(RealT) : convert(RealT, 1.245)

prim_other = SVector{7, real(equations)}(v1, 0.0, 0.0, p, 1.0, 1.0, 1.0)
prim_other = SVector(v1, 0, 0, p, 1, 1, 1)

return prim2cons(vcat(prim_other, prim_rho), equations)
end
Expand All @@ -153,8 +157,8 @@ end
v1 = rho_v1 / rho
v2 = rho_v2 / rho
v3 = rho_v3 / rho
kin_en = 0.5 * rho * (v1^2 + v2^2 + v3^2)
mag_en = 0.5 * (B1^2 + B2^2 + B3^2)
kin_en = 0.5f0 * rho * (v1^2 + v2^2 + v3^2)
mag_en = 0.5f0 * (B1^2 + B2^2 + B3^2)
gamma = totalgamma(u, equations)
p = (gamma - 1) * (rho_e - kin_en - mag_en)

Expand All @@ -164,11 +168,11 @@ end
f3 = rho_v1 * v3 - B1 * B3
f4 = (kin_en + gamma * p / (gamma - 1) + 2 * mag_en) * v1 -
B1 * (v1 * B1 + v2 * B2 + v3 * B3)
f5 = 0.0
f5 = 0
f6 = v1 * B2 - v2 * B1
f7 = v1 * B3 - v3 * B1

f_other = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7)
f_other = SVector(f1, f2, f3, f4, f5, f6, f7)

return vcat(f_other, f_rho)
end
Expand Down Expand Up @@ -198,7 +202,7 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer,
rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 7],
u_rr[i + 7])
for i in eachcomponent(equations))
rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 7] +
rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 7] +
u_rr[i + 7])
for i in eachcomponent(equations))

Expand All @@ -217,58 +221,60 @@ function flux_derigs_etal(u_ll, u_rr, orientation::Integer,
vel_dot_mag_rr = v1_rr * B1_rr + v2_rr * B2_rr + v3_rr * B3_rr

# Compute the necessary mean values needed for either direction
v1_avg = 0.5 * (v1_ll + v1_rr)
v2_avg = 0.5 * (v2_ll + v2_rr)
v3_avg = 0.5 * (v3_ll + v3_rr)
v1_avg = 0.5f0 * (v1_ll + v1_rr)
v2_avg = 0.5f0 * (v2_ll + v2_rr)
v3_avg = 0.5f0 * (v3_ll + v3_rr)
v_sum = v1_avg + v2_avg + v3_avg
B1_avg = 0.5 * (B1_ll + B1_rr)
B2_avg = 0.5 * (B2_ll + B2_rr)
B3_avg = 0.5 * (B3_ll + B3_rr)
vel_norm_avg = 0.5 * (vel_norm_ll + vel_norm_rr)
mag_norm_avg = 0.5 * (mag_norm_ll + mag_norm_rr)
vel_dot_mag_avg = 0.5 * (vel_dot_mag_ll + vel_dot_mag_rr)

enth = zero(v_sum)
help1_ll = zero(v1_ll)
help1_rr = zero(v1_rr)
B1_avg = 0.5f0 * (B1_ll + B1_rr)
B2_avg = 0.5f0 * (B2_ll + B2_rr)
B3_avg = 0.5f0 * (B3_ll + B3_rr)
vel_norm_avg = 0.5f0 * (vel_norm_ll + vel_norm_rr)
mag_norm_avg = 0.5f0 * (mag_norm_ll + mag_norm_rr)
vel_dot_mag_avg = 0.5f0 * (vel_dot_mag_ll + vel_dot_mag_rr)

RealT = eltype(u_ll)
enth = zero(RealT)
help1_ll = zero(RealT)
help1_rr = zero(RealT)

for i in eachcomponent(equations)
enth += rhok_avg[i] * gas_constants[i]
help1_ll += u_ll[i + 7] * cv[i]
help1_rr += u_rr[i + 7] * cv[i]
end

T_ll = (rho_e_ll - 0.5 * rho_ll * (vel_norm_ll) - 0.5 * mag_norm_ll) / help1_ll
T_rr = (rho_e_rr - 0.5 * rho_rr * (vel_norm_rr) - 0.5 * mag_norm_rr) / help1_rr
T = 0.5 * (1.0 / T_ll + 1.0 / T_rr)
T_log = ln_mean(1.0 / T_ll, 1.0 / T_rr)
T_ll = (rho_e_ll - 0.5f0 * rho_ll * (vel_norm_ll) - 0.5f0 * mag_norm_ll) / help1_ll
T_rr = (rho_e_rr - 0.5f0 * rho_rr * (vel_norm_rr) - 0.5f0 * mag_norm_rr) / help1_rr
T = 0.5f0 * (1 / T_ll + 1 / T_rr)
T_log = ln_mean(1 / T_ll, 1 / T_rr)

# Calculate fluxes depending on orientation with specific direction averages
help1 = zero(T_ll)
help2 = zero(T_rr)
help1 = zero(RealT)
help2 = zero(RealT)

f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg
for i in eachcomponent(equations))
for i in eachcomponent(equations)
help1 += f_rho[i] * cv[i]
help2 += f_rho[i]
end
f1 = help2 * v1_avg + enth / T + 0.5 * mag_norm_avg - B1_avg * B1_avg
f1 = help2 * v1_avg + enth / T + 0.5f0 * mag_norm_avg - B1_avg * B1_avg
f2 = help2 * v2_avg - B1_avg * B2_avg
f3 = help2 * v3_avg - B1_avg * B3_avg
f5 = 0.0
f5 = 0
f6 = v1_avg * B2_avg - v2_avg * B1_avg
f7 = v1_avg * B3_avg - v3_avg * B1_avg

# total energy flux is complicated and involves the previous eight components
v1_mag_avg = 0.5 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr)
v1_mag_avg = 0.5f0 * (v1_ll * mag_norm_ll + v1_rr * mag_norm_rr)

f4 = (help1 / T_log) - 0.5 * (vel_norm_avg) * (help2) + f1 * v1_avg + f2 * v2_avg +
f4 = (help1 / T_log) - 0.5f0 * (vel_norm_avg) * (help2) + f1 * v1_avg +
f2 * v2_avg +
f3 * v3_avg +
f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5 * v1_mag_avg +
f5 * B1_avg + f6 * B2_avg + f7 * B3_avg - 0.5f0 * v1_mag_avg +
B1_avg * vel_dot_mag_avg

f_other = SVector{7, real(equations)}(f1, f2, f3, f4, f5, f6, f7)
f_other = SVector(f1, f2, f3, f4, f5, f6, f7)

return vcat(f_other, f_rho)
end
Expand Down Expand Up @@ -309,23 +315,24 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations.
# log((ϱₗ/pₗ) / (ϱᵣ/pᵣ)) / (ϱₗ/pₗ - ϱᵣ/pᵣ)
# = pₗ pᵣ log((ϱₗ pᵣ) / (ϱᵣ pₗ)) / (ϱₗ pᵣ - ϱᵣ pₗ)
inv_rho_p_mean = p_ll * p_rr * inv_ln_mean(rho_ll * p_rr, rho_rr * p_ll)
v1_avg = 0.5 * (v1_ll + v1_rr)
v2_avg = 0.5 * (v2_ll + v2_rr)
v3_avg = 0.5 * (v3_ll + v3_rr)
p_avg = 0.5 * (p_ll + p_rr)
velocity_square_avg = 0.5 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)
magnetic_square_avg = 0.5 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr)
v1_avg = 0.5f0 * (v1_ll + v1_rr)
v2_avg = 0.5f0 * (v2_ll + v2_rr)
v3_avg = 0.5f0 * (v3_ll + v3_rr)
p_avg = 0.5f0 * (p_ll + p_rr)
velocity_square_avg = 0.5f0 * (v1_ll * v1_rr + v2_ll * v2_rr + v3_ll * v3_rr)
magnetic_square_avg = 0.5f0 * (B1_ll * B1_rr + B2_ll * B2_rr + B3_ll * B3_rr)

inv_gamma_minus_one = 1 / (totalgamma(0.5 * (u_ll + u_rr), equations) - 1)
inv_gamma_minus_one = 1 / (totalgamma(0.5f0 * (u_ll + u_rr), equations) - 1)

rhok_mean = SVector{ncomponents(equations), real(equations)}(ln_mean(u_ll[i + 7],
u_rr[i + 7])
for i in eachcomponent(equations))
rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5 * (u_ll[i + 7] +
rhok_avg = SVector{ncomponents(equations), real(equations)}(0.5f0 * (u_ll[i + 7] +
u_rr[i + 7])
for i in eachcomponent(equations))

f1 = zero(rho_ll)
RealT = eltype(u_ll)
f1 = zero(RealT)
f_rho = SVector{ncomponents(equations), real(equations)}(rhok_mean[i] * v1_avg
for i in eachcomponent(equations))
for i in eachcomponent(equations)
Expand All @@ -334,25 +341,25 @@ Hindenlang (2019), extending [`flux_ranocha`](@ref) to the MHD equations.

# Calculate fluxes depending on orientation with specific direction averages
f2 = f1 * v1_avg + p_avg + magnetic_square_avg -
0.5 * (B1_ll * B1_rr + B1_rr * B1_ll)
f3 = f1 * v2_avg - 0.5 * (B1_ll * B2_rr + B1_rr * B2_ll)
f4 = f1 * v3_avg - 0.5 * (B1_ll * B3_rr + B1_rr * B3_ll)
#f5 below
f6 = 0.0
f7 = 0.5 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr)
f8 = 0.5 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr)
0.5f0 * (B1_ll * B1_rr + B1_rr * B1_ll)
f3 = f1 * v2_avg - 0.5f0 * (B1_ll * B2_rr + B1_rr * B2_ll)
f4 = f1 * v3_avg - 0.5f0 * (B1_ll * B3_rr + B1_rr * B3_ll)
# f5 below
f6 = 0
f7 = 0.5f0 * (v1_ll * B2_ll - v2_ll * B1_ll + v1_rr * B2_rr - v2_rr * B1_rr)
f8 = 0.5f0 * (v1_ll * B3_ll - v3_ll * B1_ll + v1_rr * B3_rr - v3_rr * B1_rr)
# total energy flux is complicated and involves the previous components
f5 = (f1 * (velocity_square_avg + inv_rho_p_mean * inv_gamma_minus_one)
+
0.5 * (+p_ll * v1_rr + p_rr * v1_ll
0.5f0 * (+p_ll * v1_rr + p_rr * v1_ll
+ (v1_ll * B2_ll * B2_rr + v1_rr * B2_rr * B2_ll)
+ (v1_ll * B3_ll * B3_rr + v1_rr * B3_rr * B3_ll)
-
(v2_ll * B1_ll * B2_rr + v2_rr * B1_rr * B2_ll)
-
(v3_ll * B1_ll * B3_rr + v3_rr * B1_rr * B3_ll)))

f_other = SVector{7, real(equations)}(f2, f3, f4, f5, f6, f7, f8)
f_other = SVector(f2, f3, f4, f5, f6, f7, f8)

return vcat(f_other, f_rho)
end
Expand Down Expand Up @@ -405,8 +412,9 @@ function cons2prim(u, equations::IdealGlmMhdMulticomponentEquations1D)
gamma = totalgamma(u, equations)

p = (gamma - 1) *
(rho_e - 0.5 * rho * (v1^2 + v2^2 + v3^2) - 0.5 * (B1^2 + B2^2 + B3^2))
prim_other = SVector{7, real(equations)}(v1, v2, v3, p, B1, B2, B3)
(rho_e - 0.5f0 * rho * (v1^2 + v2^2 + v3^2) - 0.5f0 * (B1^2 + B2^2 + B3^2))
prim_other = SVector(v1, v2, v3, p, B1, B2, B3)

return vcat(prim_other, prim_rho)
end

Expand All @@ -422,20 +430,21 @@ end
v3 = rho_v3 / rho
v_square = v1^2 + v2^2 + v3^2
gamma = totalgamma(u, equations)
p = (gamma - 1) * (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2))
p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2))
s = log(p) - gamma * log(rho)
rho_p = rho / p

# Multicomponent stuff
help1 = zero(v1)
RealT = eltype(u)
help1 = zero(RealT)

for i in eachcomponent(equations)
help1 += u[i + 7] * cv[i]
end

T = (rho_e - 0.5 * rho * v_square - 0.5 * (B1^2 + B2^2 + B3^2)) / (help1)
T = (rho_e - 0.5f0 * rho * v_square - 0.5f0 * (B1^2 + B2^2 + B3^2)) / (help1)

entrop_rho = SVector{ncomponents(equations), real(equations)}(-1.0 *
entrop_rho = SVector{ncomponents(equations), real(equations)}(-1 *
(cv[i] * log(T) -
gas_constants[i] *
log(u[i + 7])) +
Expand All @@ -447,12 +456,12 @@ end
w1 = v1 / T
w2 = v2 / T
w3 = v3 / T
w4 = -1.0 / T
w4 = -1 / T
w5 = B1 / T
w6 = B2 / T
w7 = B3 / T

entrop_other = SVector{7, real(equations)}(w1, w2, w3, w4, w5, w6, w7)
entrop_other = SVector(w1, w2, w3, w4, w5, w6, w7)

return vcat(entrop_other, entrop_rho)
end
Expand All @@ -470,10 +479,10 @@ end
rho_v3 = rho * v3

gamma = totalgamma(prim, equations)
rho_e = p / (gamma - 1) + 0.5 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) +
0.5 * (B1^2 + B2^2 + B3^2)
rho_e = p / (gamma - 1) + 0.5f0 * (rho_v1 * v1 + rho_v2 * v2 + rho_v3 * v3) +
0.5f0 * (B1^2 + B2^2 + B3^2)

cons_other = SVector{7, real(equations)}(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3)
cons_other = SVector(rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3)

return vcat(cons_other, cons_rho)
end
Expand All @@ -482,9 +491,9 @@ end
rho_v1, rho_v2, rho_v3, rho_e, B1, B2, B3 = u
rho = density(u, equations)
gamma = totalgamma(u, equations)
p = (gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
p = (gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho
-
0.5 * (B1^2 + B2^2 + B3^2))
0.5f0 * (B1^2 + B2^2 + B3^2))
return rho * p
end

Expand All @@ -498,22 +507,23 @@ end
v3 = rho_v3 / rho
v_mag = sqrt(v1^2 + v2^2 + v3^2)
gamma = totalgamma(cons, equations)
p = (gamma - 1) * (rho_e - 0.5 * rho * v_mag^2 - 0.5 * (B1^2 + B2^2 + B3^2))
p = (gamma - 1) * (rho_e - 0.5f0 * rho * v_mag^2 - 0.5f0 * (B1^2 + B2^2 + B3^2))
a_square = gamma * p / rho
sqrt_rho = sqrt(rho)
b1 = B1 / sqrt_rho
b2 = B2 / sqrt_rho
b3 = B3 / sqrt_rho
b_square = b1^2 + b2^2 + b3^2

c_f = sqrt(0.5 * (a_square + b_square) +
0.5 * sqrt((a_square + b_square)^2 - 4.0 * a_square * b1^2))
c_f = sqrt(0.5f0 * (a_square + b_square) +
0.5f0 * sqrt((a_square + b_square)^2 - 4 * a_square * b1^2))

return c_f
end

@inline function density(u, equations::IdealGlmMhdMulticomponentEquations1D)
rho = zero(u[1])
RealT = eltype(u)
rho = zero(RealT)

for i in eachcomponent(equations)
rho += u[i + 7]
Expand All @@ -525,8 +535,9 @@ end
@inline function totalgamma(u, equations::IdealGlmMhdMulticomponentEquations1D)
@unpack cv, gammas = equations

help1 = zero(u[1])
help2 = zero(u[1])
RealT = eltype(u)
help1 = zero(RealT)
help2 = zero(RealT)

for i in eachcomponent(equations)
help1 += u[i + 7] * cv[i] * gammas[i]
Expand Down
Loading
Loading