doc/notebooks/06_tutorial_modifying_method_smagorinsky.ipynb

 %% Cell type:code id: tags:
 ``` python
 from lbmpy.session import *
 from lbmpy.relaxationrates import *
 ```
 %% Cell type:markdown id: tags:
 # Tutorial 06: Modifying a LBM method: Smagorinsky model
 In this demo, we show how to modify a lattice Boltzmann method. As example we are going to add a simple turbulence model, by introducing a rule that locally computes the relaxation parameter dependent on the local strain rate tensor. The Smagorinsky model is implemented directly in *lbmpy* as well, however here we take the manual approach to demonstrate how a LB method can be changed in *lbmpy*.
 ## 1) Theoretical background
 Since we have *sympy* available, we want to start out with the basic model equations and derive the concrete equations ourselves. This approach is less error prone, since the calculations are done by the computer algebra system, and oftentimes this approach is also more general and easier to understand.
 ### a) Smagorinsky model
 The basic idea of the Smagorinsky turbulence model is to safe compute time, by not resolving the smallest eddies of the flow on the grid, but model them by an artifical dissipation term.
 The energy dissipation of small scale vortices is taken into account by introducing a "turbulent viscosity". This additional viscosity depends on local flow properties, namely the local shear rates. The larger the local shear rates the higher the turbulent viscosity and the more artifical dissipation is added.
 The total viscosity is
 $$\nu_{total} = \nu_0 + \underbrace{(C_S \Delta)^2 |S|}_{\nu_{t}}$$
 where $\nu_0$ is the normal viscosity, $C_S$ is the Smagorinsky constant, not to be confused with the speed of sound! Typical values of the Smagorinsky constant are between 0.1 - 0.2. The filter length $\Delta$ is chosen as 1 in lattice coordinates.
 The quantity $|S|$ is computed from the local strain rate tensor $S$ that is given by
 $$S_{ij} = \frac{1}{2} \left( \partial_i u_j + \partial_j u_i \right)$$
 and
 $$|S| = \sqrt{2 S_{ij} S_{ij}}$$
 ### b) LBM implementation of Smagorinsky model
 To add the Smagorinsky model to a LB scheme one has to first compute the strain rate tensor $S_{ij}$ in each cell, and compute the turbulent viscosity $\nu_t$ from it. Then the local relaxation rate has to be adapted to match the total viscosity $\nu_{total}$ instead of the standard viscosity $\nu_0$.
 A fortunate property of LB methods is, that the strain rate tensor can be computed locally from the non-equilibrium part of the distribution function. This is somewhat surprising, since the strain rate tensor contains first order derivatives. The strain rate tensor can be obtained by
 $$S_{ij} = - \frac{3 \omega_s}{2 \rho_{(0)}} \Pi_{ij}^{(neq)}$$
 where $\omega_s$ is the relaxation rate that determines the viscosity, $\rho_{(0)}$ is $\rho$ in compressible models and $1$ for incompressible schemes.
 $\Pi_{ij}^{(neq)}$ is the second order moment tensor of the non-equilibrium part of the distribution functions $f^{(neq)} = f - f^{(eq)}$ and can be computed as
 $$\Pi_{ij}^{(neq)} = \sum_q c_{qi} c_{qj} \; f_q^{(neq)}$$
 %% Cell type:markdown id: tags:
 We first have to find a closed form for $S_{ij}$ since in the formula above, it depends on $\omega$, which should be adapated according to $S_{ij}$.
 So we compute $\omega$ and insert it into the formula for $S$:
 %% Cell type:code id: tags:
 ``` python
 τ_0, ρ, ω, ω_total, ω_0 = sp.symbols("tau_0 rho omega omega_total omega_0", positive=True, real=True)
 ν_0, C_S, S, Π = sp.symbols("nu_0, C_S, |S|, Pi", positive=True, real=True)
 Seq = sp.Eq(S, 3 * ω / 2 * Π)
 Seq
 ```
 %% Output
     $\displaystyle |S| = \frac{3 \Pi \omega}{2}$
           3⋅Π⋅ω
     |S| = ─────
             2
 %% Cell type:markdown id: tags:
 Note that we left of the minus, since we took the absolute value of both tensor. The absolute value is defined as above, with the factor of two inside the square root. The $\rho_{(0)}$ has been left out, remembering that $\Pi^{(neq)}$ has to be divided by $\rho$ in case of compressible models|.
 Next, we compute $\omega$ from the total viscosity as given by the Smagorinsky equation:
 %% Cell type:code id: tags:
 ``` python
 relaxation_rate_from_lattice_viscosity(ν_0 + C_S ** 2 * S)
 ```
 %% Output
     $\displaystyle \frac{2}{6 C_{S}^{2} |S| + 6 \nu_{0} + 1}$
               2
     ─────────────────────
          2
     6⋅C_S ⋅|S| + 6⋅ν₀ + 1
 %% Cell type:markdown id: tags:
 and insert it into the equation for $|S|$
 %% Cell type:code id: tags:
 ``` python
 Seq2 = Seq.subs(ω, relaxation_rate_from_lattice_viscosity(ν_0 + C_S **2 * S ))
 Seq2
 ```
 %% Output
     $\displaystyle |S| = \frac{3 \Pi}{6 C_{S}^{2} |S| + 6 \nu_{0} + 1}$
                    3⋅Π
     |S| = ─────────────────────
                2
           6⋅C_S ⋅|S| + 6⋅ν₀ + 1
 %% Cell type:markdown id: tags:
 This equation contains only known quantities, such that we can solve it for $|S|$.
 Additionally we substitute the lattice viscosity $\nu_0$ by the original relaxation time $\tau_0$. The resulting equations get simpler using relaxation times instead of rates.
 %% Cell type:code id: tags:
 ``` python
 solveRes = sp.solve(Seq2, S)
 assert len(solveRes) == 1
 SVal = solveRes[0]
 SVal = SVal.subs(ν_0, lattice_viscosity_from_relaxation_rate(1 / τ_0)).expand()
 SVal
 ```
 %% Output
     $\displaystyle - \frac{\tau_{0}}{6 C_{S}^{2}} + \frac{\sqrt{72 C_{S}^{2} \Pi + 4 \tau_{0}^{2}}}{12 C_{S}^{2}}$
                   ___________________
                  ╱       2         2
         τ₀     ╲╱  72⋅C_S ⋅Π + 4⋅τ₀
     - ────── + ──────────────────────
            2                2
       6⋅C_S           12⋅C_S
 %% Cell type:markdown id: tags:
 Knowning $|S|$ we can compute the total relaxation time using
 $$\nu_{total} = \nu_0 +C_S^2 |S|$$
 %% Cell type:code id: tags:
 ``` python
 τ_val = 1 / (relaxation_rate_from_lattice_viscosity(lattice_viscosity_from_relaxation_rate(1/τ_0) + C_S**2 * SVal)).cancel()
 τ_val
 ```
 %% Output
     $\displaystyle \frac{\tau_{0}}{2} + \frac{\sqrt{18 C_{S}^{2} \Pi + \tau_{0}^{2}}}{2}$
             _________________
            ╱       2       2
     τ₀   ╲╱  18⋅C_S ⋅Π + τ₀
     ── + ────────────────────
     2             2
 %% Cell type:markdown id: tags:
 To compute $\Pi^{(neq)}$ we use the following functions:
 %% Cell type:code id: tags:
 ``` python
 def second_order_moment_tensor(function_values, stencil):
     assert len(function_values) == len(stencil)
     dim = len(stencil[0])
     return sp.Matrix(dim, dim, lambda i, j: sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
 def frobenius_norm(matrix, factor=1):
     return sp.sqrt(sum(i*i for i in matrix) * factor)
 ```
 %% Cell type:markdown id: tags:
 In the next cell we construct equations that take an standard relaxation rate $\omega_0$ and compute a new relaxation rate $\omega_{total}$ according to the Smagorinksy model, using `τ_val` computed above
 %% Cell type:code id: tags:
 ``` python
 def smagorinsky_equations(ω_0, ω_total, method):
     f_neq = sp.Matrix(method.pre_collision_pdf_symbols) - method.get_equilibrium_terms()
     return [sp.Eq(τ_0, 1 / ω_0),
             sp.Eq(Π, frobenius_norm(second_order_moment_tensor(f_neq, method.stencil), factor=2)),
             sp.Eq(ω_total, 1 / τ_val)]
 smagorinsky_equations(ω_0, ω_total, create_lb_method())
 ```
 %% Output
     $\displaystyle \left[ \tau_{0} = \frac{1}{\omega_{0}}, \  \Pi = \sqrt{4 \left(- f_{5} + f_{6} + f_{7} - f_{8} - u_{0} u_{1}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{1} + f_{2} + f_{5} + f_{6} + f_{7} + f_{8} - u_{1}^{2}\right)^{2} + 2 \left(- \frac{\delta_{\rho}}{3} + f_{3} + f_{4} + f_{5} + f_{6} + f_{7} + f_{8} - u_{0}^{2}\right)^{2}}, \  \omega_{total} = \frac{1}{\frac{\tau_{0}}{2} + \frac{\sqrt{18 C_{S}^{2} \Pi + \tau_{0}^{2}}}{2}}\right]$
     ⎡                  ___________________________________________________________
     ⎢                 ╱
     ⎢     1          ╱                                2     ⎛  δᵨ
     ⎢τ₀ = ──, Π =   ╱   4⋅(-f₅ + f₆ + f₇ - f₈ - u₀⋅u₁)  + 2⋅⎜- ── + f₁ + f₂ + f₅ +
     ⎢     ω₀      ╲╱                                        ⎝  3
     ⎢
     ⎢
     ⎢
     ⎣
     
     ______________________________________________________________________
                         2                                               2
                       2⎞      ⎛  δᵨ                                   2⎞
      f₆ + f₇ + f₈ - u₁ ⎟  + 2⋅⎜- ── + f₃ + f₄ + f₅ + f₆ + f₇ + f₈ - u₀ ⎟  , ωₜₒₜₐₗ
                        ⎠      ⎝  3                                     ⎠
     
     
     
     
     
                                 ⎤
                                 ⎥
                    1            ⎥
      = ─────────────────────────⎥
                _________________⎥
               ╱       2       2 ⎥
        τ₀   ╲╱  18⋅C_S ⋅Π + τ₀  ⎥
        ── + ────────────────────⎥
        2             2          ⎦
 %% Cell type:markdown id: tags:
 ## 2) Application: Channel flow
 Next we modify a *lbmpy* scenario to use the Smagorinsky model.
 We create a MRT method, where we fix all relaxation rates except the relaxation rate that controls the viscosity.
 %% Cell type:code id: tags:
 ``` python
-lbm_conifg = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),
+lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, force=(1e-6, 0),
                        force_model=ForceModel.LUO, relaxation_rates=[ω, 1.9, 1.9, 1.9])
-method = create_lb_method(lbm_config=lbm_conifg)
+method = create_lb_method(lbm_config=lbm_config)
 method
 ```
 %% Output
     <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f745d5c3b20>
 %% Cell type:markdown id: tags:
 Only the collision rule has to be changed. Thus we first construct the collision rule, add the Smagorinsky equations and create a normal scenario from the modified collision rule. To avoid that the macroscopic quantity symbols in the Smagorinsky equations fall prey to optimization, we must disable simplification:
 %% Cell type:code id: tags:
 ``` python
 optimization = {'simplification' : False}
 collision_rule = create_lb_collision_rule(lb_method=method, optimization=optimization)
 collision_rule = collision_rule.new_with_substitutions({ω: ω_total})
 collision_rule.subexpressions += smagorinsky_equations(ω, ω_total, method)
 collision_rule.topological_sort(sort_subexpressions=True, sort_main_assignments=False)
 collision_rule
 ```
 %% Output
     AssignmentCollection: d_6, d_5, d_2, d_3, d_0, d_1, d_7, d_4, d_8 <- f(f_2, f_4, f_7, f_5, omega_total, f_6, f_3, f_1, f_0, f_8)
 %% Cell type:markdown id: tags:
 In the next cell the collision rule is simplified by extracting common subexpressions
 %% Cell type:code id: tags:
 ``` python
 from pystencils.simp import sympy_cse
 #collision_rule = sympy_cse(collision_rule)
 ```
 %% Cell type:markdown id: tags:
 A channel scenario can be created from a modified collision rule:
 %% Cell type:code id: tags:
 ``` python
 ch = create_channel((300, 100), force=1e-6, collision_rule=collision_rule,
                     kernel_params={"C_S": 0.12, "omega": 1.999})
 ```
 %% Cell type:code id: tags:
 ``` python
 #show_code(ch.ast)
 ```
 %% Cell type:code id: tags:
 ``` python
 ch.run(5000)
 ```
 %% Cell type:code id: tags:
 ``` python
 plt.figure(dpi=200)
 plt.vector_field(ch.velocity[:, :])
 np.max(ch.velocity[:, :])
 ```
 %% Output
     $\displaystyle 0.00504266401703371$
     0.00504266401703371
 %% Cell type:markdown id: tags:
 ## Appendix: Strain rate tensor formula from Chapman Enskog
 The connection between $S_{ij}$ and $\Pi_{ij}^{(neq)}$ can be seen using a Chapman Enskog expansion. Since *lbmpy* has a module that automatically does this expansions we can have a look at it:
 %% Cell type:code id: tags:
 ``` python
 from lbmpy.chapman_enskog import ChapmanEnskogAnalysis, CeMoment
 from lbmpy.chapman_enskog.chapman_enskog import remove_higher_order_u
 compressible_model = create_lb_method(stencil=Stencil.D2Q9, compressible=True, zero_centered=False)
 incompressible_model = create_lb_method(stencil=Stencil.D2Q9, compressible=False, zero_centered=False)
 ce_compressible = ChapmanEnskogAnalysis(compressible_model)
 ce_incompressible = ChapmanEnskogAnalysis(incompressible_model)
 ```
 %% Cell type:markdown id: tags:
 The Chapman Enskog analysis yields expresssions for the moment
 $\Pi = \Pi^{(eq)} + \epsilon \Pi^{(1)} +  \epsilon^2 \Pi^{(2)} \cdots$
 and the strain rate tensor is related to $\Pi^{(1)}$. However the best approximation we have for $\Pi^{(1)}$ is
 $\Pi^{(neq)}$. For details, see the paper "Shear stress in lattice Boltzmann simulations" by Krüger, Varnik and Raabe from 2009.
 Lets look at the values of $\Pi^{(1)}$ obtained from the Chapman enskog expansion:
 %% Cell type:code id: tags:
 ``` python
 Π_1_xy = CeMoment("\\Pi", moment_tuple=(1,1), superscript=1)
 Π_1_xx = CeMoment("\\Pi", moment_tuple=(2,0), superscript=1)
 Π_1_yy = CeMoment("\\Pi", moment_tuple=(0,2), superscript=1)
 components = (Π_1_xx, Π_1_yy, Π_1_xy)
 Π_1_xy_val = ce_compressible.higher_order_moments[Π_1_xy]
 Π_1_xy_val
 ```
 %% Output
     $\displaystyle \frac{\rho u_{0}^{2} {\partial^{(1)}_{0} u_{1}} + 2 \rho u_{0} u_{1} {\partial^{(1)}_{0} u_{0}} + 2 \rho u_{0} u_{1} {\partial^{(1)}_{1} u_{1}} + \rho u_{1}^{2} {\partial^{(1)}_{1} u_{0}} - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3} + u_{0}^{2} u_{1} {\partial^{(1)}_{0} \rho} + u_{0} u_{1}^{2} {\partial^{(1)}_{1} \rho}}{\omega}$
         2                                                    2          ρ⋅D(u_0)
     ρ⋅u₀ ⋅D(u_1) + 2⋅ρ⋅u₀⋅u₁⋅D(u_0) + 2⋅ρ⋅u₀⋅u₁⋅D(u_1) + ρ⋅u₁ ⋅D(u_0) - ──────── -
                                                                            3
     ──────────────────────────────────────────────────────────────────────────────
                                                                ω
     
      ρ⋅D(u_1)     2                  2
      ──────── + u₀ ⋅u₁⋅D(rho) + u₀⋅u₁ ⋅D(rho)
         3
     ─────────────────────────────────────────
     
 %% Cell type:markdown id: tags:
 This term has lots of higher order error terms in it. We assume that $u$ is small in lattice coordinates, so if we neglect all terms in $u$ that are quadratic or higher we get:
 %% Cell type:code id: tags:
 ``` python
 remove_higher_order_u(Π_1_xy_val.expand())
 ```
 %% Output
     $\displaystyle - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3 \omega}$
       ρ⋅D(u_0)   ρ⋅D(u_1)
     - ──────── - ────────
         3⋅ω        3⋅ω
 %% Cell type:markdown id: tags:
 Putting these steps together into a function, we can display them for the different cases quickly:
 %% Cell type:code id: tags:
 ``` python
 def get_Π_1(ce_analysis, component):
     val = ce_analysis.higher_order_moments[component]
     return remove_higher_order_u(val.expand())
 ```
 %% Cell type:markdown id: tags:
 Compressible case:
 %% Cell type:code id: tags:
 ``` python
 tuple(get_Π_1(ce_compressible, Pi) for Pi in components)
 ```
 %% Output
     $\displaystyle \left( - \frac{2 \rho {\partial^{(1)}_{0} u_{0}}}{3 \omega}, \  - \frac{2 \rho {\partial^{(1)}_{1} u_{1}}}{3 \omega}, \  - \frac{\rho {\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{\rho {\partial^{(1)}_{0} u_{1}}}{3 \omega}\right)$
     ⎛-2⋅ρ⋅D(u_0)   -2⋅ρ⋅D(u_1)     ρ⋅D(u_0)   ρ⋅D(u_1)⎞
     ⎜────────────, ────────────, - ──────── - ────────⎟
     ⎝    3⋅ω           3⋅ω           3⋅ω        3⋅ω   ⎠
 %% Cell type:markdown id: tags:
 Incompressible case:
 %% Cell type:code id: tags:
 ``` python
 tuple(get_Π_1(ce_incompressible, Pi) for Pi in components)
 ```
 %% Output
     $\displaystyle \left( \frac{2 u_{0} {\partial^{(1)}_{0} \rho}}{3 \omega} - \frac{2 {\partial^{(1)}_{0} u_{0}}}{3 \omega}, \  \frac{2 u_{1} {\partial^{(1)}_{1} \rho}}{3 \omega} - \frac{2 {\partial^{(1)}_{1} u_{1}}}{3 \omega}, \  \frac{u_{0} {\partial^{(1)}_{1} \rho}}{3 \omega} + \frac{u_{1} {\partial^{(1)}_{0} \rho}}{3 \omega} - \frac{{\partial^{(1)}_{1} u_{0}}}{3 \omega} - \frac{{\partial^{(1)}_{0} u_{1}}}{3 \omega}\right)$
     ⎛2⋅u₀⋅D(rho)   2⋅D(u_0)  2⋅u₁⋅D(rho)   2⋅D(u_1)  u₀⋅D(rho)   u₁⋅D(rho)   D(u_0
     ⎜─────────── - ────────, ─────────── - ────────, ───────── + ───────── - ─────
     ⎝    3⋅ω         3⋅ω         3⋅ω         3⋅ω        3⋅ω         3⋅ω       3⋅ω
     
     )   D(u_1)⎞
     ─ - ──────⎟
          3⋅ω  ⎠
 %% Cell type:markdown id: tags:
 In the incompressible case has some terms $\partial \rho$ which are zero, since $\rho$ is assumed constant.
 Leaving out the error terms we finally obtain:
 $$\Pi_{ij}^{(neq)} \approx \Pi_{ij}^{(1)} = -\frac{2 \rho_{(0)}}{3 \omega_s} \left( \partial_i u_j + \partial_j u_i \right)$$

doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb

 %% Cell type:markdown id: tags:
 
 # The conservative Allen-Cahn model for high Reynolds number, two phase flow with large-density and viscosity constrast
 
 %% Cell type:code id: tags:
 
 ``` python
 from pystencils.session import *
 from lbmpy.session import *
 
-from pystencils.simp import sympy_cse
 from pystencils.boundaries import BoundaryHandling
 
 from lbmpy.phasefield_allen_cahn.contact_angle import ContactAngle
 from lbmpy.phasefield_allen_cahn.kernel_equations import *
 from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti, AllenCahnParameters
 from lbmpy.advanced_streaming import LBMPeriodicityHandling
 from lbmpy.boundaries import NoSlip, LatticeBoltzmannBoundaryHandling
 ```
 
 %% Cell type:markdown id: tags:
 
-If `pycuda` is installed the simulation automatically runs on GPU
+If `cupy` is installed the simulation automatically runs on GPU
 
 %% Cell type:code id: tags:
 
 ``` python
 try:
-    import pycuda
+    import cupy
 except ImportError:
-    pycuda = None
+    cupy = None
     gpu = False
     target = ps.Target.CPU
-    print('No pycuda installed')
+    print('No cupy installed')
 
-if pycuda:
+if cupy:
     gpu = True
     target = ps.Target.GPU
 ```
 
-%% Output
-
-    No pycuda installed
-
 %% Cell type:markdown id: tags:
 
 The conservative Allen-Cahn model (CACM) for two-phase flow is based on the work of Fakhari et al. (2017) [Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). The model can be created for two-dimensional problems as well as three-dimensional problems, which have been described by Mitchell et al. (2018) [Development of a three-dimensional
 phase-field lattice Boltzmann method for the study of immiscible fluids at high density ratios](http://dx.doi.org/10.1103/PhysRevE.96.053301). Furthermore, cascaded lattice Boltzmann methods can be combined with the model which was described in [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018)
 
 
 The CACM is suitable for simulating highly complex two phase flow problems with high-density ratios and high Reynolds numbers. In this tutorial, an overview is provided on how to derive the model with lbmpy. For this, the model is defined with two LBM populations. One for the interface tracking, which we call the phase-field LB step and one for recovering the hydrodynamic properties. The latter is called the hydrodynamic LB step.
 
 %% Cell type:markdown id: tags:
 
 ## Geometry Setup
 
 First of all, the stencils for the phase-field LB step as well as the stencil for the hydrodynamic LB step are defined. According to the stencils, the simulation can be performed in either 2D- or 3D-space. For 2D simulations, only the D2Q9 stencil is supported. For 3D simulations, the D3Q15, D3Q19 and the D3Q27 stencil are supported. Note here that the cascaded LBM can not be derived for D3Q15 stencils.
 
 %% Cell type:code id: tags:
 
 ``` python
 stencil_phase = LBStencil(Stencil.D2Q9)
 stencil_hydro = LBStencil(Stencil.D2Q9)
 assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
 
 dimensions = len(stencil_phase[0])
 ```
 
 %% Cell type:markdown id: tags:
 
 Definition of the domain size
 
 %% Cell type:code id: tags:
 
 ``` python
 # domain
 L0 = 256
 domain_size = (L0, 4 * L0)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Parameter definition
 
 The next step is to calculate all parameters which are needed for the simulation. In this example, a Rayleigh-Taylor instability test case is set up. The parameter calculation for this setup is already implemented in lbmpy and can be used with the dimensionless parameters which describe the problem.
 
 %% Cell type:code id: tags:
 
 ``` python
 # time step
 timesteps = 8000
 
 # reference time
 reference_time = 4000
 
 # calculate the parameters for the RTI
 parameters = calculate_parameters_rti(reference_length=L0,
                                       reference_time=reference_time,
                                       density_heavy=1.0,
                                       capillary_number=0.44,
                                       reynolds_number=3000,
                                       atwood_number=0.998,
                                       peclet_number=1000,
                                       density_ratio=1000,
                                       viscosity_ratio=100)
 ```
 
 %% Cell type:markdown id: tags:
 
 This function returns a `AllenCahnParameters` class. It is struct like class holding all parameters for the conservative Allen Cahn model:
 
 %% Cell type:code id: tags:
 
 ``` python
 parameters
 ```
 
 %% Output
 
-    <lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x126d30cd0>
+    <lbmpy.phasefield_allen_cahn.parameter_calculation.AllenCahnParameters at 0x1329b88b0>
 
 %% Cell type:markdown id: tags:
 
 ## Fields
 
 As a next step all fields which are needed get defined. To do so, we create a `datahandling` object. More details about it can be found in the third tutorial of the [pystencils framework]( http://pycodegen.pages.walberla.net/pystencils/). This object holds all fields and manages the kernel runs.
 
 %% Cell type:code id: tags:
 
 ``` python
 # create a datahandling object
 dh = ps.create_data_handling((domain_size), periodicity=(True, False), parallel=False, default_target=target)
 
 # pdf fields. g is used for hydrodynamics and h for the interface tracking
 g = dh.add_array("g", values_per_cell=len(stencil_hydro))
 dh.fill("g", 0.0, ghost_layers=True)
 h = dh.add_array("h",values_per_cell=len(stencil_phase))
 dh.fill("h", 0.0, ghost_layers=True)
 
 g_tmp = dh.add_array("g_tmp", values_per_cell=len(stencil_hydro))
 dh.fill("g_tmp", 0.0, ghost_layers=True)
 h_tmp = dh.add_array("h_tmp",values_per_cell=len(stencil_phase))
 dh.fill("h_tmp", 0.0, ghost_layers=True)
 
 # velocity field
 u = dh.add_array("u", values_per_cell=dh.dim)
 dh.fill("u", 0.0, ghost_layers=True)
 
 # phase-field
 C = dh.add_array("C")
 dh.fill("C", 0.0, ghost_layers=True)
 C_tmp = dh.add_array("C_tmp")
 dh.fill("C_tmp", 0.0, ghost_layers=True)
 ```
 
 %% Cell type:markdown id: tags:
 
 As a next step the relaxation time is stated in a symbolic form. It is calculated via interpolation.
 
 %% Cell type:code id: tags:
 
 ``` python
 rho_L = parameters.symbolic_density_light
 rho_H = parameters.symbolic_density_heavy
 density = rho_L + C.center * (rho_H - rho_L)
 
 body_force = [0, 0, 0]
 body_force[1] = parameters.symbolic_gravitational_acceleration * density
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the lattice Boltzmann methods
 
 %% Cell type:markdown id: tags:
 
 For both LB steps, a weighted orthogonal MRT (WMRT) method is used. It is also possible to change the method to a simpler SRT scheme or a more complicated CLBM scheme. The CLBM scheme can be obtained with `Method.CENTRAL_MOMENT`. Note here that the hydrodynamic LB step is formulated as an incompressible velocity-based LBM. Thus, the velocity terms can not be removed from the equilibrium in the central moment space.
 
 %% Cell type:code id: tags:
 
 ``` python
 w_c = parameters.symbolic_omega_phi
 
 config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
                          delta_equilibrium=False,
                          force=sp.symbols("F_:2"), velocity_input=u,
                          weighted=True, relaxation_rates=[0, w_c, w_c, 1, 1, 1, 1, 1, 1],
-                         output={'density': C_tmp}, kernel_type='stream_pull_collide')
+                         output={'density': C_tmp})
 
 method_phase = create_lb_method(lbm_config=config_phase)
 
 method_phase
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x126d44ee0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x1346b22e0>
 
 %% Cell type:code id: tags:
 
 ``` python
 omega = parameters.omega(C)
 config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
                          weighted=True, relaxation_rates=[omega, 1, 1, 1],
-                         force=sp.symbols("F_:2"),
-                         output={'velocity': u}, kernel_type='collide_stream_push')
+                         force=sp.symbols("F_:2"), output={'velocity': u})
 
 method_hydro = create_lb_method(lbm_config=config_hydro)
 
 method_hydro
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x126d22eb0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x137772a30>
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:markdown id: tags:
 
 The probability distribution functions (pdfs) are initialised with the equilibrium distribution for the LB methods.
 
 %% Cell type:code id: tags:
 
 ``` python
 h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
 g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
 
 h_init = ps.create_kernel(h_updates, target=dh.default_target, cpu_openmp=True).compile()
 g_init = ps.create_kernel(g_updates, target=dh.default_target, cpu_openmp=True).compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 Following this, the phase field is initialised directly in python.
 
 %% Cell type:code id: tags:
 
 ``` python
 # initialize the domain
 def Initialize_distributions():
     Nx = domain_size[0]
     Ny = domain_size[1]
 
     for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
         x = np.zeros_like(block.midpoint_arrays[0])
         x[:, :] = block.midpoint_arrays[0]
 
         y = np.zeros_like(block.midpoint_arrays[1])
         y[:, :] = block.midpoint_arrays[1]
 
         y -= 2 * L0
         tmp = 0.1 * Nx * np.cos((2 * np.pi * x) / Nx)
         init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (parameters.interface_thickness / 2))
         block["C"][:, :] = init_values
         block["C_tmp"][:, :] = init_values
 
     if gpu:
         dh.all_to_gpu()
 
     dh.run_kernel(h_init, **parameters.symbolic_to_numeric_map)
     dh.run_kernel(g_init)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 Initialize_distributions()
 plt.scalar_field(dh.gather_array(C.name))
 ```
 
 %% Output
 
-    <matplotlib.image.AxesImage at 0x1417f86d0>
+    <matplotlib.image.AxesImage at 0x137ad6bb0>
 
-
+
 
 %% Cell type:markdown id: tags:
 
 ## Source Terms
 
 %% Cell type:markdown id: tags:
 
 For the Allen-Cahn LB step, the Allen-Cahn equation needs to be applied as a source term. Here, a simple forcing model is used which is directly applied in the moment space:
 $$
 F_i^\phi (\boldsymbol{x}, t) = \Delta t \frac{\left[1 - 4 \left(\phi - \phi_0\right)^2\right]}{\xi} w_i \boldsymbol{c}_i \cdot \frac{\nabla \phi}{|{\nabla \phi}|},
 $$
 where $\phi$ is the phase-field, $\phi_0$ is the interface location, $\Delta t$ it the timestep size $\xi$ is the interface width, $\boldsymbol{c}_i$ is the discrete direction from stencil_phase and $w_i$ are the weights. Furthermore, the equilibrium needs to be shifted:
 
 $$
 \bar{h}^{eq}_\alpha = h^{eq}_\alpha - \frac{1}{2} F^\phi_\alpha
 $$
 
 The hydrodynamic force is given by:
 $$
 F_i (\boldsymbol{x}, t) = \Delta t w_i \frac{\boldsymbol{c}_i \boldsymbol{F}}{\rho c_s^2},
 $$
 where $\rho$ is the interpolated density and $\boldsymbol{F}$ is the source term which consists of the pressure force
 $$
 \boldsymbol{F}_p = -p^* c_s^2 \nabla \rho,
 $$
 the surface tension force:
 $$
 \boldsymbol{F}_s = \mu_\phi \nabla \phi
 $$
 and the viscous force term:
 $$
 F_{\mu, i}^{\mathrm{MRT}} = - \frac{\nu}{c_s^2 \Delta t} \left[\sum_{\beta} c_{\beta i} c_{\beta j} \times \sum_{\alpha} \Omega_{\beta \alpha}(g_\alpha - g_\alpha^{\mathrm{eq}})\right] \frac{\partial \rho}{\partial x_j}.
 $$
 
 In the above equations $p^*$ is the normalised pressure which can be obtained from the zeroth order moment of the hydrodynamic distribution function $g$. The lattice speed of sound is given with $c_s$ and the chemical potential is $\mu_\phi$. Furthermore, the viscosity is $\nu$ and $\Omega$ is the moment-based collision operator. Note here that the hydrodynamic equilibrium is also adjusted as shown above for the phase-field distribution functions.
 
 
 For CLBM methods the forcing is applied directly in the central moment space. This is done with the `CentralMomentMultiphaseForceModel`. Furthermore, the GUO force model is applied here to be consistent with [A cascaded phase-field lattice Boltzmann model for the simulation of incompressible, immiscible fluids with high density contrast](http://dx.doi.org/10.1016/j.camwa.2019.08.018). Here we refer to equation D.7 which can be derived for 3D stencils automatically with lbmpy.
 
 %% Cell type:code id: tags:
 
 ``` python
 force_h = interface_tracking_force(C, stencil_phase, parameters)
