Skip to content

Commit 26ff522

Browse files
DanielDoehringranochaSimonCan
authored
Nontrivial Coupling Simulation: LBM + Polytropic Euler (#2421)
* First try coupling LBM Euler * simplify flux * seemingly working version * real coupling example * better examples * test vals * comment * fix * Update examples/structured_2d_dgsem/elixir_lbm_eulerpolytropic_coupled.jl Co-authored-by: Hendrik Ranocha <[email protected]> * Update examples/structured_2d_dgsem/elixir_lbm_lid_driven_cavity.jl * comments + docstrings * comment * simplify * Update examples/structured_2d_dgsem/elixir_lbm_lid_driven_cavity.jl * revisit * fmt * Apply suggestions from code review --------- Co-authored-by: Hendrik Ranocha <[email protected]> Co-authored-by: Simon Candelaresi <[email protected]>
1 parent 8e057e7 commit 26ff522

File tree

6 files changed

+452
-7
lines changed

6 files changed

+452
-7
lines changed

‎examples/structured_2d_dgsem/elixir_lbm_eulerpolytropic_coupled.jl‎

Lines changed: 231 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,231 @@
1+
using OrdinaryDiffEqSSPRK
2+
using Trixi
3+
using LinearAlgebra: norm
4+
5+
###############################################################################
6+
# Semidiscretizations of the polytropic Euler equations and Lattice-Boltzmann method (LBM)
7+
# coupled using converter functions across their respective domains to generate a periodic system.
8+
#
9+
# In this elixir, we have a rectangular domain that is divided into a left and right half.
10+
# On each half of the domain, an independent SemidiscretizationHyperbolic is created for each set of equations.
11+
# The two systems are coupled in the x-direction and are periodic in the y-direction.
12+
# For a high-level overview, see also the figure below:
13+
#
14+
# (-2, 1) ( 2, 1)
15+
# ┌────────────────────┬────────────────────┐
16+
# │ ↑ periodic ↑ │ ↑ periodic ↑ │
17+
# │ │ │
18+
# │ ========= │ ========= │
19+
# │ system #1 │ system #2 │
20+
# | Euler | LBM |
21+
# │ ========= │ ========= │
22+
# │ │ │
23+
# │<-- coupled │<-- coupled │
24+
# │ coupled -->│ coupled -->│
25+
# │ │ │
26+
# │ ↓ periodic ↓ │ ↓ periodic ↓ │
27+
# └────────────────────┴────────────────────┘
28+
# (-2, -1) ( 2, -1)
29+
30+
polydeg = 2
31+
cells_per_dim_per_section = (16, 8)
32+
33+
###########
34+
# system #1
35+
# Euler
36+
###########
37+
38+
### Setup taken from "elixir_eulerpolytropic_isothermal_wave.jl" ###
39+
40+
gamma = 1.0 # Isothermal gas
41+
kappa = 1.0 # Scaling factor for the pressure, must fit to LBM `c_s`
42+
eqs_euler = PolytropicEulerEquations2D(gamma, kappa)
43+
44+
volume_flux = flux_winters_etal
45+
solver_euler = DGSEM(polydeg = polydeg, surface_flux = flux_hll,
46+
volume_integral = VolumeIntegralFluxDifferencing(volume_flux))
47+
48+
# Linear pressure wave/Gaussian bump moving in the positive x-direction.
49+
function initial_condition_pressure_bump(x, t, equations::PolytropicEulerEquations2D)
50+
rho = ((1.0 + 0.01 * exp(-(x[1] - 1)^2 / 0.1)) / equations.kappa)^(1 / equations.gamma)
51+
v1 = ((0.01 * exp(-(x[1] - 1)^2 / 0.1)) / equations.kappa)
52+
v2 = 0.0
53+
54+
return prim2cons(SVector(rho, v1, v2), equations)
55+
end
56+
initial_condition_euler = initial_condition_pressure_bump
57+
58+
coords_min_euler = (-2.0, -1.0)
59+
coords_max_euler = (0.0, 1.0)
60+
mesh_euler = StructuredMesh(cells_per_dim_per_section,
61+
coords_min_euler, coords_max_euler,
62+
periodicity = (false, true))
63+
64+
# Use macroscopic variables derived from LBM populations
65+
# as boundary values for the Euler equations
66+
function coupling_function_LBM2Euler(x, u, equations_other, equations_own)
67+
rho, v1, v2, _ = cons2macroscopic(u, equations_other)
68+
return prim2cons(SVector(rho, v1, v2), equations_own)
69+
end
70+
71+
boundary_conditions_euler = (x_neg = BoundaryConditionCoupled(2, (:end, :i_forward),
72+
Float64,
73+
coupling_function_LBM2Euler),
74+
x_pos = BoundaryConditionCoupled(2, (:begin, :i_forward),
75+
Float64,
76+
coupling_function_LBM2Euler),
77+
y_neg = boundary_condition_periodic,
78+
y_pos = boundary_condition_periodic)
79+
80+
semi_euler = SemidiscretizationHyperbolic(mesh_euler, eqs_euler, initial_condition_euler,
81+
solver_euler;
82+
boundary_conditions = boundary_conditions_euler)
83+
84+
###########
85+
# system #2
86+
# LBM
87+
###########
88+
89+
# Results in c_s = c/sqrt(3) = 1.
90+
# This in turn implies that also in the LBM, p = c_s^2 * rho = 1 * rho = kappa * rho holds
91+
# This is absolutely essential for the correct coupling between the two systems.
92+
c = sqrt(3)
93+
94+
# Reference values `rho0, u0` correspond to the initial condition of the Euler equations.
95+
# The gas should be inviscid (Re = Inf) to be consistent with the inviscid Euler equations.
96+
# The Mach number `Ma` is computed internally from the speed of sound `c_s = c / sqrt(3)` and `u0`.
97+
eqs_lbm = LatticeBoltzmannEquations2D(c = c, Re = Inf, rho0 = 1.0, u0 = 0.0, Ma = nothing)
98+
99+
# Quick & dirty implementation of the `flux_godunov` for Cartesian, yet structured meshes.
100+
@inline function Trixi.flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
101+
equations::LatticeBoltzmannEquations2D)
102+
RealT = eltype(normal_direction)
103+
if isapprox(normal_direction[2], zero(RealT), atol = 2 * eps(RealT))
104+
v_alpha = equations.v_alpha1 * abs(normal_direction[1])
105+
elseif isapprox(normal_direction[1], zero(RealT), atol = 2 * eps(RealT))
106+
v_alpha = equations.v_alpha2 * abs(normal_direction[2])
107+
else
108+
error("Invalid normal direction for flux_godunov: $normal_direction")
109+
end
110+
return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
111+
end
112+
113+
solver_lbm = DGSEM(polydeg = 2, surface_flux = flux_godunov)
114+
115+
function initial_condition_lbm(x, t, equations::LatticeBoltzmannEquations2D)
116+
rho = (1.0 + 0.01 * exp(-(x[1] - 1)^2 / 0.1))
117+
v1 = 0.01 * exp(-(x[1] - 1)^2 / 0.1)
118+
119+
v2 = 0.0
120+
121+
return equilibrium_distribution(rho, v1, v2, equations)
122+
end
123+
124+
coords_min_lbm = (0.0, -1.0)
125+
coords_max_lbm = (2.0, 1.0)
126+
mesh_lbm = StructuredMesh(cells_per_dim_per_section,
127+
coords_min_lbm, coords_max_lbm,
128+
periodicity = (false, true))
129+
130+
# Supply equilibrium (Maxwellian) distribution function computed
131+
# from the Euler-variables as boundary values for the LBM equations
132+
function coupling_function_Euler2LBM(x, u, equations_other, equations_own)
133+
u_prim_euler = cons2prim(u, equations_other)
134+
rho = u_prim_euler[1]
135+
v1 = u_prim_euler[2]
136+
v2 = u_prim_euler[3]
137+
138+
return equilibrium_distribution(rho, v1, v2, equations_own)
139+
end
140+
141+
boundary_conditions_lbm = (x_neg = BoundaryConditionCoupled(1, (:end, :i_forward),
142+
Float64,
143+
coupling_function_Euler2LBM),
144+
x_pos = BoundaryConditionCoupled(1, (:begin, :i_forward),
145+
Float64,
146+
coupling_function_Euler2LBM),
147+
y_neg = boundary_condition_periodic,
148+
y_pos = boundary_condition_periodic)
149+
150+
semi_lbm = SemidiscretizationHyperbolic(mesh_lbm, eqs_lbm, initial_condition_lbm,
151+
solver_lbm;
152+
boundary_conditions = boundary_conditions_lbm)
153+
154+
# Create a semidiscretization that bundles the two semidiscretizations
155+
semi = SemidiscretizationCoupled(semi_euler, semi_lbm)
156+
157+
###############################################################################
158+
# ODE solvers, callbacks etc.
159+
160+
tspan = (0.0, 10.0)
161+
ode = semidiscretize(semi, tspan)
162+
163+
summary_callback = SummaryCallback()
164+
165+
analysis_interval = 100
166+
analysis_callback_euler = AnalysisCallback(semi_euler, interval = analysis_interval)
167+
analysis_callback_lbm = AnalysisCallback(semi_lbm, interval = analysis_interval)
168+
analysis_callback = AnalysisCallbackCoupled(semi,
169+
analysis_callback_euler,
170+
analysis_callback_lbm)
171+
172+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
173+
174+
# Need to implement `cons2macroscopic` for `PolytropicEulerEquations2D`
175+
# in order to be able to use this in the `SaveSolutionCallback` below
176+
@inline function Trixi.cons2macroscopic(u, equations::PolytropicEulerEquations2D)
177+
u_prim = cons2prim(u, equations)
178+
p = pressure(u, equations)
179+
return SVector(u_prim[1], u_prim[2], u_prim[3], p)
180+
end
181+
function Trixi.varnames(::typeof(cons2macroscopic), ::PolytropicEulerEquations2D)
182+
("rho", "v1", "v2", "p")
183+
end
184+
185+
save_solution = SaveSolutionCallback(interval = 50,
186+
save_initial_solution = true,
187+
save_final_solution = true,
188+
solution_variables = cons2macroscopic)
189+
190+
cfl = 2.0
191+
stepsize_callback = StepsizeCallback(cfl = cfl)
192+
193+
# Need special version of the LBM collision callback for a `SemidiscretizationCoupled`
194+
@inline function Trixi.lbm_collision_callback(integrator)
195+
dt = get_proposed_dt(integrator)
196+
semi_coupled = integrator.p # Here `p` is the `SemidiscretizationCoupled`
197+
u_ode_full = integrator.u # ODE Vector for the entire coupled system
198+
for (semi_index, semi_i) in enumerate(semi_coupled.semis)
199+
mesh, equations, solver, cache = Trixi.mesh_equations_solver_cache(semi_i)
200+
if equations isa LatticeBoltzmannEquations2D
201+
@unpack collision_op = equations
202+
203+
u_ode_i = Trixi.get_system_u_ode(u_ode_full, semi_index, semi_coupled)
204+
u = Trixi.wrap_array(u_ode_i, mesh, equations, solver, cache)
205+
206+
Trixi.@trixi_timeit Trixi.timer() "LBM collision" Trixi.apply_collision!(u, dt,
207+
collision_op,
208+
mesh,
209+
equations,
210+
solver,
211+
cache)
212+
end
213+
end
214+
215+
return nothing
216+
end
217+
218+
collision_callback = LBMCollisionCallback()
219+
220+
callbacks = CallbackSet(summary_callback,
221+
analysis_callback, alive_callback,
222+
save_solution,
223+
stepsize_callback,
224+
collision_callback)
225+
226+
###############################################################################
227+
# run the simulation
228+
229+
sol = solve(ode, SSPRK83();
230+
dt = 0.01, # solve needs some value here but it will be overwritten by the stepsize_callback
231+
ode_default_options()..., callback = callbacks);

