CHANGELOG.md

+# Change Log
+
+## Unreleased
+
+### Removed
+* Removing OpenCL support because it is not supported by pystencils anymore

doc/notebooks/03_tutorial_lbm_formulation.ipynb

 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.session import *
 ```
 
 %% Cell type:markdown id: tags:
 
 # Tutorial 03: Defining LB methods in *lbmpy*
 
 
 ## A) General Form
 
 %% Cell type:markdown id: tags:
 
 The lattice Boltzmann equation in its most general form is:
 
 $$f_q(\mathbf{x} + \mathbf{c}_q \delta t, t+\delta t) = K\left( f_q(\mathbf{x}, t) \right)$$
 
 with a discrete velocity set $\mathbf{c}_q$ (stencil) and a generic collision operator $K$.
 
 So a lattice Boltzmann method can be fully defined by picking a stencil and a collision operator.
 The collision operator $K$ has the following structure:
 - Transformation of particle distribution function $f$ into collision space. This transformation has to be invertible and may be nonlinear.
 - The collision operation is an convex combination of the pdf representation in collision space $c$ and some equilibrium vector $c^{(eq)}$. This equilibrium can also be defined in physical space, then $c^{(eq)} = C( f^{(eq)} ) $. The convex combination is done elementwise using a diagonal matrix $S$ where the diagonal entries are the relaxation rates.
 - After collision, the collided state $c'$ is transformed back into physical space
 
 ![](../img/collision.svg)
 
 
 The full collision operator is:
 
 $$K(f) = C^{-1}\left( (I-S)C(f) + SC(f^{(eq}) \right)$$
 
 or
 
 $$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$$
 
 %% Cell type:markdown id: tags:
 
 ## B) Moment-based relaxation
 
 The most commonly used LBM collision operator is the multi relaxation time (MRT) collision.
 In MRT methods the collision space is spanned by moments of the distribution function. This is a very natural approach, since the pdf moments are the quantities that go into the Chapman Enskog analysis that is used to show that LB methods can solve the Navier Stokes equations. Also the lower order moments correspond to the macroscopic quantities of interest (density/pressure, velocity, shear rates, heat flux). Furthermore the transformation to collision space is linear in this case, simplifying the collision equations:
 
 $$K(f) = C^{-1}\left( C(f) - S (C(f) - C(f^{(eq})) \right)$$
 
 $$K(f) = f - \underbrace{ C^{-1}SC}_{A}(f - f^{(eq)})$$
 
 in *lbmpy* the following formulation is used, since it is more natural to define the equilibrium in moment-space instead of physical space:
 
 $$K(f) = f -  C^{-1}S(Cf - c^{(eq)})$$
 
 %% Cell type:markdown id: tags:
 
 ### Use a pre-defined method
 
 Lets create a moment-based method in *lbmpy* and see how the moment transformation $C$ and the relaxation rates that comprise the diagonal matrix $S$ can be defined. We start with a function that creates a basic MRT model.
 Don't use this for real simulations, there orthogonalized MRT methods should be used, as discussed in the next section.
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.creationfunctions import create_lb_method
 lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT_RAW)
 
 method = create_lb_method(lbm_config=lbm_config)
 # check also method='srt', 'trt', 'mrt'
 method
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x107ac7a30>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x117fd1d90>
 
 %% Cell type:markdown id: tags:
 
 The first column labeled "Moment" defines the collision space and thus the transformation matrix $C$.
 The remaining columns specify the equilibrium vector in moment space $c^{(eq)}$ and the corresponding relaxation rate.
 
 Each row of the "Moment" column defines one row of $C$. In the next cells this matrix and the discrete velocity set (stencil) of our method are shown. Check for example the second last row of the table $x^2 y$: In the corresponding second last row of the moment matrix $C$ where each column stands for a lattice velocity (for ordering visualized stencil below) and each entry is the expression $x^2 y$ where $x$ and $y$ are the components of the lattice velocity.
 
 In general the transformation matrix $C_{iq}$ is defined as;
 
-$$c_i = C_{iq} f_q = \sum_q  m_i(c_q)$$
+$$c_i = C_{iq} f_q = \sum_q  m_i(\mathbf{c}_q)$$
 
-where $m_i(c_q)$ is the $i$'th moment polynomial where $x$ and $y$ are substituted with the components of the $q$'th lattice velocity
+where $m_i(\mathbf{c}_q)$ is the $i$'th moment polynomial where $x$ and $y$ are substituted with the components of the $q$'th lattice velocity
 
 %% Cell type:code id: tags:
 
 ``` python
 # Transformation matrix C
 method.moment_matrix
 ```
 
 %% Output
 
-
+
     $\displaystyle \left[\begin{matrix}1 & 1 & 1 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 0 & 0 & -1 & 1 & -1 & 1 & -1 & 1\\0 & 1 & -1 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 1\\0 & 1 & 1 & 0 & 0 & 1 & 1 & 1 & 1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & 1 & -1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & -1 & -1\\0 & 0 & 0 & 0 & 0 & -1 & 1 & -1 & 1\\0 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & 1\end{matrix}\right]$
     ⎡1  1  1   1   1  1   1  1   1 ⎤
     ⎢                              ⎥
     ⎢0  0  0   -1  1  -1  1  -1  1 ⎥
     ⎢                              ⎥
     ⎢0  1  -1  0   0  1   1  -1  -1⎥
     ⎢                              ⎥
     ⎢0  0  0   1   1  1   1  1   1 ⎥
     ⎢                              ⎥
     ⎢0  1  1   0   0  1   1  1   1 ⎥
     ⎢                              ⎥
     ⎢0  0  0   0   0  -1  1  1   -1⎥
     ⎢                              ⎥
     ⎢0  0  0   0   0  1   1  -1  -1⎥
     ⎢                              ⎥
     ⎢0  0  0   0   0  -1  1  -1  1 ⎥
     ⎢                              ⎥
     ⎣0  0  0   0   0  1   1  1   1 ⎦
 
 %% Cell type:code id: tags:
 
 ``` python
 method.stencil.plot()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 method.stencil
 ```
 
 %% Output
 
-    <lbmpy.stencils.LBStencil at 0x13ae6b670>
+    <lbmpy.stencils.LBStencil at 0x117c0bdf0>
 
 %% Cell type:markdown id: tags:
 
 ### Orthogonal MRTs
 
 For a real MRT method, the moments should be orthogonalized.
 One can either orthogonalize using the standard scalar product or a scalar product that is weighted with the lattice weights. If unsure, use the weighted version.
 
 The next cell shows how to get both orthogonalized MRT versions in lbmpy.
 
 %% Cell type:code id: tags:
 
 ``` python
 rr = [sp.Symbol('omega_shear'), sp.Symbol('omega_bulk'), sp.Symbol('omega_3'), sp.Symbol('omega_4')]
 
 lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=True,
                        relaxation_rates=rr)
 weighted_ortho_mrt = create_lb_method(lbm_config=lbm_config)
 weighted_ortho_mrt
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13ae91df0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x117e9e5b0>
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, weighted=False,
                        relaxation_rates=rr)
 ortho_mrt = create_lb_method(lbm_config=lbm_config)
 ortho_mrt
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13accd610>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x118216b20>
 
 %% Cell type:markdown id: tags:
 
 One can check if a method is orthogonalized:
 
 %% Cell type:code id: tags:
 
 ``` python
 ortho_mrt.is_orthogonal, weighted_ortho_mrt.is_weighted_orthogonal
 ```
 
 %% Output
 
     (True, True)
 
 %% Cell type:markdown id: tags:
 
 Note here how the relaxation rate for each moment is set. By default, a single relaxation rate per moment group is assigned. With a moment group, we mean here moment polynomials of the same order. However, the second-order moments are split between the shear and bulk moments (defining the shear and bulk viscosity). Thus the first two relaxation rates in the `rr`-list are used for these.
 
 If only a single relaxation rate would be defined, it is applied to the shear moments, and all other higher-order moments are set to one. This means the PDFs are directly fixed to the equilibrium in the collision process. Furthermore, the relaxation rate for conserved moments can be chosen freely as they are conserved in the collision. Thus by default, we set it to zero.
 
 It is also possible to define an individual relaxation rate per moment. In this case, `Q` relaxation rates need to be defined, where `Q` is the number of moments.
 
 %% Cell type:markdown id: tags:
 
 ### Central moment lattice Boltzmann methods
 
 %% Cell type:markdown id: tags:
 
 Another popular method is the cascaded lattice Boltzmann method. The cascaded LBM increases the numerical stability by shifting the collision step to the central moment space. Thus it is applied in the non-moving frame and achieves a better Galilean invariance. Typically the central moment collision operator is derived for compressible LB methods, and a higher-order equilibrium is used. Although incompressible LB methods with a second-order equilibrium can be derived with lbmpy, it violates the Galilean invariance and thus reduces the advantages of the method.
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CENTRAL_MOMENT, equilibrium_order=4,
                        compressible=True, relaxation_rates=rr)
 
 central_moment_method = create_lb_method(lbm_config)
 central_moment_method
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x13ad7e460>
+    <lbmpy.methods.momentbased.centralmomentbasedmethod.CentralMomentBasedLbMethod at 0x117e7d3d0>
 
 %% Cell type:markdown id: tags:
 
 The shift to the central moment space is done by applying a so-called shift matrix. Usually, this introduces a high numerical overhead. This problem is solved with lbmpy because each transformation stage can be specifically optimised individually. Therefore, it is possible to derive a CLBM with only a little numerical overhead.
 
 %% Cell type:code id: tags:
 
 ``` python
 central_moment_method.shift_matrix
 ```
 
 %% Output
 
