Skip to content

Commit

Permalink
complete the first
Browse files Browse the repository at this point in the history
  • Loading branch information
huiyuxie committed May 28, 2024
1 parent 7e068bc commit 0b9589e
Showing 1 changed file with 17 additions and 17 deletions.
34 changes: 17 additions & 17 deletions src/equations/compressible_euler_quasi_1d.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,8 +81,8 @@ function initial_condition_convergence_test(x, t,
RealT = eltype(x)
c = 2
A = convert(RealT, 0.1)
L = convert(RealT, 2)
f = 1 / L
L = 2
f = 1.0f0 / L
ω = 2 * convert(RealT, pi) * f
ini = c + A * sin* (x[1] - t))

Expand Down Expand Up @@ -110,18 +110,18 @@ as defined in [`initial_condition_convergence_test`](@ref).
# Same settings as in `initial_condition_convergence_test`.
# Derivatives calculated with ForwardDiff.jl
c = 2
A = 0.1
A = convert(RealT, 0.1)
L = 2
f = 1 / L
ω = 2 * pi * f
f = 1.0f0 / L
ω = 2 * convert(RealT, pi) * f
x1, = x
ini(x1, t) = c + A * sin* (x1 - t))

rho(x1, t) = ini(x1, t)
v1(x1, t) = 1.0
v1(x1, t) = 1
e(x1, t) = ini(x1, t)^2 / rho(x1, t)
p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5 * rho(x1, t) * v1(x1, t)^2)
a(x1, t) = 1.5 - 0.5 * cos(x1 * pi)
p1(x1, t) = (equations.gamma - 1) * (e(x1, t) - 0.5f0 * rho(x1, t) * v1(x1, t)^2)
a(x1, t) = 1.5f0 - 0.5f0 * cos(x1 * pi)

arho(x1, t) = a(x1, t) * rho(x1, t)
arhou(x1, t) = arho(x1, t) * v1(x1, t)
Expand All @@ -143,7 +143,7 @@ as defined in [`initial_condition_convergence_test`](@ref).
du2 = darhou_dt(x1, t) + darhouu_dx(x1, t) + a(x1, t) * dp1_dx(x1, t)
du3 = daE_dt(x1, t) + dauEp_dx(x1, t)

return SVector(du1, du2, du3, 0.0)
return SVector(du1, du2, du3, 0)
end

# Calculate 1D flux for a single point
Expand Down Expand Up @@ -240,17 +240,17 @@ Further details are available in the paper:
# 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)
a_v1_avg = 0.5 * (a_ll * v1_ll + a_rr * v1_rr)
p_avg = 0.5 * (p_ll + p_rr)
velocity_square_avg = 0.5 * (v1_ll * v1_rr)
v1_avg = 0.5f0 * (v1_ll + v1_rr)
a_v1_avg = 0.5f0 * (a_ll * v1_ll + a_rr * v1_rr)
p_avg = 0.5f0 * (p_ll + p_rr)
velocity_square_avg = 0.5f0 * (v1_ll * v1_rr)

# Calculate fluxes
# Ignore orientation since it is always "1" in 1D
f1 = rho_mean * a_v1_avg
f2 = rho_mean * a_v1_avg * v1_avg
f3 = f1 * (velocity_square_avg + inv_rho_p_mean * equations.inv_gamma_minus_one) +
0.5 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll)
0.5f0 * (p_ll * a_rr * v1_rr + p_rr * a_ll * v1_ll)

return SVector(f1, f2, f3, zero(eltype(u_ll)))
end
Expand All @@ -277,13 +277,13 @@ end
e_ll = a_e_ll / a_ll
v1_ll = a_rho_v1_ll / a_rho_ll
v_mag_ll = abs(v1_ll)
p_ll = (equations.gamma - 1) * (e_ll - 0.5 * rho_ll * v_mag_ll^2)
p_ll = (equations.gamma - 1) * (e_ll - 0.5f0 * rho_ll * v_mag_ll^2)
c_ll = sqrt(equations.gamma * p_ll / rho_ll)
rho_rr = a_rho_rr / a_rr
e_rr = a_e_rr / a_rr
v1_rr = a_rho_v1_rr / a_rho_rr
v_mag_rr = abs(v1_rr)
p_rr = (equations.gamma - 1) * (e_rr - 0.5 * rho_rr * v_mag_rr^2)
p_rr = (equations.gamma - 1) * (e_rr - 0.5f0 * rho_rr * v_mag_rr^2)
c_rr = sqrt(equations.gamma * p_rr / rho_rr)

λ_max = max(v_mag_ll, v_mag_rr) + max(c_ll, c_rr)
Expand All @@ -294,7 +294,7 @@ end
rho = a_rho / a
v1 = a_rho_v1 / a_rho
e = a_e / a
p = (equations.gamma - 1) * (e - 0.5 * rho * v1^2)
p = (equations.gamma - 1) * (e - 0.5f0 * rho * v1^2)
c = sqrt(equations.gamma * p / rho)

return (abs(v1) + c,)
Expand Down

0 comments on commit 0b9589e

Please sign in to comment.