[1]:
from lbmpy.session import *

# Tutorial 04: The Cumulant Lattice Boltzmann Method in lbmpy¶

## A) Principles of the centered cumulant collision operator¶

Recently, an advanced Lattice Boltzmann collision operator based on relaxation in cumulant space has gained popularity. Similar to moments, cumulants are statistical quantities inherent to a probability distribution. A significant advantage of the cumulants is that they are statistically independent by construction. Moments can be defined by using the so-called moment generating function, which for the discrete particle distribution of the LB equation is stated as

\begin{align} M( \mathbf{X} ) = \sum_i f_i \exp \left( \mathbf{c}_i \cdot \mathbf{X} \right) \end{align}

The raw moments $$m_{\alpha \beta \gamma}$$ can be expressed as its derivatives, evaluated at zero:

\begin{align} m_{\alpha \beta \gamma} = \left. \frac{\partial^{\alpha}}{\partial X^{\alpha}} \frac{\partial^{\beta}}{\partial Y^{\beta}} \frac{\partial^{\gamma}}{\partial Z^{\gamma}} M(X, Y, Z) \right\vert_{\mathbf{X}=0} \end{align}

The cumulant-generating function is defined as the natural logarithm of this moment-generating function, and the cumulants $$c_{\alpha \beta \gamma}$$ are defined as its derivatives evaluated at zero:

\begin{split}\begin{align} C(\mathbf{X}) :=& \, \log ( M(\mathbf{X}) ) \\ c_{\alpha \beta \gamma} =& \, \left. \frac{\partial^{\alpha}}{\partial X^{\alpha}} \frac{\partial^{\beta}}{\partial Y^{\beta}} \frac{\partial^{\gamma}}{\partial Z^{\gamma}} C(X, Y, Z) \right\vert_{\mathbf{X}=0} \end{align}\end{split}

Other than with moments, there is no straightforward way to express cumulants in terms of the populations. However, their generating functions can be related to allowing the computation of cumulants from both raw and central moments, computed from populations. In lbmpy, the transformations from populations to cumulants and back are implemented using central moments as intermediaries. All cumulants of orders 2 and 3 are equal to their corresponding central moments, up to the density $$\rho$$ as a proportionality factor.

The central moment-generating function $$K$$ can be related to the moment-generating function through $$K(\mathbf{X}) = \exp( - \mathbf{X} \cdot \mathbf{u} ) M(\mathbf{X})$$. It is possible to recombine the equation with the definition of the cumulant-generating function

\begin{align} C( \mathbf{X} ) = \mathbf{X} \cdot \mathbf{u} + \log( K( \mathbf{X} ) ). \end{align}

Derivatives of $$C$$ can thus be expressed in terms of derivatives of $$K$$, directly yielding equations of the cumulants in terms of central moments.

In the cumulant lattice Boltzmann method, forces are applied symmetrically in central moment space.

## B) Method Creation¶

Using the create_lb_method interface, creation of a cumulant-based lattice Boltzmann method in lbmpy is straightforward. Cumulants can either be relaxed in their raw (monomial forms) or in polynomial combinations. Both variants are available as predefined setups. Attention: Cumulant-based methods are only available for the compressible equilibrium.

### Relaxation of Monomial Cumulants¶

Monomial cumulant relaxation is available through the method="monomial_cumulant" parameter setting. This variant requires fewer transformation steps and is a little faster than polynomial relaxation, but it does not allow separation of bulk and shear viscosity. Default monomial cumulant sets are available for the D2Q9, D3Q19 and D3Q27 stencils. It is also possible to define a custom set of monomial cumulants.

When creating a monomial cumulant method, one option is to specify only a single relaxation rate which will then be assigned to all cumulants related to the shear viscosity. In this case, all other non-conserved cumulants will be relaxed to equilibrium. Alternatively, individual relaxation rates for all non-conserved cumulants can be specified. The conserved cumulants are set to zero by default to save computational cost. They can be adjusted with set_zeroth_moment_relaxation_rate, set_first_moment_relaxation_rate and set_conserved_moments_relaxation_rate.

