Shan-Chen Two-Phase Single-Component Lattice Boltzmann

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

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

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]\).

def psi(dens):
    return rho0 * (1. - sp.exp(-dens / rho0));

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

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

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

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

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

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

init()
plot_ρs()

Run the simulation until converged

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

assert np.isfinite(dh.gather_array(ρ.name)).all()