# 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 = get_stencil("D2Q9")
weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))


## Data structures¶

[3]:

dim = len(stencil[0])

dh = ps.create_data_handling((N,)*dim, periodicity=True, default_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] * dh.dim)

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]:

collision = create_lb_update_rule(stencil=stencil,
relaxation_rate=omega_a,
compressible=True,
force_model='guo',
force=force,
kernel_type='collide_only',
optimization={'symbolic_field': src})

stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})

opts = {'cpu_openmp': False,
'target': dh.default_target}

stream_kernel = ps.create_kernel(stream, **opts).compile()
collision_kernel = ps.create_kernel(collision, **opts).compile()


## Initialization¶

[7]:

method_without_force = create_lb_method(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).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.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()


### Check the first time step against reference data¶

The reference data was obtained with the sample code after making the following changes:

const int nsteps = 1000;
const int noutput = 1;


Remove the next cell if you changed the parameters at the beginning of this notebook.

[12]:

init()
time_loop(1)
ref = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.136756, 0.220324, 1.2382, 2.26247, 2.26183, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.26183, 2.26247, 1.2382, 0.220324, 0.136756, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])

assert np.allclose(dh.gather_array(ρ.name)[N//2], ref)


### Run the simulation until converged¶

[13]:

init()
time_loop(1000)
plot_ρs()