-
+
     $\displaystyle \left[\begin{matrix}1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{0} & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\- u_{1} & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0\\u_{0} u_{1} & - u_{1} & - u_{0} & 1 & 0 & 0 & 0 & 0 & 0\\u_{0}^{2} - u_{1}^{2} & - 2 u_{0} & 2 u_{1} & 0 & 1 & 0 & 0 & 0 & 0\\u_{0}^{2} + u_{1}^{2} & - 2 u_{0} & - 2 u_{1} & 0 & 0 & 1 & 0 & 0 & 0\\- u_{0}^{2} u_{1} & 2 u_{0} u_{1} & u_{0}^{2} & - 2 u_{0} & - \frac{u_{1}}{2} & - \frac{u_{1}}{2} & 1 & 0 & 0\\- u_{0} u_{1}^{2} & u_{1}^{2} & 2 u_{0} u_{1} & - 2 u_{1} & \frac{u_{0}}{2} & - \frac{u_{0}}{2} & 0 & 1 & 0\\u_{0}^{2} u_{1}^{2} & - 2 u_{0} u_{1}^{2} & - 2 u_{0}^{2} u_{1} & 4 u_{0} u_{1} & - \frac{u_{0}^{2}}{2} + \frac{u_{1}^{2}}{2} & \frac{u_{0}^{2}}{2} + \frac{u_{1}^{2}}{2} & - 2 u_{1} & - 2 u_{0} & 1\end{matrix}\right]$
     ⎡    1          0          0         0          0           0        0      0    0⎤
     ⎢                                                                                 ⎥
     ⎢   -u₀         1          0         0          0           0        0      0    0⎥
     ⎢                                                                                 ⎥
     ⎢   -u₁         0          1         0          0           0        0      0    0⎥
     ⎢                                                                                 ⎥
     ⎢  u₀⋅u₁       -u₁        -u₀        1          0           0        0      0    0⎥
     ⎢                                                                                 ⎥
     ⎢  2     2                                                                        ⎥
     ⎢u₀  - u₁     -2⋅u₀      2⋅u₁        0          1           0        0      0    0⎥
     ⎢                                                                                 ⎥
     ⎢  2     2                                                                        ⎥
     ⎢u₀  + u₁     -2⋅u₀      -2⋅u₁       0          0           1        0      0    0⎥
     ⎢                                                                                 ⎥
     ⎢    2                      2                 -u₁         -u₁                     ⎥
     ⎢ -u₀ ⋅u₁    2⋅u₀⋅u₁      u₀       -2⋅u₀      ────        ────       1      0    0⎥
     ⎢                                              2           2                      ⎥
     ⎢                                                                                 ⎥
     ⎢       2        2                             u₀         -u₀                     ⎥
     ⎢ -u₀⋅u₁       u₁       2⋅u₀⋅u₁    -2⋅u₁       ──         ────       0      1    0⎥
     ⎢                                              2           2                      ⎥
     ⎢                                                                                 ⎥
     ⎢                                              2     2    2     2                 ⎥
     ⎢   2   2           2       2                u₀    u₁   u₀    u₁                  ⎥
     ⎢ u₀ ⋅u₁    -2⋅u₀⋅u₁   -2⋅u₀ ⋅u₁  4⋅u₀⋅u₁  - ─── + ───  ─── + ───  -2⋅u₁  -2⋅u₀  1⎥
     ⎣                                             2     2    2     2                  ⎦
 
 %% Cell type:markdown id: tags:
 
 ### Define custom MRT method
 
 To choose custom values for the left moment column one can pass a nested list of moments.
 Moments that should be relaxed with the same paramter are grouped together.
 
 *lbmpy* also comes with a few templates for this list taken from literature:
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.methods import mrt_orthogonal_modes_literature
 from lbmpy.stencils import get_stencil
 from lbmpy.moments import MOMENT_SYMBOLS
 
 x, y, z = MOMENT_SYMBOLS
 
 moments = mrt_orthogonal_modes_literature(LBStencil(Stencil.D2Q9), is_weighted=True)
 moments
 ```
 
 %% Output
 
-
+
     $\displaystyle \left[ \left[ 1\right], \  \left[ x, \  y\right], \  \left[ 3 x^{2} + 3 y^{2} - 2, \  x^{2} - y^{2}, \  x y\right], \  \left[ x \left(3 x^{2} + 3 y^{2} - 4\right), \  y \left(3 x^{2} + 3 y^{2} - 4\right)\right], \  \left[ - 15 x^{2} - 15 y^{2} + 9 \left(x^{2} + y^{2}\right)^{2} + 2\right]\right]$
     ⎡                                                                                          ⎡                             2
     ⎢             ⎡   2      2       2    2     ⎤  ⎡  ⎛   2      2    ⎞    ⎛   2      2    ⎞⎤  ⎢      2       2     ⎛ 2    2⎞
     ⎣[1], [x, y], ⎣3⋅x  + 3⋅y  - 2, x  - y , x⋅y⎦, ⎣x⋅⎝3⋅x  + 3⋅y  - 4⎠, y⋅⎝3⋅x  + 3⋅y  - 4⎠⎦, ⎣- 15⋅x  - 15⋅y  + 9⋅⎝x  + y ⎠
     
         ⎤⎤
         ⎥⎥
      + 2⎦⎦
 
 %% Cell type:markdown id: tags:
 
 This nested moment list can be passed to `create_lb_method`:
 
 %% Cell type:code id: tags:
 
 ``` python
 create_lb_method(LBMConfig(stencil=Stencil.D2Q9, method=Method.MRT, nested_moments=moments,
                            relaxation_rates=rr))
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13acdf850>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x117f86f70>
 
 %% Cell type:markdown id: tags:
 
 If one needs to also specify custom equilibrium moments the following approach can be used
 
 %% Cell type:code id: tags:
 
 ``` python
 rho = sp.symbols("rho")
 u = sp.symbols("u_:3")
 omega = sp.symbols("omega_:4")
 
 method_table = [
   # Conserved moments
   (1,    rho,  0 ),
   (x,    u[0], 0 ),
   (y,    u[1], 0 ),
 
   # Shear moments
   (x*y,       u[0]*u[1],                   omega[0]),
   (x**2-y**2, u[0]**2 - u[1]**2,           omega[0]),
   (x**2+y**2, 2*rho/3 + u[0]**2 + u[1]**2, omega[1]),
 
   # Higher order
   (x * y**2,    u[0]/3,                        omega[2]),
   (x**2 * y,    u[1]/3,                        omega[2]),
   (x**2 * y**2, rho/9 + u[0]**2/3 + u[1]**2/3, omega[3]),
 ]
 method = create_generic_mrt(LBStencil(Stencil.D2Q9), method_table)
 method
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13b22f5b0>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x117edb4f0>
 
 %% Cell type:markdown id: tags:
 
 Instead of manually defining all entries in the method table, *lbmpy* has functions to fill the table according to a specific pattern. For example:
 - for a full stencil (D2Q9, D3Q27) there exist exactly 9 or 27 linearly independent moments. These can either be taken as they are, or orthogonalized using Gram-Schmidt, weighted Gram-Schmidt or a Hermite approach
 - equilibrium values can be computed from the standard discrete equilibrium of order 1,2 or 3. Alternatively they can also be computed as continuous moments of a Maxwellian distribution
 
 One option is to start with one of *lbmpy*'s built-in methods and modify it with `create_lb_method_from_existing`.
 In the next cell we fix the fourth order relaxation rate to a constant, by writing a function that defines how to alter each row of the collision table. This is for demonstration only, of course we could have done it right away when passing in the collision table.
 
 %% Cell type:code id: tags:
 
 ``` python
 def modification_func(moment, eq, rate):
     if rate == omega[3]:
         return moment, eq, 1.0
     return moment, eq, rate
 
 
 method = create_lb_method_from_existing(method, modification_func)
 method
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x13b212460>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x1182b41c0>
 
 %% Cell type:markdown id: tags:
 
 Our customized method can be directly passed into one of the scenarios. We can for example set up a channel flow with it. Since we used symbols as relaxation rates, we have to pass them in as `kernel_params`.
 
 %% Cell type:code id: tags:
 
 ``` python
 ch = create_channel(domain_size=(100, 30), lb_method=method, u_max=0.05,
                     kernel_params={'omega_0': 1.8, 'omega_1': 1.4, 'omega_2': 1.5})
 ch.run(500)
 plt.figure(dpi=200)
 plt.vector_field(ch.velocity[:, :]);
 ```
 
 %% Output
 
 
 %% Cell type:markdown id: tags:
 
 #### Bonus: Automatic analysis
 Above we have created a non-orthogonal MRT, where the shear viscosity and bulk viscosity can be independently controlled. For moment-based methods, *lbmpy* also offers an automatic Chapman Enskog analysis that can find the relation between viscosity and relaxation rate(s):
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.chapman_enskog import ChapmanEnskogAnalysis
 analysis = ChapmanEnskogAnalysis(method)
 analysis.get_dynamic_viscosity()
 ```
 
 %% Output
 
-
+
     $\displaystyle - \frac{\omega_{0} - 2}{6 \omega_{0}}$
     -(ω₀ - 2)
     ──────────
        6⋅ω₀
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis.get_bulk_viscosity()
 ```
 
 %% Output
 
