Skip to content

Commit

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

* complete equations

* revert

* complete 1D

* complete 2D/3D
  • Loading branch information
@huiyuxie
huiyuxie authored Jul 17, 2024
1 parent d0d704d commit 33aee2a
Show file tree
Hide file tree
Showing 4 changed files with 148 additions and 33 deletions.
21 changes: 14 additions & 7 deletions src/equations/lattice_boltzmann_2d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,12 +101,16 @@ function LatticeBoltzmannEquations2D(; Ma, Re, collision_op = collision_bgk,
# The relation between the isothermal speed of sound `c_s` and the mean thermal molecular velocity
# `c` depends on the used phase space discretization, and is valid for D2Q9 (and others). For
# details, see, e.g., [3] in the docstring above.
c_s = c / sqrt(3)
# c_s = c / sqrt(3)

# Calculate missing quantities
if isnothing(Ma)
RealT = eltype(u0)
c_s = c / sqrt(convert(RealT, 3))
Ma = u0 / c_s
elseif isnothing(u0)
RealT = eltype(Ma)
c_s = c / sqrt(convert(RealT, 3))
u0 = Ma * c_s
end
if isnothing(Re)
Expand All @@ -119,9 +123,10 @@ function LatticeBoltzmannEquations2D(; Ma, Re, collision_op = collision_bgk,
Ma, Re, c, L, rho0, u0, nu = promote(Ma, Re, c, L, rho0, u0, nu)

# Source for weights and speeds: [4] in the docstring above
weights = SVector(1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36, 1 / 36, 4 / 9)
v_alpha1 = SVector(c, 0, -c, 0, c, -c, -c, c, 0)
v_alpha2 = SVector(0, c, 0, -c, c, c, -c, -c, 0)
weights = SVector{9, RealT}(1 / 9, 1 / 9, 1 / 9, 1 / 9, 1 / 36, 1 / 36, 1 / 36,
1 / 36, 4 / 9)
v_alpha1 = SVector{9, RealT}(c, 0, -c, 0, c, -c, -c, c, 0)
v_alpha2 = SVector{9, RealT}(0, c, 0, -c, c, c, -c, -c, 0)

LatticeBoltzmannEquations2D(c, c_s, rho0, Ma, u0, Re, L, nu,
weights, v_alpha1, v_alpha2,
Expand Down Expand Up @@ -154,7 +159,9 @@ A constant initial condition to test free-stream preservation.
"""
function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations2D)
@unpack u0 = equations
rho = pi

RealT = eltype(x)
rho = convert(RealT, pi)
v1 = u0
v2 = u0

Expand Down Expand Up @@ -253,7 +260,7 @@ end
else
v_alpha = equations.v_alpha2
end
return 0.5 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
end

"""
Expand Down Expand Up @@ -358,7 +365,7 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
@inline function collision_bgk(u, dt, equations::LatticeBoltzmannEquations2D)
@unpack c_s, nu = equations
tau = nu / (c_s^2 * dt)
return -(u - equilibrium_distribution(u, equations)) / (tau + 1 / 2)
return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
end

@inline have_constant_speed(::LatticeBoltzmannEquations2D) = True()
Expand Down
49 changes: 29 additions & 20 deletions src/equations/lattice_boltzmann_3d.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -141,12 +141,16 @@ function LatticeBoltzmannEquations3D(; Ma, Re, collision_op = collision_bgk,
# The relation between the isothermal speed of sound `c_s` and the mean thermal molecular velocity
# `c` depends on the used phase space discretization, and is valid for D3Q27 (and others). For
# details, see, e.g., [3] in the docstring above.
c_s = c / sqrt(3)
# c_s = c / sqrt(3)

# Calculate missing quantities
if isnothing(Ma)
RealT = eltype(u0)
c_s = c / sqrt(convert(RealT, 3))
Ma = u0 / c_s
elseif isnothing(u0)
RealT = eltype(Ma)
c_s = c / sqrt(convert(RealT, 3))
u0 = Ma * c_s
end
if isnothing(Re)
Expand All @@ -159,21 +163,24 @@ function LatticeBoltzmannEquations3D(; Ma, Re, collision_op = collision_bgk,
Ma, Re, c, L, rho0, u0, nu = promote(Ma, Re, c, L, rho0, u0, nu)

# Source for weights and speeds: [4] in docstring above
weights = SVector(2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 1 / 54, 1 / 54,
1 / 54,
1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54,
1 / 54,
1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216,
1 / 216, 8 / 27)
v_alpha1 = SVector(c, -c, 0, 0, 0, 0, c, -c, c,
-c, 0, 0, c, -c, c, -c, 0, 0,
c, -c, c, -c, c, -c, -c, c, 0)
v_alpha2 = SVector(0, 0, c, -c, 0, 0, c, -c, 0,
0, c, -c, -c, c, 0, 0, c, -c,
c, -c, c, -c, -c, c, c, -c, 0)
v_alpha3 = SVector(0, 0, 0, 0, c, -c, 0, 0, c,
-c, c, -c, 0, 0, -c, c, -c, c,
c, -c, -c, c, c, -c, c, -c, 0)
weights = SVector{27, RealT}(2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 2 / 27, 1 / 54,
1 / 54,
1 / 54,
1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54, 1 / 54,
1 / 54,
1 / 54,
1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216, 1 / 216,
1 / 216,
1 / 216, 8 / 27)
v_alpha1 = SVector{27, RealT}(c, -c, 0, 0, 0, 0, c, -c, c,
-c, 0, 0, c, -c, c, -c, 0, 0,
c, -c, c, -c, c, -c, -c, c, 0)
v_alpha2 = SVector{27, RealT}(0, 0, c, -c, 0, 0, c, -c, 0,
0, c, -c, -c, c, 0, 0, c, -c,
c, -c, c, -c, -c, c, c, -c, 0)
v_alpha3 = SVector{27, RealT}(0, 0, 0, 0, c, -c, 0, 0, c,
-c, c, -c, 0, 0, -c, c, -c, c,
c, -c, -c, c, c, -c, c, -c, 0)

LatticeBoltzmannEquations3D(c, c_s, rho0, Ma, u0, Re, L, nu,
weights, v_alpha1, v_alpha2, v_alpha3,
Expand Down Expand Up @@ -206,7 +213,9 @@ A constant initial condition to test free-stream preservation.
"""
function initial_condition_constant(x, t, equations::LatticeBoltzmannEquations3D)
@unpack u0 = equations
rho = pi

RealT = eltype(x)
rho = convert(RealT, pi)
v1 = u0
v2 = u0
v3 = u0
Expand Down Expand Up @@ -243,7 +252,7 @@ end
else # z-direction
v_alpha = equations.v_alpha3
end
return 0.5 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
end

"""
Expand Down Expand Up @@ -369,7 +378,7 @@ Collision operator for the Bhatnagar, Gross, and Krook (BGK) model.
@inline function collision_bgk(u, dt, equations::LatticeBoltzmannEquations3D)
@unpack c_s, nu = equations
tau = nu / (c_s^2 * dt)
return -(u - equilibrium_distribution(u, equations)) / (tau + 1 / 2)
return -(u - equilibrium_distribution(u, equations)) / (tau + 0.5f0)
end

@inline have_constant_speed(::LatticeBoltzmannEquations3D) = True()
Expand All @@ -391,7 +400,7 @@ end
rho = density(u, equations)
v1, v2, v3 = velocity(u, equations)

return 0.5 * (v1^2 + v2^2 + v3^2) / rho / equations.rho0
return 0.5f0 * (v1^2 + v2^2 + v3^2) / rho / equations.rho0
end

# Calculate nondimensionalized kinetic energy for a conservative state `u`
Expand Down
93 changes: 93 additions & 0 deletions test/test_type.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1517,6 +1517,99 @@ isdir(outdir) && rm(outdir, recursive = true)
end
end

@timed_testset "Lattice Boltzmann 2D" begin
for RealT in (Float32, Float64)
equations = @inferred LatticeBoltzmannEquations2D(Ma = RealT(0.1), Re = 1000)

x = SVector(zero(RealT), zero(RealT))
t = zero(RealT)
u = u_ll = u_rr = u_inner = SVector(one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT))
orientations = [1, 2]
directions = [1, 2, 3, 4]
p = rho = dt = one(RealT)

surface_flux_function = flux_godunov

@test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT

for orientation in orientations
for direction in directions
@test eltype(@inferred boundary_condition_noslip_wall(u_inner,
orientation,
direction, x, t,
surface_flux_function,
equations)) ==
RealT
end

@test eltype(@inferred flux(u, orientation, equations)) == RealT
@test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) ==
RealT

@test typeof(@inferred velocity(u, orientation, equations)) == RealT
end

@test typeof(@inferred density(p, equations)) == RealT
@test typeof(@inferred density(u, equations)) == RealT
@test eltype(@inferred velocity(u, equations)) == RealT
@test typeof(@inferred pressure(rho, equations)) == RealT
@test typeof(@inferred pressure(u, equations)) == RealT

@test eltype(@inferred Trixi.collision_bgk(u, dt, equations)) == RealT

@test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT
@test eltype(@inferred cons2prim(u, equations)) == RealT
@test eltype(@inferred cons2entropy(u, equations)) == RealT
end
end

@timed_testset "Lattice Boltzmann 3D" begin
for RealT in (Float32, Float64)
equations = @inferred LatticeBoltzmannEquations3D(Ma = RealT(0.1), Re = 1000)

x = SVector(zero(RealT), zero(RealT), zero(RealT))
t = zero(RealT)
u = u_ll = u_rr = SVector(one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT),
one(RealT), one(RealT), one(RealT))
orientations = [1, 2, 3]
p = rho = dt = one(RealT)

@test eltype(@inferred initial_condition_constant(x, t, equations)) == RealT

for orientation in orientations
@test eltype(@inferred flux(u, orientation, equations)) == RealT
@test eltype(@inferred flux_godunov(u_ll, u_rr, orientation, equations)) ==
RealT

@test typeof(@inferred velocity(u, orientation, equations)) == RealT
end

@test typeof(@inferred density(p, equations)) == RealT
@test typeof(@inferred density(u, equations)) == RealT
@test eltype(@inferred velocity(u, equations)) == RealT
@test typeof(@inferred pressure(rho, equations)) == RealT
@test typeof(@inferred pressure(u, equations)) == RealT

@test eltype(@inferred Trixi.collision_bgk(u, dt, equations)) == RealT

@test eltype(@inferred Trixi.max_abs_speeds(equations)) == RealT
@test eltype(@inferred cons2prim(u, equations)) == RealT
@test eltype(@inferred cons2entropy(u, equations)) == RealT
@test typeof(@inferred energy_kinetic(u, equations)) == RealT
@test typeof(@inferred Trixi.energy_kinetic_nondimensional(u, equations)) ==
RealT
end
end

@timed_testset "Linear Scalar Advection 1D" begin
for RealT in (Float32, Float64)
equations = @inferred LinearScalarAdvectionEquation1D(RealT(1))
Expand Down
18 changes: 12 additions & 6 deletions test/test_unit.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -439,29 +439,35 @@ end
# Neither Mach number nor velocity set
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = nothing, Re = 1000)
# Both Mach number and velocity set
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, u0 = 1)
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000,
u0 = 1.0)
# Neither Reynolds number nor viscosity set
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = nothing)
# Both Reynolds number and viscosity set
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000, nu = 1)
@test_throws ErrorException LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000,
nu = 1.0)

# No non-dimensional values set
@test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa
@test LatticeBoltzmannEquations2D(Ma = nothing, Re = nothing, u0 = 1.0,
nu = 1.0) isa
LatticeBoltzmannEquations2D
end

@timed_testset "LBM 3D constructor" begin
# Neither Mach number nor velocity set
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = nothing, Re = 1000)
# Both Mach number and velocity set
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, u0 = 1)
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000,
u0 = 1.0)
# Neither Reynolds number nor viscosity set
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = nothing)
# Both Reynolds number and viscosity set
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000, nu = 1)
@test_throws ErrorException LatticeBoltzmannEquations3D(Ma = 0.1, Re = 1000,
nu = 1.0)

# No non-dimensional values set
@test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1, nu = 1) isa
@test LatticeBoltzmannEquations3D(Ma = nothing, Re = nothing, u0 = 1.0,
nu = 1.0) isa
LatticeBoltzmannEquations3D
end

Expand Down

0 comments on commit 33aee2a

Please sign in to comment.