[2]:
lbm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, relaxation_rate=sp.Symbol('omega_v'), compressible=True)
method_monomial = create_lb_method(lbm_config=lbm_config)
method_monomial
[2]:
 Cumulant-Based Method Stencil: D2Q9 Zero-Centered Storage: ✓ Force Model: None
 Continuous Hydrodynamic Maxwellian Equilibrium $f (\rho, \left( u_{0}, \ u_{1}\right), \left( v_{0}, \ v_{1}\right)) = \frac{3 \rho e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$ Compressible: ✓ Deviation Only: ✗ Order: 2
Relaxation Info
Cumulant Eq. Value Relaxation Rate
$1$ $\rho \log{\left(\rho \right)}$ $0.0$
$x$ $\rho u_{0}$ $0.0$
$y$ $\rho u_{1}$ $0.0$
$x^{2}$ $\frac{\rho}{3}$ $\omega_{v}$
$y^{2}$ $\frac{\rho}{3}$ $\omega_{v}$
$x y$ $0$ $\omega_{v}$
$x^{2} y$ $0$ $1.0$
$x y^{2}$ $0$ $1.0$
$x^{2} y^{2}$ $0$ $1.0$

### Relaxation of Polynomial Cumulants¶

By setting method="cumulant", a set of default polynomial cumulants is chosen to be relaxed. Those cumulants are taken from literature and assembled into groups selected to enforce rotational invariance (see: lbmpy.methods.centeredcumulant.centered_cumulants). Default polynomial groups are available for the D2Q9, D3Q19 and D3Q27 stencils. As before it is possible to specify only a single relaxation rate assigned to the moments governing the shear viscosity. All other relaxation rates are then automatically set to one.

[3]:
lbm_config = LBMConfig(method=Method.CUMULANT, relaxation_rate=sp.Symbol('omega_v'), compressible=True)
method_polynomial = create_lb_method(lbm_config=lbm_config)
method_polynomial
[3]:
 Cumulant-Based Method Stencil: D2Q9 Zero-Centered Storage: ✓ Force Model: None
 Continuous Hydrodynamic Maxwellian Equilibrium $f (\rho, \left( u_{0}, \ u_{1}\right), \left( v_{0}, \ v_{1}\right)) = \frac{3 \rho e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$ Compressible: ✓ Deviation Only: ✗ Order: 2
Relaxation Info
Cumulant Eq. Value Relaxation Rate
$1$ $\rho \log{\left(\rho \right)}$ $0.0$
$x$ $\rho u_{0}$ $0.0$
$y$ $\rho u_{1}$ $0.0$
$x y$ $0$ $\omega_{v}$
$x^{2} - y^{2}$ $0$ $\omega_{v}$
$x^{2} + y^{2}$ $\frac{2 \rho}{3}$ $1.0$
$x^{2} y$ $0$ $1.0$
$x y^{2}$ $0$ $1.0$
$x^{2} y^{2}$ $0$ $1.0$

## C) Exemplary simulation: flow around sphere¶

To end this tutorial, we show an example simulation of a channel flow with a spherical obstacle. This example is shown in 2D with the D2Q9 stencil but can be adapted easily for a 3D simulation.

[4]:
from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments

To define the channel flow with dimensionless parameters, we first define the reference length in lattice cells. The reference length will be the diameter of the spherical obstacle. Furthermore, we define a maximal velocity which will be set for the inflow later. The Reynolds number is set relatively high to 100000, which will cause a turbulent flow in the channel.

From these definitions, we can calculate the relaxation rate omega for the lattice Boltzmann method.

[5]:
reference_length = 30
maximal_velocity = 0.05
reynolds_number = 100000
kinematic_vicosity = (reference_length * maximal_velocity) / reynolds_number
initial_velocity=(maximal_velocity, 0)

omega = relaxation_rate_from_lattice_viscosity(kinematic_vicosity)

As a next step, we define the domain size of our set up.

[6]:
stencil = LBStencil(Stencil.D2Q9)
domain_size = (reference_length * 12, reference_length * 4)
dim = len(domain_size)

Now the data for the simulation is allocated. We allocate a field src for the PDFs and field dst used as a temporary field to implement the two grid pull pattern. Additionally, we allocate a velocity field velField

[7]:
dh = ps.create_data_handling(domain_size=domain_size, periodicity=(False, False))

dh.fill('src', 0.0, ghost_layers=True)
dh.fill('dst', 0.0, ghost_layers=True)

dh.fill('velField', 0.0, ghost_layers=True)

We choose a cumulant lattice Boltzmann method, as described above. Here the second-order cumulants are relaxed with the relaxation rate calculated above. All higher-order cumulants are relaxed with one, which means we set them to the equilibrium.

[8]:
lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CUMULANT, relaxation_rate=omega,
compressible=True,
output={'velocity': velField}, kernel_type='stream_pull_collide')

method = create_lb_method(lbm_config=lbm_config)
method
[8]:
 Cumulant-Based Method Stencil: D2Q9 Zero-Centered Storage: ✓ Force Model: None
 Continuous Hydrodynamic Maxwellian Equilibrium $f (\rho, \left( u_{0}, \ u_{1}\right), \left( v_{0}, \ v_{1}\right)) = \frac{3 \rho e^{- \frac{3 \left(- u_{0} + v_{0}\right)^{2}}{2} - \frac{3 \left(- u_{1} + v_{1}\right)^{2}}{2}}}{2 \pi}$ Compressible: ✓ Deviation Only: ✗ Order: 2
Relaxation Info
Cumulant Eq. Value Relaxation Rate
$1$ $\rho \log{\left(\rho \right)}$ $0.0$
$x$ $\rho u_{0}$ $0.0$
$y$ $\rho u_{1}$ $0.0$
$x y$ $0$ $1.9998200161985422$
$x^{2} - y^{2}$ $0$ $1.9998200161985422$
$x^{2} + y^{2}$ $\frac{2 \rho}{3}$ $1.0$
$x^{2} y$ $0$ $1.0$
$x y^{2}$ $0$ $1.0$
$x^{2} y^{2}$ $0$ $1.0$

### Initialization with equilibrium¶

[9]:
init = pdf_initialization_assignments(method, 1.0, initial_velocity, src.center_vector)

ast_init = ps.create_kernel(init, target=dh.default_target)
kernel_init = ast_init.compile()

dh.run_kernel(kernel_init)
---------------------------------------------------------------------------
KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'zero'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'zero'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'zero'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'extended_negative'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'extended_real'

During handling of the above exception, another exception occurred:

KeyError                                  Traceback (most recent call last)
~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
453         try:
--> 454             return self._assumptions[fact]
455         except KeyError:

KeyError: 'finite'

During handling of the above exception, another exception occurred:

KeyboardInterrupt                         Traceback (most recent call last)
<ipython-input-9-75e28407fdf7> in <module>
----> 1 init = pdf_initialization_assignments(method, 1.0, initial_velocity, src.center_vector)
2
3 ast_init = ps.create_kernel(init, target=dh.default_target)
4 kernel_init = ast_init.compile()
5

~/LSS/lssgit/python_includes/lbmpy/macroscopic_value_kernels.py in pdf_initialization_assignments(lb_method, density, velocity, pdfs, streaming_pattern, previous_timestep, set_pre_collision_pdfs)
33     cqc = lb_method.conserved_quantity_computation
34     inp_eqs = cqc.equilibrium_input_equations_from_init_values(density, velocity, force_substitution=False)
---> 35     setter_eqs = lb_method.get_equilibrium(conserved_quantity_equations=inp_eqs)
36     setter_eqs = setter_eqs.new_with_substitutions({sym: field_accesses[i]
37                                                     for i, sym in enumerate(lb_method.post_collision_pdf_symbols)})

~/LSS/lssgit/python_includes/lbmpy/methods/cumulantbased/cumulantbasedmethod.py in get_equilibrium(self, conserved_quantity_equations, subexpressions, pre_simplification, keep_cqc_subexpressions, include_force_terms)
229         r_info_dict = {c: RelaxationInfo(info.equilibrium_value, sp.Integer(1))
230                        for c, info in self.relaxation_info_dict.items()}
--> 231         ac = self._centered_cumulant_collision_rule(
232             r_info_dict, conserved_quantity_equations, pre_simplification,
233             include_force_terms=include_force_terms, symbolic_relaxation_rates=False)

~/LSS/lssgit/python_includes/lbmpy/methods/cumulantbased/cumulantbasedmethod.py in _centered_cumulant_collision_rule(self, cumulant_to_relaxation_info_dict, conserved_quantity_equations, pre_simplification, include_force_terms, symbolic_relaxation_rates)
334         k_to_c_transform = self._cumulant_transform_class(stencil, polynomial_cumulants, density, velocity)
335         k_to_c_eqs = k_to_c_transform.forward_transform(simplification=pre_simplification)
--> 336         c_post_to_k_post_eqs = k_to_c_transform.backward_transform(simplification=pre_simplification)
337
338         C_pre = k_to_c_transform.pre_collision_symbols

~/LSS/lssgit/python_includes/lbmpy/moment_transforms/cumulanttransforms.py in backward_transform(self, central_moment_base, simplification, subexpression_base, start_from_monomials)
143         main_assignments = []
144         for exp in self.central_moment_exponents:
--> 145             eq = self.central_moment_from_cumulants(exp, self.mono_base_post)
146             k_symbol = statistical_quantity_symbol(central_moment_base, exp)
147             main_assignments.append(Assignment(k_symbol, eq))

~/LSS/lssgit/python_includes/lbmpy/moment_transforms/cumulanttransforms.py in central_moment_from_cumulants(self, moment_exponents, cumulant_symbol_base)
208
209         c_indices = [(0,) * dim] + [moment_index_from_derivative(d, wave_numbers) for d in derivatives]
--> 210         moment = moment.subs([(statistical_quantity_symbol('c', idx),
211                                statistical_quantity_symbol(cumulant_symbol_base, idx) / rho)
212                               for idx in c_indices])

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/basic.py in subs(self, *args, **kwargs)
976             rv = self
977             for old, new in sequence:
--> 978                 rv = rv._subs(old, new, **kwargs)
979                 if not isinstance(rv, Basic):
980                     break

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/basic.py in _subs(self, old, new, **hints)
1090         rv = self._eval_subs(old, new)
1091         if rv is None:
-> 1092             rv = fallback(self, old, new)
1093         return rv
1094

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/basic.py in fallback(self, old, new)
1067                     args[i] = arg
1068             if hit:
-> 1069                 rv = self.func(*args)
1070                 hack2 = hints.get('hack2', False)
1071                 if hack2 and self.is_Mul and not rv.is_Mul:  # 2-arg hack

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/operations.py in __new__(cls, *args, **options)
50             return args[0]
51
---> 52         c_part, nc_part, order_symbols = cls.flatten(args)
53         is_commutative = not nc_part
54         obj = cls._from_args(c_part + nc_part, is_commutative)

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/mul.py in flatten(cls, seq)
195                 seq = [a, b]
196             assert not a is S.One
--> 197             if not a.is_zero and a.is_Rational:
198                 r, b = b.as_coeff_Mul()

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

499         pass
500     else:
--> 501         a = evaluate(obj)
502         if a is not None:
503             assumptions.deduce_all_facts(((fact, a),))

649         pos = nonneg = nonpos = unknown_sign = False
650         saw_INF = set()
--> 651         args = [a for a in self.args if not a.is_zero]
652         if not args:
653             return False

649         pos = nonneg = nonpos = unknown_sign = False
650         saw_INF = set()
--> 651         args = [a for a in self.args if not a.is_zero]
652         if not args:
653             return False

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

499         pass
500     else:
--> 501         a = evaluate(obj)
502         if a is not None:
503             assumptions.deduce_all_facts(((fact, a),))

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/power.py in _eval_is_algebraic(self)
1307                 return False
1308
-> 1309         if self.base.is_zero or _is_one(self.base):
1310             return True
1311         elif self.exp.is_rational:

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/power.py in _is_one(expr)
1302         def _is_one(expr):
1303             try:
-> 1304                 return (expr - 1).is_zero
1305             except ValueError:
1306                 # when the operation is not allowed

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

499         pass
500     else:
--> 501         a = evaluate(obj)
502         if a is not None:
503             assumptions.deduce_all_facts(((fact, a),))

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/expr.py in _eval_is_negative(self)
869         if finite is False:
870             return False
--> 871         extended_negative = self.is_extended_negative
872         if finite is True:
873             return extended_negative

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

499         pass
500     else:
--> 501         a = evaluate(obj)
502         if a is not None:
503             assumptions.deduce_all_facts(((fact, a),))

723         if not c.is_zero:
--> 724             v = _monotonic_sign(a)
725             if v is not None:
726                 s = v + c

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/exprtools.py in _monotonic_sign(self)
64     >>> F(nn - 1)  # could be negative, zero or positive
65     """
---> 66     if not self.is_extended_real:
67         return
68

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

511             continue
512         if pk in handler_map:
514
515             # we might have found the value of fact

499         pass
500     else:
--> 501         a = evaluate(obj)
502         if a is not None:
503             assumptions.deduce_all_facts(((fact, a),))

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/expr.py in _eval_is_positive(self)
856
857     def _eval_is_positive(self):
--> 858         finite = self.is_finite
859         if finite is False:
860             return False

~/anaconda3/envs/walberla_dev/lib/python3.8/site-packages/sympy/core/assumptions.py in getit(self)
456             if self._assumptions is self.default_assumptions:
457                 self._assumptions = self.default_assumptions.copy()
459
460     getit.func_name = as_property(fact)

506     # Try assumption's prerequisites
507     prereq = list(_assume_rules.prereq[fact])
--> 508     shuffle(prereq)
509     for pk in prereq:
510         if pk in assumptions:

~/anaconda3/envs/walberla_dev/lib/python3.8/random.py in shuffle(self, x, random)
304             for i in reversed(range(1, len(x))):
305                 # pick an element in x[:i+1] with which to exchange x[i]
--> 306                 j = randbelow(i+1)
307                 x[i], x[j] = x[j], x[i]
308         else:

~/anaconda3/envs/walberla_dev/lib/python3.8/random.py in _randbelow_with_getrandbits(self, n)
253         getrandbits = self.getrandbits
254         k = n.bit_length()  # don't use (n-1) here because n can be 1
--> 255         r = getrandbits(k)          # 0 <= r < 2**k
256         while r >= n:
257             r = getrandbits(k)

KeyboardInterrupt:

### Definition of the update rule¶

[ ]:
lbm_optimisation = LBMOptimisation(symbolic_field=src, symbolic_temporary_field=dst)
update = create_lb_update_rule(lb_method=method,
lbm_config=lbm_config,
lbm_optimisation=lbm_optimisation)

ast_kernel = ps.create_kernel(update, target=dh.default_target, cpu_openmp=True)
kernel = ast_kernel.compile()

### Definition of the boundary set up¶

[ ]:
def set_sphere(x, y, *_):
mid = (domain_size[0] // 3, domain_size[1] // 2)
return (x-mid[0])**2 + (y-mid[1])**2 < radius**2
[ ]:
bh = LatticeBoltzmannBoundaryHandling(method, dh, 'src', name="bh")

inflow = UBB(initial_velocity)
outflow = ExtrapolationOutflow(stencil[4], method)
wall = NoSlip("wall")

bh.set_boundary(inflow, slice_from_direction('W', dim))
bh.set_boundary(outflow, slice_from_direction('E', dim))
for direction in ('N', 'S'):
bh.set_boundary(wall, slice_from_direction(direction, dim))

plt.figure(dpi=200)
plt.boundary_handling(bh)
[ ]:
def timeloop(timeSteps):
for i in range(timeSteps):
bh()
dh.run_kernel(kernel)
dh.swap("src", "dst")

### Run the simulation¶

[ ]:
mask = np.fromfunction(set_sphere, (domain_size[0], domain_size[1], len(domain_size)))
if 'is_test_run' not in globals():
timeloop(50000)  # initial steps

def run():
timeloop(100)