-
+
     $\displaystyle - \frac{1}{9} - \frac{1}{3 \omega_{1}} + \frac{5}{9 \omega_{0}}$
       1    1      5
     - ─ - ──── + ────
       9   3⋅ω₁   9⋅ω₀

doc/notebooks/04_tutorial_cumulant_LBM.ipynb

 %% Cell type:code id: tags:
 ``` python
 from lbmpy.session import *
 ```
 %% Cell type:markdown id: tags:
 # 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( \vec{X} ) =
+    M( \mathbf{X} ) =
     \sum_i f_i
     \exp \left(
-        \vec{c}_i \cdot \vec{X}
+        \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_{\vec{X}=0}
+    \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{align}
-    C(\vec{X}) :=& \, \log ( M(\vec{X}) ) \\
+    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_{\vec{X}=0}
+    \right\vert_{\mathbf{X}=0}
 \end{align}
 $$
 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. This is done for two primary reasons:
  1. All cumulants of orders 2 and 3 are equal to their corresponding central moments, up to the density $\rho$ as a proportionality factor.
  2. The conserved modes of the first order, which correspond to momentum, are relaxed in central moment space to allow for a more efficient implicit forcing scheme.
-The central moment-generating function $K$ can be related to the moment-generating function through $K(\vec{X}) = \exp( - \vec{X} \cdot \vec{u} ) M(\vec{X})$. It is possible to recombine the equation with the definition of the cumulant-generating function
+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( \vec{X} ) = \vec{X} \cdot \vec{u} + \log( K( \vec{X} ) ).
+    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 collision operator, forces can be applied through a simple forcing scheme which in lbmpy is called *implicit forcing*. Ony the force contributions to the first-order (momentum) modes are considered. When the first-order central moments are taken in the frame of reference shifted by $F/2$, they trail behind zero by just that much. Forcing can be applied by first adding half the force *before* collision, and then adding the remaining half *after* collision. Due to this, the first-order central moments entering the cumulant equations are always zero. By applying the force already in central-moment space, the cumulant equations are simplified significantly as no forces need to be taken into account. The pre- and post-collision momentum central moments take the form:
 $$
 \begin{align}
     \kappa_{100} &= - \frac{F_x}{2}, \\
     \kappa^{\ast}_{100} &= \kappa_{100} + F_x = \frac{F_x}{2}
 \end{align}
 $$
 In total, through forcing, the first central moments change sign. This is equivalent to relaxation with a relaxation rate $\omega = 2$. For this reason, lbmpy's implementation of the cumulant LBM calculates the collision of the momentum modes in central moment space, and the default force model overrides their relaxation rate by setting it to 2.
 %% Cell type:markdown id: tags:
 ## 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.
 ### 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`.
 %% Cell type:code id: tags:
 ``` python
 lbm_config = LBMConfig(method=Method.MONOMIAL_CUMULANT, relaxation_rate=sp.Symbol('omega_v'))
 method_monomial = create_lb_method(lbm_config=lbm_config)
 method_monomial
 ```
 %% Output
     <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x127652c40>
 %% Cell type:markdown id: tags:
 ### 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.
 %% Cell type:code id: tags:
 ``` python
 lbm_config = LBMConfig(method=Method.CUMULANT, relaxation_rate=sp.Symbol('omega_v'))
 method_polynomial = create_lb_method(lbm_config=lbm_config)
 method_polynomial
 ```
 %% Output
     <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x12745b5e0>
 %% Cell type:markdown id: tags:
 ### Central Moments and Forcing
 The conserved modes are marked with the note *(central moment)* in the table above. It highlights the fact that these modes are relaxed in central moment space, other than cumulant space. As described in section A, this is done to enable the implicit forcing scheme. When a force is specified, the momentum modes' relaxation rates are overridden by the force model. In the following cell, a symbolic force is specified. Further, a full list of relaxation rates is passed to allow the specification of bulk viscosity.
 %% Cell type:code id: tags:
 ``` python
 lbm_config = LBMConfig(method=Method.CUMULANT, force=sp.symbols('F_:3'),
                        relaxation_rates= [sp.Symbol('omega_shear'),
                                           sp.Symbol('omega_bulk'),
                                           sp.Symbol('omega_3'),
                                           sp.Symbol('omega_4')])
 method_with_force = create_lb_method(lbm_config=lbm_config)
 method_with_force
 ```
 %% Output
     <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x127652520>
 %% Cell type:markdown id: tags:
 ## 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.
 %% Cell type:code id: tags:
 ``` python
 from lbmpy.relaxationrates import relaxation_rate_from_lattice_viscosity
 from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
 ```
 %% Cell type:markdown id: tags:
 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.
 %% Cell type:code id: tags:
 ``` python
 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)
 ```
 %% Cell type:markdown id: tags:
 As a next step, we define the domain size of our set up.
 %% Cell type:code id: tags:
 ``` python
 stencil = LBStencil(Stencil.D2Q9)
 domain_size = (reference_length * 12, reference_length * 4)
 dim = len(domain_size)
 ```
 %% Cell type:markdown id: tags:
 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`
 %% Cell type:code id: tags:
 ``` python
 dh = ps.create_data_handling(domain_size=domain_size, periodicity=(False, False))
 src = dh.add_array('src', values_per_cell=len(stencil), alignment=True)
 dh.fill('src', 0.0, ghost_layers=True)
 dst = dh.add_array('dst', values_per_cell=len(stencil), alignment=True)
 dh.fill('dst', 0.0, ghost_layers=True)
 velField = dh.add_array('velField', values_per_cell=dh.dim, alignment=True)
 dh.fill('velField', 0.0, ghost_layers=True)
 ```
 %% Cell type:markdown id: tags:
 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.
 %% Cell type:code id: tags:
 ``` python
 lbm_config = LBMConfig(stencil=Stencil.D2Q9, method=Method.CUMULANT, relaxation_rate=omega,
                        output={'velocity': velField}, kernel_type='stream_pull_collide')
 method = create_lb_method(lbm_config=lbm_config)
 method
 ```
 %% Output
     <lbmpy.methods.centeredcumulant.centeredcumulantmethod.CenteredCumulantBasedLbMethod at 0x103b6d610>
 %% Cell type:markdown id: tags:
 ### Initialisation with equilibrium
 %% Cell type:code id: tags:
 ``` python
 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)
 ```
 %% Cell type:markdown id: tags:
 ### Definition of the update rule
 %% Cell type:code id: tags:
 ``` python
 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()
 ```
 %% Cell type:markdown id: tags:
 ### Definition of the boundary set up
 %% Cell type:code id: tags:
 ``` python
 def set_sphere(x, y, *_):
     mid = (domain_size[0] // 3, domain_size[1] // 2)
     radius = reference_length // 2
     return (x-mid[0])**2 + (y-mid[1])**2 < radius**2
 ```
 %% Cell type:code id: tags:
 ``` python
 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))
 bh.set_boundary(NoSlip("obstacle"), mask_callback=set_sphere)
 plt.figure(dpi=200)
 plt.boundary_handling(bh)
 ```
 %% Output
 %% Cell type:code id: tags:
 ``` python
 def timeloop(timeSteps):
     for i in range(timeSteps):
         bh()
         dh.run_kernel(kernel)
         dh.swap("src", "dst")
 ```
 %% Cell type:markdown id: tags:
 ### Run the simulation
 %% Cell type:code id: tags:
 ``` python
 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)
         return np.ma.array(dh.gather_array('velField'), mask=mask)
     animation = plt.vector_field_magnitude_animation(run, frames=600, rescale=True)
     set_display_mode('video')
     res = display_animation(animation)
 else:
     timeloop(10)
     res = None
 res
 ```
 %% Output
     <IPython.core.display.HTML object>

doc/notebooks/08_tutorial_shanchen_twophase.ipynb

 %% Cell type:markdown id: tags:
 
 # Shan-Chen Two-Phase Single-Component Lattice Boltzmann
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.session import *
 from lbmpy.updatekernels import create_stream_pull_with_output_kernel
 from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
 from lbmpy.maxwellian_equilibrium import get_weights
 ```
 
 %% Cell type:markdown id: tags:
 
 This is based on section 9.3.2 of Krüger et al.'s "The Lattice Boltzmann Method", Springer 2017 (http://www.lbmbook.com).
 Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp).
 
 %% Cell type:markdown id: tags:
 
 ## Parameters
 
 %% Cell type:code id: tags:
 
 ``` python
 N = 64
 omega_a = 1.
 g_aa = -4.7
 rho0 = 1.
 
 stencil = LBStencil(Stencil.D2Q9)
 weights = get_weights(stencil, c_s_sq=sp.Rational(1,3))
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Data structures
 
 %% Cell type:code id: tags:
 
 ``` python
