# 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)



## 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()