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 (navier) #1968

Merged
merged 19 commits into from
Jun 25, 2024
Merged

Add numerical support of other real types (navier) #1968

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
19 commits
Select commit Hold shift + click to select a range
24d9fe8
start
huiyuxie Jun 6, 2024
7b98442
Merge branch 'main' into main-navier
huiyuxie Jun 7, 2024
663de65
complete equations
huiyuxie Jun 7, 2024
38db407
Merge branch 'main' into main-navier
huiyuxie Jun 9, 2024
58de8d6
Merge branch 'main' into main-navier
huiyuxie Jun 11, 2024
4979b27
fix comment
huiyuxie Jun 12, 2024
1400f97
Merge branch 'main' into main-navier
huiyuxie Jun 12, 2024
e653e19
Merge branch 'main' into main-navier
huiyuxie Jun 16, 2024
1efba1e
Merge branch 'main' into main-navier
huiyuxie Jun 16, 2024
782b8a5
part of unit tests
huiyuxie Jun 16, 2024
6c8976d
Merge branch 'main-navier' of https://github.com/huiyuxie/Trixi.jl in…
huiyuxie Jun 16, 2024
f182f89
bc tests
huiyuxie Jun 17, 2024
6006510
apply code change request
huiyuxie Jun 18, 2024
33510d3
Merge branch 'main' into main-navier
huiyuxie Jun 20, 2024
335baf9
complete with issue open
huiyuxie Jun 20, 2024
9b37b10
fix typos
huiyuxie Jun 20, 2024
f6bbe8c
minor fix
huiyuxie Jun 20, 2024
1736f50
Merge branch 'main' into main-navier
huiyuxie Jun 25, 2024
948914e
complete
huiyuxie Jun 25, 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
4 changes: 2 additions & 2 deletions src/equations/compressible_navier_stokes_1d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -171,7 +171,7 @@ function flux(u, gradients, orientation::Integer,
mu = dynamic_viscosity(u, equations)

# viscous flux components in the x-direction
f1 = zero(rho)
f1 = 0
f2 = tau_11 * mu
f3 = (v1 * tau_11 + q1) * mu

Expand Down Expand Up @@ -252,7 +252,7 @@ end
@inline function temperature(u, equations::CompressibleNavierStokesDiffusion1D)
rho, rho_v1, rho_e = u

p = (equations.gamma - 1) * (rho_e - 0.5 * rho_v1^2 / rho)
p = (equations.gamma - 1) * (rho_e - 0.5f0 * rho_v1^2 / rho)
T = p / rho
return T
end
Expand Down
14 changes: 7 additions & 7 deletions src/equations/compressible_navier_stokes_2d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -157,13 +157,13 @@ function flux(u, gradients, orientation::Integer,

# Components of viscous stress tensor

# (4/3 * (v1)_x - 2/3 * (v2)_y)
tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * dv2dy
# (4 * (v1)_x / 3 - 2 * (v2)_y / 3)
tau_11 = 4 * dv1dx / 3 - 2 * dv2dy / 3
# ((v1)_y + (v2)_x)
# stress tensor is symmetric
tau_12 = dv1dy + dv2dx # = tau_21
# (4/3 * (v2)_y - 2/3 * (v1)_x)
tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * dv1dx
tau_22 = 4 * dv2dy / 3 - 2 * dv1dx / 3

# Fick's law q = -kappa * grad(T) = -kappa * grad(p / (R rho))
# with thermal diffusivity constant kappa = gamma μ R / ((gamma-1) Pr)
Expand All @@ -181,7 +181,7 @@ function flux(u, gradients, orientation::Integer,

if orientation == 1
# viscous flux components in the x-direction
f1 = zero(rho)
f1 = 0
f2 = tau_11 * mu
f3 = tau_12 * mu
f4 = (v1 * tau_11 + v2 * tau_12 + q1) * mu
Expand All @@ -190,7 +190,7 @@ function flux(u, gradients, orientation::Integer,
else # if orientation == 2
# viscous flux components in the y-direction
# Note, symmetry is exploited for tau_12 = tau_21
g1 = zero(rho)
g1 = 0
g2 = tau_12 * mu # tau_21 * mu
g3 = tau_22 * mu
g4 = (v1 * tau_12 + v2 * tau_22 + q2) * mu
Expand Down Expand Up @@ -276,7 +276,7 @@ end
@inline function temperature(u, equations::CompressibleNavierStokesDiffusion2D)
rho, rho_v1, rho_v2, rho_e = u

p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2) / rho)
p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2) / rho)
T = p / rho
return T
end
Expand All @@ -285,7 +285,7 @@ end
# Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v