-dim = len(stencil[0])
+dh = ps.create_data_handling((N,) * stencil.D, periodicity=True, default_target=ps.Target.CPU)
 
-dh = ps.create_data_handling((N,)*dim, periodicity=True, default_target=ps.Target.CPU)
-
-src = dh.add_array('src', values_per_cell=len(stencil))
+src = dh.add_array('src', values_per_cell=stencil.Q)
 dst = dh.add_array_like('dst', 'src')
 
 ρ = dh.add_array('rho')
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Force & combined velocity
 
 %% Cell type:markdown id: tags:
 
 The force on the fluid is
 $\vec{F}_A(\vec{x})=-\psi(\rho_A(\vec{x}))g_{AA}\sum\limits_{i=1}^{q}w_i\psi(\rho_A(\vec{x}+\vec{c}_i))\vec{c}_i$
 with
 $\psi(\rho)=\rho_0\left[1-\exp(-\rho/\rho_0)\right]$.
 
 %% Cell type:code id: tags:
 
 ``` python
 def psi(dens):
     return rho0 * (1. - sp.exp(-dens / rho0));
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
-zero_vec = sp.Matrix([0] * dh.dim)
+zero_vec = sp.Matrix([0] * stencil.D)
 
 force = sum((psi(ρ[d]) * w_d * sp.Matrix(d)
             for d, w_d in zip(stencil, weights)), zero_vec) * psi(ρ.center) * -1 * g_aa
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Kernels
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_config = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
                        force_model=ForceModel.GUO, force=force, kernel_type='collide_only')
 
 collision = create_lb_update_rule(lbm_config=lbm_config,
                                   optimization={'symbolic_field': src})
 
 stream = create_stream_pull_with_output_kernel(collision.method, src, dst, {'density': ρ})
 
 
-opts = {'cpu_openmp': False,
-        'target': dh.default_target}
+config = ps.CreateKernelConfig(target=dh.default_target, cpu_openmp=False)
 
-stream_kernel = ps.create_kernel(stream, **opts).compile()
-collision_kernel = ps.create_kernel(collision, **opts).compile()
+stream_kernel = ps.create_kernel(stream, config=config).compile()
+collision_kernel = ps.create_kernel(collision, config=config).compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:code id: tags:
 
 ``` python
 method_without_force = create_lb_method(LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True))
 init_assignments = macroscopic_values_setter(method_without_force, velocity=(0, 0),
                                              pdfs=src.center_vector, density=ρ.center)
 
 
-init_kernel = ps.create_kernel(init_assignments, ghost_layers=0).compile()
+init_kernel = ps.create_kernel(init_assignments, ghost_layers=0, config=config).compile()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def init():
     for x in range(N):
         for y in range(N):
             if (x-N/2)**2 + (y-N/2)**2 <= 15**2:
                 dh.fill(ρ.name, 2.1, slice_obj=[x,y])
             else:
                 dh.fill(ρ.name, 0.15, slice_obj=[x,y])
 
     dh.run_kernel(init_kernel)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Timeloop
 
 %% Cell type:code id: tags:
 
 ``` python
 sync_pdfs = dh.synchronization_function([src.name])
 sync_ρs = dh.synchronization_function([ρ.name])
 
 def time_loop(steps):
     dh.all_to_gpu()
     for i in range(steps):
         sync_ρs()
         dh.run_kernel(collision_kernel)
 
         sync_pdfs()
         dh.run_kernel(stream_kernel)
 
         dh.swap(src.name, dst.name)
     dh.all_to_cpu()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def plot_ρs():
     plt.figure(dpi=200)
     plt.title("$\\rho$")
     plt.scalar_field(dh.gather_array(ρ.name), vmin=0, vmax=2.5)
     plt.colorbar()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Run the simulation
 ### Initial state
 
 %% Cell type:code id: tags:
 
 ``` python
 init()
 plot_ρs()
 ```
 
 %% Output
 
-
+
 
 %% Cell type:markdown id: tags:
 
-### Check the first time step against reference data
-
-The reference data was obtained with the [sample code](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp) after making the following changes:
-```c++
-const int nsteps = 1000;
-const int noutput = 1;
-```
-
-Remove the next cell if you changed the parameters at the beginning of this notebook.
+### Run the simulation until converged
 
 %% Cell type:code id: tags:
 
 ``` python
 init()
