From 8e1cca58d72c3288e56bfca97714c6ab163406d9 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 24 Jun 2024 21:02:02 -1000 Subject: [PATCH 1/5] start --- src/equations/lattice_boltzmann_2d.jl | 2 +- src/equations/lattice_boltzmann_3d.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl index 272dd897ce3..3ff23f357a4 100644 --- a/src/equations/lattice_boltzmann_2d.jl +++ b/src/equations/lattice_boltzmann_2d.jl @@ -253,7 +253,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 """ diff --git a/src/equations/lattice_boltzmann_3d.jl b/src/equations/lattice_boltzmann_3d.jl index d3eada15f56..6a29a0a7310 100644 --- a/src/equations/lattice_boltzmann_3d.jl +++ b/src/equations/lattice_boltzmann_3d.jl @@ -243,7 +243,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 """ @@ -391,7 +391,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` From aa21f335272c71d0afc8de3ef0acc65a0c9042cd Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Mon, 15 Jul 2024 11:02:56 -1000 Subject: [PATCH 2/5] complete equations --- src/equations/lattice_boltzmann_2d.jl | 19 +++++++---- src/equations/lattice_boltzmann_3d.jl | 45 ++++++++++++++++----------- 2 files changed, 40 insertions(+), 24 deletions(-) diff --git a/src/equations/lattice_boltzmann_2d.jl b/src/equations/lattice_boltzmann_2d.jl index 3ff23f357a4..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 @@ -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 6a29a0a7310..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 @@ -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() From 8de1337e0e2f04fe801306c65c23f6228c37349e Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 16 Jul 2024 22:23:56 -1000 Subject: [PATCH 3/5] revert --- test/test_unit.jl | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) 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 From aff062b42780a2b5dbc591eb58d1bea9e6e98c61 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Tue, 16 Jul 2024 23:51:27 -1000 Subject: [PATCH 4/5] complete 1D --- test/test_type.jl | 48 +++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 48 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index c45a750c07f..e3bea53ab9d 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1517,6 +1517,54 @@ 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 "Linear Scalar Advection 1D" begin for RealT in (Float32, Float64) equations = @inferred LinearScalarAdvectionEquation1D(RealT(1)) From afc30168a4fc61aaf10c1f3c4b7e36eed17d21d9 Mon Sep 17 00:00:00 2001 From: huiyuxie Date: Wed, 17 Jul 2024 00:30:18 -1000 Subject: [PATCH 5/5] complete 2D/3D --- test/test_type.jl | 45 +++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 45 insertions(+) diff --git a/test/test_type.jl b/test/test_type.jl index e3bea53ab9d..24c005d1408 100644 --- a/test/test_type.jl +++ b/test/test_type.jl @@ -1565,6 +1565,51 @@ isdir(outdir) && rm(outdir, recursive = true) 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))