omega = vorticity(u, gradients, equations)
return 0.5 * u[1] * omega^2
return 0.5f0 * u[1] * omega^2
end

@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion2D)
Expand Down
22 changes: 11 additions & 11 deletions src/equations/compressible_navier_stokes_3d.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -162,12 +162,12 @@ function flux(u, gradients, orientation::Integer,
# Components of viscous stress tensor

# Diagonal parts
# (4/3 * (v1)_x - 2/3 * ((v2)_y + (v3)_z)
tau_11 = 4.0 / 3.0 * dv1dx - 2.0 / 3.0 * (dv2dy + dv3dz)
# (4/3 * (v2)_y - 2/3 * ((v1)_x + (v3)_z)
tau_22 = 4.0 / 3.0 * dv2dy - 2.0 / 3.0 * (dv1dx + dv3dz)
# (4/3 * (v3)_z - 2/3 * ((v1)_x + (v2)_y)
tau_33 = 4.0 / 3.0 * dv3dz - 2.0 / 3.0 * (dv1dx + dv2dy)
# (4 * (v1)_x / 3 - 2 * ((v2)_y + (v3)_z)) / 3)
tau_11 = 4 * dv1dx / 3 - 2 * (dv2dy + dv3dz) / 3
# (4 * (v2)_y / 3 - 2 * ((v1)_x + (v3)_z) / 3)
tau_22 = 4 * dv2dy / 3 - 2 * (dv1dx + dv3dz) / 3
# (4 * (v3)_z / 3 - 2 * ((v1)_x + (v2)_y) / 3)
tau_33 = 4 * dv3dz / 3 - 2 * (dv1dx + dv2dy) / 3

# Off diagonal parts, exploit that stress tensor is symmetric
# ((v1)_y + (v2)_x)
Expand All @@ -194,7 +194,7 @@ function flux(u, gradients, orientation::Integer,

if orientation == 1
# viscous flux components in the x-direction
f1 = zero(rho)
f1 = 0
f2 = tau_11 * mu
f3 = tau_12 * mu
f4 = tau_13 * mu
Expand All @@ -204,7 +204,7 @@ function flux(u, gradients, orientation::Integer,
elseif orientation == 2
# viscous flux components in the y-direction
# Note, symmetry is exploited for tau_12 = tau_21
g1 = zero(rho)
g1 = 0
g2 = tau_12 * mu # tau_21 * mu
g3 = tau_22 * mu
g4 = tau_23 * mu
Expand All @@ -214,7 +214,7 @@ function flux(u, gradients, orientation::Integer,
else # if orientation == 3
# viscous flux components in the z-direction
# Note, symmetry is exploited for tau_13 = tau_31, tau_23 = tau_32
h1 = zero(rho)
h1 = 0
h2 = tau_13 * mu # tau_31 * mu
h3 = tau_23 * mu # tau_32 * mu
h4 = tau_33 * mu
Expand Down Expand Up @@ -304,7 +304,7 @@ end
@inline function temperature(u, equations::CompressibleNavierStokesDiffusion3D)
rho, rho_v1, rho_v2, rho_v3, rho_e = u

p = (equations.gamma - 1) * (rho_e - 0.5 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)
p = (equations.gamma - 1) * (rho_e - 0.5f0 * (rho_v1^2 + rho_v2^2 + rho_v3^2) / rho)
T = p / rho
return T
end
Expand All @@ -313,7 +313,7 @@ end
# Enstrophy is 0.5 rho ω⋅ω where ω = ∇ × v

omega = vorticity(u, gradients, equations)
return 0.5 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2)
return 0.5f0 * u[1] * (omega[1]^2 + omega[2]^2 + omega[3]^2)
end

@inline function vorticity(u, gradients, equations::CompressibleNavierStokesDiffusion3D)
Expand Down
71 changes: 70 additions & 1 deletion test/test_type.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -109,7 +109,7 @@ isdir(outdir) && rm(outdir, recursive = true)

x = SVector(zero(RealT))
t = zero(RealT)
u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT))
u = u_ll = u_rr = u_inner = SVector(one(RealT), zero(RealT), zero(RealT))
huiyuxie marked this conversation as resolved.
Show resolved Hide resolved
orientation = 1
directions = [1, 2]
cons = SVector(one(RealT), one(RealT), one(RealT))
Expand Down Expand Up @@ -686,6 +686,75 @@ isdir(outdir) && rm(outdir, recursive = true)
end
end

@timed_testset "Compressible Navier Stokes Diffusion 1D" begin
for RealT in (Float32, Float64)
equations = @inferred CompressibleEulerEquations1D(RealT(1.4))
prandtl_number() = RealT(0.72)
mu() = RealT(0.01)
equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion1D(equations,
mu = mu(),
Prandtl = prandtl_number(),
gradient_variables = GradientVariablesPrimitive())
equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion1D(equations,
mu = mu(),
Prandtl = prandtl_number(),
gradient_variables = GradientVariablesEntropy())
huiyuxie marked this conversation as resolved.
Show resolved Hide resolved

x = SVector(zero(RealT))
t = zero(RealT)
u = u_inner = u_transformed = flux_inner = SVector(one(RealT), zero(RealT),
zero(RealT))
orientation = 1
directions = [1, 2]
gradients = SVector(RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1))