-time_loop(1)
-ref = np.array([0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.136756, 0.220324, 1.2382, 2.26247, 2.26183, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.1, 2.26183, 2.26247, 1.2382, 0.220324, 0.136756, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15, 0.15])
-
-assert np.allclose(dh.gather_array(ρ.name)[N//2], ref)
+time_loop(1000)
+plot_ρs()
 ```
 
-%% Cell type:markdown id: tags:
+%% Output
 
-### Run the simulation until converged
+
 
 %% Cell type:code id: tags:
 
 ``` python
-init()
-time_loop(1000)
-plot_ρs()
+assert np.isfinite(dh.gather_array(ρ.name)).all()
 ```
-
-%% Output
-
-

doc/notebooks/09_tutorial_shanchen_twocomponent.ipynb

 %% Cell type:markdown id: tags:
 
 # Shan-Chen Two-Component Lattice Boltzmann
 
 %% Cell type:code id: tags:
 
 ``` python
 from lbmpy.session import *
 from lbmpy.updatekernels import create_stream_pull_with_output_kernel
 from lbmpy.macroscopic_value_kernels import macroscopic_values_getter, macroscopic_values_setter
 from lbmpy.maxwellian_equilibrium import get_weights
 ```
 
 %% Cell type:markdown id: tags:
 
 This is based on section 9.3.3 of Krüger et al.'s "The Lattice Boltzmann Method", Springer 2017 (http://www.lbmbook.com).
 Sample code is available at [https://github.com/lbm-principles-practice/code/](https://github.com/lbm-principles-practice/code/blob/master/chapter9/shanchen.cpp).
 
 %% Cell type:markdown id: tags:
 
 ## Parameters
 
 %% Cell type:code id: tags:
 
 ``` python
 N = 64       # domain size
 omega_a = 1. # relaxation rate of first component
 omega_b = 1. # relaxation rate of second component
 
 # interaction strength
 g_aa = 0.
 g_ab = g_ba = 6.
 g_bb = 0.
 
 rho0 = 1.
 
 stencil = LBStencil(Stencil.D2Q9)
 weights = get_weights(stencil, c_s_sq=sp.Rational(1, 3))
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Data structures
 
 We allocate two sets of PDF's, one for each phase. Additionally, for each phase there is one field to store its density and velocity.
 
 To run the simulation on GPU, change the `default_target` to gpu
 
 %% Cell type:code id: tags:
 
 ``` python
 dim = stencil.D
 dh = ps.create_data_handling((N, ) * dim, periodicity=True, default_target=ps.Target.CPU)
 
 src_a = dh.add_array('src_a', values_per_cell=stencil.Q)
 dst_a = dh.add_array_like('dst_a', 'src_a')
 
 src_b = dh.add_array('src_b', values_per_cell=stencil.Q)
 dst_b = dh.add_array_like('dst_b', 'src_b')
 
 ρ_a = dh.add_array('rho_a')
 ρ_b = dh.add_array('rho_b')
 u_a = dh.add_array('u_a', values_per_cell=stencil.D)
 u_b = dh.add_array('u_b', values_per_cell=stencil.D)
 u_bary = dh.add_array_like('u_bary', u_a.name)
 
 f_a = dh.add_array('f_a', values_per_cell=stencil.D)
 f_b = dh.add_array_like('f_b', f_a.name)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Force & combined velocity
 
 The two LB methods are coupled using a force term. Its symbolic representation is created in the next cells.
 The force value is not written to a field, but directly evaluated inside the collision kernel.
 
 %% Cell type:markdown id: tags:
 
 The force between the two components is
-$\vec{F}_k(\vec{x})=-\psi(\rho_k(\vec{x}))\sum\limits_{k^\prime\in\{A,B\}}g_{kk^\prime}\sum\limits_{i=1}^{q}w_i\psi(\rho_{k^\prime}(\vec{x}+\vec{c}_i))\vec{c}_i$
+$\mathbf{F}_k(\mathbf{x})=-\psi(\rho_k(\mathbf{x}))\sum\limits_{k^\prime\in\{A,B\}}g_{kk^\prime}\sum\limits_{i=1}^{q}w_i\psi(\rho_{k^\prime}(\mathbf{x}+\mathbf{c}_i))\mathbf{c}_i$
 for $k\in\{A,B\}$
 and with
 $\psi(\rho)=\rho_0\left[1-\exp(-\rho/\rho_0)\right]$.
 
 %% Cell type:code id: tags:
 
 ``` python
 def psi(dens):
     return rho0 * (1. - sp.exp(-dens / rho0));
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 zero_vec = sp.Matrix([0] * dh.dim)
 
 force_a = zero_vec
 for factor, ρ in zip([g_aa, g_ab], [ρ_a, ρ_b]):
     force_a += sum((psi(ρ[d]) * w_d * sp.Matrix(d)
                     for d, w_d in zip(stencil, weights)),
                    zero_vec) * psi(ρ_a.center) * -1 * factor
 
 force_b = zero_vec
 for factor, ρ in zip([g_ba, g_bb], [ρ_a, ρ_b]):
     force_b += sum((psi(ρ[d]) * w_d * sp.Matrix(d)
                     for d, w_d in zip(stencil, weights)),
                    zero_vec) * psi(ρ_b.center) * -1 * factor
 
 f_expressions = ps.AssignmentCollection([
     ps.Assignment(f_a.center_vector, force_a),
     ps.Assignment(f_b.center_vector, force_b)
 ])
 ```
 
 %% Cell type:markdown id: tags:
 
 The barycentric velocity, which is used in place of the individual components' velocities in the equilibrium distribution and Guo force term, is
 $\vec{u}=\frac{1}{\rho_a+\rho_b}\left(\rho_a\vec{u}_a+\frac{1}{2}\vec{F}_a+\rho_b\vec{u}_b+\frac{1}{2}\vec{F}_b\right)$.
 
 %% Cell type:code id: tags:
 
 ``` python
 u_full = ps.Assignment(u_bary.center_vector,
                        (ρ_a.center * u_a.center_vector + ρ_b.center * u_b.center_vector) / (ρ_a.center + ρ_b.center))
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Kernels
 
 %% Cell type:code id: tags:
 
 ``` python
 lbm_config_a = LBMConfig(stencil=stencil, relaxation_rate=omega_a, compressible=True,
                          velocity_input=u_bary, density_input=ρ_a, force_model=ForceModel.GUO,
                          force=f_a, kernel_type='collide_only')
 
 lbm_config_b = LBMConfig(stencil=stencil, relaxation_rate=omega_b, compressible=True,
                          velocity_input=u_bary, density_input=ρ_b, force_model=ForceModel.GUO,
                          force=f_b, kernel_type='collide_only')
 
 
 
 collision_a = create_lb_update_rule(lbm_config=lbm_config_a,
                                     optimization={'symbolic_field': src_a})
 
 collision_b = create_lb_update_rule(lbm_config=lbm_config_b,
                                     optimization={'symbolic_field': src_b})
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 stream_a = create_stream_pull_with_output_kernel(collision_a.method, src_a, dst_a,
                                                  {'density': ρ_a, 'velocity': u_a})
 stream_b = create_stream_pull_with_output_kernel(collision_b.method, src_b, dst_b,
                                                  {'density': ρ_b, 'velocity': u_b})
 
 config = ps.CreateKernelConfig(target=dh.default_target)
 
 stream_a_kernel = ps.create_kernel(stream_a, config=config).compile()
 stream_b_kernel = ps.create_kernel(stream_b, config=config).compile()
 collision_a_kernel = ps.create_kernel(collision_a, config=config).compile()
 collision_b_kernel = ps.create_kernel(collision_b, config=config).compile()
 
 force_kernel = ps.create_kernel(f_expressions, config=config).compile()
 u_kernel = ps.create_kernel(u_full, config=config).compile()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Initialization
 
 %% Cell type:code id: tags:
 
 ``` python
 init_a = macroscopic_values_setter(collision_a.method, velocity=(0, 0),
                                    pdfs=src_a.center_vector, density=ρ_a.center)
 init_b = macroscopic_values_setter(collision_b.method, velocity=(0, 0),
                                    pdfs=src_b.center_vector, density=ρ_b.center)
 init_a_kernel = ps.create_kernel(init_a, ghost_layers=0).compile()
 init_b_kernel = ps.create_kernel(init_b, ghost_layers=0).compile()
 
 
 dh.run_kernel(init_a_kernel)
 dh.run_kernel(init_b_kernel)
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def init():
     dh.fill(ρ_a.name, 0.1, slice_obj=ps.make_slice[:, :0.5])
     dh.fill(ρ_a.name, 0.9, slice_obj=ps.make_slice[:, 0.5:])
 
     dh.fill(ρ_b.name, 0.9, slice_obj=ps.make_slice[:, :0.5])
     dh.fill(ρ_b.name, 0.1, slice_obj=ps.make_slice[:, 0.5:])
 
+    dh.fill(f_a.name, 0.0)
+    dh.fill(f_b.name, 0.0)
+
     dh.run_kernel(init_a_kernel)
     dh.run_kernel(init_b_kernel)
 
     dh.fill(u_a.name, 0.0)
     dh.fill(u_b.name, 0.0)
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Timeloop
 
 %% Cell type:code id: tags:
 
 ``` python
 sync_pdfs = dh.synchronization_function([src_a.name, src_b.name])
 sync_ρs = dh.synchronization_function([ρ_a.name, ρ_b.name])
 
 def time_loop(steps):
     dh.all_to_gpu()
     for i in range(steps):
         sync_ρs()  # force values depend on neighboring ρ's
         dh.run_kernel(force_kernel)
         dh.run_kernel(u_kernel)
         dh.run_kernel(collision_a_kernel)
         dh.run_kernel(collision_b_kernel)
 
         sync_pdfs()
         dh.run_kernel(stream_a_kernel)
         dh.run_kernel(stream_b_kernel)
 
         dh.swap(src_a.name, dst_a.name)
         dh.swap(src_b.name, dst_b.name)
     dh.all_to_cpu()
 ```
 
 %% Cell type:code id: tags:
 
 ``` python
 def plot_ρs():
     plt.figure(dpi=200)
     plt.subplot(1,2,1)
     plt.title("$\\rho_A$")
     plt.scalar_field(dh.gather_array(ρ_a.name), vmin=0, vmax=2)
     plt.colorbar()
     plt.subplot(1,2,2)
     plt.title("$\\rho_B$")
     plt.scalar_field(dh.gather_array(ρ_b.name), vmin=0, vmax=2)
     plt.colorbar()
 ```
 
 %% Cell type:markdown id: tags:
 
 ## Run the simulation
 ### Initial state
 
 %% Cell type:code id: tags:
 
 ``` python
 init()
 plot_ρs()
 ```
 
 %% Output
 
 
 %% Cell type:markdown id: tags:
 
 ### Run the simulation until converged
 
 %% Cell type:code id: tags:
 
 ``` python
 init()
 time_loop(10000)
 plot_ρs()
 ```
 
 %% Output
 
 
 %% Cell type:code id: tags:
 
 ``` python
 assert np.isfinite(dh.gather_array(ρ_a.name)).all()
 assert np.isfinite(dh.gather_array(ρ_b.name)).all()
 ```

lbmpy/advanced_streaming/communication.py

@@ -104,8 +104,7 @@ def get_communication_slices(
 
 
 def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
-                             domain_size=None, target=Target.GPU,
-                             opencl_queue=None, opencl_ctx=None):
+                             domain_size=None, target=Target.GPU):
     """Copies a rectangular array slice onto another non-overlapping array slice"""
     from pystencils.gpucuda.kernelcreation import create_cuda_kernel
 
@@ -136,9 +135,6 @@ def periodic_pdf_copy_kernel(pdf_field, src_slice, dst_slice,
     if target == Target.GPU:
         from pystencils.gpucuda import make_python_function
         return make_python_function(ast)
-    elif target == Target.OPENCL:
-        from pystencils.opencl import make_python_function
-        return make_python_function(ast, opencl_queue, opencl_ctx)
     else:
         raise ValueError('Invalid target:', target)
 
@@ -147,22 +143,17 @@ class LBMPeriodicityHandling:
 
     def __init__(self, stencil, data_handling, pdf_field_name,
                  streaming_pattern='pull', ghost_layers=1,
-                 opencl_queue=None, opencl_ctx=None,
                  pycuda_direct_copy=True):
         """
             Periodicity Handling for Lattice Boltzmann Streaming.
 
-            **On the usage with cuda/opencl:** 
+            **On the usage with cuda:**
             - pycuda allows the copying of sliced arrays within device memory using the numpy syntax,
             e.g. `dst[:,0] = src[:,-1]`. In this implementation, this is the default for periodicity
             handling. Alternatively, if you set `pycuda_direct_copy=False`, GPU kernels are generated and
             compiled. The compiled kernels are almost twice as fast in execution as pycuda array copying,
             but especially for large stencils like D3Q27, their compilation can take up to 20 seconds. 
             Choose your weapon depending on your use case.
-
-            - pyopencl does not support copying of non-contiguous sliced arrays, so the usage of compiled
-            copy kernels is forced on us. On the positive side, compilation of the OpenCL kernels appears
-            to be about four times faster.
         """
         if not isinstance(data_handling, SerialDataHandling):
             raise ValueError('Only serial data handling is supported!')
@@ -172,7 +163,7 @@ class LBMPeriodicityHandling:
         self.dh = data_handling
 
         target = data_handling.default_target
-        assert target in [Target.CPU, Target.GPU, Target.OPENCL]
+        assert target in [Target.CPU, Target.GPU]
 
         self.pdf_field_name = pdf_field_name
         self.ghost_layers = ghost_layers
@@ -180,8 +171,6 @@ class LBMPeriodicityHandling:
         self.inplace_pattern = is_inplace(streaming_pattern)
         self.target = target
         self.cpu = target == Target.CPU
-        self.opencl_queue = opencl_queue
-        self.opencl_ctx = opencl_ctx
         self.pycuda_direct_copy = target == Target.GPU and pycuda_direct_copy
 
         def is_copy_direction(direction):
@@ -205,7 +194,7 @@ class LBMPeriodicityHandling:
                                                            ghost_layers=ghost_layers)
             self.comm_slices.append(list(chain.from_iterable(v for k, v in slices_per_comm_dir.items())))
 
-        if target == Target.OPENCL or (target == Target.GPU and not pycuda_direct_copy):
+        if target == Target.GPU and not pycuda_direct_copy:
             self.device_copy_kernels = []
             for timestep in timesteps:
                 self.device_copy_kernels.append(self._compile_copy_kernels(timestep))
@@ -227,9 +216,7 @@ class LBMPeriodicityHandling:
         kernels = []
         for src, dst in self.comm_slices[timestep.idx]:
             kernels.append(
-                periodic_pdf_copy_kernel(
-                    pdf_field, src, dst, target=self.target,
-                    opencl_queue=self.opencl_queue, opencl_ctx=self.opencl_ctx))
+                periodic_pdf_copy_kernel(pdf_field, src, dst, target=self.target))
         return kernels
 
     def _periodicity_handling_gpu(self, prev_timestep):

