Skip to content

Commit

Permalink
Merge branch 'main' into bug_fix_PERK2
Browse files Browse the repository at this point in the history
  • Loading branch information
@warisa-r
warisa-r authored Jul 7, 2024
2 parents 372590f + 2d0aa80 commit 78e2585
Show file tree
Hide file tree
Showing 6 changed files with 400 additions and 209 deletions.
2 changes: 1 addition & 1 deletion src/equations/compressible_euler_multicomponent_1d.jl
View file
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) :
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
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) :
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
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

0 comments on commit 78e2585

Please sign in to comment.