-hydro_force = hydrodynamic_force(g, C, method_hydro, parameters, body_force)
+hydro_force = hydrodynamic_force(C, method_hydro, parameters, body_force)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Definition of the LB update rules
 
 %% Cell type:markdown id: tags:
 
 The update rule for the phase-field LB step is defined as:
 
 $$
 h_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = h_i(\boldsymbol{x}, t) + \Omega_{ij}^h(\bar{h_j}^{eq} - h_j)|_{(\boldsymbol{x}, t)} + F_i^\phi(\boldsymbol{x}, t).
 $$
 In our framework the pull scheme is applied as streaming step. Furthermore, the update of the phase-field is directly integrated into the kernel. As a result of this, a second temporary phase-field is needed.
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_optimisation = LBMOptimisation(symbolic_field=h, symbolic_temporary_field=h_tmp)
 allen_cahn_update_rule = create_lb_update_rule(lbm_config=config_phase,
                                                lbm_optimisation=lbm_optimisation)
 
 allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
 
 ast_kernel = ps.create_kernel(allen_cahn_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_allen_cahn_lb = ast_kernel.compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 The update rule for the hydrodynmaic LB step is defined as:
 
 $$
 g_i (\boldsymbol{x} + \boldsymbol{c}_i \Delta t, t + \Delta t) = g_i(\boldsymbol{x}, t) + \Omega_{ij}^g(\bar{g_j}^{eq} - g_j)|_{(\boldsymbol{x}, t)} + F_i(\boldsymbol{x}, t).
 $$
 
 Here, the push scheme is applied which is easier due to the data access required for the viscous force term. Furthermore, the velocity update is directly done in the kernel.
 
 %% Cell type:code id: tags:
 
 ``` python
-force_Assignments = hydrodynamic_force_assignments(g, u, C, method_hydro, parameters, body_force)
+force_Assignments = hydrodynamic_force_assignments(u, C, method_hydro, parameters, body_force)
 
 lbm_optimisation = LBMOptimisation(symbolic_field=g, symbolic_temporary_field=g_tmp)
 hydro_lb_update_rule = create_lb_update_rule(lbm_config=config_hydro,
                                              lbm_optimisation=lbm_optimisation)
 
-hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
+hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters,
+                                              config_hydro)
 
 ast_kernel = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
 kernel_hydro_lb = ast_kernel.compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Boundary Conditions
 
 %% Cell type:markdown id: tags:
 
 As a last step suitable boundary conditions are applied
 
 %% Cell type:code id: tags:
 
 ``` python
 # periodic Boundarys for g, h and C
 periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
 
 periodic_BC_g = LBMPeriodicityHandling(stencil=stencil_hydro, data_handling=dh, pdf_field_name=g.name,
-                                       streaming_pattern='push')
+                                       streaming_pattern=config_hydro.streaming_pattern)
 periodic_BC_h = LBMPeriodicityHandling(stencil=stencil_phase, data_handling=dh, pdf_field_name=h.name,
-                                       streaming_pattern='pull')
+                                       streaming_pattern=config_phase.streaming_pattern)
 
 # No slip boundary for the phasefield lbm
 bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
                                                  target=dh.default_target, name='boundary_handling_h',