lbmpy/boundaries/boundaryconditions.py

@@ -115,8 +115,8 @@ class NoSlip(LbBoundary):
 
 class FreeSlip(LbBoundary):
     """
-    Free-Slip boundary condition, which enforces a zero normal fluid velocity $u_n = 0$ but places no restrictions
-    on the tangential fluid velocity $u_t$.
+    Free-Slip boundary condition, which enforces a zero normal fluid velocity :math:`u_n = 0` but places no restrictions
+    on the tangential fluid velocity :math:`u_t`.
 
     Args:
         stencil: LBM stencil which is used for the simulation
@@ -283,7 +283,7 @@ class UBB(LbBoundary):
         Returns:
             list containing LbmWeightInfo and NeighbourOffsetArrays
         """
-        return [LbmWeightInfo(lb_method), NeighbourOffsetArrays(lb_method.stencil)]
+        return [LbmWeightInfo(lb_method, self.data_type), NeighbourOffsetArrays(lb_method.stencil)]
 
     @property
     def velocity_is_callable(self):
@@ -312,7 +312,8 @@ class UBB(LbBoundary):
             velocity = [eq.rhs for eq in shifted_vel_eqs.new_filtered(cqc.first_order_moment_symbols).main_assignments]
 
         c_s_sq = sp.Rational(1, 3)
-        weight_of_direction = LbmWeightInfo.weight_of_direction
+        weight_info = LbmWeightInfo(lb_method, data_type=self.data_type)
+        weight_of_direction = weight_info.weight_of_direction
         vel_term = 2 / c_s_sq * sum([d_i * v_i for d_i, v_i in zip(neighbor_offset, velocity)]) * weight_of_direction(
             dir_symbol, lb_method)
 
@@ -595,10 +596,11 @@ class DiffusionDirichlet(LbBoundary):
         name: optional name of the boundary.
     """
 
-    def __init__(self, concentration, name=None):
+    def __init__(self, concentration, name=None, data_type='double'):
         if name is None:
             name = "Diffusion Dirichlet " + str(concentration)
         self.concentration = concentration
+        self._data_type = data_type
 
         super(DiffusionDirichlet, self).__init__(name)
 
@@ -611,10 +613,11 @@ class DiffusionDirichlet(LbBoundary):
         Returns:
             list containing LbmWeightInfo
         """
-        return [LbmWeightInfo(lb_method)]
+        return [LbmWeightInfo(lb_method, self._data_type)]
 
     def __call__(self, f_out, f_in, dir_symbol, inv_dir, lb_method, index_field):
-        w_dir = LbmWeightInfo.weight_of_direction(dir_symbol, lb_method)
+        weight_info = LbmWeightInfo(lb_method, self._data_type)
+        w_dir = weight_info.weight_of_direction(dir_symbol, lb_method)
         return [Assignment(f_in(inv_dir[dir_symbol]),
                            2 * w_dir * self.concentration - f_out(dir_symbol))]
 

lbmpy/boundaries/boundaryhandling.py

@@ -38,7 +38,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
         self._prev_timestep = None
 
     def add_fixed_steps(self, fixed_loop, **kwargs):
-        if self._inplace:   # Fixed Loop can't do timestep selection
+        if self._inplace:  # Fixed Loop can't do timestep selection
             raise NotImplementedError("Adding to fixed loop is currently not supported for inplace kernels")
         super(LatticeBoltzmannBoundaryHandling, self).add_fixed_steps(fixed_loop, **kwargs)
 
@@ -52,10 +52,12 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
         if boundary_obj not in self._boundary_object_to_boundary_info:
             sym_index_field = Field.create_generic('indexField', spatial_dimensions=1,
                                                    dtype=numpy_data_type_for_boundary_object(boundary_obj, self.dim))
-            kernels = [self._create_boundary_kernel(
-                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.EVEN).compile(),
-                self._create_boundary_kernel(
-                self._data_handling.fields[self._field_name], sym_index_field, boundary_obj, Timestep.ODD).compile()]
+
+            ast_even = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field,
+                                                    boundary_obj, Timestep.EVEN)
+            ast_odd = self._create_boundary_kernel(self._data_handling.fields[self._field_name], sym_index_field,
+                                                   boundary_obj, Timestep.ODD)
+            kernels = [ast_even.compile(), ast_odd.compile()]
             if flag is None:
                 flag = self.flag_interface.reserve_next_flag()
             boundary_info = self.InplaceStreamingBoundaryInfo(self, boundary_obj, flag, kernels)
@@ -84,6 +86,7 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
             self.boundary_object = boundary_obj
             self.flag = flag
             self._kernels = kernels
+
     #   end class InplaceStreamingBoundaryInfo
 
     # ------------------------------ Force On Boundary ------------------------------------------------------------
@@ -148,29 +151,32 @@ class LatticeBoltzmannBoundaryHandling(BoundaryHandling):
 
         return dh.reduce_float_sequence(list(result), 'sum')
 
+
 # end class LatticeBoltzmannBoundaryHandling
 
 
 class LbmWeightInfo(CustomCodeNode):
