Shan-Chen Two-Phase Single-Component Lattice Boltzmann¶
[1]:
from lbmpy.session import *
from lbmpy.updatekernels import create_stream_pull_with_output_kernel
from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
from lbmpy.maxwellian_equilibrium import get_weights
This is based on section 9.3.2 of Krüger et al.’s “The Lattice Boltzmann Method”, Springer 2017 (http://www.lbmbook.com). Sample code is available at https://github.com/lbm-principles-practice/code/.
Parameters¶
[2]:
N = 64
omega_a = 1.
g_aa = -4.7
rho0 = 1.
stencil = LBStencil(Stencil.D2Q9)
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
Data structures¶
[3]:
dh = ps.create_data_handling((N,) * stencil.D, periodicity=True, default_target=ps.Target.CPU)
src = dh.add_array('src', values_per_cell=stencil.Q)
dst = dh.add_array_like('dst', 'src')
ρ = dh.add_array('rho')
Force & combined velocity¶
The force on the fluid is \(\vec{F}_A(\vec{x})=-\psi(\rho_A(\vec{x}))g_{AA}\sum\limits_{i=1}^{q}w_i\psi(\rho_A(\vec{x}+\vec{c}_i))\vec{c}_i\) with \(\psi(\rho)=\rho_0\left[1-\exp(-\rho/\rho_0)\right]\).
[4]:
def psi(dens):
return rho0 * (1. - sp.exp(-dens / rho0));
[5]:
zero_vec = sp.Matrix([0] * stencil.D)
force = sum((psi(ρ[d]) * w_d * sp.Matrix(d)
for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa
Kernels¶
[6]:
lbm_config = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
force_model=ForceModel.GUO, force=force, kernel_type='collide_only')
collision = create_lb_update_rule(lbm_config=lbm_config,
optimization={'symbolic_field': src})
stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})
config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=False)
stream_kernel = ps.create_kernel(stream, config=config).compile()
collision_kernel = ps.create_kernel(collision, config=config).compile()
Initialization¶
[7]:
method_without_force = create_lb_method(LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True))
init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0),
pdfs=src.center_vector, density=ρ.center)
init_kernel = ps.create_kernel(init_assignments, ghost_layers=0, config=config).compile()
[8]:
def init():
for x in range(N):
for y in range(N):
if (x-N/2)**2 + (y-N/2)**2 <= 15**2:
dh.fill(ρ.name, 2.1, slice_obj=[x,y])
else:
dh.fill(ρ.name, 0.15, slice_obj=[x,y])
dh.run_kernel(init_kernel)
Timeloop¶
[9]:
sync_pdfs = dh.synchronization_function([src.name])
sync_ρs = dh.synchronization_function([ρ.name])
def time_loop(steps):
dh.all_to_gpu()
for i in range(steps):
sync_ρs()
dh.run_kernel(collision_kernel)
sync_pdfs()
dh.run_kernel(stream_kernel)
dh.swap(src.name, dst.name)
dh.all_to_cpu()
[10]:
def plot_ρs():
plt.figure(dpi=200)
plt.title("$\\rho$")
plt.scalar_field(dh.gather_array(ρ.name), vmin=0, vmax=2.5)
plt.colorbar()
Run the simulation¶
Initial state¶
[11]:
init()
plot_ρs()

Run the simulation until converged¶
[12]:
init()
time_loop(1000)
plot_ρs()

[13]:
assert np.isfinite(dh.gather_array(ρ.name)).all()