-                                                 streaming_pattern='pull')
+                                                 streaming_pattern=config_phase.streaming_pattern)
 
 # No slip boundary for the velocityfield lbm
 bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
                                             target=dh.default_target, name='boundary_handling_g',
-                                            streaming_pattern='push')
+                                            streaming_pattern=config_hydro.streaming_pattern)
 
 contact_angle = BoundaryHandling(dh, C.name, stencil_hydro, target=dh.default_target)
 contact = ContactAngle(90, parameters.interface_thickness)
 
 wall = NoSlip()
 if dimensions == 2:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0])
     bh_hydro.set_boundary(wall, make_slice[:, -1])
 
     contact_angle.set_boundary(contact, make_slice[:, 0])
     contact_angle.set_boundary(contact, make_slice[:, -1])
 else:
     bh_allen_cahn.set_boundary(wall, make_slice[:, 0, :])
     bh_allen_cahn.set_boundary(wall, make_slice[:, -1, :])
 
     bh_hydro.set_boundary(wall, make_slice[:, 0, :])
     bh_hydro.set_boundary(wall, make_slice[:, -1, :])
 
     contact_angle.set_boundary(contact, make_slice[:, 0, :])
     contact_angle.set_boundary(contact, make_slice[:, -1, :])
 
 
 bh_allen_cahn.prepare()
 bh_hydro.prepare()
 contact_angle.prepare()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Full timestep
 
 %% Cell type:code id: tags:
 
 ``` python
 # definition of the timestep for the immiscible fluids model
 def timeloop():
     # Solve the interface tracking LB step with boundary conditions
     periodic_BC_h()
     bh_allen_cahn()
     dh.run_kernel(kernel_allen_cahn_lb, **parameters.symbolic_to_numeric_map)
-    dh.swap("C", "C_tmp")
+    # Solve the hydro LB step with boundary conditions
+    periodic_BC_g()
+    bh_hydro()
+    dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
 
+    dh.swap("C", "C_tmp")
     # apply the three phase-phase contact angle
     contact_angle()
     # periodic BC of the phase-field
     periodic_BC_C()
 
-    # solve the hydro LB step with boundary conditions
-    dh.run_kernel(kernel_hydro_lb, **parameters.symbolic_to_numeric_map)
-    periodic_BC_g()
-    bh_hydro()
-
-
     # field swaps
     dh.swap("h", "h_tmp")
     dh.swap("g", "g_tmp")
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 Initialize_distributions()
 
 frames = 300
 steps_per_frame = (timesteps//frames) + 1
 
 if 'is_test_run' not in globals():
     def run():
         for i in range(steps_per_frame):
             timeloop()
 
         if gpu:
             dh.to_cpu("C")
         return dh.gather_array(C.name)
 
     animation = plt.scalar_field_animation(run, frames=frames, rescale=True)
     set_display_mode('video')
     res = display_animation(animation)
 else:
     timeloop()
     res = None
 res
 ```
 
 %% Output
 
     <IPython.core.display.HTML object>
 
 %% Cell type:markdown id: tags:
 
 Note that the video is played for 10 seconds while the simulation time is only 2 seconds!