+    def __init__(self, lb_method, data_type='double'):
+        self.weights_symbol = TypedSymbol("weights", data_type)
+        data_type_string = "double" if self.weights_symbol.dtype.numpy_dtype == np.float64 else "float"
+
+        weights = [str(w.evalf(17)) for w in lb_method.weights]
+        if data_type_string == "float":
+            weights = "f, ".join(weights)
+            weights += "f"  # suffix for the last element
+        else:
+            weights = ", ".join(weights)
+        w_sym = self.weights_symbol
+        code = f"const {data_type_string} {w_sym.name} [] = {{{ weights }}};\n"
+        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
 
-    # --------------------------- Functions to be used by boundaries --------------------------
-
-    @staticmethod
-    def weight_of_direction(dir_idx, lb_method=None):
+    def weight_of_direction(self, dir_idx, lb_method=None):
         if isinstance(sp.sympify(dir_idx), sp.Integer):
-            return lb_method.weights[dir_idx].evalf()
+            return lb_method.weights[dir_idx].evalf(17)
         else:
-            return sp.IndexedBase(LbmWeightInfo.WEIGHTS_SYMBOL, shape=(1,))[dir_idx]
+            return sp.IndexedBase(self.weights_symbol, shape=(1,))[dir_idx]
 
-    # ---------------------------------- Internal ---------------------------------------------
 
-    WEIGHTS_SYMBOL = TypedSymbol("weights", "double")
-
-    def __init__(self, lb_method):
-        weights = [str(w.evalf()) for w in lb_method.weights]
-        w_sym = LbmWeightInfo.WEIGHTS_SYMBOL
-        code = "const double %s [] = { %s };\n" % (w_sym.name, ",".join(weights))
-        super(LbmWeightInfo, self).__init__(code, symbols_read=set(), symbols_defined={w_sym})
 # end class LbmWeightInfo
 
 

lbmpy/continuous_distribution_measures.py