for equations_parabolic in (equations_parabolic_primitive,
equations_parabolic_entropy)
@test eltype(@inferred flux(u, gradients, orientation, equations_parabolic)) ==
RealT

@test eltype(@inferred cons2prim(u, equations_parabolic)) == RealT
@test eltype(@inferred prim2cons(u, equations_parabolic)) == RealT
@test eltype(@inferred cons2entropy(u, equations_parabolic)) == RealT
@test eltype(@inferred entropy2cons(u, equations_parabolic)) == RealT
@test typeof(@inferred Trixi.temperature(u, equations_parabolic)) == RealT

@test eltype(@inferred Trixi.convert_transformed_to_primitive(u_transformed,
equations_parabolic)) ==
RealT
@test eltype(@inferred Trixi.convert_derivative_to_primitive(u, gradients,
equations_parabolic)) ==
RealT
end

# test for GradientVariablesPrimitive
huiyuxie marked this conversation as resolved.
Show resolved Hide resolved
velocity_bc_left_right = NoSlip((x, t, equations) -> initial_condition_navier_stokes_convergence_test(x,
t,
equations)[2])
heat_bc_left = Isothermal((x, t, equations) -> Trixi.temperature(initial_condition_navier_stokes_convergence_test(x,
t,
equations),
equations_parabolic))
heat_bc_right = Adiabatic((x, t, equations) -> zero(RealT))
boundary_condition_left = BoundaryConditionNavierStokesWall(velocity_bc_left_right,
heat_bc_left)
boundary_condition_right = BoundaryConditionNavierStokesWall(velocity_bc_left_right,
heat_bc_right)
gradient = Trixi.Gradient()

for direction in directions
@test eltype(@inferred boundary_condition_left(flux_inner, u_inner,
orientation, direction, x, t,
gradient,
equations_parabolic_primitive)) ==
RealT
end

# test for GradientVariablesEntropy
huiyuxie marked this conversation as resolved.
Show resolved Hide resolved

end
end

@timed_testset "Linear Scalar Advection 1D" begin
for RealT in (Float32, Float64)
equations = @inferred LinearScalarAdvectionEquation1D(RealT(1))
Expand Down
Loading