‎examples/structured_2d_dgsem/elixir_lbm_lid_driven_cavity.jl‎

Lines changed: 128 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,128 @@
1+
using OrdinaryDiffEqLowStorageRK
2+
using Trixi
3+
4+
###############################################################################
5+
# semidiscretization of the Lattice-Boltzmann equations for the D2Q9 scheme
6+
7+
equations = LatticeBoltzmannEquations2D(Ma = 0.1, Re = 1000)
8+
9+
"""
10+
initial_condition_lid_driven_cavity(x, t, equations::LatticeBoltzmannEquations2D)
11+
12+
Initial state for a lid-driven cavity flow setup. To be used in combination with
13+
[`boundary_condition_lid_driven_cavity`](@ref) and [`boundary_condition_noslip_wall`](@ref).
14+
"""
15+
function initial_condition_lid_driven_cavity(x, t, equations::LatticeBoltzmannEquations2D)
16+
@unpack L, u0, nu = equations
17+
18+
rho = 1
19+
v1 = 0
20+
v2 = 0
21+
22+
return equilibrium_distribution(rho, v1, v2, equations)
23+
end
24+
initial_condition = initial_condition_lid_driven_cavity
25+
26+
"""
27+
boundary_condition_lid_driven_cavity(u_inner, orientation, direction, x, t,
28+
surface_flux_function,
29+
equations::LatticeBoltzmannEquations2D)
30+
31+
Boundary condition for a lid-driven cavity flow setup, where the top lid (+y boundary) is a moving
32+
no-slip wall. To be used in combination with [`initial_condition_lid_driven_cavity`](@ref).
33+
"""
34+
function boundary_condition_lid_driven_cavity(u_inner, orientation, direction, x, t,
35+
surface_flux_function,
36+
equations::LatticeBoltzmannEquations2D)
37+
return boundary_condition_moving_wall_ypos(u_inner, orientation, direction, x, t,
38+
surface_flux_function, equations)
39+
end
40+
41+
function boundary_condition_moving_wall_ypos(u_inner, orientation, direction, x, t,
42+
surface_flux_function,
43+
equations::LatticeBoltzmannEquations2D)
44+
@assert direction==4 "moving wall assumed in +y direction"
45+
46+
@unpack rho0, u0, weights, c_s = equations
47+
cs_squared = c_s^2
48+
49+
pdf1 = u_inner[3] + 2 * weights[1] * rho0 * u0 / cs_squared
50+
pdf2 = u_inner[2] # outgoing
51+
pdf3 = u_inner[1] + 2 * weights[3] * rho0 * (-u0) / cs_squared
52+
pdf4 = u_inner[2]
53+
pdf5 = u_inner[5] # outgoing
54+
pdf6 = u_inner[6] # outgoing
55+
pdf7 = u_inner[5] + 2 * weights[7] * rho0 * (-u0) / cs_squared
56+
pdf8 = u_inner[6] + 2 * weights[8] * rho0 * u0 / cs_squared
57+
pdf9 = u_inner[9]
58+
59+
u_boundary = SVector(pdf1, pdf2, pdf3, pdf4, pdf5, pdf6, pdf7, pdf8, pdf9)
60+
61+
# Calculate boundary flux (u_inner is "left" of boundary, u_boundary is "right" of boundary)
62+
return surface_flux_function(u_inner, u_boundary, orientation, equations)
63+
end
64+
boundary_conditions = (x_neg = boundary_condition_noslip_wall,
65+
x_pos = boundary_condition_noslip_wall,
66+
y_neg = boundary_condition_noslip_wall,
67+
y_pos = boundary_condition_lid_driven_cavity)
68+
69+
# Quick & dirty implementation of the `flux_godunov` for Cartesian, yet structured meshes.
70+
@inline function Trixi.flux_godunov(u_ll, u_rr, normal_direction::AbstractVector,
71+
equations::LatticeBoltzmannEquations2D)
72+
RealT = eltype(normal_direction)
73+
if isapprox(normal_direction[2], zero(RealT), atol = 10 * eps(RealT))
74+
v_alpha = equations.v_alpha1 * abs(normal_direction[1])
75+
elseif isapprox(normal_direction[1], zero(RealT), atol = 10 * eps(RealT))
76+
v_alpha = equations.v_alpha2 * abs(normal_direction[2])
77+
else
78+
error("Invalid normal direction for flux_godunov: $normal_direction")
79+
end
80+
return 0.5f0 * (v_alpha .* (u_ll + u_rr) - abs.(v_alpha) .* (u_rr - u_ll))
81+
end
82+
83+
solver = DGSEM(polydeg = 5, surface_flux = flux_godunov)
84+
85+
cells_per_dimension = (16, 16)
86+
coordinates_min = (0.0, 0.0)
87+
coordinates_max = (1.0, 1.0)
88+
mesh = StructuredMesh(cells_per_dimension,
89+
coordinates_min, coordinates_max,
90+
periodicity = false)
91+
92+
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver,
93+
boundary_conditions = boundary_conditions)
94+
95+
###############################################################################
96+
# ODE solvers, callbacks etc.
97+
98+
tspan = (0.0, 1.0)
99+
ode = semidiscretize(semi, tspan)
100+
101+
summary_callback = SummaryCallback()
102+
103+
analysis_interval = 100_000
104+
analysis_callback = AnalysisCallback(semi, interval = analysis_interval)
105+
106+
alive_callback = AliveCallback(analysis_interval = analysis_interval)
107+
108+
save_solution = SaveSolutionCallback(interval = 5000,
109+
save_initial_solution = true,
110+
save_final_solution = true,
111+
solution_variables = cons2macroscopic)
112+
113+
stepsize_callback = StepsizeCallback(cfl = 1.0)
114+
115+
collision_callback = LBMCollisionCallback()
116+
117+
callbacks = CallbackSet(summary_callback,
118+
analysis_callback, alive_callback,
119+
save_solution,
120+
stepsize_callback,
121+
collision_callback)
122+
123+
###############################################################################
124+
# run the simulation
125+
126+
sol = solve(ode, CarpenterKennedy2N54(williamson_condition = false);
127+
dt = 1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
128+
ode_default_options()..., callback = callbacks);

‎src/equations/lattice_boltzmann_2d.jl‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -140,7 +140,12 @@ function varnames(::typeof(cons2prim), equations::LatticeBoltzmannEquations2D)
140140
varnames(cons2cons, equations)
141141
end
142142

143-
# Convert conservative variables to macroscopic
143+
"""
144+
cons2macroscopic(u, equations::LatticeBoltzmannEquations2D)
145+
146+
Convert the conservative variables `u` (the particle distribution functions)
147+
to the macroscopic variables (density, velocity_1, velocity_2, pressure).
148+
"""
144149
@inline function cons2macroscopic(u, equations::LatticeBoltzmannEquations2D)
145150
rho = density(u, equations)
146151
v1, v2 = velocity(u, equations)
@@ -308,7 +313,7 @@ Calculate the macroscopic pressure from the density `rho` or the particle distr
308313
`u`.
309314
"""
310315
@inline function pressure(rho::Real, equations::LatticeBoltzmannEquations2D)
311-
rho * equations.c_s^2
316+
return rho * equations.c_s^2
312317
end
313318
@inline function pressure(u, equations::LatticeBoltzmannEquations2D)
314319
pressure(density(u, equations), equations)

0 commit comments

Comments
 (0)