@@ -15,11 +15,11 @@ def moment_generating_function(generating_function, symbols, symbols_in_result,
     Computes the moment generating function of a probability distribution. It is defined as:
 
     .. math ::
-        F[f(\mathbf{x})](\mathbf{t}) = \int e^{<\mathbf{x}, \mathbf{t}>} f(x)\; dx
+        F[f(\mathbf{x})](t) = \int e^{<\mathbf{x}, t>} f(\mathbf{x})\; dx
 
     Args:
         generating_function: sympy expression
-        symbols: a sequence of symbols forming the vector x
+        symbols: a sequence of symbols forming the vector :math:`\mathbf{x}`
         symbols_in_result: a sequence forming the vector t
         velocity: if the generating function generates central moments, the velocity needs to be substracted. Thus the
                   velocity symbols need to be passed. All generating functions need to have the same parameters.
@@ -62,7 +62,7 @@ def central_moment_generating_function(func, symbols, symbols_in_result, velocit
     Computes central moment generating func, which is defined as:
 
     .. math ::
-        K( \vec{\Xi} ) = \exp ( - \vec{\Xi} \cdot \vec{u} ) M( \vec{\Xi}.
+        K( \mathbf{\Xi} ) = \exp ( - \mathbf{\Xi} \cdot \mathbf{u} ) M( \mathbf{\Xi} ).
 
     For parameter description see :func:`moment_generating_function`.
     """
@@ -76,7 +76,7 @@ def cumulant_generating_function(func, symbols, symbols_in_result, velocity=None
     Computes cumulant generating func, which is the logarithm of the moment generating func:
 
     .. math ::
-        C(\vec{\Xi}) = \log M(\vec{\Xi})
+        C(\mathbf{\Xi}) = \log M(\mathbf{\Xi})
 
     For parameter description see :func:`moment_generating_function`.
     """

lbmpy/creationfunctions.py

@@ -247,21 +247,21 @@ class LBMConfig:
 
     kernel_type: Union[str, Type[PdfFieldAccessor]] = 'default_stream_collide'
     """
-    Supported values: 'default_stream_collide' (default), 'collide_only', 'stream_pull_only'. 
-    With 'default_stream_collide', streaming pattern and even/odd time-step (for in-place patterns) can be specified
+    Supported values: ``'default_stream_collide'`` (default), ``'collide_only'``, ``'stream_pull_only'``. 
+    With ``'default_stream_collide'``, streaming pattern and even/odd time-step (for in-place patterns) can be specified
     by the ``streaming_pattern`` and ``timestep`` arguments. For backwards compatibility, ``kernel_type`` also accepts
-    'stream_pull_collide', 'collide_stream_push', 'esotwist_even', 'esotwist_odd', 'aa_even' and 'aa_odd' for selection 
-    of the streaming pattern. 
+    ``'stream_pull_collide'``, ``'collide_stream_push'``, ``'esotwist_even'``, ``'esotwist_odd'``, ``'aa_even'``
+    and ``'aa_odd'`` for selection of the streaming pattern. 
     """
     streaming_pattern: str = 'pull'
     """
-    The streaming pattern to be used with a 'default_stream_collide' kernel. Accepted values are
-    'pull', 'push', 'aa' and 'esotwist'.
+    The streaming pattern to be used with a ``'default_stream_collide'`` kernel. Accepted values are
+    ``'pull'``, ``'push'``, ``'aa'`` and ``'esotwist'``.
     """
     timestep: Timestep = Timestep.BOTH
     """
     Timestep modulus for the streaming pattern. For two-fields patterns, this argument is irrelevant and
-    by default set to ``Timestep.BOTH``. For in-place patterns, ``Timestep.EVEN`` or ``Timestep.ODD`` must be speficied.
+    by default set to ``Timestep.BOTH``. For in-place patterns, ``Timestep.EVEN`` or ``Timestep.ODD`` must be specified.
     """
 
     field_name: str = 'src'

lbmpy/enums.py

@@ -21,6 +21,10 @@ class Stencil(Enum):
     """
     A two dimensional stencil using 37 discrete velocities. (long range stencil).
     """
+    D3Q7 = auto()
+    """
+    A three dimensional stencil using 7 discrete velocities.
+    """
     D3Q15 = auto()
     """
     A three dimensional stencil using 15 discrete velocities.

lbmpy/forcemodels.py

@@ -17,16 +17,16 @@ Force models add a term :math:`C_F` to the collision equation:
 
 .. math ::
 
-    f(\pmb{x} + c_q \Delta t, t + \Delta t) - f(\pmb{x},t) = \Omega(f, f^{(\mathrm{eq})})
+    f(\mathbf{x} + c_q \Delta t, t + \Delta t) - f(\mathbf{x},t) = \Omega(f, f^{(\mathrm{eq})})
                                                             + \underbrace{F_q}_{\mbox{forcing term}}
 
 The form of this term depends on the concrete force model: the first moment of this forcing term is equal
-to the acceleration :math:`\pmb{a}` for all force models. Here :math:`\mathbf{F}` is the D dimensional force vector,
+to the acceleration :math:`\mathbf{a}` for all force models. Here :math:`\mathbf{F}` is the D dimensional force vector,
 which defines the force for each spatial dircetion.
 
 .. math ::
 
-    \sum_q \pmb{c}_q \mathbf{F} = \pmb{a}
+    \sum_q \mathbf{c}_q \mathbf{F} = \mathbf{a}
 
 
 The second order moment is different for the forcing models - if it is zero the model is suited for
@@ -57,7 +57,7 @@ For all force models the computation of the macroscopic velocity has to be adapt
     
     .. math ::
     
-        \pmb{u} &= \sum_q \pmb{c}_q f_q + S_{\mathrm{macro}}
+        \mathbf{u} &= \sum_q \mathbf{c}_q f_q + S_{\mathrm{macro}}
         
         S_{\mathrm{macro}} &= \frac{\Delta t}{2 \cdot \rho} \sum_q F_q
 
@@ -296,7 +296,7 @@ class He(AbstractForceModel):
                 F_x m^{\mathrm{eq}}_{\alpha+1,\beta,\gamma} 
                 + F_y m^{\mathrm{eq}}_{\alpha,\beta+1,\gamma} 
                 + F_z m^{\mathrm{eq}}_{\alpha,\beta,\gamma+1} 
-                - m^{eq}_{\alpha\beta\gamma} ( \mathbf{F} \cdot \vec{u} ) 
+                - m^{eq}_{\alpha\beta\gamma} ( \mathbf{F} \cdot \mathbf{u} )
             \right)
     """
 

lbmpy/lbstep.py

@@ -74,7 +74,7 @@ class LatticeBoltzmannStep:
         self.density_data_name = name + "_density" if density_data_name is None else density_data_name
         self.density_data_index = density_data_index
 
-        self._gpu = target == Target.GPU or target == Target.OPENCL
+        self._gpu = target == Target.GPU
         layout = lbm_optimisation.field_layout
 
         alignment = False

lbmpy/macroscopic_value_kernels.py

@@ -8,12 +8,16 @@ from lbmpy.advanced_streaming.utility import get_accessor, Timestep
 
 
 def pdf_initialization_assignments(lb_method, density, velocity, pdfs,
-                                   streaming_pattern='pull', previous_timestep=Timestep.BOTH):
+                                   streaming_pattern='pull', previous_timestep=Timestep.BOTH,
+                                   set_pre_collision_pdfs=False):
     """Assignments to initialize the pdf field with equilibrium"""
     if isinstance(pdfs, Field):
-        previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
-        field_accesses = previous_step_accessor.write(pdfs, lb_method.stencil)
-    elif streaming_pattern == 'pull':
+        accessor = get_accessor(streaming_pattern, previous_timestep)
+        if set_pre_collision_pdfs:
+            field_accesses = accessor.read(pdfs, lb_method.stencil)
+        else:
+            field_accesses = accessor.write(pdfs, lb_method.stencil)
+    elif streaming_pattern == 'pull' and not set_pre_collision_pdfs:
         field_accesses = pdfs
     else:
         raise ValueError("Invalid value of pdfs: A PDF field reference is required to derive "
@@ -28,11 +32,15 @@ def pdf_initialization_assignments(lb_method, density, velocity, pdfs,
 
 
 def macroscopic_values_getter(lb_method, density, velocity, pdfs,
-                              streaming_pattern='pull', previous_timestep=Timestep.BOTH):
+                              streaming_pattern='pull', previous_timestep=Timestep.BOTH,
+                              use_pre_collision_pdfs=False):
     if isinstance(pdfs, Field):
-        previous_step_accessor = get_accessor(streaming_pattern, previous_timestep)
-        field_accesses = previous_step_accessor.write(pdfs, lb_method.stencil)
-    elif streaming_pattern == 'pull':
+        accessor = get_accessor(streaming_pattern, previous_timestep)
+        if use_pre_collision_pdfs:
+            field_accesses = accessor.read(pdfs, lb_method.stencil)
+        else:
+            field_accesses = accessor.write(pdfs, lb_method.stencil)
+    elif streaming_pattern == 'pull' and not use_pre_collision_pdfs:
         field_accesses = pdfs
     else:
         raise ValueError("Invalid value of pdfs: A PDF field reference is required to derive "

lbmpy/maxwellian_equilibrium.py

@@ -31,6 +31,10 @@ get_weights.weights = {
         1: sp.Rational(1, 9),
         2: sp.Rational(1, 36),
     },
+    7: {
+        0: sp.simplify(0.0),
+        1: sp.Rational(1, 6),
+    },
     15: {
         0: sp.Rational(2, 9),
         1: sp.Rational(1, 9),

lbmpy/methods/abstractlbmethod.py

@@ -2,6 +2,7 @@ import abc
 from collections import namedtuple
 
 import sympy as sp
+from sympy.core.numbers import Zero
 
 from pystencils import Assignment, AssignmentCollection
 
@@ -52,6 +53,7 @@ class AbstractLbMethod(abc.ABC):
         """Returns a qxq diagonal matrix which contains the relaxation rate for each moment on the diagonal"""
         d = sp.zeros(len(self.relaxation_rates))
         for i in range(0, len(self.relaxation_rates)):
+            # note that 0.0 is converted to sp.Zero here. It is not possible to prevent this.
             d[i, i] = self.relaxation_rates[i]
         return d
 
@@ -104,6 +106,9 @@ class AbstractLbMethod(abc.ABC):
         for relaxation_rate in rr:
             if relaxation_rate not in unique_relaxation_rates:
                 relaxation_rate = sp.sympify(relaxation_rate)
+                # special treatment for zero, sp.Zero would be an integer ..
+                if isinstance(relaxation_rate, Zero):
+                    relaxation_rate = 0.0
                 if not isinstance(relaxation_rate, sp.Symbol):
                     rt_symbol = sp.Symbol(f"rr_{len(subexpressions)}")
                     subexpressions[relaxation_rate] = rt_symbol

lbmpy/methods/conservedquantitycomputation.py

@@ -286,7 +286,7 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs, symb
 
     dim = stencil.D
 
-    subexpressions = []
+    subexpressions = list()
     pdf_sum = sum(symbolic_pdfs)
     u = [0] * dim
     for f, offset in zip(symbolic_pdfs, stencil):
@@ -308,12 +308,16 @@ def get_equations_for_zeroth_and_first_order_moment(stencil, symbolic_pdfs, symb
             rhs -= plus_terms[j]
         eq = Assignment(velo_terms[i], sum(rhs))
         subexpressions.append(eq)
+        if len(rhs) == 0:  # if one of the substitutions is not found the simplification can not be applied
+            subexpressions = []
+            break
 
     for subexpression in subexpressions:
         pdf_sum = pdf_sum.subs(subexpression.rhs, subexpression.lhs)
 
-    for i in range(dim):
-        u[i] = u[i].subs(subexpressions[i].rhs, subexpressions[i].lhs)
+    if len(subexpressions) > 0:
+        for i in range(dim):
+            u[i] = u[i].subs(subexpressions[i].rhs, subexpressions[i].lhs)
 
     equations = []
     equations += [Assignment(symbolic_zeroth_moment, pdf_sum)]

lbmpy/methods/creationfunctions.py

@@ -423,8 +423,10 @@ def create_mrt_orthogonal(stencil, relaxation_rates, maxwellian_moments=False, w
             diagonal_viscous_moments = [x ** 2 + y ** 2, x ** 2]
         else:
             diagonal_viscous_moments = [x ** 2 + y ** 2 + z ** 2, x ** 2, y ** 2 - z ** 2]
+
         for i, d in enumerate(MOMENT_SYMBOLS[:stencil.D]):
-            moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
+            if d ** 2 in moments:
+                moments[moments.index(d ** 2)] = diagonal_viscous_moments[i]
         orthogonal_moments = gram_schmidt(moments, stencil, weights)
         orthogonal_moments_scaled = [e * common_denominator(e) for e in orthogonal_moments]
         nested_moments = list(sort_moments_into_groups_of_same_order(orthogonal_moments_scaled).values())
@@ -603,11 +605,11 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
         for group in nested_moments:
             for moment in group:
                 if get_order(moment) <= 1:
-                    result[moment] = 0
+                    result[moment] = 0.0
                 elif is_shear_moment(moment, dim):
                     result[moment] = relaxation_rates[0]
                 else:
-                    result[moment] = 1
+                    result[moment] = 1.0
 
     # if relaxation rate for each moment is specified they are all used
     if len(relaxation_rates) == number_of_moments:
@@ -632,7 +634,7 @@ def _get_relaxation_info_dict(relaxation_rates, nested_moments, dim):
                 next_rr = False
                 for moment in group:
                     if get_order(moment) <= 1:
-                        result[moment] = 0
+                        result[moment] = 0.0
                     elif is_shear_moment(moment, dim):
                         result[moment] = shear_rr
                     elif is_bulk_moment(moment, dim):

lbmpy/methods/default_moment_sets.py

@@ -18,7 +18,7 @@ def cascaded_moment_sets_literature(stencil):
     - Remaining groups do not govern hydrodynamic properties
 
     Args:
-        stencil: instance of :class:`lbmpy.stencils.LBStencil`. Can be D2Q9, D3Q15, D3Q19 or D3Q27
+        stencil: instance of :class:`lbmpy.stencils.LBStencil`. Can be D2Q9, D3Q7, D3Q15, D3Q19 or D3Q27
     """
     x, y, z = MOMENT_SYMBOLS
     if have_same_entries(stencil, LBStencil(Stencil.D2Q9)):
@@ -42,6 +42,18 @@ def cascaded_moment_sets_literature(stencil):
             [x ** 2 * y ** 2]
         ]
 
+    elif have_same_entries(stencil, LBStencil(Stencil.D3Q7)):
+        # D3Q7 moments: https://arxiv.org/ftp/arxiv/papers/1611/1611.03329.pdf
+        return [
+            [sp.sympify(1)],  # density is conserved
+            [x, y, z],  # momentum might be affected by forcing
+
+            [x ** 2 - y ** 2,
+             x ** 2 - z ** 2],  # shear
+
+            [x ** 2 + y ** 2 + z ** 2],  # bulk
+        ]
+
     elif have_same_entries(stencil, LBStencil(Stencil.D3Q15)):
         #   D3Q15 central moments by Premnath et al. https://arxiv.org/pdf/1202.6081.pdf.
         return [