doc/notebooks/demo_shallow_water_lbm.ipynb

@@ -541,7 +541,7 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -555,7 +555,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.2"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,

doc/sphinx/api.rst

@@ -8,13 +8,10 @@ API Reference
    enums.rst
    stencils.rst
    kernelcreation.rst
-   methods.rst
    equilibrium.rst
    moment_transforms.rst
    maxwellian_equilibrium.rst
    continuous_distribution_measures.rst
    moments.rst
    cumulants.rst
-   boundary_conditions.rst
-   forcemodels.rst
    zbibliography.rst

doc/sphinx/lbmpy.bib

@@ -110,6 +110,24 @@ issn = {0898-1221},
 doi = {10.1016/j.camwa.2015.05.001},
 }
 
+@article{geier2017,
+author = {Geier, Martin and Pasquali, Andrea and Sch{\"{o}}nherr, Martin},
+title = {Parametrization of the Cumulant Lattice Boltzmann Method for Fourth Order Accurate Diffusion Part I},
+year = {2017},
+issue_date = {November 2017},
+publisher = {Academic Press Professional, Inc.},
+address = {USA},
+volume = {348},
+number = {C},
+issn = {0021-9991},
+url = {https://doi.org/10.1016/j.jcp.2017.05.040},
+doi = {10.1016/j.jcp.2017.05.040},
+journal = {J. Comput. Phys.},
+month = {nov},
+pages = {862–888},
+numpages = {27}
+}
+
 @Article{Coreixas2019,
   title = {Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations},
   author = {Coreixas, Christophe and Chopard, Bastien and Latt, Jonas},
@@ -248,4 +266,52 @@ journal = {Communications in Computational Physics}
   publisher = {Springer Link},
 }
 
+@article{BouzidiBC,
+    author = {Bouzidi, M’hamed and Firdaouss, Mouaouia and Lallemand, Pierre},
+    title = "{Momentum transfer of a Boltzmann-lattice fluid with boundaries}",
+    journal = {Physics of Fluids},
+    year = {2001},
+    month = {11},
+    doi = {10.1063/1.1399290},
+    url = {https://doi.org/10.1063/1.1399290},
+}
+
+@Article{rozema15,
+  doi       = {10.1063/1.4928700},
+  year      = {2015},
+  month     = {aug},
+  publisher = {AIP Publishing},
+  volume    = {27},
+  number    = {8},
+  author    = {Wybe Rozema and Hyun J. Bae and Parviz Moin and Roel Verstappen},
+  title     = {Minimum-dissipation models for large-eddy simulation},
+  journal   = {Physics of Fluids}
+}
+
+@article{Han2021,
+   doi = {10.1088/1873-7005/ac1782},
+   url = {https://dx.doi.org/10.1088/1873-7005/ac1782}                    ,
+   year = {2021},
+   month = {aug},
+   publisher = {IOP Publishing},
+   volume = {53},
+   number = {4},
+   pages = {045506},
+   author = {Mengtao Han and Ryozo Ooka and Hideki Kikumoto},
+   title = {A wall function approach in lattice Boltzmann method: algorithm and validation using turbulent channel flow},
+   journal = {Fluid Dynamics Research}
+}
+
+@article{Maronga2020,
+ author = {Maronga, Bj{\"o}rn and Knigge, Christoph and Raasch, Siegfried},
+ year = {2020},
+ title = {{An Improved Surface Boundary Condition for Large-Eddy Simulations Based on Monin--Obukhov Similarity Theory: Evaluation and Consequences for Grid Convergence in Neutral and Stable Conditions}},
+ pages = {297--325},
+ volume = {174},
+ number = {2},
+ issn = {0006-8314},
+ journal = {{Boundary-layer meteorology}},
+ doi = {10.1007/s10546-019-00485-w}
+}
+
 @Comment{jabref-meta: databaseType:bibtex;}

doc/sphinx/methods.rst

-*******************************************
-Methods and Method Creation (lbmpy.methods)
-*******************************************
+********************************
+Collision models (lbmpy.methods)
+********************************
 
 This module defines the classes defining all types of lattice Boltzmann methods available in *lbmpy*,
 together with a set of factory functions used to create their instances. The factory functions are
@@ -22,21 +22,6 @@ Abstract LB Method Base Class
 .. autoclass:: lbmpy.methods.AbstractLbMethod
     :members:
 
-
-Conserved Quantity Computation
-==============================
-
-The classes of the conserved quantity computation (CQC) submodule define an LB Method's conserved quantities and
-the equations to compute them. For hydrodynamic methods, :class:`lbmpy.methods.DensityVelocityComputation` is
-the typical choice. For custom methods, however, a custom CQC class might have to be created.
-
-.. autoclass:: lbmpy.methods.AbstractConservedQuantityComputation
-    :members:
-
-.. autoclass:: lbmpy.methods.DensityVelocityComputation
-    :members:
-
-
 .. _methods_rawmomentbased:
 
 Raw Moment-based methods
@@ -163,3 +148,17 @@ Low-Level Creation Functions
 .. autofunction:: lbmpy.methods.creationfunctions.create_with_continuous_maxwellian_equilibrium
 
 .. autofunction:: lbmpy.methods.creationfunctions.create_from_equilibrium
+
+
+Conserved Quantity Computation
+==============================
+
+The classes of the conserved quantity computation (CQC) submodule define an LB Method's conserved quantities and
+the equations to compute them. For hydrodynamic methods, :class:`lbmpy.methods.DensityVelocityComputation` is
+the typical choice. For custom methods, however, a custom CQC class might have to be created.
+
+.. autoclass:: lbmpy.methods.AbstractConservedQuantityComputation
+    :members:
+
+.. autoclass:: lbmpy.methods.DensityVelocityComputation
+    :members:
\ No newline at end of file

doc/sphinx/tutorials.rst

@@ -4,6 +4,10 @@ Tutorials
 All tutorials are automatically created by Jupyter Notebooks.
 You can open the notebooks directly to play around with the code examples.
 
+=================
+Basics
+=================
+
 .. toctree::
     :maxdepth: 1
 
@@ -13,17 +17,68 @@ You can open the notebooks directly to play around with the code examples.
     /notebooks/03_tutorial_lbm_formulation.ipynb
     /notebooks/04_tutorial_cumulant_LBM.ipynb
     /notebooks/05_tutorial_nondimensionalization_and_scaling.ipynb
+
+===================
+Turbulence modeling
+===================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/06_tutorial_modifying_method_smagorinsky.ipynb
+
+=================
+Thermal flows
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/07_tutorial_thermal_lbm.ipynb
+    /notebooks/demo_thermalized_lbm.ipynb
+
+=================
+Multiphase flows
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/08_tutorial_shanchen_twophase.ipynb
     /notebooks/09_tutorial_shanchen_twocomponent.ipynb
     /notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb
+
+========================
+Thermocapillary flows
+========================
+
+.. toctree::
+    :maxdepth: 1
+
+    /notebooks/12_Thermocapillary_flows_heated_channel.ipynb
+    /notebooks/13_Thermocapillary_flows_droplet_motion.ipynb
+
+===================
+Non Newtonian flow
+===================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/11_tutorial_Non_Newtonian_Flow.ipynb
+
+=================
+Diverse
+=================
+
+.. toctree::
+    :maxdepth: 1
+
     /notebooks/demo_stencils.ipynb
     /notebooks/demo_streaming_patterns.ipynb
     /notebooks/demo_create_method_from_scratch.ipynb
     /notebooks/demo_moments_cumulants_and_maxwellian_equilibrium.ipynb
     /notebooks/demo_automatic_chapman_enskog_analysis.ipynb
-    /notebooks/demo_thermalized_lbm.ipynb
+    /notebooks/demo_interpolation_boundary_conditions.ipynb
     /notebooks/demo_shallow_water_lbm.ipynb
     /notebooks/demo_theoretical_background_generic_equilibrium_construction.ipynb

lbmpy/__init__.py

-from .creationfunctions import create_lb_ast, create_lb_collision_rule, create_lb_function,\
-    create_lb_method, create_lb_update_rule, LBMConfig, LBMOptimisation
-from .enums import Stencil, Method, ForceModel, CollisionSpace
-from .lbstep import LatticeBoltzmannStep
-from .macroscopic_value_kernels import pdf_initialization_assignments, macroscopic_values_getter,\
-    compile_macroscopic_values_getter, compile_macroscopic_values_setter, create_advanced_velocity_setter_collision_rule
-from .maxwellian_equilibrium import get_weights
-from .relaxationrates import relaxation_rate_from_lattice_viscosity, lattice_viscosity_from_relaxation_rate,\
-    relaxation_rate_from_magic_number
-from .scenarios import create_lid_driven_cavity, create_fully_periodic_flow
-from .stencils import LBStencil
-
-
-__all__ = ['create_lb_ast', 'create_lb_collision_rule', 'create_lb_function', 'create_lb_method',
-           'create_lb_method_from_existing', 'create_lb_update_rule', 'LBMConfig', 'LBMOptimisation',
-           'Stencil', 'Method', 'ForceModel', 'CollisionSpace',
-           'LatticeBoltzmannStep',
-           'pdf_initialization_assignments', 'macroscopic_values_getter', 'compile_macroscopic_values_getter',
-           'compile_macroscopic_values_setter', 'create_advanced_velocity_setter_collision_rule',
-           'get_weights',
-           'relaxation_rate_from_lattice_viscosity', 'lattice_viscosity_from_relaxation_rate',
-           'relaxation_rate_from_magic_number',
-           'create_lid_driven_cavity', 'create_fully_periodic_flow',
-           'LBStencil']
-
-
-from ._version import get_versions
-__version__ = get_versions()['version']
-del get_versions

lbmpy/boundaries/__init__.py

-from lbmpy.boundaries.boundaryconditions import (
-    UBB, FixedDensity, DiffusionDirichlet, SimpleExtrapolationOutflow,
-    ExtrapolationOutflow, NeumannByCopy, NoSlip, StreamInConstant, FreeSlip)
-from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-
-__all__ = ['NoSlip', 'FreeSlip', 'UBB', 'SimpleExtrapolationOutflow', 'ExtrapolationOutflow',
-           'FixedDensity', 'DiffusionDirichlet', 'NeumannByCopy',
-           'LatticeBoltzmannBoundaryHandling', 'StreamInConstant']

lbmpy/phasefield_allen_cahn/analytical.py

-
-def analytic_rising_speed(gravitational_acceleration, bubble_diameter, viscosity_gas):
-    r"""
-    Calculated the analytical rising speed of a bubble. This is the expected end rising speed.
-    Args:
-        gravitational_acceleration: the gravitational acceleration acting in the simulation scenario. Usually it gets
-                                    calculated based on dimensionless parameters which describe the scenario
-        bubble_diameter: the diameter of the bubble at the beginning of the simulation
-        viscosity_gas: the viscosity of the fluid inside the bubble
-    """
-    result = -(gravitational_acceleration * bubble_diameter * bubble_diameter) / (12.0 * viscosity_gas)
-    return result

lbmpy/phasefield_allen_cahn/phasefield_simplifications.py

-import sympy as sp
-
-from pystencils.simp import (
-    SimplificationStrategy, apply_to_all_assignments,
-    insert_aliases, insert_zeros, insert_constants)
-
-
-def create_phasefield_simplification_strategy(lb_method):
-    s = SimplificationStrategy()
-    expand = apply_to_all_assignments(sp.expand)
-
-    s.add(expand)
-
-    s.add(insert_zeros)
-    s.add(insert_aliases)
-    s.add(insert_constants)
-    s.add(lambda ac: ac.new_without_unused_subexpressions())
-
-    return s

lbmpy_tests/cumulantmethod/test_periodic_pipe_flow.ipynb

-%% Cell type:code id: tags:
-
-``` python
-%load_ext autoreload
-%autoreload 2
-```
-
-%% Cell type:code id: tags:
-
-``` python
-from dataclasses import replace
-
-import sympy as sp
-import numpy as np
-from numpy.testing import assert_allclose
-
-from pystencils.session import *
-from pystencils import Target
-from lbmpy.session import *
-
-from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-    import pycuda.gpuarray as gpuarray
-except ImportError:
-    pycuda = None
-    target = ps.Target.CPU
-    print('No pycuda installed')
-
-if pycuda:
-    target = ps.Target.GPU
-```
-
-%% Output
-
-    No pycuda installed
-
-%% Cell type:markdown id: tags:
-
-## Pipe Flow Scenario
-
-%% Cell type:code id: tags:
-
-``` python
-class PeriodicPipeFlow:
-    def __init__(self, lbm_config, lbm_optimisation, config):
-
-        wall_boundary = NoSlip()
-
-        self.stencil = lbm_config.stencil
-        self.streaming_pattern = lbm_config.streaming_pattern
-        self.target = config.target
-
-        self.gpu = self.target == Target.GPU
-
-        #   Stencil
-        self.q = stencil.Q
-        self.dim = stencil.D
-
-        #   Streaming
-        self.inplace = is_inplace(self.streaming_pattern)
-        self.timesteps = get_timesteps(self.streaming_pattern)
-        self.zeroth_timestep = self.timesteps[0]
-
-        #   Domain, Data Handling and PDF fields
-        self.pipe_length = 60
-        self.pipe_radius = 15
-        self.domain_size = (self.pipe_length, ) + (2 * self.pipe_radius,) * (self.dim - 1)
-        self.periodicity = (True, ) + (False, ) * (self.dim - 1)
-        force = (0.0001, ) + (0.0,) * (self.dim - 1)
-        self.force = lbm_config.force
-
-        self.dh = create_data_handling(domain_size=self.domain_size,
-                                       periodicity=self.periodicity, default_target=self.target)
-
-        self.pdfs = self.dh.add_array('pdfs', self.q)
-        if not self.inplace:
-            self.pdfs_tmp = self.dh.add_array_like('pdfs_tmp', self.pdfs.name)
-
-        lbm_optimisation = replace(lbm_optimisation, symbolic_field=self.pdfs)
-
-        if not self.inplace:
-            lbm_optimisation = replace(lbm_optimisation, symbolic_temporary_field=self.pdfs_tmp)
-
-        self.lb_collision = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_optimisation)
-        self.lb_method = self.lb_collision.method
-
-        self.lb_kernels = []
-        for t in self.timesteps:
-            lbm_config = replace(lbm_config, timestep=t, collision_rule=self.lb_collision)
-            self.lb_kernels.append(create_lb_function(lbm_config=lbm_config,
-                                                      lbm_optimisation=lbm_optimisation,
-                                                      config=config))
-
-        #   Macroscopic Values
-        self.density = 1.0
-        self.density_field = self.dh.add_array('rho', 1)
-        u_x = 0.0
-        self.velocity = (u_x,) * self.dim
-        self.velocity_field = self.dh.add_array('u', self.dim)
-
-        single_gl_config = replace(config, ghost_layers=1)
-        setter = macroscopic_values_setter(
-            self.lb_method, self.density, self.velocity, self.pdfs,
-            streaming_pattern=self.streaming_pattern, previous_timestep=self.zeroth_timestep)
-        self.init_kernel = create_kernel(setter, config=single_gl_config).compile()
-
-        self.getter_kernels = []
-        for t in self.timesteps:
-            getter = macroscopic_values_getter(
-                self.lb_method, self.density_field, self.velocity_field, self.pdfs,
-                streaming_pattern=self.streaming_pattern, previous_timestep=t)
-            self.getter_kernels.append(create_kernel(getter, config=single_gl_config).compile())
-
-        #   Periodicity
-        self.periodicity_handler = LBMPeriodicityHandling(
-            self.stencil, self.dh, self.pdfs.name, streaming_pattern=self.streaming_pattern)
-
-        #   Boundary Handling
-        self.wall = wall_boundary
-        self.bh = LatticeBoltzmannBoundaryHandling(
-            self.lb_method, self.dh, self.pdfs.name,
-            streaming_pattern=self.streaming_pattern, target=self.target)
-
-        self.bh.set_boundary(boundary_obj=self.wall, mask_callback=self.mask_callback)
-
-        self.current_timestep = self.zeroth_timestep
-
-    def mask_callback(self, x, y, z=None):
-        y = y - self.pipe_radius
-        z = z - self.pipe_radius if z is not None else 0
-        return np.sqrt(y**2 + z**2) >= self.pipe_radius
-
-    def init(self):
-        self.current_timestep = self.zeroth_timestep
-        self.dh.run_kernel(self.init_kernel)
-
-    def step(self):
-        #   Order matters! First communicate, then boundaries, otherwise
-        #   periodicity handling overwrites reflected populations
-        # Periodicty
-        self.periodicity_handler(self.current_timestep)
-
-        # Boundaries
-        self.bh(prev_timestep=self.current_timestep)
-
-        # Here, the next time step begins
-        self.current_timestep = self.current_timestep.next()
-
-        # LBM Step
-        self.dh.run_kernel(self.lb_kernels[self.current_timestep.idx])
-
-        # Field Swaps
-        if not self.inplace:
-            self.dh.swap(self.pdfs.name, self.pdfs_tmp.name)
-
-        # Macroscopic Values
-        self.dh.run_kernel(self.getter_kernels[self.current_timestep.idx])
-
-    def run(self, iterations):
-        for _ in range(iterations):
-            self.step()
-
-    @property
-    def velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        return self.dh.gather_array(self.velocity_field.name)
-
-    def get_trimmed_velocity_array(self):
-        if self.gpu:
-            self.dh.to_cpu(self.velocity_field.name)
-        u = np.copy(self.dh.gather_array(self.velocity_field.name))
-        mask = self.bh.get_mask(None, self.wall)
-        for idx in np.ndindex(u.shape[:-1]):
-            if mask[idx] != 0:
-                u[idx] = np.full((self.dim, ), np.nan)
-        return u
-```
-
-%% Cell type:markdown id: tags:
-
-## General Setup
-
-%% Cell type:code id: tags:
-
-``` python
-stencil = LBStencil(Stencil.D3Q19)
-dim = stencil.D
-streaming_pattern = 'aa'
-
-config = ps.CreateKernelConfig(target=target)
-
-viscous_rr = 1.1
-force = (0.0001, ) + (0.0,) * (dim - 1)
-```
-
-%% Cell type:markdown id: tags:
-
-## 1. Reference: SRT Method
-
-%% Cell type:code id: tags:
-
-``` python
-srt_config = LBMConfig(stencil=stencil, method=Method.SRT, relaxation_rate=viscous_rr,
-                       force_model=ForceModel.SIMPLE, force=force, streaming_pattern=streaming_pattern)
-
-srt_flow = PeriodicPipeFlow(srt_config, LBMOptimisation(), config)
-srt_flow.init()
-srt_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-srt_u = srt_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(srt_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x127d43bb0>
-
-
-
-%% Cell type:markdown id: tags:
-
-## 2. Centered Cumulant Method with Central Moment-Space Forcing
-
-%% Cell type:code id: tags:
-
-``` python
-cm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, stencil=stencil, compressible=True,
-                      relaxation_rate=viscous_rr,
-                      force=force, streaming_pattern=streaming_pattern)
-
-lbm_opt = LBMOptimisation(pre_simplification=True)
-
-cm_impl_f_flow = PeriodicPipeFlow(cm_config, lbm_opt, config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_impl_f_flow.lb_method
-assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_flow.init()
-cm_impl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_impl_f_u = cm_impl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_impl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x127d1df10>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_impl_f_u, srt_u, rtol=1, atol=1e-4)
-```
-
-%% Cell type:markdown id: tags:
-
-## 3. Centered Cumulant Method with Simple force model
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.forcemodels import Simple
-
-cm_expl_force_config = LBMConfig(method=Method.CUMULANT, stencil=stencil, compressible=True,
-                      relaxation_rate=viscous_rr, force_model=Simple(force),
-                      force=force, streaming_pattern=streaming_pattern)
-
-lbm_opt = LBMOptimisation(pre_simplification=True)
-
-cm_expl_f_flow = PeriodicPipeFlow(cm_expl_force_config, lbm_opt, config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-lb_method = cm_expl_f_flow.lb_method
-assert all(rr == 0 for rr in lb_method.relaxation_rates[1:4])
-assert all(rr == viscous_rr for rr in lb_method.relaxation_rates[4:9])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_flow.init()
-cm_expl_f_flow.run(400)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-cm_expl_f_u = cm_expl_f_flow.get_trimmed_velocity_array()
-ps.plot.vector_field_magnitude(cm_expl_f_u[30,:,:,:])
-ps.plot.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x1272c5760>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert_allclose(cm_expl_f_u, srt_u, rtol=1, atol=1e-4)
-assert_allclose(cm_expl_f_u, cm_impl_f_u, rtol=1, atol=1e-4)
-```

lbmpy_tests/test_compiled_in_boundaries.ipynb

-%% Cell type:code id: tags:
-
-``` python
-import pystencils as ps
-from lbmpy.session import *
-from lbmpy.boundaries.boundaries_in_kernel import update_rule_with_push_boundaries
-from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
-from collections import OrderedDict
-from time import perf_counter
-```
-
-%% Cell type:markdown id: tags:
-
-# Version 1: compile-in boundaries
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (32, 32, 32)
-relaxation_rate = 1.8
-time_steps = 100
-lid_velocity = 0.05
-stencil = LBStencil(Stencil.D3Q19)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, default_target=ps.Target.CPU)
-pdfs = dh.add_array('pdfs', values_per_cell=stencil.Q)
-u = dh.add_array('u', values_per_cell=stencil.D)
-streaming_pattern = 'aa'
-```
-
-%% Cell type:code id: tags:
-
-``` python
-boundaries = OrderedDict((
-    ((0, 1, 0), UBB([lid_velocity, 0, 0])),
-    ((1, 0, 0), NoSlip()),
-    ((-1, 0, 0), NoSlip()),
-    ((0, -1, 0), NoSlip()),
-    ((0, 0, 1), NoSlip()),
-    ((0, 0, -1), NoSlip()),
-))
-
-lbm_opt = LBMOptimisation(symbolic_field=pdfs, cse_global=False, cse_pdfs=False)
-lbm_config = LBMConfig(stencil=stencil, relaxation_rate=relaxation_rate, compressible=False)
-
-cr_even = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-cr_odd = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
-update_rule_aa_even = update_rule_with_push_boundaries(cr_even, pdfs, boundaries, streaming_pattern, Timestep.EVEN)
-update_rule_aa_odd = update_rule_with_push_boundaries(cr_odd, pdfs, boundaries, streaming_pattern, Timestep.ODD)
-
-getter_assignments = macroscopic_values_getter(update_rule_aa_even.method, velocity=u.center_vector,
-                                               pdfs=pdfs, density=None,
-                                               streaming_pattern=streaming_pattern,
-                                               previous_timestep=Timestep.EVEN)
-
-config = ps.CreateKernelConfig(target=dh.default_target)
-getter_kernel = ps.create_kernel(getter_assignments, config=config).compile()
-even_kernel = ps.create_kernel(update_rule_aa_even, config=config).compile()
-odd_kernel = ps.create_kernel(update_rule_aa_odd, config=config).compile()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:code id: tags:
-
-``` python
-def init():
-    dh.fill(pdfs.name, 0, ghost_layers=True)
-
-def aa_time_loop(steps=100):
-    assert steps % 2 == 0, "Works only for an even number of time steps"
-    dh.all_to_gpu()
-    for i in range(steps // 2):
-        dh.run_kernel(odd_kernel)
-        dh.run_kernel(even_kernel)
-    dh.run_kernel(getter_kernel)
-    dh.all_to_cpu()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-aa_time_loop(time_steps)
-vel_version1 = dh.gather_array(u.name, ghost_layers=False).copy()
-plt.vector_field_magnitude(vel_version1[:, :, domain_size[2]//2, :])
-plt.colorbar()
-```
-
-%% Output
-
-    <matplotlib.colorbar.Colorbar at 0x13b12c820>
-
-
-
-%% Cell type:markdown id: tags:
-
-# Version 2: Normal boundary handling
-
-%% Cell type:code id: tags:
-
-``` python
-ldc = create_lid_driven_cavity(domain_size, relaxation_rate=relaxation_rate, lid_velocity=lid_velocity)
-ldc.run(time_steps)
-vel_version2 = ldc.velocity[:, :, :, :]
-
-plt.vector_field_magnitude(vel_version2[:, :, domain_size[2]//2, :])
-plt.colorbar()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-
-    <matplotlib.colorbar.Colorbar at 0x149a4a190>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-np.testing.assert_almost_equal(vel_version1[1:-1, 1:-1, :], vel_version2[1:-1, 1:-1, :], decimal=2)
-```

lbmpy_tests/test_float_kernel.py

-import pytest
-
-import pystencils as ps
-
-from lbmpy.creationfunctions import create_lb_function, LBMConfig
-from lbmpy.enums import Method
-from lbmpy.scenarios import create_lid_driven_cavity
-
-
-@pytest.mark.parametrize('double_precision', [False, True])
-@pytest.mark.parametrize('method_enum', [Method.SRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
-def test_creation(double_precision, method_enum):
-    """Simple test that makes sure that only float variables are created"""
-    lbm_config = LBMConfig(method=method_enum, relaxation_rate=1.5, compressible=True)
-    config = ps.CreateKernelConfig(data_type="float64" if double_precision else "float32",
-                                   default_number_float="float64" if double_precision else "float32")
-    func = create_lb_function(lbm_config=lbm_config, config=config)
-    code = ps.get_code_str(func)
-
-    if double_precision:
-        assert 'float' not in code
-        assert 'double' in code
-    else:
-        assert 'double' not in code
-        assert 'float' in code
-
-
-@pytest.mark.parametrize('double_precision', [False, True])
-@pytest.mark.parametrize('method_enum', [Method.SRT, Method.CENTRAL_MOMENT, Method.CUMULANT])
-def test_scenario(method_enum, double_precision):
-    lbm_config = LBMConfig(method=method_enum, relaxation_rate=1.45, compressible=True)
-    config = ps.CreateKernelConfig(data_type="float64" if double_precision else "float32",
-                                   default_number_float="float64" if double_precision else "float32")
-    sc = create_lid_driven_cavity((16, 16, 8), lbm_config=lbm_config, config=config)
-    sc.run(1)
-    code = ps.get_code_str(sc.ast)
-
-    if double_precision:
-        assert 'float' not in code
-        assert 'double' in code
-    else:
-        assert 'double' not in code
-        assert 'float' in code

lbmpy_tests/test_sparse_lbm.ipynb

-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.sparse import *
-from pystencils.field import FieldType
-
-import pystencils as ps
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-    import pycuda.gpuarray as gpuarray
-except ImportError:
-    pycuda = None
-    target = ps.Target.CPU
-    print('No pycuda installed')
-
-if pycuda:
-    target = ps.Target.GPU
-```
-
-%% Output
-
-    No pycuda installed
-
-%% Cell type:markdown id: tags:
-
-# Sparse (list based) LBM
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (20, 20)
-omega = 1.8
-target = ps.Target.CPU
-
-ghost_layers = 1
-arr_size = tuple(ds + 2 * ghost_layers for ds in domain_size)
-lid_velocity = 0.01
-force = 1e-6
-
-channel = True
-if channel:
-    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), relaxation_rate=omega,
-                           compressible=False, force=(force, 0))
-else:
-    lbm_config = LBMConfig(stencil=LBStencil(Stencil.D2Q9), relaxation_rate=omega, compressible=False)
-
-
-method = create_lb_method(lbm_config=lbm_config)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-ubb = UBB( (lid_velocity, 0) )
-noslip = NoSlip()
-flags = {
-    'fluid': 1,
-    noslip: 2,
-    ubb: 4,
-}
-flag_arr = np.zeros(arr_size, dtype=np.uint16)
-flag_arr.fill(flags['fluid'])
-
-if channel:
-    flag_arr[0, :] = 0
-    flag_arr[-1, :] = 0
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[:, -1] = flags[noslip]
-else:
-    flag_arr[:, -1] = flags[ubb]
-    flag_arr[:, 0] = flags[noslip]
-    flag_arr[0, :] = flags[noslip]
-    flag_arr[-1, :] = flags[noslip]
-
-plt.scalar_field(flag_arr)
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Mappings
-
-%% Cell type:code id: tags:
-
-``` python
-mapping = SparseLbMapper(method.stencil, flag_arr, flags['fluid'], flags[noslip], flags[ubb])
-index_arr = mapping.create_index_array(ghost_layers)
-
-# Arrays
-#index_arr = index_arr_linear.reshape([len(method.stencil), mapping.num_fluid_cells])
-#index_arr = index_arr.swapaxes(0, 1)
-
-pdf_arr = np.empty((len(mapping), len(method.stencil)), order='f')
-pdf_arr_tmp = np.empty_like(pdf_arr)
-
-vel_arr = np.ones([mapping.num_fluid_cells, method.dim], order='f')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-pdf_field, pdf_field_tmp, vel_field = ps.fields("f(9), d(9), u(2): [1D]",
-                                                #f=pdf_arr[:mapping.num_fluid_cells],
-                                                #d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                                                #u=vel_arr
-                                               )
-pdf_field.field_type = FieldType.CUSTOM
-pdf_field.pdf_field_tmp = FieldType.CUSTOM
-```
-
-%% Cell type:markdown id: tags:
-
-### Macroscopic quantities
-
-%% Cell type:code id: tags:
-
-``` python
-cqc = method.conserved_quantity_computation
-inp_eqs = cqc.equilibrium_input_equations_from_init_values(force_substitution=False)
-setter_eqs = method.get_equilibrium(conserved_quantity_equations=inp_eqs)
-setter_eqs = setter_eqs.new_with_substitutions({sym: pdf_field(i)
-                                                for i, sym in enumerate(method.post_collision_pdf_symbols)})
-
-config = ps.CreateKernelConfig(ghost_layers=((0, 0),))
-kernel_initialize = ps.create_kernel(setter_eqs, config=config).compile()
-
-def init():
-    kernel_initialize(f=pdf_arr[:mapping.num_fluid_cells])
-init()
-
-
-getter_eqs = cqc.output_equations_from_pdfs(pdf_field.center_vector,
-                                            {'velocity': vel_field})
-kernel_compute_u = ps.create_kernel(getter_eqs, config=config).compile()
-
-def get_velocity(arr=pdf_arr):
-    fluid_cell_arr = mapping.coordinates
-    kernel_compute_u(f=pdf_arr[:mapping.num_fluid_cells], u=vel_arr)
-    full_velocity_arr = np.zeros(flag_arr.shape + (2,), dtype=np.float64)
-    full_velocity_arr.fill(np.nan)
-    arr = fluid_cell_arr[:mapping.num_fluid_cells]
-    full_velocity_arr[arr['x'], arr['y']] = vel_arr
-    return full_velocity_arr[1:-1, 1:-1]
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:markdown id: tags:
-
-### Stream Collide Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-#index_field = ps.Field.create_from_numpy_array("idx", index_arr, index_dimensions=1)
-index_field = ps.Field.create_generic("idx", spatial_dimensions=1, index_dimensions=1, dtype=index_arr.dtype)
-
-collision_rule = method.get_collision_rule()
-
-Q = len(method.stencil)
-symbol_subs = {sym: pdf_field.absolute_access((index_field(i-1),),())
-               for i, sym in enumerate(method.pre_collision_pdf_symbols)}
-
-symbol_subs.update({sym: pdf_field_tmp(i) for i, sym in enumerate(method.post_collision_pdf_symbols)})
-
-symbol_subs[method.pre_collision_pdf_symbols[0]] = pdf_field(0) # special case for center
-symbol_subs
-```
-
-%% Output
-
-
-    $\displaystyle \left\{ d_{0} : {d}_{0}^{0}, \  d_{1} : {d}_{0}^{1}, \  d_{2} : {d}_{0}^{2}, \  d_{3} : {d}_{0}^{3}, \  d_{4} : {d}_{0}^{4}, \  d_{5} : {d}_{0}^{5}, \  d_{6} : {d}_{0}^{6}, \  d_{7} : {d}_{0}^{7}, \  d_{8} : {d}_{0}^{8}, \  f_{0} : {f}_{0}^{0}, \  f_{1} : {f}_{\mathbf{idx}_{0}^{0}}^{0}, \  f_{2} : {f}_{\mathbf{idx}_{0}^{1}}^{0}, \  f_{3} : {f}_{\mathbf{idx}_{0}^{2}}^{0}, \  f_{4} : {f}_{\mathbf{idx}_{0}^{3}}^{0}, \  f_{5} : {f}_{\mathbf{idx}_{0}^{4}}^{0}, \  f_{6} : {f}_{\mathbf{idx}_{0}^{5}}^{0}, \  f_{7} : {f}_{\mathbf{idx}_{0}^{6}}^{0}, \  f_{8} : {f}_{\mathbf{idx}_{0}^{7}}^{0}\right\}$
-    {d₀: d_C__0, d₁: d_C__1, d₂: d_C__2, d₃: d_C__3, d₄: d_C__4, d₅: d_C__5, d₆: d_C__6, d₇: d_C__7, d₈: d_C__8, f₀: f_C__0, f₁: f_26b74c363b15, f₂: f_e111196926c6, f₃: f_e3f72afe7d66, f₄: f_0b929cb3f1da, f₅: f_3f75acc2de6d,
-    f₆: f_3e20adce708c, f₇: f_3a33f411da6b, f₈: f_68da5b60e7d8}
-
-%% Cell type:code id: tags:
-
-``` python
-collision_rule = method.get_collision_rule()
-update_rule = collision_rule.new_with_substitutions(symbol_subs)
-
-config = ps.CreateKernelConfig(ghost_layers=((0, 0),), target=target)
-kernel_stream_collide = ps.create_kernel(update_rule, config=config).compile()
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-
-%% Cell type:markdown id: tags:
-
-### Boundary Kernels
-
-%% Cell type:code id: tags:
-
-``` python
-if not channel:
-    if target == ps.Target.GPU:
-        raise NotImplementedError("UBB on GPU not working yet")
-
-    ubb_mapper = SparseLbBoundaryMapper(ubb, method, pdf_field)
-    #TODO the following line is wrong: kernel contains accesses to index_field and pdf_field which have
-    #different size: a correct kernel comes out when by change the shape is taken from index field,
-    # when taken from pdf field, a wrong kernel is generated
-    ubb_kernel = ps.create_kernel(ubb_mapper.assignments(), ghost_layers=0).compile()
-    ubb_idx_arr = ubb_mapper.create_index_arr(mapping, flags[ubb])
-    ps.show_code(ubb_kernel.ast)
-    def handle_ubb():
-        ubb_kernel(indexField=ubb_idx_arr, f=pdf_arr[:mapping.num_fluid_cells])
-else:
-    def handle_ubb():
-        pass
-```
-
-%% Cell type:markdown id: tags:
-
-### Time Loop
-
-%% Cell type:code id: tags:
-
-``` python
-def time_step():
-    global pdf_arr, pdf_arr_tmp, index_arr
-    handle_ubb()
-    if target == ps.Target.GPU:
-        gpu_pdf_arr = gpuarray.to_gpu(pdf_arr)
-        gpu_pdf_arr_tmp = gpuarray.to_gpu(pdf_arr_tmp)
-        gpu_index_arr = gpuarray.to_gpu(index_arr)
-
-        kernel_stream_collide(f=gpu_pdf_arr[:mapping.num_fluid_cells],
-                          d=gpu_pdf_arr_tmp[:mapping.num_fluid_cells],
-                          idx=gpu_index_arr)
-
-        pdf_arr = gpu_pdf_arr.get()
-        pdf_arr_tmp = gpu_pdf_arr_tmp.get()
-        index_arr = gpu_index_arr.get()
-    else:
-        kernel_stream_collide(f=pdf_arr[:mapping.num_fluid_cells],
-                              d=pdf_arr_tmp[:mapping.num_fluid_cells],
-                              idx=index_arr)
-    pdf_arr_tmp, pdf_arr = pdf_arr, pdf_arr_tmp
-
-def run(steps=100):
-    for t in range(steps):
-        time_step()
-    return get_velocity()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-result = run(100)
-plt.vector_field(result, step=1);
-```
-
-%% Output
-
-
-
-%% Cell type:markdown id: tags:
-
-### Check against reference
-
-%% Cell type:code id: tags:
-
-``` python
-if channel:
-    reference = create_channel(domain_size, force=force, lb_method=method)
-else:
-    reference = create_lid_driven_cavity(domain_size, relaxation_rate=omega, lid_velocity=lid_velocity,
-                                         compressible=False)
-reference.run(100)
-```
-
-%% Output
-
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Using Nodes is experimental and not fully tested. Double check your generated code!
-    WARNING:root:Lhs"dir of type "int64_t" is assigned with a different datatype rhs: "indexField[0](dir)" of type "int32_t".
-
-%% Cell type:code id: tags:
-
-``` python
-np.testing.assert_almost_equal(reference.velocity[:, :], result)
-```