From 33aee2aca1289f6169f7e72db841e0b7e081bde2 Mon Sep 17 00:00:00 2001 From: Huiyu Xie Date: Wed, 17 Jul 2024 04:43:22 -1000 Subject: [PATCH] Add numerical support of other real types (`lattice`) (#1991) * start * complete equations * revert * complete 1D * complete 2D/3D --- src/equations/lattice_boltzmann_2d.jl | 21 ++++-- src/equations/lattice_boltzmann_3d.jl | 49 ++++++++------ test/test_type.jl | 93 +++++++++++++++++++++++++++ test/test_unit.jl | 18 ++++-- 4 files changed, 148 insertions(+), 33 deletions(-) diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl index 272dd897ce3..bad417d9c84 100644 --- a/src/equations/lattice_boltzmann_2d.jl +++ b/src/equations/lattice_boltzmann_2d.jl @@ -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) @@ -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, @@ -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 @@ -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 """ @@ -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() diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl index d3eada15f56..bf4c365e0fe 100644 --- a/src/equations/lattice_boltzmann_3d.jl +++ b/src/equations/lattice_boltzmann_3d.jl @@ -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) @@ -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, @@ -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 @@ -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 """ @@ -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() @@ -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` diff --git a/test/test_type.jl b/test/test_type.jl index c45a750c07f..24c005d1408 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -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)) diff --git a/test/test_unit.jl b/test/test_unit.jl index 09cf506576c..aa66ce2556f 100644 --- a/test/test_unit.jl +++ b/test/test_unit.jl @@ -439,14 +439,17 @@ 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 @@ -454,14 +457,17 @@ end # 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