Skip to content

Commit

Permalink
Add numerical support of other real types (navier) (#1968)
Browse files Browse the repository at this point in the history 
* start

* complete equations

* fix comment

* part of unit tests

* bc tests

* apply code change request

Co-authored-by: Hendrik Ranocha <[email protected]>

* complete with issue open

* fix typos

* minor fix

* complete

---------

Co-authored-by: Hendrik Ranocha <[email protected]>
  • Loading branch information
@huiyuxie @ranocha
huiyuxie and ranocha authored Jun 25, 2024
1 parent e114d01 commit 0358146
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 20 deletions.
4 changes: 2 additions & 2 deletions src/equations/compressible_navier_stokes_1d.jl
View file
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
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
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
142 changes: 142 additions & 0 deletions test/test_type.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -686,6 +686,148 @@ 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())

u = u_transformed = SVector(one(RealT), zero(RealT),
zero(RealT))
orientation = 1
gradients = SVector(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

# TODO: BC tests for GradientVariablesPrimitive
# TODO: BC tests for GradientVariablesEntropy
end
end

@timed_testset "Compressible Navier Stokes Diffusion 2D" begin
for RealT in (Float32, Float64)
equations = @inferred CompressibleEulerEquations2D(RealT(1.4))
prandtl_number = RealT(0.72)
mu = RealT(0.01)
equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion2D(equations,
mu = mu,
Prandtl = prandtl_number,
gradient_variables = GradientVariablesPrimitive())
equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion2D(equations,
mu = mu,
Prandtl = prandtl_number,
gradient_variables = GradientVariablesEntropy())

u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT))
orientations = [1, 2]
gradient = SVector(RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1))
gradients = SVector(gradient, gradient)

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

@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 typeof(@inferred Trixi.enstrophy(u, gradients, equations_parabolic)) ==
RealT
@test typeof(@inferred Trixi.vorticity(u, gradients, 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, gradient,
equations_parabolic)) ==
RealT
end

# TODO: BC tests for GradientVariablesPrimitive
# TODO: BC tests for GradientVariablesEntropy
end
end

@timed_testset "Compressible Navier Stokes Diffusion 3D" begin
for RealT in (Float32, Float64)
equations = @inferred CompressibleEulerEquations3D(RealT(1.4))
prandtl_number = RealT(0.72)
mu = RealT(0.01)
equations_parabolic_primitive = @inferred CompressibleNavierStokesDiffusion3D(equations,
mu = mu,
Prandtl = prandtl_number,
gradient_variables = GradientVariablesPrimitive())
equations_parabolic_entropy = @inferred CompressibleNavierStokesDiffusion3D(equations,
mu = mu,
Prandtl = prandtl_number,
gradient_variables = GradientVariablesEntropy())

u = u_transformed = SVector(one(RealT), zero(RealT), zero(RealT), zero(RealT),
zero(RealT))
orientations = [1, 2, 3]
gradient = SVector(RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1), RealT(0.1))
gradients = SVector(gradient, gradient, gradient)

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

@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 typeof(@inferred Trixi.enstrophy(u, gradients, equations_parabolic)) ==
RealT
@test eltype(@inferred Trixi.vorticity(u, gradients, 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, gradient,
equations_parabolic)) ==
RealT
end

# TODO: BC tests for GradientVariablesPrimitive
# TODO: BC tests for GradientVariablesEntropy
end
end

@timed_testset "Hyperbolic Diffusion 1D" begin
for RealT in (Float32, Float64)
nu = one(RealT)
Expand Down

0 comments on commit 0358146

Please sign in to comment.