Compare revisions

107c35da · 107c35da · 107c35da · 107c35da · 107c35da · 107c35da
--- a/lbmpy_tests/phasefield/test_analytical_n_phase.ipynb
+++ b/lbmpy_tests/phasefield/test_analytical_n_phase.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from functools import partial
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for N-Phase model
-
-%% Cell type:markdown id: tags:
-
-### Formulation of free energy
-
-%% Cell type:markdown id: tags:
-
-In the next cell you can inspect the formulation of the free energy functional.
-Bulk and interface terms can be viewed separately.
-
-%% Cell type:code id: tags:
-
-``` python
-num_phases = 3
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-f2 = partial(n_phases_correction_function, beta=1)
-F = free_energy_functional_n_phases(order_parameters=order_params,
-                                    include_interface=False, f2=f2,
-                                    include_bulk=True )
-F
-```
-
-%% Output
-
-
-    $$\frac{3}{2 \alpha} \left(\tau_{0 1} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{1} + 1\right)^{2} & \text{for}\: - \phi_{1} + 1 < 0 \\- \phi_{1}^{2} & \text{for}\: - \phi_{1} + 1 > 1 \\\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{0} + 1\right)^{2} & \text{for}\: - \phi_{0} + 1 < 0 \\- \phi_{0}^{2} & \text{for}\: - \phi_{0} + 1 > 1 \\\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} & \text{otherwise} \end{cases}\right)\right)$$
-      ⎛     ⎛                                  ⎛⎧                 2
-      ⎜     ⎜                                  ⎜⎪       -(φ₀ + φ₁)           for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜  2          2     2          2   ⎜⎪                    2
-    3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(-φ₀ + 1)  + φ₁ ⋅(-φ₁ + 1)  - ⎜⎨     -(-φ₀ - φ₁ + 1)        for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜                                  ⎜⎪         2               2
-      ⎜     ⎜                                  ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)      ot
-      ⎝     ⎝                                  ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-              ⎞⎞        ⎛                                              ⎛⎧
-    ₀ + φ₁ < 0⎟⎟        ⎜                                              ⎜⎪ -(-φ₁ +
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜  2          2            2               2   ⎜⎪        2
-    ₀ + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(-φ₀ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  - ⎜⎨     -φ₁
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜                                              ⎜⎪  2
-    herwise   ⎟⎟        ⎜                                              ⎜⎪φ₁ ⋅(-φ₁
-              ⎠⎠        ⎝                                              ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-                                                            2⋅α
-    
-      2                   ⎞⎞        ⎛
-    1)     for -φ₁ + 1 < 0⎟⎟        ⎜
-                          ⎟⎟        ⎜
-                          ⎟⎟        ⎜  2          2            2               2
-           for -φ₁ + 1 > 1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(-φ₁ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  -
-                          ⎟⎟        ⎜
-        2                 ⎟⎟        ⎜
-    + 1)      otherwise   ⎟⎟        ⎜
-                          ⎠⎠        ⎝
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-     ⎛⎧           2                   ⎞⎞⎞
-     ⎜⎪ -(-φ₀ + 1)     for -φ₀ + 1 < 0⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪        2                      ⎟⎟⎟
-     ⎜⎨     -φ₀        for -φ₀ + 1 > 1⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪  2          2                 ⎟⎟⎟
-     ⎜⎪φ₀ ⋅(-φ₀ + 1)      otherwise   ⎟⎟⎟
-     ⎝⎩                               ⎠⎠⎠
-    ─────────────────────────────────────
-    
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-First we define the order parameters and free energy, including bulk and interface terms:
-
-%% Cell type:code id: tags:
-
-``` python
-F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
-```
-
-%% Cell type:markdown id: tags:
-
-Then we automatically derive the differential equations for the chemial potential $\mu$
-
-%% Cell type:code id: tags:
-
-``` python
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
-
-# there is one equation less than phases
-assert len(mu_diff_eq) == num_phases - 1
-
-# show the first one
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3}}{\alpha} \tau_{0 2} - \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 1} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 2} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{1 2} - \frac{18 \phi_{0}^{2}}{\alpha} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 1} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{1 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 1} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{1 2} + \frac{6 \phi_{0}}{\alpha} \tau_{0 2} - \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 1} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 2} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{1 2} + \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 1} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 2} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{1 2} - \frac{3 \phi_{1}}{\alpha} \tau_{0 1} + \frac{3 \phi_{1}}{\alpha} \tau_{0 2} + \frac{3 \phi_{1}}{\alpha} \tau_{1 2}$$
-    
-    
-    3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
-    
-    
-                       3             2                2                2
-                  12⋅φ₀ ⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₀ ₁   18⋅φ₀ ⋅φ₁⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₁ ₂
-    D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
-                       α              α                α                α
-    
-          2                2                2                2
-     18⋅φ₀ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₀ ₁   18⋅φ₀⋅φ₁ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₁ ₂   18⋅φ₀⋅φ₁⋅τ₀
-     ─────────── - ────────────── + ────────────── + ────────────── + ────────────
-          α              α                α                α                α
-    
-                                                        3            3
-    ₁   18⋅φ₀⋅φ₁⋅τ₀ ₂   18⋅φ₀⋅φ₁⋅τ₁ ₂   6⋅φ₀⋅τ₀ ₂   6⋅φ₁ ⋅τ₀ ₁   6⋅φ₁ ⋅τ₀ ₂   6⋅φ₁
-    ─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
-              α               α             α           α            α
-    
-    3            2            2            2
-     ⋅τ₁ ₂   9⋅φ₁ ⋅τ₀ ₁   9⋅φ₁ ⋅τ₀ ₂   9⋅φ₁ ⋅τ₁ ₂   3⋅φ₁⋅τ₀ ₁   3⋅φ₁⋅τ₀ ₂   3⋅φ₁⋅τ
-    ────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
-    α            α            α            α            α           α           α
-    
-    
-    ₁ ₂
-    ───
-    
-
-%% Cell type:markdown id: tags:
-
-Next we check, that the $tanh$ is indeed a solution of this equation...
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expected_profile = analytic_interface_profile(x)
-expected_profile
-```
-
-%% Output
-
-
-    $$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
-        ⎛ x ⎞
-    tanh⎜───⎟
-        ⎝2⋅α⎠   1
-    ───────── + ─
-        2       2
-
-%% Cell type:markdown id: tags:
-
-... by inserting it for $\phi_0$, and setting all other order parameters to zero
-
-%% Cell type:code id: tags:
-
-``` python
-# zero other parameters
-diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
-
-# insert analytical solution
-diff_eq = diff_eq.subs(order_params[0], expected_profile)
-diff_eq.factor()
-```
-
-%% Output
-
-    $$- \frac{3 \tau_{0 2}}{2 \alpha} \left(4 \alpha^{2} {\partial {\partial (\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}) }} - \tanh^{3}{\left (\frac{x}{2 \alpha} \right )} + \tanh{\left (\frac{x}{2 \alpha} \right )}\right)$$
-            ⎛   2                                       3⎛ x ⎞       ⎛ x ⎞⎞
-    -3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
-            ⎝                                            ⎝2⋅α⎠       ⎝2⋅α⎠⎠
-    ────────────────────────────────────────────────────────────────────────
-                                      2⋅α
-
-%% Cell type:markdown id: tags:
-
-finally the differentials have to be evaluated...
-
-%% Cell type:code id: tags:
-
-``` python
-from pystencils.fd import evaluate_diffs
-diff_eq = evaluate_diffs(diff_eq, x)
-diff_eq
-```
-
-%% Output
-
-
-    $$\frac{12}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{3} - \frac{18}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{2} + \frac{6}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right) + \frac{3 \tau_{0 2}}{2 \alpha} \left(- \tanh^{2}{\left (\frac{x}{2 \alpha} \right )} + 1\right) \tanh{\left (\frac{x}{2 \alpha} \right )}$$
-                           3                          2
-            ⎛    ⎛ x ⎞    ⎞            ⎛    ⎛ x ⎞    ⎞           ⎛    ⎛ x ⎞    ⎞
-            ⎜tanh⎜───⎟    ⎟            ⎜tanh⎜───⎟    ⎟           ⎜tanh⎜───⎟    ⎟
-            ⎜    ⎝2⋅α⎠   1⎟            ⎜    ⎝2⋅α⎠   1⎟           ⎜    ⎝2⋅α⎠   1⎟
-    12⋅τ₀ ₂⋅⎜───────── + ─⎟    18⋅τ₀ ₂⋅⎜───────── + ─⎟    6⋅τ₀ ₂⋅⎜───────── + ─⎟
-            ⎝    2       2⎠            ⎝    2       2⎠           ⎝    2       2⎠
-    ──────────────────────── - ──────────────────────── + ────────────────────── +
-               α                          α                         α
-    
-    
-    
-    
-            ⎛      2⎛ x ⎞    ⎞     ⎛ x ⎞
-     3⋅τ₀ ₂⋅⎜- tanh ⎜───⎟ + 1⎟⋅tanh⎜───⎟
-            ⎝       ⎝2⋅α⎠    ⎠     ⎝2⋅α⎠
-     ───────────────────────────────────
-                     2⋅α
-
-%% Cell type:markdown id: tags:
-
-.. and the result simplified...
-
-%% Cell type:code id: tags:
-
-``` python
-assert diff_eq.expand() == 0
-diff_eq.expand()
-```
-
-%% Output
-
-
-    $$0$$
-    0
-
-%% Cell type:markdown id: tags:
-
-...and indeed the expected tanh profile satisfies this differential equation.
-
-Next lets check for the interface between phase 0 and phase 1:
-
-%% Cell type:code id: tags:
-
-``` python
-for diff_eq in mu_diff_eq:
-    eq = diff_eq.subs({order_params[0]: expected_profile,
-                       order_params[1]: 1 - expected_profile})
-    assert evaluate_diffs(eq, x).expand() == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-Computing the excess free energy per unit area of an interface between two phases.
-This should be exactly the surface tension parameter.
-
-%% Cell type:code id: tags:
-
-``` python
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-F = free_energy_functional_n_phases(order_parameters=order_params)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-two_phase_free_energy = F.subs({order_params[0]: expected_profile,
-                                order_params[1]: 1 - expected_profile})
-
-# Evaluate differentials and simplification
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-excess_free_energy = sp.Integral(two_phase_free_energy, x)
-excess_free_energy
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}\, dx$$
-    ⌠
-    ⎮     3⋅τ₀ ₁
-    ⎮ ────────────── dx
-    ⎮         4⎛ x ⎞
-    ⎮ 8⋅α⋅cosh ⎜───⎟
-    ⎮          ⎝2⋅α⎠
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
-
-%% Cell type:code id: tags:
-
-``` python
-coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
-transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
-transformed_int
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left (u \right )}}\, du$$
-    ⌠
-    ⎮   3⋅τ₀ ₁
-    ⎮ ────────── du
-    ⎮       4
-    ⎮ 4⋅cosh (u)
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Now the integral can be done:
-
-%% Cell type:code id: tags:
-
-``` python
-result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
-assert result == symmetric_symbolic_surface_tension(0,1)
-result
-```
-
-%% Output
-
-
-    $$\tau_{0 1}$$
-    τ₀ ₁
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from functools import partial
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for N-Phase model
-
-%% Cell type:markdown id: tags:
-
-### Formulation of free energy
-
-%% Cell type:markdown id: tags:
-
-In the next cell you can inspect the formulation of the free energy functional.
-Bulk and interface terms can be viewed separately.
-
-%% Cell type:code id: tags:
-
-``` python
-num_phases = 3
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-f2 = partial(n_phases_correction_function, beta=1)
-F = free_energy_functional_n_phases(order_parameters=order_params,
-                                    include_interface=False, f2=f2,
-                                    include_bulk=True )
-F
-```
-
-%% Output
-
-
-    $$\frac{3}{2 \alpha} \left(\tau_{0 1} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} - \begin{cases} - \left(\phi_{0} + \phi_{1}\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} < 0 \\- \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{for}\: \phi_{0} + \phi_{1} > 1 \\\left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{0 2} \left(\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{1} + 1\right)^{2} & \text{for}\: - \phi_{1} + 1 < 0 \\- \phi_{1}^{2} & \text{for}\: - \phi_{1} + 1 > 1 \\\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} & \text{otherwise} \end{cases}\right) + \tau_{1 2} \left(\phi_{1}^{2} \left(- \phi_{1} + 1\right)^{2} + \left(\phi_{0} + \phi_{1}\right)^{2} \left(- \phi_{0} - \phi_{1} + 1\right)^{2} - \begin{cases} - \left(- \phi_{0} + 1\right)^{2} & \text{for}\: - \phi_{0} + 1 < 0 \\- \phi_{0}^{2} & \text{for}\: - \phi_{0} + 1 > 1 \\\phi_{0}^{2} \left(- \phi_{0} + 1\right)^{2} & \text{otherwise} \end{cases}\right)\right)$$
-      ⎛     ⎛                                  ⎛⎧                 2
-      ⎜     ⎜                                  ⎜⎪       -(φ₀ + φ₁)           for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜  2          2     2          2   ⎜⎪                    2
-    3⋅⎜τ₀ ₁⋅⎜φ₀ ⋅(-φ₀ + 1)  + φ₁ ⋅(-φ₁ + 1)  - ⎜⎨     -(-φ₀ - φ₁ + 1)        for φ
-      ⎜     ⎜                                  ⎜⎪
-      ⎜     ⎜                                  ⎜⎪         2               2
-      ⎜     ⎜                                  ⎜⎪(φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)      ot
-      ⎝     ⎝                                  ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-              ⎞⎞        ⎛                                              ⎛⎧
-    ₀ + φ₁ < 0⎟⎟        ⎜                                              ⎜⎪ -(-φ₁ +
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜  2          2            2               2   ⎜⎪        2
-    ₀ + φ₁ > 1⎟⎟ + τ₀ ₂⋅⎜φ₀ ⋅(-φ₀ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  - ⎜⎨     -φ₁
-              ⎟⎟        ⎜                                              ⎜⎪
-              ⎟⎟        ⎜                                              ⎜⎪  2
-    herwise   ⎟⎟        ⎜                                              ⎜⎪φ₁ ⋅(-φ₁
-              ⎠⎠        ⎝                                              ⎝⎩
-    ──────────────────────────────────────────────────────────────────────────────
-                                                            2⋅α
-    
-      2                   ⎞⎞        ⎛
-    1)     for -φ₁ + 1 < 0⎟⎟        ⎜
-                          ⎟⎟        ⎜
-                          ⎟⎟        ⎜  2          2            2               2
-           for -φ₁ + 1 > 1⎟⎟ + τ₁ ₂⋅⎜φ₁ ⋅(-φ₁ + 1)  + (φ₀ + φ₁) ⋅(-φ₀ - φ₁ + 1)  -
-                          ⎟⎟        ⎜
-        2                 ⎟⎟        ⎜
-    + 1)      otherwise   ⎟⎟        ⎜
-                          ⎠⎠        ⎝
-    ──────────────────────────────────────────────────────────────────────────────
-    
-    
-     ⎛⎧           2                   ⎞⎞⎞
-     ⎜⎪ -(-φ₀ + 1)     for -φ₀ + 1 < 0⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪        2                      ⎟⎟⎟
-     ⎜⎨     -φ₀        for -φ₀ + 1 > 1⎟⎟⎟
-     ⎜⎪                               ⎟⎟⎟
-     ⎜⎪  2          2                 ⎟⎟⎟
-     ⎜⎪φ₀ ⋅(-φ₀ + 1)      otherwise   ⎟⎟⎟
-     ⎝⎩                               ⎠⎠⎠
-    ─────────────────────────────────────
-    
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-First we define the order parameters and free energy, including bulk and interface terms:
-
-%% Cell type:code id: tags:
-
-``` python
-F = free_energy_functional_n_phases(order_parameters=symbolic_order_parameters(num_symbols=num_phases-1))
-```
-
-%% Cell type:markdown id: tags:
-
-Then we automatically derive the differential equations for the chemial potential $\mu$
-
-%% Cell type:code id: tags:
-
-``` python
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_params)
-
-# there is one equation less than phases
-assert len(mu_diff_eq) == num_phases - 1
-
-# show the first one
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$3 \alpha \tau_{0 1} {\partial {\partial \phi_{1}}} - 6 \alpha \tau_{0 2} {\partial {\partial \phi_{0}}} - 3 \alpha \tau_{0 2} {\partial {\partial \phi_{1}}} - 3 \alpha \tau_{1 2} {\partial {\partial \phi_{1}}} + \frac{12 \phi_{0}^{3}}{\alpha} \tau_{0 2} - \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 1} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{0 2} + \frac{18 \phi_{1}}{\alpha} \phi_{0}^{2} \tau_{1 2} - \frac{18 \phi_{0}^{2}}{\alpha} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 1} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{0 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1}^{2} \tau_{1 2} + \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 1} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{0 2} - \frac{18 \phi_{0}}{\alpha} \phi_{1} \tau_{1 2} + \frac{6 \phi_{0}}{\alpha} \tau_{0 2} - \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 1} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{0 2} + \frac{6 \phi_{1}^{3}}{\alpha} \tau_{1 2} + \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 1} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{0 2} - \frac{9 \phi_{1}^{2}}{\alpha} \tau_{1 2} - \frac{3 \phi_{1}}{\alpha} \tau_{0 1} + \frac{3 \phi_{1}}{\alpha} \tau_{0 2} + \frac{3 \phi_{1}}{\alpha} \tau_{1 2}$$
-    
-    
-    3⋅α⋅τ₀ ₁⋅D(D(phi_1)) - 6⋅α⋅τ₀ ₂⋅D(D(phi_0)) - 3⋅α⋅τ₀ ₂⋅D(D(phi_1)) - 3⋅α⋅τ₁ ₂⋅
-    
-    
-                       3             2                2                2
-                  12⋅φ₀ ⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₀ ₁   18⋅φ₀ ⋅φ₁⋅τ₀ ₂   18⋅φ₀ ⋅φ₁⋅τ₁ ₂
-    D(D(phi_1)) + ─────────── - ────────────── + ────────────── + ────────────── -
-                       α              α                α                α
-    
-          2                2                2                2
-     18⋅φ₀ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₀ ₁   18⋅φ₀⋅φ₁ ⋅τ₀ ₂   18⋅φ₀⋅φ₁ ⋅τ₁ ₂   18⋅φ₀⋅φ₁⋅τ₀
-     ─────────── - ────────────── + ────────────── + ────────────── + ────────────
-          α              α                α                α                α
-    
-                                                        3            3
-    ₁   18⋅φ₀⋅φ₁⋅τ₀ ₂   18⋅φ₀⋅φ₁⋅τ₁ ₂   6⋅φ₀⋅τ₀ ₂   6⋅φ₁ ⋅τ₀ ₁   6⋅φ₁ ⋅τ₀ ₂   6⋅φ₁
-    ─ - ───────────── - ───────────── + ───────── - ────────── + ────────── + ────
-              α               α             α           α            α
-    
-    3            2            2            2
-     ⋅τ₁ ₂   9⋅φ₁ ⋅τ₀ ₁   9⋅φ₁ ⋅τ₀ ₂   9⋅φ₁ ⋅τ₁ ₂   3⋅φ₁⋅τ₀ ₁   3⋅φ₁⋅τ₀ ₂   3⋅φ₁⋅τ
-    ────── + ────────── - ────────── - ────────── - ───────── + ───────── + ──────
-    α            α            α            α            α           α           α
-    
-    
-    ₁ ₂
-    ───
-    
-
-%% Cell type:markdown id: tags:
-
-Next we check, that the $tanh$ is indeed a solution of this equation...
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expected_profile = analytic_interface_profile(x)
-expected_profile
-```
-
-%% Output
-
-
-    $$\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}$$
-        ⎛ x ⎞
-    tanh⎜───⎟
-        ⎝2⋅α⎠   1
-    ───────── + ─
-        2       2
-
-%% Cell type:markdown id: tags:
-
-... by inserting it for $\phi_0$, and setting all other order parameters to zero
-
-%% Cell type:code id: tags:
-
-``` python
-# zero other parameters
-diff_eq = mu_diff_eq[0].subs({p: 0 for p in order_params[1:]})
-
-# insert analytical solution
-diff_eq = diff_eq.subs(order_params[0], expected_profile)
-diff_eq.factor()
-```
-
-%% Output
-
-    $$- \frac{3 \tau_{0 2}}{2 \alpha} \left(4 \alpha^{2} {\partial {\partial (\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}) }} - \tanh^{3}{\left (\frac{x}{2 \alpha} \right )} + \tanh{\left (\frac{x}{2 \alpha} \right )}\right)$$
-            ⎛   2                                       3⎛ x ⎞       ⎛ x ⎞⎞
-    -3⋅τ₀ ₂⋅⎜4⋅α ⋅D(D(tanh(x/(2*alpha))/2 + 1/2)) - tanh ⎜───⎟ + tanh⎜───⎟⎟
-            ⎝                                            ⎝2⋅α⎠       ⎝2⋅α⎠⎠
-    ────────────────────────────────────────────────────────────────────────
-                                      2⋅α
-
-%% Cell type:markdown id: tags:
-
-finally the differentials have to be evaluated...
-
-%% Cell type:code id: tags:
-
-``` python
-from pystencils.fd import evaluate_diffs
-diff_eq = evaluate_diffs(diff_eq, x)
-diff_eq
-```
-
-%% Output
-
-
-    $$\frac{12}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{3} - \frac{18}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right)^{2} + \frac{6}{\alpha} \tau_{0 2} \left(\frac{1}{2} \tanh{\left (\frac{x}{2 \alpha} \right )} + \frac{1}{2}\right) + \frac{3 \tau_{0 2}}{2 \alpha} \left(- \tanh^{2}{\left (\frac{x}{2 \alpha} \right )} + 1\right) \tanh{\left (\frac{x}{2 \alpha} \right )}$$
-                           3                          2
-            ⎛    ⎛ x ⎞    ⎞            ⎛    ⎛ x ⎞    ⎞           ⎛    ⎛ x ⎞    ⎞
-            ⎜tanh⎜───⎟    ⎟            ⎜tanh⎜───⎟    ⎟           ⎜tanh⎜───⎟    ⎟
-            ⎜    ⎝2⋅α⎠   1⎟            ⎜    ⎝2⋅α⎠   1⎟           ⎜    ⎝2⋅α⎠   1⎟
-    12⋅τ₀ ₂⋅⎜───────── + ─⎟    18⋅τ₀ ₂⋅⎜───────── + ─⎟    6⋅τ₀ ₂⋅⎜───────── + ─⎟
-            ⎝    2       2⎠            ⎝    2       2⎠           ⎝    2       2⎠
-    ──────────────────────── - ──────────────────────── + ────────────────────── +
-               α                          α                         α
-    
-    
-    
-    
-            ⎛      2⎛ x ⎞    ⎞     ⎛ x ⎞
-     3⋅τ₀ ₂⋅⎜- tanh ⎜───⎟ + 1⎟⋅tanh⎜───⎟
-            ⎝       ⎝2⋅α⎠    ⎠     ⎝2⋅α⎠
-     ───────────────────────────────────
-                     2⋅α
-
-%% Cell type:markdown id: tags:
-
-.. and the result simplified...
-
-%% Cell type:code id: tags:
-
-``` python
-assert diff_eq.expand() == 0
-diff_eq.expand()
-```
-
-%% Output
-
-
-    $$0$$
-    0
-
-%% Cell type:markdown id: tags:
-
-...and indeed the expected tanh profile satisfies this differential equation.
-
-Next lets check for the interface between phase 0 and phase 1:
-
-%% Cell type:code id: tags:
-
-``` python
-for diff_eq in mu_diff_eq:
-    eq = diff_eq.subs({order_params[0]: expected_profile,
-                       order_params[1]: 1 - expected_profile})
-    assert evaluate_diffs(eq, x).expand() == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-Computing the excess free energy per unit area of an interface between two phases.
-This should be exactly the surface tension parameter.
-
-%% Cell type:code id: tags:
-
-``` python
-order_params = symbolic_order_parameters(num_symbols=num_phases-1)
-F = free_energy_functional_n_phases(order_parameters=order_params)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-two_phase_free_energy = F.subs({order_params[0]: expected_profile,
-                                order_params[1]: 1 - expected_profile})
-
-# Evaluate differentials and simplification
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-excess_free_energy = sp.Integral(two_phase_free_energy, x)
-excess_free_energy
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{8 \alpha \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}\, dx$$
-    ⌠
-    ⎮     3⋅τ₀ ₁
-    ⎮ ────────────── dx
-    ⎮         4⎛ x ⎞
-    ⎮ 8⋅α⋅cosh ⎜───⎟
-    ⎮          ⎝2⋅α⎠
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Sympy cannot integrate this automatically - help with a manual substitution $\frac{x}{2\alpha} \rightarrow u$
-
-%% Cell type:code id: tags:
-
-``` python
-coshTerm = list(excess_free_energy.atoms(sp.cosh))[0]
-transformed_int = excess_free_energy.transform(coshTerm.args[0], sp.Symbol("u", real=True))
-transformed_int
-```
-
-%% Output
-
-
-    $$\int \frac{3 \tau_{0 1}}{4 \cosh^{4}{\left (u \right )}}\, du$$
-    ⌠
-    ⎮   3⋅τ₀ ₁
-    ⎮ ────────── du
-    ⎮       4
-    ⎮ 4⋅cosh (u)
-    ⌡
-
-%% Cell type:markdown id: tags:
-
-Now the integral can be done:
-
-%% Cell type:code id: tags:
-
-``` python
-result = sp.integrate(transformed_int.args[0], (transformed_int.args[1][0], -sp.oo, sp.oo))
-assert result == symmetric_symbolic_surface_tension(0,1)
-result
-```
-
-%% Output
-
-
-    $$\tau_{0 1}$$
-    τ₀ ₁
--- a/lbmpy_tests/phasefield/test_liquid_lens_angles.ipynb
+++ b/lbmpy_tests/phasefield/test_liquid_lens_angles.ipynb
--- a/lbmpy_tests/phasefield/test_n_phase_paper_boyer_lbm.ipynb
+++ b/lbmpy_tests/phasefield/test_n_phase_paper_boyer_lbm.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-
-from pystencils.fd.spatial import discretize_spatial
-from pystencils.simp import sympy_cse_on_assignment_list
-from lbmpy.phasefield.cahn_hilliard_lbm import *
-from pystencils.fd.derivation import *
-
-one = sp.sympify(1)
-```
-
-%% Cell type:markdown id: tags:
-
-# Chemical potential
-
-Current state:
-
- not stable (yet)
- without LB coupling the model is stable
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (100, 100)
-
-n = 4
-c = sp.symbols(f"c_:{n}")
-simple_potential = False
-omega_h = 1.3
-
-if simple_potential:
-    f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
-else:
-    ε = one * 4
-    mobility = one * 2 / 1000
-    κ = (one,  one/2, one/3, one/4)
-    sigma_factor = one / 15
-    σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-
-    α, _ = diffusion_coefficients(σ)
-    f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
-                                                              assume_nonnegative=True, zero_threshold=1e-10)
-```
-
-%% Cell type:markdown id: tags:
-
-# Data Setup
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
-p_linearization = symmetric_tensor_linearization(dh.dim)
-
-c_field = dh.add_array("c", values_per_cell=n)
-rho_field = dh.add_array("rho")
-mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
-pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
-force_field = dh.add_array("F", values_per_cell=dh.dim)
-u_field = dh.add_array("u", values_per_cell=dh.dim)
-
-# Distribution functions for each order parameter
-pdf_field = []
-pdf_dst_field = []
-for i in range(n):
-    pdf_field_local = dh.add_array("f%d" % i, values_per_cell=9) # 9 for D2Q9
-    pdf_dst_field_local = dh.add_array("f%d_dst"%i, values_per_cell=9, latex_name="f%d_{dst}" % i)
-    pdf_field.append(pdf_field_local)
-    pdf_dst_field.append(pdf_dst_field_local)
-
-# Distribution functions for the hydrodynamics
-pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
-pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
-```
-
-%% Cell type:markdown id: tags:
-
-### μ-Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
-else:
-    mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
-    mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
-    mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
-
-    mu_assignments = [Assignment(mu_field(i),
-                                 mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
-
-    mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
-
-μ_kernel = create_kernel(mu_assignments).compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-second_neighbor_stencil = [(i, j)
-                           for i in (-2, -1, 0, 1, 2)
-                           for j in (-2, -1, 0, 1, 2)
-                          ]
-x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
-x_diff.set_weight((2, 0), sp.Rational(1, 10))
-x_diff.assume_symmetric(0, anti_symmetric=True)
-x_diff.assume_symmetric(1)
-x_diff_stencil = x_diff.get_stencil(isotropic=True)
-
-y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
-y_diff.set_weight((0, 2), sp.Rational(1, 10))
-y_diff.assume_symmetric(1, anti_symmetric=True)
-y_diff.assume_symmetric(0)
-y_diff_stencil = y_diff.get_stencil(isotropic=True)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ = mu_field
-μ_discretization = {}
-
-for i in range(n):
-    μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
-                             Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
-force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
-force_rhs = force_rhs.subs(μ_discretization)
-force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
-
-force_kernel = create_kernel(force_assignments).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-## Lattice Boltzmann kernels
-
-For each order parameter a Cahn-Hilliard LB method computes the time evolution.
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_alpha = mu_field.center_vector
-else:
-    mu_alpha = mobility * α * mu_field.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# Defining the Cahn-Hilliard Collision assignments
-ch_kernels = []
-
-for i in range(n):
-    ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu_alpha[i],
-                                        relaxation_rate=1.0, gamma=1.0)
-    kernel = create_lb_function(lb_method=ch_method,
-                                velocity_input=u_field.center_vector,
-                                compressible=True,
-                                output={'density': c_field(i)},
-                                optimization={"symbolic_field":pdf_field[i],
-                                              "symbolic_temporary_field": pdf_dst_field[i]})
-
-    ch_kernels.append(kernel)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-hydro_lbm = create_lb_function(relaxation_rate=omega_h, force=force_field,
-                               compressible=True,
-                               optimization={"symbolic_field": pdf_hydro_field,
-                                             "symbolic_temporary_field": pdf_hydro_dst_field},
-                               output={'velocity': u_field}
-                               )
-#hydro_lbm.update_rule
-```
-
-%% Cell type:markdown id: tags:
-
-# Initialization
-
-%% Cell type:code id: tags:
-
-``` python
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c_field.name]
-        arr[..., : ] = values
-
-def init():
-    dh.fill(u_field.name, 0)
-    dh.fill(mu_field.name, 0)
-    dh.fill(force_field.name, 0)
-
-    set_c(make_slice[:, :], [0, 0, 0, 0])
-    set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-    set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-    set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-
-pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
-hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
-c_sync_fn = dh.synchronization_function([c_field.name])
-mu_sync_fn = dh.synchronization_function([mu_field.name])
-
-def time_loop(steps):
-    for i in range(steps):
-        c_sync_fn()
-        dh.run_kernel(μ_kernel)
-
-        mu_sync_fn()
-        dh.run_kernel(force_kernel)
-
-        hydro_sync_fn()
-        dh.run_kernel(hydro_lbm)
-        dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
-
-        pdf_sync_fns()
-        for i in range(n):
-            dh.run_kernel(ch_kernels[i])
-            dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-plt.phase_plot(dh.gather_array(c_field.name))
-```
-
-%% Output
-
-    /local/bauer/miniconda3/envs/pystencils_dev/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range.
-      warnings.warn("No contour levels were found"
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-time_loop(10)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
-plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
-plt.colorbar();
-```
-
-%% Output
-
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-
-from pystencils.fd.spatial import discretize_spatial
-from pystencils.simp import sympy_cse_on_assignment_list
-from lbmpy.phasefield.cahn_hilliard_lbm import *
-from pystencils.fd.derivation import *
-
-one = sp.sympify(1)
-```
-
-%% Cell type:markdown id: tags:
-
-# Chemical potential
-
-Current state:
-
- not stable (yet)
- without LB coupling the model is stable
-
-%% Cell type:code id: tags:
-
-``` python
-domain_size = (100, 100)
-
-n = 4
-c = sp.symbols(f"c_:{n}")
-simple_potential = False
-omega_h = 1.3
-
-if simple_potential:
-    f = free_energy_functional_n_phases_penalty_term(c, interface_width=1, kappa=(0.05, 0.05/2, 0.05/4))
-else:
-    ε = one * 4
-    mobility = one * 2 / 1000
-    κ = (one,  one/2, one/3, one/4)
-    sigma_factor = one / 15
-    σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-
-    α, _ = diffusion_coefficients(σ)
-    f_b, f_if, mu_b, mu_if = chemical_potential_n_phase_boyer(c, ε, σ, one,
-                                                              assume_nonnegative=True, zero_threshold=1e-10)
-```
-
-%% Cell type:markdown id: tags:
-
-# Data Setup
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_ghost_layers=2)
-p_linearization = symmetric_tensor_linearization(dh.dim)
-
-c_field = dh.add_array("c", values_per_cell=n)
-rho_field = dh.add_array("rho")
-mu_field = dh.add_array("mu", values_per_cell=n, latex_name="\\mu")
-pressure_field = dh.add_array("P", values_per_cell=len(p_linearization))
-force_field = dh.add_array("F", values_per_cell=dh.dim)
-u_field = dh.add_array("u", values_per_cell=dh.dim)
-
-# Distribution functions for each order parameter
-pdf_field = []
-pdf_dst_field = []
-for i in range(n):
-    pdf_field_local = dh.add_array("f%d" % i, values_per_cell=9) # 9 for D2Q9
-    pdf_dst_field_local = dh.add_array("f%d_dst"%i, values_per_cell=9, latex_name="f%d_{dst}" % i)
-    pdf_field.append(pdf_field_local)
-    pdf_dst_field.append(pdf_dst_field_local)
-
-# Distribution functions for the hydrodynamics
-pdf_hydro_field = dh.add_array("fh", values_per_cell=9)
-pdf_hydro_dst_field = dh.add_array("fh_dst", values_per_cell=9, latex_name="fh_{dst}")
-```
-
-%% Cell type:markdown id: tags:
-
-### μ-Kernel
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_assignments = mu_kernel(f, c, c_field, mu_field, discretization='isotropic')
-else:
-    mu_subs = {a: b for a, b in zip(c, c_field.center_vector)}
-    mu_if_discrete = [discretize_spatial(e.subs(mu_subs), dx=1, stencil='isotropic') for e in mu_if]
-    mu_b_discrete = [e.subs(mu_subs) for e in mu_b]
-
-    mu_assignments = [Assignment(mu_field(i),
-                                 mu_if_discrete[i] + mu_b_discrete[i]) for i in range(n)]
-
-    mu_assignments = sympy_cse_on_assignment_list(mu_assignments)
-
-μ_kernel = create_kernel(mu_assignments).compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-second_neighbor_stencil = [(i, j)
-                           for i in (-2, -1, 0, 1, 2)
-                           for j in (-2, -1, 0, 1, 2)
-                          ]
-x_diff = FiniteDifferenceStencilDerivation((0,), second_neighbor_stencil)
-x_diff.set_weight((2, 0), sp.Rational(1, 10))
-x_diff.assume_symmetric(0, anti_symmetric=True)
-x_diff.assume_symmetric(1)
-x_diff_stencil = x_diff.get_stencil(isotropic=True)
-
-y_diff = FiniteDifferenceStencilDerivation((1,), second_neighbor_stencil)
-y_diff.set_weight((0, 2), sp.Rational(1, 10))
-y_diff.assume_symmetric(1, anti_symmetric=True)
-y_diff.assume_symmetric(0)
-y_diff_stencil = y_diff.get_stencil(isotropic=True)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ = mu_field
-μ_discretization = {}
-
-for i in range(n):
-    μ_discretization.update({Diff(μ(i), 0): x_diff_stencil.apply(μ(i)),
-                             Diff(μ(i), 1): y_diff_stencil.apply(μ(i))})
-force_rhs = force_from_phi_and_mu(order_parameters=c_field.center_vector, dim=dh.dim, mu=mu_field.center_vector)
-force_rhs = force_rhs.subs(μ_discretization)
-force_assignments = [Assignment(force_field(i), force_rhs[i]) for i in range(dh.dim)]
-
-force_kernel = create_kernel(force_assignments).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-## Lattice Boltzmann kernels
-
-For each order parameter a Cahn-Hilliard LB method computes the time evolution.
-
-%% Cell type:code id: tags:
-
-``` python
-if simple_potential:
-    mu_alpha = mu_field.center_vector
-else:
-    mu_alpha = mobility * α * mu_field.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# Defining the Cahn-Hilliard Collision assignments
-ch_kernels = []
-
-for i in range(n):
-    ch_method = cahn_hilliard_lb_method(get_stencil("D2Q9"), mu_alpha[i],
-                                        relaxation_rate=1.0, gamma=1.0)
-    kernel = create_lb_function(lb_method=ch_method,
-                                velocity_input=u_field.center_vector,
-                                compressible=True,
-                                output={'density': c_field(i)},
-                                optimization={"symbolic_field":pdf_field[i],
-                                              "symbolic_temporary_field": pdf_dst_field[i]})
-
-    ch_kernels.append(kernel)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-hydro_lbm = create_lb_function(relaxation_rate=omega_h, force=force_field,
-                               compressible=True,
-                               optimization={"symbolic_field": pdf_hydro_field,
-                                             "symbolic_temporary_field": pdf_hydro_dst_field},
-                               output={'velocity': u_field}
-                               )
-#hydro_lbm.update_rule
-```
-
-%% Cell type:markdown id: tags:
-
-# Initialization
-
-%% Cell type:code id: tags:
-
-``` python
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c_field.name]
-        arr[..., : ] = values
-
-def init():
-    dh.fill(u_field.name, 0)
-    dh.fill(mu_field.name, 0)
-    dh.fill(force_field.name, 0)
-
-    set_c(make_slice[:, :], [0, 0, 0, 0])
-    set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-    set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-    set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-
-pdf_sync_fns = dh.synchronization_function([f.name for f in pdf_field])
-hydro_sync_fn=dh.synchronization_function([pdf_hydro_field.name])
-c_sync_fn = dh.synchronization_function([c_field.name])
-mu_sync_fn = dh.synchronization_function([mu_field.name])
-
-def time_loop(steps):
-    for i in range(steps):
-        c_sync_fn()
-        dh.run_kernel(μ_kernel)
-
-        mu_sync_fn()
-        dh.run_kernel(force_kernel)
-
-        hydro_sync_fn()
-        dh.run_kernel(hydro_lbm)
-        dh.swap(pdf_hydro_field.name, pdf_hydro_dst_field.name)
-
-        pdf_sync_fns()
-        for i in range(n):
-            dh.run_kernel(ch_kernels[i])
-            dh.swap(pdf_field[i].name,pdf_dst_field[i].name)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-init()
-plt.phase_plot(dh.gather_array(c_field.name))
-```
-
-%% Output
-
-    /local/bauer/miniconda3/envs/pystencils_dev/lib/python3.7/site-packages/matplotlib/contour.py:1243: UserWarning: No contour levels were found within the data range.
-      warnings.warn("No contour levels were found"
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-time_loop(10)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#plt.scalar_field(dh.cpu_arrays[mu_field.name][..., 0])
-plt.scalar_field(dh.cpu_arrays[u_field.name][..., 0])
-plt.colorbar();
-```
-
-%% Output
-
-
--- a/lbmpy_tests/phasefield/test_n_phase_penaltyterm_model.ipynb
+++ b/lbmpy_tests/phasefield/test_n_phase_penaltyterm_model.ipynb
--- a/lbmpy_tests/phasefield/test_neumann_evaluation_circle_fitting.ipynb
+++ b/lbmpy_tests/phasefield/test_neumann_evaluation_circle_fitting.ipynb
-%% Cell type:code id: tags:
-
-``` python
-import numpy as np
-from lbmpy.session import *
-from lbmpy.phasefield.experiments2D import *
-from lbmpy.phasefield.post_processing import *
-from lbmpy.phasefield.contact_angle_circle_fitting import *
-```
-
-%% Cell type:markdown id: tags:
-
-# Testing Neumann angle evaluation based on circle fitting
-
-Set up a 3 phase model to have example data
-
-%% Cell type:code id: tags:
-
-``` python
-kappa3 = 0.03
-alpha = 1
-
-sc = liquid_lens_setup(domain_size=(150, 60), optimization={'target': 'cpu'},
-                       kappas=(0.01, 0.02, kappa3),
-                       cahn_hilliard_relaxation_rates=[np.nan, 1, 3/2],
-                       cahn_hilliard_gammas=[1, 1, 1/3],
-                       alpha=alpha)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-sc.run(10000)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot_for_step(sc)
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-angles = liquid_lens_neumann_angles(sc.concentration[:, :, :])
-assert sum(angles) == 360
-angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 97.07526088545221, \  120.64387551356494, \  142.28086360098288\right]$
-    [97.07526088545221, 120.64387551356494, 142.28086360098288]
-
-%% Cell type:code id: tags:
-
-``` python
-analytic_angles = analytic_neumann_angles([0.01, 0.02, kappa3])
-analytic_angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 89.99999999999999, \  126.86989764584402, \  143.13010235415598\right]$
-    [89.99999999999999, 126.86989764584402, 143.13010235415598]
-
-%% Cell type:code id: tags:
-
-``` python
-for ref, simulated in zip(analytic_angles, angles):
-    assert np.abs(ref - simulated) < 8
-```
-%% Cell type:code id: tags:
-
-``` python
-import numpy as np
-from lbmpy.session import *
-from lbmpy.phasefield.experiments2D import *
-from lbmpy.phasefield.post_processing import *
-from lbmpy.phasefield.contact_angle_circle_fitting import *
-```
-
-%% Cell type:markdown id: tags:
-
-# Testing Neumann angle evaluation based on circle fitting
-
-Set up a 3 phase model to have example data
-
-%% Cell type:code id: tags:
-
-``` python
-kappa3 = 0.03
-alpha = 1
-
-sc = liquid_lens_setup(domain_size=(150, 60), optimization={'target': 'cpu'},
-                       kappas=(0.01, 0.02, kappa3),
-                       cahn_hilliard_relaxation_rates=[np.nan, 1, 3/2],
-                       cahn_hilliard_gammas=[1, 1, 1/3],
-                       alpha=alpha)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-sc.run(10000)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot_for_step(sc)
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-angles = liquid_lens_neumann_angles(sc.concentration[:, :, :])
-assert sum(angles) == 360
-angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 97.07526088545221, \  120.64387551356494, \  142.28086360098288\right]$
-    [97.07526088545221, 120.64387551356494, 142.28086360098288]
-
-%% Cell type:code id: tags:
-
-``` python
-analytic_angles = analytic_neumann_angles([0.01, 0.02, kappa3])
-analytic_angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 89.99999999999999, \  126.86989764584402, \  143.13010235415598\right]$
-    [89.99999999999999, 126.86989764584402, 143.13010235415598]
-
-%% Cell type:code id: tags:
-
-``` python
-for ref, simulated in zip(analytic_angles, angles):
-    assert np.abs(ref - simulated) < 8
-```
--- a/lbmpy_tests/phasefield/test_numerical_1D_3phase_model.ipynb
+++ b/lbmpy_tests/phasefield/test_numerical_1D_3phase_model.ipynb
--- a/lbmpy_tests/phasefield/test_numerical_1D_nphase_model.ipynb
+++ b/lbmpy_tests/phasefield/test_numerical_1D_nphase_model.ipynb
--- a/lbmpy_tests/phasefield/test_phasefieldstep_direct.ipynb
+++ b/lbmpy_tests/phasefield/test_phasefieldstep_direct.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
-from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
-from pystencils.fd import Diff
-```
-
-%% Cell type:markdown id: tags:
-
-# Test of phase field with 4th order FD
-
-%% Cell type:markdown id: tags:
-
-### Free energy definition:
-
-%% Cell type:code id: tags:
-
-``` python
-num_phases = 3
-kappa = [0.01, 0.01, 0.01]
-penalty_factor = 0.01
-domain_size = (40, 40)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-c = sp.symbols("c_:{}".format(num_phases))
-
-def f(c):
-    return c**2 * (1 - c)**2
-
-free_energy = sum((kappa[i] / 2) * f( c[i] ) + kappa[i]/2 * (Diff(c[i]))**2
-                  for i in range(num_phases))
-free_energy += penalty_factor*(1-sum(c[i] for i in range(num_phases)))**2
-
-free_energy
-```
-
-%% Output
-
-
-    $\displaystyle 0.005 c_{0}^{2} \left(1 - c_{0}\right)^{2} + 0.005 c_{1}^{2} \left(1 - c_{1}\right)^{2} + 0.005 c_{2}^{2} \left(1 - c_{2}\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1\right)^{2} + 0.005 {\partial c_{0}}^{2} + 0.005 {\partial c_{1}}^{2} + 0.005 {\partial c_{2}}^{2}$
-            2         2           2         2           2         2
-    0.005⋅c₀ ⋅(1 - c₀)  + 0.005⋅c₁ ⋅(1 - c₁)  + 0.005⋅c₂ ⋅(1 - c₂)  + 0.01⋅(-c₀ -
-    
-                2               2               2               2
-    c₁ - c₂ + 1)  + 0.005⋅D(c_0)  + 0.005⋅D(c_1)  + 0.005⋅D(c_2)
-
-%% Cell type:markdown id: tags:
-
-### Simulation:
-
-%% Cell type:code id: tags:
-
-``` python
-step = PhaseFieldStepDirect(free_energy, c, domain_size)
-
-# geometric setup
-step.set_concentration(make_slice[:, :], [0, 0, 0])
-step.set_single_concentration(make_slice[:, 0.5:], phase_idx=0)
-step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
-step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
-
-step.set_pdf_fields_from_macroscopic_values()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-for i in range(1500):
-    step.time_step()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-angles = liquid_lens_neumann_angles(step.phi[:, :])
-assert angles[0] > 107
-angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 110.18261758183783, \  104.07972456756643, \  145.73765785059575\right]$
-    [110.18261758183783, 104.07972456756643, 145.73765785059575]
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
-from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
-from pystencils.fd import Diff
-```
-
-%% Cell type:markdown id: tags:
-
-# Test of phase field with 4th order FD
-
-%% Cell type:markdown id: tags:
-
-### Free energy definition:
-
-%% Cell type:code id: tags:
-
-``` python
-num_phases = 3
-kappa = [0.01, 0.01, 0.01]
-penalty_factor = 0.01
-domain_size = (40, 40)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-c = sp.symbols("c_:{}".format(num_phases))
-
-def f(c):
-    return c**2 * (1 - c)**2
-
-free_energy = sum((kappa[i] / 2) * f( c[i] ) + kappa[i]/2 * (Diff(c[i]))**2
-                  for i in range(num_phases))
-free_energy += penalty_factor*(1-sum(c[i] for i in range(num_phases)))**2
-
-free_energy
-```
-
-%% Output
-
-
-    $\displaystyle 0.005 c_{0}^{2} \left(1 - c_{0}\right)^{2} + 0.005 c_{1}^{2} \left(1 - c_{1}\right)^{2} + 0.005 c_{2}^{2} \left(1 - c_{2}\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1\right)^{2} + 0.005 {\partial c_{0}}^{2} + 0.005 {\partial c_{1}}^{2} + 0.005 {\partial c_{2}}^{2}$
-            2         2           2         2           2         2
-    0.005⋅c₀ ⋅(1 - c₀)  + 0.005⋅c₁ ⋅(1 - c₁)  + 0.005⋅c₂ ⋅(1 - c₂)  + 0.01⋅(-c₀ -
-    
-                2               2               2               2
-    c₁ - c₂ + 1)  + 0.005⋅D(c_0)  + 0.005⋅D(c_1)  + 0.005⋅D(c_2)
-
-%% Cell type:markdown id: tags:
-
-### Simulation:
-
-%% Cell type:code id: tags:
-
-``` python
-step = PhaseFieldStepDirect(free_energy, c, domain_size)
-
-# geometric setup
-step.set_concentration(make_slice[:, :], [0, 0, 0])
-step.set_single_concentration(make_slice[:, 0.5:], phase_idx=0)
-step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
-step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
-
-step.set_pdf_fields_from_macroscopic_values()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-for i in range(1500):
-    step.time_step()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-angles = liquid_lens_neumann_angles(step.phi[:, :])
-assert angles[0] > 107
-angles
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 110.18261758183783, \  104.07972456756643, \  145.73765785059575\right]$
-    [110.18261758183783, 104.07972456756643, 145.73765785059575]
--- a/lbmpy_tests/phasefield/test_phasefieldstep_direct_boyer.ipynb
+++ b/lbmpy_tests/phasefield/test_phasefieldstep_direct_boyer.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
-from lbmpy.phasefield.n_phase_boyer import *
-
-from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
-from pystencils.fd import Diff
-
-one = sp.sympify(1)
-
-import pyximport
-pyximport.install(language_level=3)
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-```
-
-%% Cell type:markdown id: tags:
-
-# Test of phase field with 4th order FD
-
-%% Cell type:markdown id: tags:
-
-### Free energy definition:
-
-%% Cell type:code id: tags:
-
-``` python
-n = 3
-domain_size = (150, 150)
-
-ε = one * 4
-penalty_factor = 0.01
-stabilization_factor = 10
-
-#κ = (one,  one/2, one/3, one/4)
-κ = (one, one, one, one)
-sigma_factor = one / 4 * 67 * 100
-σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-c = sp.symbols("c_:{}".format(n))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-α, _ = diffusion_coefficients(σ)
-
-f = lambda c: c**2 * ( 1 - c ) **2
-a, b = compute_ab(f)
-
-capital_f = capital_f0(c, σ) + correction_g(c, σ) + stabilization_factor * stabilization_term(c, α)
-
-f_bulk = free_energy_bulk(capital_f, b, ε)
-f_if = free_energy_interfacial(c, σ, a, ε)
-free_energy = (f_bulk + f_if) / 20000 / 100 + penalty_factor * (one - sum(c))**2
-```
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy.evalf()
-```
-
-%% Output
-
-    $$1.5 \cdot 10^{-5} c_{0}^{2} c_{1}^{2} c_{2}^{2} + 0.005025 c_{0}^{2} \left(- c_{0} + 1.0\right)^{2} + 0.005025 c_{1}^{2} \left(- c_{1} + 1.0\right)^{2} + 0.005025 c_{2}^{2} \left(- c_{2} + 1.0\right)^{2} - 0.0025125 \left(c_{0} + c_{1}\right)^{2} \left(- c_{0} - c_{1} + 1.0\right)^{2} - 0.0025125 \left(c_{0} + c_{2}\right)^{2} \left(- c_{0} - c_{2} + 1.0\right)^{2} - 0.0025125 \left(c_{1} + c_{2}\right)^{2} \left(- c_{1} - c_{2} + 1.0\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1.0\right)^{2} - 0.005025 {\partial c_{0}} {\partial c_{1}} - 0.005025 {\partial c_{0}} {\partial c_{2}} - 0.005025 {\partial c_{1}} {\partial c_{2}}$$
-             2   2   2              2            2              2            2
-    1.5e-5⋅c₀ ⋅c₁ ⋅c₂  + 0.005025⋅c₀ ⋅(-c₀ + 1.0)  + 0.005025⋅c₁ ⋅(-c₁ + 1.0)  + 0
-    
-              2            2                      2                 2
-    .005025⋅c₂ ⋅(-c₂ + 1.0)  - 0.0025125⋅(c₀ + c₁) ⋅(-c₀ - c₁ + 1.0)  - 0.0025125⋅
-    
-             2                 2                      2                 2
-    (c₀ + c₂) ⋅(-c₀ - c₂ + 1.0)  - 0.0025125⋅(c₁ + c₂) ⋅(-c₁ - c₂ + 1.0)  + 0.01⋅(
-    
-                        2
-    -c₀ - c₁ - c₂ + 1.0)  - 0.005025⋅D(c_0)⋅D(c_1) - 0.005025⋅D(c_0)⋅D(c_2) - 0.00
-    
-    
-    5025⋅D(c_1)⋅D(c_2)
-
-%% Cell type:markdown id: tags:
-
-### Simulation:
-
-%% Cell type:code id: tags:
-
-``` python
-step = PhaseFieldStepDirect(free_energy, c, domain_size)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# geometric setup
-step.set_concentration(make_slice[:, :], [0, 0, 0])
-step.set_single_concentration(make_slice[:, :], phase_idx=0)
-step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
-step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
-#step.smooth(4)
-step.set_pdf_fields_from_macroscopic_values()
-
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-for i in range(10):
-    step.time_step()
-    #simplex_projection_2d(step.data_handling.cpu_arrays[step.phi_field.name])
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-plt.plot(step.phi[25, :])
-```
-
-%% Output
-
-    [<matplotlib.lines.Line2D at 0x7f1106c0eda0>,
-     <matplotlib.lines.Line2D at 0x7f1106c0eef0>,
-     <matplotlib.lines.Line2D at 0x7f1106b98080>]
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-plt.plot(step.phi[15, :])
-```
-
-%% Output
-
-    [<matplotlib.lines.Line2D at 0x7f1106c3bef0>,
-     <matplotlib.lines.Line2D at 0x7f11077832e8>,
-     <matplotlib.lines.Line2D at 0x7f1106c27048>]
-
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.phasefieldstep_direct import PhaseFieldStepDirect
-from lbmpy.phasefield.n_phase_boyer import *
-
-from lbmpy.phasefield.contact_angle_circle_fitting import liquid_lens_neumann_angles
-from pystencils.fd import Diff
-
-one = sp.sympify(1)
-
-import pyximport
-pyximport.install(language_level=3)
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-```
-
-%% Cell type:markdown id: tags:
-
-# Test of phase field with 4th order FD
-
-%% Cell type:markdown id: tags:
-
-### Free energy definition:
-
-%% Cell type:code id: tags:
-
-``` python
-n = 3
-domain_size = (150, 150)
-
-ε = one * 4
-penalty_factor = 0.01
-stabilization_factor = 10
-
-#κ = (one,  one/2, one/3, one/4)
-κ = (one, one, one, one)
-sigma_factor = one / 4 * 67 * 100
-σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-c = sp.symbols("c_:{}".format(n))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-α, _ = diffusion_coefficients(σ)
-
-f = lambda c: c**2 * ( 1 - c ) **2
-a, b = compute_ab(f)
-
-capital_f = capital_f0(c, σ) + correction_g(c, σ) + stabilization_factor * stabilization_term(c, α)
-
-f_bulk = free_energy_bulk(capital_f, b, ε)
-f_if = free_energy_interfacial(c, σ, a, ε)
-free_energy = (f_bulk + f_if) / 20000 / 100 + penalty_factor * (one - sum(c))**2
-```
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy.evalf()
-```
-
-%% Output
-
-    $$1.5 \cdot 10^{-5} c_{0}^{2} c_{1}^{2} c_{2}^{2} + 0.005025 c_{0}^{2} \left(- c_{0} + 1.0\right)^{2} + 0.005025 c_{1}^{2} \left(- c_{1} + 1.0\right)^{2} + 0.005025 c_{2}^{2} \left(- c_{2} + 1.0\right)^{2} - 0.0025125 \left(c_{0} + c_{1}\right)^{2} \left(- c_{0} - c_{1} + 1.0\right)^{2} - 0.0025125 \left(c_{0} + c_{2}\right)^{2} \left(- c_{0} - c_{2} + 1.0\right)^{2} - 0.0025125 \left(c_{1} + c_{2}\right)^{2} \left(- c_{1} - c_{2} + 1.0\right)^{2} + 0.01 \left(- c_{0} - c_{1} - c_{2} + 1.0\right)^{2} - 0.005025 {\partial c_{0}} {\partial c_{1}} - 0.005025 {\partial c_{0}} {\partial c_{2}} - 0.005025 {\partial c_{1}} {\partial c_{2}}$$
-             2   2   2              2            2              2            2
-    1.5e-5⋅c₀ ⋅c₁ ⋅c₂  + 0.005025⋅c₀ ⋅(-c₀ + 1.0)  + 0.005025⋅c₁ ⋅(-c₁ + 1.0)  + 0
-    
-              2            2                      2                 2
-    .005025⋅c₂ ⋅(-c₂ + 1.0)  - 0.0025125⋅(c₀ + c₁) ⋅(-c₀ - c₁ + 1.0)  - 0.0025125⋅
-    
-             2                 2                      2                 2
-    (c₀ + c₂) ⋅(-c₀ - c₂ + 1.0)  - 0.0025125⋅(c₁ + c₂) ⋅(-c₁ - c₂ + 1.0)  + 0.01⋅(
-    
-                        2
-    -c₀ - c₁ - c₂ + 1.0)  - 0.005025⋅D(c_0)⋅D(c_1) - 0.005025⋅D(c_0)⋅D(c_2) - 0.00
-    
-    
-    5025⋅D(c_1)⋅D(c_2)
-
-%% Cell type:markdown id: tags:
-
-### Simulation:
-
-%% Cell type:code id: tags:
-
-``` python
-step = PhaseFieldStepDirect(free_energy, c, domain_size)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# geometric setup
-step.set_concentration(make_slice[:, :], [0, 0, 0])
-step.set_single_concentration(make_slice[:, :], phase_idx=0)
-step.set_single_concentration(make_slice[:, :0.5], phase_idx=1)
-step.set_single_concentration(make_slice[0.25:0.75, 0.25:0.75], phase_idx=2)
-#step.smooth(4)
-step.set_pdf_fields_from_macroscopic_values()
-
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-for i in range(10):
-    step.time_step()
-    #simplex_projection_2d(step.data_handling.cpu_arrays[step.phi_field.name])
-plt.phase_plot(step.phi[:, :])
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-plt.plot(step.phi[25, :])
-```
-
-%% Output
-
-    [<matplotlib.lines.Line2D at 0x7f1106c0eda0>,
-     <matplotlib.lines.Line2D at 0x7f1106c0eef0>,
-     <matplotlib.lines.Line2D at 0x7f1106b98080>]
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-plt.plot(step.phi[15, :])
-```
-
-%% Output
-
-    [<matplotlib.lines.Line2D at 0x7f1106c3bef0>,
-     <matplotlib.lines.Line2D at 0x7f11077832e8>,
-     <matplotlib.lines.Line2D at 0x7f1106c27048>]
-
-
--- a/lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
+++ b/lbmpy_tests/phasefield_allen_cahn/test_codegen_3d.py
-import numpy as np
-
-from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
-from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
-from lbmpy.phasefield_allen_cahn.kernel_equations import (
-    get_collision_assignments_hydro, hydrodynamic_force, initializer_kernel_hydro_lb, initializer_kernel_phase_field_lb,
-    interface_tracking_force)
-from lbmpy.phasefield_allen_cahn.parameter_calculation import (
-    calculate_dimensionless_rising_bubble, calculate_parameters_rti)
-from lbmpy.stencils import get_stencil
-from pystencils import AssignmentCollection, fields
-
-
-def test_codegen_3d():
-    stencil_phase = get_stencil("D3Q15")
-    stencil_hydro = get_stencil("D3Q27")
-    assert (len(stencil_phase[0]) == len(stencil_hydro[0]))
-    dimensions = len(stencil_hydro[0])
-
-    parameters = calculate_dimensionless_rising_bubble(reference_time=18000,
-                                                       density_heavy=1.0,
-                                                       bubble_radius=16,
-                                                       bond_number=30,
-                                                       reynolds_number=420,
-                                                       density_ratio=1000,
-                                                       viscosity_ratio=100)
-
-    np.isclose(parameters["density_light"], 0.001, rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters["gravitational_acceleration"], -9.876543209876543e-08, rtol=1e-05, atol=1e-08, equal_nan=False)
-
-    parameters = calculate_parameters_rti(reference_length=128,
-                                          reference_time=18000,
-                                          density_heavy=1.0,
-                                          capillary_number=9.1,
-                                          reynolds_number=128,
-                                          atwood_number=1.0,
-                                          peclet_number=744,
-                                          density_ratio=3,
-                                          viscosity_ratio=3)
-
-    np.isclose(parameters["density_light"], 1/3, rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters["gravitational_acceleration"], -3.9506172839506174e-07, rtol=1e-05, atol=1e-08, equal_nan=False)
-    np.isclose(parameters["mobility"], 0.0012234169653524492, rtol=1e-05, atol=1e-08, equal_nan=False)
-
-    rs = analytic_rising_speed(1-6, 20, 0.01)
-    np.isclose(rs, 16666.666666666668, rtol=1e-05, atol=1e-08, equal_nan=False)
-
-    density_liquid = 1.0
-    density_gas = 0.001
-    surface_tension = 0.0001
-    M = 0.02
-
-    # phase-field parameter
-    drho3 = (density_liquid - density_gas) / 3
-    # interface thickness
-    W = 5
-    # coefficient related to surface tension
-    beta = 12.0 * (surface_tension / W)
-    # coefficient related to surface tension
-    kappa = 1.5 * surface_tension * W
-    # relaxation rate allen cahn (h)
-    w_c = 1.0 / (0.5 + (3.0 * M))
-
-    # fields
-    u = fields("vel_field(" + str(dimensions) + "): [" + str(dimensions) + "D]", layout='fzyx')
-    C = fields("phase_field: [" + str(dimensions) + "D]", layout='fzyx')
-
-    h = fields("lb_phase_field(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
-    h_tmp = fields("lb_phase_field_tmp(" + str(len(stencil_phase)) + "): [" + str(dimensions) + "D]", layout='fzyx')
-
-    g = fields("lb_velocity_field(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')
-    g_tmp = fields("lb_velocity_field_tmp(" + str(len(stencil_hydro)) + "): [" + str(dimensions) + "D]", layout='fzyx')
-
-    # calculate the relaxation rate for the hydro lb as well as the body force
-    density = density_gas + C.center * (density_liquid - density_gas)
-    # force acting on all phases of the model
-    body_force = np.array([0, 1e-6, 0])
-
-    relaxation_time = 0.03 + 0.5
-    relaxation_rate = 1.0 / relaxation_time
-
-    method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible=True)
-
-    method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
-                                    relaxation_rates=[relaxation_rate, 1, 1, 1, 1, 1],
-                                    maxwellian_moments=True, entropic=False)
-
-    # create the kernels for the initialization of the g and h field
-    h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
-    g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-
-    force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
-    force_model_h = MultiphaseForceModel(force=force_h)
-
-    force_g = hydrodynamic_force(g, C, method_hydro,
-                                 relaxation_time, density_liquid, density_gas, kappa, beta, body_force)
-    force_model_g = MultiphaseForceModel(force=force_g, rho=density)
-
-    h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-    sum_h = np.sum(h_tmp_symbol_list[:])
-
-    method_phase.set_force_model(force_model_h)
-    
-    allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
-                                          velocity_input=u,
-                                          compressible=True,
-                                          optimization={"symbolic_field": h,
-                                                        "symbolic_temporary_field": h_tmp},
-                                          kernel_type='stream_pull_collide')
-
-    allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
-    allen_cahn_update_rule = AssignmentCollection(main_assignments=allen_cahn_lb.main_assignments,
-                                                  subexpressions=allen_cahn_lb.subexpressions)
-    # ---------------------------------------------------------------------------------------------------------
-
-    hydro_lb_update_rule_normal = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                                  density=density,
-                                                                  velocity_input=u,
-                                                                  force=force_g,
-                                                                  sub_iterations=2,
-                                                                  symbolic_fields={"symbolic_field": g,
-                                                                                   "symbolic_temporary_field": g_tmp},
-                                                                  kernel_type='collide_only')
-
-    hydro_lb_update_rule_push = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                                density=density,
-                                                                velocity_input=u,
-                                                                force=force_g,
-                                                                sub_iterations=2,
-                                                                symbolic_fields={"symbolic_field": g,
-                                                                                 "symbolic_temporary_field": g_tmp},
-                                                                kernel_type='collide_stream_push')
--- a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
+++ b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn_2d.ipynb
-%% Cell type:code id: tags:
-
-``` python
-import math
-
-from lbmpy.session import *
-from pystencils.session import *
-
-from lbmpy.moments import MOMENT_SYMBOLS
-from collections import OrderedDict
-
-from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
-
-from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
-from lbmpy.phasefield_allen_cahn.kernel_equations import *
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-except ImportError:
-    pycuda = None
-    gpu = False
-    target = 'cpu'
-    print('No pycuda installed')
-
-if pycuda:
-    gpu = True
-    target = 'gpu'
-```
-
-%% Cell type:code id: tags:
-
-``` python
-stencil_phase = get_stencil("D2Q9")
-stencil_hydro = get_stencil("D2Q9")
-assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
-
-dimensions = len(stencil_phase[0])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# domain
-L0 = 256
-domain_size = (L0, 4 * L0)
-# time step
-timesteps = 1000
-
-# reference time
-reference_time = 8000
-# density of the heavier fluid
-rho_H = 1.0
-
-# calculate the parameters for the RTI
-parameters = calculate_parameters_rti(reference_length=L0,
-                                      reference_time=reference_time,
-                                      density_heavy=rho_H,
-                                      capillary_number=0.44,
-                                      reynolds_number=3000,
-                                      atwood_number=0.998,
-                                      peclet_number=1000,
-                                      density_ratio=1000,
-                                      viscosity_ratio=100)
-# get the parameters
-rho_L = parameters.get("density_light")
-
-mu_H = parameters.get("dynamic_viscosity_heavy")
-mu_L = parameters.get("dynamic_viscosity_light")
-
-tau_H = parameters.get("relaxation_time_heavy")
-tau_L = parameters.get("relaxation_time_light")
-
-sigma = parameters.get("surface_tension")
-M = parameters.get("mobility")
-gravitational_acceleration = parameters.get("gravitational_acceleration")
-
-
-drho3 = (rho_H - rho_L)/3
-# interface thickness
-W = 5
-# coeffcient related to surface tension
-beta = 12.0 * (sigma/W)
-# coeffcient related to surface tension
-kappa = 1.5 * sigma*W
-# relaxation rate allen cahn (h)
-w_c = 1.0/(0.5 + (3.0 * M))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = ps.create_data_handling((domain_size), periodicity = (True, False), parallel=False, default_target=target)
-
-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)
-
-u = dh.add_array("u", values_per_cell=dh.dim)
-dh.fill("u", 0.0, ghost_layers=True)
-
-C = dh.add_array("C")
-dh.fill("C", 0.0, ghost_layers=True)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)
-s8 = 1/(tau)
-
-rho = rho_L + (C.center) * (rho_H - rho_L)
-
-body_force = [0, 0, 0]
-body_force[1] = gravitational_acceleration * rho
-```
-
-%% Cell type:code id: tags:
-
-``` python
-method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible = True)
-
-method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
-                                relaxation_rates=[s8, 1, 1, 1, 1, 1], maxwellian_moments=True, entropic=False)
-```
-
-%% 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 * math.pi * x) / Nx)
-        init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))
-        block["C"][:, :] = init_values
-
-    if gpu:
-        dh.all_to_gpu()
-
-    dh.run_kernel(h_init)
-    dh.run_kernel(g_init)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
-g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-
-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:code id: tags:
-
-``` python
-force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
-force_model_h = MultiphaseForceModel(force=force_h)
-
-force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
-
-method_phase.set_force_model(force_model_h)
-
-allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
-                                      velocity_input = u,
-                                      compressible = True,
-                                      optimization = {"symbolic_field": h,
-                                                      "symbolic_temporary_field": h_tmp},
-                                      kernel_type = 'stream_pull_collide')
-
-allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
-allen_cahn_lb = sympy_cse(allen_cahn_lb)
-
-ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)
-kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                       density=rho,
-                                                       velocity_input=u,
-                                                       force = force_g,
-                                                       sub_iterations=2,
-                                                       symbolic_fields={"symbolic_field": g,
-                                                                        "symbolic_temporary_field": g_tmp},
-                                                       kernel_type='collide_only')
-
-hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)
-
-ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
-kernel_hydro_lb = ast_hydro_lb.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# periodic Boundarys for g, h and C
-periodic_BC_g = dh.synchronization_function(g.name, target=dh.default_target, optimization = {"openmp": True})
-periodic_BC_h = dh.synchronization_function(h.name, target=dh.default_target, optimization = {"openmp": True})
-periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
-
-# No slip boundary for the phasefield lbm
-bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
-                                                 target=dh.default_target, name='boundary_handling_h')
-
-# No slip boundary for the velocityfield lbm
-bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
-                                            target=dh.default_target, name='boundary_handling_g')
-
-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])
-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, :])
-
-
-bh_allen_cahn.prepare()
-bh_hydro.prepare()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-ac_g = create_lb_update_rule(stencil = stencil_hydro,
-                             optimization={"symbolic_field": g,
-                                           "symbolic_temporary_field": g_tmp},
-                             kernel_type='stream_pull_only')
-ast_g = ps.create_kernel(ac_g, target=dh.default_target, cpu_openmp=True)
-stream_g = ast_g.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# definition of the timestep for the immiscible fluids model
-def Improved_PhaseField_h_g():
-    bh_allen_cahn()
-    periodic_BC_h()
-
-    dh.run_kernel(kernel_allen_cahn_lb)
-    periodic_BC_C()
-    dh.run_kernel(kernel_hydro_lb)
-
-    bh_hydro()
-    periodic_BC_g()
-
-    dh.run_kernel(stream_g)
-    dh.swap("h", "h_tmp")
-    dh.swap("g", "g_tmp")
-```
-
-%% Cell type:code id: tags:
-
-``` python
-Initialize_distributions()
-
-for i in range(0, timesteps):
-    Improved_PhaseField_h_g()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-if gpu:
-    dh.to_cpu("C")
-
-Ny = domain_size[1]
-
-pos_liquid_front = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
-pos_bubble_front = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
-```
-
-%% Cell type:code id: tags:
-
-``` python
-assert(np.isclose(pos_liquid_front, -1e-1, atol=0.01))
-assert(np.isclose(pos_bubble_front, 9e-2, atol=0.01))
-```
-%% Cell type:code id: tags:
-
-``` python
-import math
-
-from lbmpy.session import *
-from pystencils.session import *
-
-from lbmpy.moments import MOMENT_SYMBOLS
-from collections import OrderedDict
-
-from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
-
-from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_parameters_rti
-from lbmpy.phasefield_allen_cahn.kernel_equations import *
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
-```
-
-%% Cell type:code id: tags:
-
-``` python
-try:
-    import pycuda
-except ImportError:
-    pycuda = None
-    gpu = False
-    target = 'cpu'
-    print('No pycuda installed')
-
-if pycuda:
-    gpu = True
-    target = 'gpu'
-```
-
-%% Cell type:code id: tags:
-
-``` python
-stencil_phase = get_stencil("D2Q9")
-stencil_hydro = get_stencil("D2Q9")
-assert(len(stencil_phase[0]) == len(stencil_hydro[0]))
-
-dimensions = len(stencil_phase[0])
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# domain
-L0 = 256
-domain_size = (L0, 4 * L0)
-# time step
-timesteps = 1000
-
-# reference time
-reference_time = 8000
-# density of the heavier fluid
-rho_H = 1.0
-
-# calculate the parameters for the RTI
-parameters = calculate_parameters_rti(reference_length=L0,
-                                      reference_time=reference_time,
-                                      density_heavy=rho_H,
-                                      capillary_number=0.44,
-                                      reynolds_number=3000,
-                                      atwood_number=0.998,
-                                      peclet_number=1000,
-                                      density_ratio=1000,
-                                      viscosity_ratio=100)
-# get the parameters
-rho_L = parameters.get("density_light")
-
-mu_H = parameters.get("dynamic_viscosity_heavy")
-mu_L = parameters.get("dynamic_viscosity_light")
-
-tau_H = parameters.get("relaxation_time_heavy")
-tau_L = parameters.get("relaxation_time_light")
-
-sigma = parameters.get("surface_tension")
-M = parameters.get("mobility")
-gravitational_acceleration = parameters.get("gravitational_acceleration")
-
-
-drho3 = (rho_H - rho_L)/3
-# interface thickness
-W = 5
-# coeffcient related to surface tension
-beta = 12.0 * (sigma/W)
-# coeffcient related to surface tension
-kappa = 1.5 * sigma*W
-# relaxation rate allen cahn (h)
-w_c = 1.0/(0.5 + (3.0 * M))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = ps.create_data_handling((domain_size), periodicity = (True, False), parallel=False, default_target=target)
-
-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)
-
-u = dh.add_array("u", values_per_cell=dh.dim)
-dh.fill("u", 0.0, ghost_layers=True)
-
-C = dh.add_array("C")
-dh.fill("C", 0.0, ghost_layers=True)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-tau = 0.5 + tau_L + (C.center) * (tau_H - tau_L)
-s8 = 1/(tau)
-
-rho = rho_L + (C.center) * (rho_H - rho_L)
-
-body_force = [0, 0, 0]
-body_force[1] = gravitational_acceleration * rho
-```
-
-%% Cell type:code id: tags:
-
-``` python
-method_phase = create_lb_method(stencil=stencil_phase, method='srt', relaxation_rate=w_c, compressible = True)
-
-method_hydro = create_lb_method(stencil=stencil_hydro, method="mrt", weighted=True,
-                                relaxation_rates=[s8, 1, 1, 1, 1, 1], maxwellian_moments=True, entropic=False)
-```
-
-%% 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 * math.pi * x) / Nx)
-        init_values = 0.5 + 0.5 * np.tanh((y - tmp) / (W / 2))
-        block["C"][:, :] = init_values
-
-    if gpu:
-        dh.all_to_gpu()
-
-    dh.run_kernel(h_init)
-    dh.run_kernel(g_init)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-h_updates = initializer_kernel_phase_field_lb(h, C, u, method_phase, W)
-g_updates = initializer_kernel_hydro_lb(g, u, method_hydro)
-
-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:code id: tags:
-
-``` python
-force_h = [f / 3 for f in interface_tracking_force(C, stencil_phase, W)]
-force_model_h = MultiphaseForceModel(force=force_h)
-
-force_g = hydrodynamic_force(g, C, method_hydro, tau, rho_H, rho_L, kappa, beta, body_force)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
-
-method_phase.set_force_model(force_model_h)
-
-allen_cahn_lb = create_lb_update_rule(lb_method=method_phase,
-                                      velocity_input = u,
-                                      compressible = True,
-                                      optimization = {"symbolic_field": h,
-                                                      "symbolic_temporary_field": h_tmp},
-                                      kernel_type = 'stream_pull_collide')
-
-allen_cahn_lb.set_main_assignments_from_dict({**allen_cahn_lb.main_assignments_dict, **{C.center: sum_h}})
-allen_cahn_lb = sympy_cse(allen_cahn_lb)
-
-ast_allen_cahn_lb = ps.create_kernel(allen_cahn_lb, target=dh.default_target, cpu_openmp=True)
-kernel_allen_cahn_lb = ast_allen_cahn_lb.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                       density=rho,
-                                                       velocity_input=u,
-                                                       force = force_g,
-                                                       sub_iterations=2,
-                                                       symbolic_fields={"symbolic_field": g,
-                                                                        "symbolic_temporary_field": g_tmp},
-                                                       kernel_type='collide_only')
-
-hydro_lb_update_rule = sympy_cse(hydro_lb_update_rule)
-
-ast_hydro_lb = ps.create_kernel(hydro_lb_update_rule, target=dh.default_target, cpu_openmp=True)
-kernel_hydro_lb = ast_hydro_lb.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# periodic Boundarys for g, h and C
-periodic_BC_g = dh.synchronization_function(g.name, target=dh.default_target, optimization = {"openmp": True})
-periodic_BC_h = dh.synchronization_function(h.name, target=dh.default_target, optimization = {"openmp": True})
-periodic_BC_C = dh.synchronization_function(C.name, target=dh.default_target, optimization = {"openmp": True})
-
-# No slip boundary for the phasefield lbm
-bh_allen_cahn = LatticeBoltzmannBoundaryHandling(method_phase, dh, 'h',
-                                                 target=dh.default_target, name='boundary_handling_h')
-
-# No slip boundary for the velocityfield lbm
-bh_hydro = LatticeBoltzmannBoundaryHandling(method_hydro, dh, 'g' ,
-                                            target=dh.default_target, name='boundary_handling_g')
-
-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])
-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, :])
-
-
-bh_allen_cahn.prepare()
-bh_hydro.prepare()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-ac_g = create_lb_update_rule(stencil = stencil_hydro,
-                             optimization={"symbolic_field": g,
-                                           "symbolic_temporary_field": g_tmp},
-                             kernel_type='stream_pull_only')
-ast_g = ps.create_kernel(ac_g, target=dh.default_target, cpu_openmp=True)
-stream_g = ast_g.compile()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# definition of the timestep for the immiscible fluids model
-def Improved_PhaseField_h_g():
-    bh_allen_cahn()
-    periodic_BC_h()
-
-    dh.run_kernel(kernel_allen_cahn_lb)
-    periodic_BC_C()
-    dh.run_kernel(kernel_hydro_lb)
-
-    bh_hydro()
-    periodic_BC_g()
-
-    dh.run_kernel(stream_g)
-    dh.swap("h", "h_tmp")
-    dh.swap("g", "g_tmp")
-```
-
-%% Cell type:code id: tags:
-
-``` python
-Initialize_distributions()
-
-for i in range(0, timesteps):
-    Improved_PhaseField_h_g()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-if gpu:
-    dh.to_cpu("C")
-
-Ny = domain_size[1]
-
-pos_liquid_front = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
-pos_bubble_front = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
-```
-
-%% Cell type:code id: tags:
-
-``` python
-assert(np.isclose(pos_liquid_front, -1e-1, atol=0.01))
-assert(np.isclose(pos_bubble_front, 9e-2, atol=0.01))
-```
--- a/lbmpy_tests/test_boundary_handling.py
+++ b/lbmpy_tests/test_boundary_handling.py
-import numpy as np
-import pytest
-
-from lbmpy.boundaries import UBB, NeumannByCopy, NoSlip, StreamInConstant
-from lbmpy.boundaries.boundaryhandling import LatticeBoltzmannBoundaryHandling
-from lbmpy.creationfunctions import create_lb_function
-from lbmpy.geometry import add_box_boundary
-from lbmpy.lbstep import LatticeBoltzmannStep
-from pystencils import create_data_handling, make_slice
-
-
-@pytest.mark.parametrize("target", ['cpu', 'gpu', 'opencl'])
-def test_simple(target):
-    if target == 'gpu':
-        import pytest
-        pytest.importorskip('pycuda')
-    elif target == 'opencl':
-        import pytest
-        pytest.importorskip('pyopencl')
-        import pystencils.opencl.autoinit
-
-    dh = create_data_handling((4, 4), parallel=False, default_target=target)
-    dh.add_array('pdfs', values_per_cell=9, cpu=True, gpu=target != 'cpu')
-    for i in range(9):
-        dh.fill("pdfs", i, value_idx=i, ghost_layers=True)
-
-    if target == 'gpu' or target == 'opencl':
-        dh.all_to_gpu()
-
-    lb_func = create_lb_function(stencil='D2Q9',
-                                 compressible=False,
-                                 relaxation_rate=1.8,
-                                 optimization={'target': target})
-
-    bh = LatticeBoltzmannBoundaryHandling(lb_func.method, dh, 'pdfs', target=target)
-
-    wall = NoSlip()
-    moving_wall = UBB((1, 0))
-    bh.set_boundary(wall, make_slice[0, :])
-    bh.set_boundary(wall, make_slice[-1, :])
-    bh.set_boundary(wall, make_slice[:, 0])
-    bh.set_boundary(moving_wall, make_slice[:, -1])
-
-    bh.prepare()
-    bh()
-
-    if target == 'gpu' or target == 'opencl':
-        dh.all_to_cpu()
-    # left lower corner
-    assert (dh.cpu_arrays['pdfs'][0, 0, 6] == 7)
-
-    assert (dh.cpu_arrays['pdfs'][0, 1, 4] == 3)
-    assert (dh.cpu_arrays['pdfs'][0, 1, 6] == 7)
-
-    assert (dh.cpu_arrays['pdfs'][1, 0, 1] == 2)
-    assert (dh.cpu_arrays['pdfs'][1, 0, 6] == 7)
-
-    # left side
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 4] == 3))
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 6] == 7))
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 5] == 5))
-
-    # left upper corner
-    assert (dh.cpu_arrays['pdfs'][0, 4, 4] == 3)
-    assert (dh.cpu_arrays['pdfs'][0, 4, 8] == 5)
-
-    assert (dh.cpu_arrays['pdfs'][0, 5, 8] == 5 + 6 / 36)
-
-    assert (dh.cpu_arrays['pdfs'][1, 5, 8] == 5 + 6 / 36)
-    assert (dh.cpu_arrays['pdfs'][1, 5, 2] == 1)
-
-    # top side
-    assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 2] == 1))
-    assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 7] == 6 - 6 / 36))
-    assert (all(dh.cpu_arrays['pdfs'][2:4, 5, 8] == 5 + 6 / 36))
-
-    # right upper corner
-    assert (dh.cpu_arrays['pdfs'][4, 5, 2] == 1)
-    assert (dh.cpu_arrays['pdfs'][4, 5, 7] == 6 - 6 / 36)
-
-    assert (dh.cpu_arrays['pdfs'][5, 5, 7] == 6 - 6 / 36)
-
-    assert (dh.cpu_arrays['pdfs'][5, 4, 3] == 4)
-    assert (dh.cpu_arrays['pdfs'][5, 4, 7] == 6)
-
-    # right side
-    assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 3] == 4))
-    assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 5] == 8))
-    assert (all(dh.cpu_arrays['pdfs'][5, 2:4, 7] == 6))
-
-    # right lower corner
-    assert (dh.cpu_arrays['pdfs'][5, 1, 3] == 4)
-    assert (dh.cpu_arrays['pdfs'][5, 1, 5] == 8)
-
-    assert (dh.cpu_arrays['pdfs'][5, 0, 5] == 8)
-
-    assert (dh.cpu_arrays['pdfs'][4, 0, 1] == 2)
-    assert (dh.cpu_arrays['pdfs'][4, 0, 5] == 8)
-
-    # lower side
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 4] == 3))
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 6] == 7))
-    assert (all(dh.cpu_arrays['pdfs'][0, 2:4, 8] == 5))
-
-
-def test_exotic_boundaries():
-    step = LatticeBoltzmannStep((50, 50), relaxation_rate=1.8, compressible=False, periodicity=False)
-    add_box_boundary(step.boundary_handling, NeumannByCopy())
-    step.boundary_handling.set_boundary(StreamInConstant(0), make_slice[0, :])
-    step.run(100)
-    assert np.max(step.velocity[:, :, :]) < 1e-13
--- a/lbmpy_tests/test_code_hashequivalence.py
+++ b/lbmpy_tests/test_code_hashequivalence.py
-from hashlib import sha256
-
-from lbmpy.creationfunctions import create_lb_ast
-
-def test_hash_equivalence_llvm():
-
-    import pytest
-    pytest.importorskip("llvmlite")
-    from pystencils.llvm.llvmjit import generate_llvm
-
-
-    ref_value = "6db6ed9e2cbd05edae8fcaeb8168e3178dd578c2681133f3ae9228b23d2be432"
-    ast = create_lb_ast(stencil='D3Q19', method='srt', optimization={'target': 'llvm'})
-    code = generate_llvm(ast)
-    hash_value = sha256(str(code).encode()).hexdigest()
-    assert hash_value == ref_value
--- a/lbmpy_tests/test_compiled_in_boundaries.ipynb
+++ b/lbmpy_tests/test_compiled_in_boundaries.ipynb
-%% Cell type:code id: tags:
-
-``` python
-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
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, default_target='cpu')
-pdfs = dh.add_array('pdfs', values_per_cell=19)
-u = dh.add_array('u', values_per_cell=len(domain_size))
-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()),
-))
-opt = {'symbolic_field': pdfs, 'cse_global': False, 'cse_pdfs': True}
-cr_even = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)
-cr_odd = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False,  optimization=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)
-
-getter_kernel = ps.create_kernel(getter_assignments, target=dh.default_target).compile()
-even_kernel = ps.create_kernel(update_rule_aa_even, target=dh.default_target, ghost_layers=1).compile()
-odd_kernel = ps.create_kernel(update_rule_aa_odd, target=dh.default_target, ghost_layers=1).compile()
-```
-
-%% 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 0x7f9d260364c0>
-
-
-
-
-%% 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
-
-    <matplotlib.colorbar.Colorbar at 0x7f9d2158a520>
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-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
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, default_target='cpu')
-pdfs = dh.add_array('pdfs', values_per_cell=19)
-u = dh.add_array('u', values_per_cell=len(domain_size))
-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()),
-))
-opt = {'symbolic_field': pdfs, 'cse_global': False, 'cse_pdfs': True}
-cr_even = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False, optimization=opt)
-cr_odd = create_lb_collision_rule(stencil="D3Q19", relaxation_rate=relaxation_rate, compressible=False,  optimization=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)
-
-getter_kernel = ps.create_kernel(getter_assignments, target=dh.default_target).compile()
-even_kernel = ps.create_kernel(update_rule_aa_even, target=dh.default_target, ghost_layers=1).compile()
-odd_kernel = ps.create_kernel(update_rule_aa_odd, target=dh.default_target, ghost_layers=1).compile()
-```
-
-%% 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 0x7f9d260364c0>
-
-
-
-
-%% 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
-
-    <matplotlib.colorbar.Colorbar at 0x7f9d2158a520>
-
-
-
--- a/lbmpy_tests/test_cpu_gpu_equivalence.py
+++ b/lbmpy_tests/test_cpu_gpu_equivalence.py
-from copy import deepcopy
-
-import numpy as np
-import pytest
-
-from lbmpy.scenarios import create_channel
-
-
-def run_equivalence_test(scenario_creator, time_steps=13, **params):
-    print("Scenario", params)
-    params['optimization']['target'] = 'cpu'
-    cpu_scenario = scenario_creator(**params)
-    params['optimization']['target'] = 'gpu'
-    gpu_scenario = scenario_creator(**params)
-
-    cpu_scenario.run(time_steps)
-    gpu_scenario.run(time_steps)
-
-    max_vel_error = np.max(np.abs(cpu_scenario.velocity_slice() - gpu_scenario.velocity_slice()))
-    max_rho_error = np.max(np.abs(cpu_scenario.density_slice() - gpu_scenario.density_slice()))
-
-    np.testing.assert_allclose(max_vel_error, 0, atol=1e-14)
-    np.testing.assert_allclose(max_rho_error, 0, atol=1e-14)
-
-
-def test_force_driven_channel_short():
-    default = {
-        'scenario_creator': create_channel,
-        'domain_size': (32, 32),
-        'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
-        'method': 'mrt3',
-        'pressure_difference': 0.001,
-        'optimization': {},
-    }
-    scenarios = []
-
-    # Different methods
-    for ds, method, compressible, block_size, field_layout in [((17, 12), 'srt', False, (12, 4), 'reverse_numpy'),
-                                                               ((18, 20), 'mrt', True, (4, 2), 'zyxf'),
-                                                               ((7, 11, 18), 'trt', False, False, 'numpy')]:
-        params = deepcopy(default)
-        if block_size is not False:
-            params['optimization'].update({
-                'gpu_indexing_params': {'block_size': block_size}
-            })
-        else:
-            params['optimization']['gpu_indexing'] = 'line'
-
-        params['domain_size'] = ds
-        params['method'] = method
-        params['compressible'] = compressible
-        params['optimization']['field_layout'] = field_layout
-        scenarios.append(params)
-
-    for scenario in scenarios:
-        run_equivalence_test(**scenario)
-
-
-@pytest.mark.longrun
-def test_force_driven_channel():
-    default = {
-        'scenario_creator': create_channel,
-        'domain_size': (32, 32),
-        'relaxation_rates': [1.95, 1.9, 1.92, 1.92],
-        'method': 'mrt',
-        'pressure_difference': 0.001,
-        'optimization': {},
-    }
-
-    scenarios = []
-
-    # Different methods
-    for method in ('srt', 'mrt'):
-        for compressible in (True, False):
-            params = deepcopy(default)
-            params['optimization'].update({
-                'gpu_indexing_params': {'block_size': (16, 16)}
-            })
-            params['method'] = method
-            params['compressible'] = compressible
-            scenarios.append(params)
-
-    # Blocked indexing with different block sizes
-    for block_size in ((16, 16), (8, 16), (4, 2)):
-        params = deepcopy(default)
-        params['method'] = 'mrt'
-        params['compressible'] = True
-        params['optimization'].update({
-            'gpu_indexing': 'block',
-            'gpu_indexing_params': {'block_size': block_size}
-        })
-        scenarios.append(params)
-
-    # Line wise indexing
-    params = deepcopy(default)
-    params['optimization']['gpu_indexing'] = 'line'
-    scenarios.append(params)
-
-    # Different field layouts
-    for field_layout in ('numpy', 'reverse_numpy', 'zyxf'):
-        for fixed_size in (False, True):
-            params = deepcopy(default)
-            params['optimization'].update({
-                'gpu_indexing_params': {'block_size': (16, 16)}
-            })
-            if fixed_size:
-                params['optimization']['field_size'] = params['domain_size']
-            else:
-                params['optimization']['field_size'] = None
-
-            params['optimization']['field_layout'] = field_layout
-            scenarios.append(params)
-
-    print("Testing %d scenarios" % (len(scenarios),))
-    for scenario in scenarios:
-        run_equivalence_test(**scenario)
--- a/lbmpy_tests/test_float_kernel.py
+++ b/lbmpy_tests/test_float_kernel.py
-from lbmpy.creationfunctions import create_lb_function
-from lbmpy.scenarios import create_lid_driven_cavity
-from pystencils import show_code
-
-
-def test_creation():
-    """Simple test that makes sure that only float variables are created"""
-
-    func = create_lb_function(method='srt', relaxation_rate=1.5,
-                              optimization={'double_precision': False})
-    code = str(show_code(func.ast))
-    assert 'double' not in code
-
-
-def test_scenario():
-    sc = create_lid_driven_cavity((16, 16, 8), relaxation_rate=1.5,
-                                  optimization={'double_precision': False})
-    sc.run(1)
-    code_str = str(show_code(sc.ast))
-    assert 'double' not in code_str
--- a/lbmpy_tests/test_fluctuating_lb.py
+++ b/lbmpy_tests/test_fluctuating_lb.py
-"""Tests velocity and stress fluctuations for thermalized LB"""
-
-
-import pystencils as ps
-from lbmpy.creationfunctions import *
-from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
-import numpy as np
-from lbmpy.moments import is_bulk_moment, is_shear_moment, get_order
-
-
-def single_component_maxwell(x1, x2, kT, mass):
-    """Integrate the probability density from x1 to x2 using the trapezoidal rule"""
-    x = np.linspace(x1, x2, 1000)
-    return np.trapz(np.exp(-mass * x**2 / (2. * kT)), x) / np.sqrt(2. * np.pi * kT/mass)
-
-
-def rr_getter(moment_group):
-    """Maps a group of moments to a relaxation rate (shear, bulk, even, odd)
-    in the 4 relaxation time thermalized LB model
-    """
-    is_shear = [is_shear_moment(m, 3) for m in moment_group]
-    is_bulk = [is_bulk_moment(m, 3) for m in moment_group]
-    order = [get_order(m) for m in moment_group]
-    assert min(order) == max(order)
-    order = order[0]
-
-    if order < 2:
-        return 0
-    elif any(is_bulk):
-        assert all(is_bulk)
-        return sp.Symbol("omega_bulk")
-    elif any(is_shear):
-        assert all(is_shear)
-        return sp.Symbol("omega_shear")
-    elif order % 2 == 0:
-        assert order > 2
-        return sp.Symbol("omega_even")
-    else:
-        return sp.Symbol("omega_odd")
-
-
-def second_order_moment_tensor_assignments(function_values, stencil, output_field):
-    """Assignments for calculating the pressure tensor"""
-    assert len(function_values) == len(stencil)
-    dim = len(stencil[0])
-    return [ps.Assignment(output_field(i, j),
-                          sum(c[i] * c[j] * f for f, c in zip(function_values, stencil)))
-            for i in range(dim) for j in range(dim)]
-
-
-def add_pressure_output_to_collision_rule(collision_rule, pressure_field):
-    pressure_ouput = second_order_moment_tensor_assignments(collision_rule.method.pre_collision_pdf_symbols,
-                                                            collision_rule.method.stencil, pressure_field)
-    collision_rule.main_assignments = collision_rule.main_assignments + pressure_ouput
-
-
-def get_fluctuating_lb(size=None, kT=None, omega_shear=None, omega_bulk=None, omega_odd=None, omega_even=None, rho_0=None, target=None):
-
-    # Parameters
-    stencil = get_stencil('D3Q19')
-
-    # Setup data handling
-    dh = ps.create_data_handling(
-        [size]*3, periodicity=True, default_target=target)
-    src = dh.add_array('src', values_per_cell=len(stencil), layout='f')
-    dst = dh.add_array_like('dst', 'src')
-    rho = dh.add_array('rho', layout='f', latex_name='\\rho')
-    u = dh.add_array('u', values_per_cell=dh.dim, layout='f')
-    pressure_field = dh.add_array('pressure', values_per_cell=(
-        3, 3), layout='f', gpu=target == 'gpu')
-    force_field = dh.add_array(
-        'force', values_per_cell=3, layout='f', gpu=target == 'gpu')
-
-    # Method setup
-    method = create_mrt_orthogonal(
-        stencil=get_stencil('D3Q19'),
-        compressible=True,
-        weighted=True,
-        relaxation_rate_getter=rr_getter,
-        force_model=force_model_from_string('schiller', force_field.center_vector))
-    collision_rule = create_lb_collision_rule(
-        method,
-        fluctuating={
-            'temperature': kT
-        },
-        optimization={'cse_global': True}
-    )
-
-    add_pressure_output_to_collision_rule(collision_rule, pressure_field)
-
-    collision = create_lb_update_rule(collision_rule=collision_rule,
-                                      stencil=stencil,
-                                      method=method,
-                                      compressible=True,
-                                      kernel_type='collide_only',
-                                      optimization={'symbolic_field': src})
-    stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
-                                                   {'density': rho, 'velocity': u})
-
-    opts = {'cpu_openmp': True,
-            'cpu_vectorize_info': None,
-            'target': dh.default_target}
-
-    # Compile kernels
-    stream_kernel = ps.create_kernel(stream, **opts).compile()
-    collision_kernel = ps.create_kernel(collision, **opts).compile()
-
-    sync_pdfs = dh.synchronization_function([src.name])
-
-    # Initialization
-    init = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim,
-                                     pdfs=src.center_vector, density=rho.center)
-    init_kernel = ps.create_kernel(init, ghost_layers=0).compile()
-
-    dh.fill(rho.name, rho_0)
-    dh.fill(u.name, np.nan, ghost_layers=True, inner_ghost_layers=True)
-    dh.fill(u.name, 0)
-    dh.fill(force_field.name, np.nan,
-            ghost_layers=True, inner_ghost_layers=True)
-    dh.fill(force_field.name, 0)
-    dh.run_kernel(init_kernel)
-
-    # time loop
-    def time_loop(start, steps):
-        dh.all_to_gpu()
-        for i in range(start, start+steps):
-            dh.run_kernel(collision_kernel, omega_shear=omega_shear, omega_bulk=omega_bulk,
-                          omega_odd=omega_odd, omega_even=omega_even, seed=42, time_step=i)
-
-            sync_pdfs()
-            dh.run_kernel(stream_kernel)
-
-            dh.swap(src.name, dst.name)
-        return start+steps
-
-    return dh, time_loop
-
-
-def test_resting_fluid(target="cpu"):
-    rho_0 = 0.86
-    kT = 4E-4
-    L = [60]*3
-    dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
-                                       rho_0=rho_0, kT=kT,
-                                       omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
-
-    # Test
-    t = 0
-    # warm up
-    t = time_loop(t, 10)
-
-    # Measurement
-    for i in range(10):
-        t = time_loop(t, 5)
-
-        res_u = dh.gather_array("u").reshape((-1, 3))
-        res_rho = dh.gather_array("rho").reshape((-1,))
-
-        # mass conservation
-        np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
-
-        # momentum conservation
-        momentum = np.dot(res_rho, res_u)
-        np.testing.assert_allclose(momentum, [0, 0, 0], atol=1E-10)
-
-        # temperature
-        kinetic_energy = 1/2*np.dot(res_rho, res_u*res_u)/np.product(L)
-        np.testing.assert_allclose(
-            kinetic_energy, [kT/2]*3, atol=kT*0.01)
-
-        # velocity distribution
-        v_hist, v_bins = np.histogram(
-            res_u, bins=11, range=(-.075, .075), density=True)
-
-        # Calculate expected values from single
-        v_expected = []
-        for j in range(len(v_hist)):
-            # Maxwell distribution
-            res = 1./(v_bins[j+1]-v_bins[j]) * \
-                single_component_maxwell(
-                    v_bins[j], v_bins[j+1], kT, rho_0)
-            v_expected.append(res)
-        v_expected = np.array(v_expected)
-
-        # 10% accuracy on the entire histogram
-        np.testing.assert_allclose(v_hist, v_expected, rtol=0.1)
-        # 1% accuracy on the middle part
-        remove = 3
-        np.testing.assert_allclose(
-            v_hist[remove:-remove], v_expected[remove:-remove], rtol=0.01)
-
-        # pressure tensor against expressions from
-        # Duenweg, Schiller, Ladd, https://arxiv.org/abs/0707.1581
-
-        res_pressure = dh.gather_array("pressure").reshape((-1, 3, 3))
-
-        c_s = np.sqrt(1/3)  # speed of sound
-
-        # average of pressure tensor
-        # Diagonal elements are rho c_s^22 +<u,u>. When the fluid is
-        # thermalized, the expectation value of <u,u> = kT due to the
-        # equi-partition theorem.
-        p_av_expected = np.diag([rho_0*c_s**2 + kT]*3)
-        np.testing.assert_allclose(
-            np.mean(res_pressure, axis=0), p_av_expected, atol=c_s**2/2000)
-
-
-def test_point_force(target="cpu"):
-    """Test momentum balance for thermalized fluid with applied poitn forces"""
-    rho_0 = 0.86
-    kT = 4E-4
-    L = [8]*3
-    dh, time_loop = get_fluctuating_lb(size=L[0], target=target,
-                                       rho_0=rho_0, kT=kT,
-                                       omega_shear=0.8, omega_bulk=0.5, omega_even=.04, omega_odd=0.3)
-
-    # Test
-    t = 0
-    # warm up
-    t = time_loop(t, 100)
-
-    introduced_momentum = np.zeros(3)
-    for i in range(100):
-        point_force = 1E-5*(np.random.random(3) - .5)
-        introduced_momentum += point_force
-
-        # Note that ghost layers are included in the indexing
-        force_pos = np.random.randint(1, L[0]-2, size=3)
-
-        dh.cpu_arrays["force"][force_pos[0],
-                               force_pos[1], force_pos[2]] = point_force
-        t = time_loop(t, 1)
-        res_u = dh.gather_array("u").reshape((-1, 3))
-        res_rho = dh.gather_array("rho").reshape((-1,))
-
-        # mass conservation
-        np.testing.assert_allclose(np.mean(res_rho), rho_0, atol=3E-12)
-
-        # momentum conservation
-        momentum = np.dot(res_rho, res_u)
-        np.testing.assert_allclose(
-            momentum, introduced_momentum + 0.5 * point_force, atol=1E-10)
-        dh.cpu_arrays["force"][force_pos[0],
-                               force_pos[1], force_pos[2]] = np.zeros(3)
--- a/lbmpy_tests/test_force.py
+++ b/lbmpy_tests/test_force.py
-from pystencils.session import *
-from lbmpy.session import *
-from lbmpy.macroscopic_value_kernels import macroscopic_values_setter
-import lbmpy.forcemodels
-from lbmpy.moments import is_bulk_moment
-
-import pytest
-from contextlib import ExitStack as does_not_raise
-
-
-force_models = [fm.lower() for fm in dir(lbmpy.forcemodels) if fm[0].isupper()]
-
-def test_force_models_available():
-    assert 'guo' in force_models
-    assert 'luo' in force_models
-
-@pytest.mark.parametrize("method", ["srt", "trt"])
-@pytest.mark.parametrize("force_model", force_models)
-@pytest.mark.parametrize("omega", [0.5, 1.5])
-def test_total_momentum(method, force_model, omega):
-    L = (16, 16)
-    stencil = get_stencil("D2Q9")
-    F = [2e-4, -3e-4]
-
-    dh = ps.create_data_handling(L, periodicity=True, default_target='cpu')
-    src = dh.add_array('src', values_per_cell=len(stencil))
-    dst = dh.add_array_like('dst', 'src')
-    ρ = dh.add_array('rho')
-    u = dh.add_array('u', values_per_cell=dh.dim)
-
-    expectation = does_not_raise()
-    skip = False
-    if force_model in ['guo', 'buick'] and method != 'srt':
-        expectation = pytest.raises(AssertionError)
-        skip = True
-    with expectation:
-        collision = create_lb_update_rule(method=method,
-                                          stencil=stencil,
-                                          relaxation_rate=omega, 
-                                          compressible=True,
-                                          force_model=force_model, 
-                                          force=F,
-                                          kernel_type='collide_only',
-                                          optimization={'symbolic_field': src})
-    if skip:
-        return
-
-    stream = create_stream_pull_with_output_kernel(collision.method, src, dst,
-                                                   {'density': ρ, 'velocity': u})
-
-    opts = {'cpu_openmp': True, 
-            'cpu_vectorize_info': None,
-            'target': dh.default_target}
-
-    stream_kernel = ps.create_kernel(stream, **opts).compile()
-    collision_kernel = ps.create_kernel(collision, **opts).compile()
-
-    def init():
-        dh.fill(ρ.name, 1)
-        dh.fill(u.name, 0)
-
-        setter = macroscopic_values_setter(collision.method, velocity=(0,)*dh.dim, 
-                                       pdfs=src.center_vector, density=ρ.center)
-        kernel = ps.create_kernel(setter, ghost_layers=0).compile()
-        dh.run_kernel(kernel)
-
-    sync_pdfs = dh.synchronization_function([src.name])
-    def time_loop(steps):
-        dh.all_to_gpu()
-        for i in range(steps):
-            dh.run_kernel(collision_kernel)
-            sync_pdfs()
-            dh.run_kernel(stream_kernel)
-            dh.swap(src.name, dst.name)
-        dh.all_to_cpu()
-
-    t = 20
-    init()
-    time_loop(t)
-    total = np.sum(dh.gather_array(u.name), axis=(0,1))
-    assert np.allclose(total/np.prod(L)/F/t, 1)
-
-
-@pytest.mark.parametrize("stencil", ["D2Q9", "D3Q15", "D3Q19", "D3Q27"])
-@pytest.mark.parametrize("force_model", ["simple", "schiller"])
-def test_modes(stencil, force_model):
-    """check Schiller's force term in mode space"""
-    stencil = get_stencil(stencil)
-    dim = len(stencil[0])
-    
-    omega_s = sp.Symbol("omega_s")
-    omega_b = sp.Symbol("omega_b")
-    omega_o = sp.Symbol("omega_o")
-    omega_e = sp.Symbol("omega_e")
-    
-    F = [sp.Symbol(f"F_{i}") for i in range(dim)]
-    
-    method = create_lb_method(method="mrt", weighted=True,
-                              stencil=stencil,
-                              relaxation_rates=[omega_s, omega_b, omega_o, omega_e, omega_o, omega_e], 
-                              compressible=True,
-                              force_model=force_model,
-                              force=F)
-    
-    force_moments = sp.simplify(method.moment_matrix * sp.Matrix(method.force_model(method)))
-    
-    # The mass mode should be zero
-    assert force_moments[0] == 0
-    
-    # The momentum moments should contain the force
-    assert list(force_moments[1:dim+1]) == F
-    
-    if force_model == "schiller":
-        num_stresses = (dim*dim-dim)//2+dim
-        lambda_s, lambda_b = -omega_s, -omega_b
-
-        # The stress moments should match eq. 47 from https://doi.org/10.1023/A:1010414013942
-        u = method.first_order_equilibrium_moment_symbols
-        def traceless(m):
-            tr = sp.simplify(sum([m[i,i] for i in range(dim)]))
-            return m - tr/m.shape[0]*sp.eye(m.shape[0])
-        C = sp.Rational(1,2) * (2 + lambda_s) * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
-                                                traceless(sp.Matrix(F) * sp.Matrix(u).transpose())) + \
-            sp.Rational(1,method.dim) * (2 + lambda_b) * sp.Matrix(u).dot(F) * sp.eye(method.dim)
-
-        subs = {sp.Symbol(chr(ord("x")+i)) * sp.Symbol(chr(ord("x")+j)) : C[i,j]
-                for i in range(dim) for j in range(dim)}
-        for force_moment, moment in zip(force_moments[dim+1:dim+1+num_stresses],
-                                        method.moments[dim+1:dim+1+num_stresses]):
-            ref = moment.subs(subs)
-            diff = sp.simplify(ref - force_moment)
-            if is_bulk_moment(moment, dim):
-                assert diff == 0 or isinstance(diff, sp.Rational) # difference should be zero or a constant
-            else:
-                assert diff == 0 # difference should be zero
-
-        ff = sp.Matrix(method.force_model(method))
-        # Check eq. 4.53a from schiller2008thermal
-        assert sp.simplify(sum(ff)) == 0
-        # Check eq. 4.53b from schiller2008thermal
-        assert [sp.simplify(sum(ff[i] * stencil[i][j] for i in range(len(stencil)))) for j in range(dim)] == F
-        # Check eq. 4.61a from schiller2008thermal
-        ref = (2 + lambda_s)/2 * (traceless(sp.Matrix(u) * sp.Matrix(F).transpose()) + \
-                                 traceless(sp.Matrix(F) * sp.Matrix(u).transpose()))
-        s = sp.zeros(dim)
-        for i in range(0, len(stencil)):
-            s += ff[i] * traceless(sp.Matrix(stencil[i]) * sp.Matrix(stencil[i]).transpose())
-        assert sp.simplify(s-ref) == sp.zeros(dim)
-        # Check eq. 4.61b from schiller2008thermal
-        assert sp.simplify(sum(ff[i] * stencil[i][a]**2 for i in range(len(stencil)) for a in range(dim))
-                           - (2+lambda_b)*sp.Matrix(u).dot(F)) == 0
-
-        # All other moments should be zero
-        assert list(force_moments[dim+1+num_stresses:]) == [0] * (len(stencil)-(dim+1+num_stresses))
-    elif force_model == "simple":
-        # All other moments should be zero
-        assert list(force_moments[dim+1:]) == [0] * (len(stencil)-(dim+1))
--- a/lbmpy_tests/test_maxwellian_equilibrium.py
+++ b/lbmpy_tests/test_maxwellian_equilibrium.py
-from lbmpy.cumulants import raw_moment_as_function_of_cumulants
-from lbmpy.maxwellian_equilibrium import *
-from lbmpy.moments import (
-    MOMENT_SYMBOLS, exponents_to_polynomial_representations, moment_matrix,
-    moments_up_to_component_order, moments_up_to_order)
-from lbmpy.stencils import get_stencil
-from pystencils.sympyextensions import remove_higher_order_terms
-
-
-def test_maxwellian_moments():
-    """Check moments of continuous Maxwellian"""
-    rho = sp.Symbol("rho")
-    u = sp.symbols("u_0 u_1 u_2")
-    c_s = sp.Symbol("c_s")
-    eq_moments = get_moments_of_continuous_maxwellian_equilibrium(((0, 0, 0), (0, 0, 1)),
-                                                                  dim=3, rho=rho, u=u, c_s_sq=c_s ** 2)
-    assert eq_moments[0] == rho
-    assert eq_moments[1] == rho * u[2]
-
-    x, y, z = MOMENT_SYMBOLS
-    one = sp.Rational(1, 1)
-    eq_moments = get_moments_of_continuous_maxwellian_equilibrium((one, x, x ** 2, x * y),
-                                                                  dim=2, rho=rho, u=u[:2], c_s_sq=c_s ** 2)
-    assert eq_moments[0] == rho
-    assert eq_moments[1] == rho * u[0]
-    assert eq_moments[2] == rho * (c_s ** 2 + u[0] ** 2)
-    assert eq_moments[3] == rho * u[0] * u[1]
-
-
-def test_continuous_discrete_moment_equivalence():
-    """Check that moments up to order 3 agree with moments of the continuous Maxwellian"""
-    for stencil in [get_stencil(n) for n in ["D2Q9", "D3Q15", "D3Q19", "D3Q27"]]:
-        dim = len(stencil[0])
-        c_s_sq = sp.Rational(1, 3)
-        moments = tuple(moments_up_to_order(3, dim=dim, include_permutations=False))
-        cm = sp.Matrix(get_moments_of_continuous_maxwellian_equilibrium(moments, order=2, dim=dim, c_s_sq=c_s_sq))
-        dm = sp.Matrix(get_moments_of_discrete_maxwellian_equilibrium(stencil, moments, order=2,
-                                                                      compressible=True, c_s_sq=c_s_sq))
-
-        diff = sp.simplify(cm - dm)
-        for d in diff:
-            assert d == 0
-
-
-def test_moment_cumulant_continuous_equivalence():
-    """Test that discrete equilibrium is the same up to order 3 when obtained with following methods
-
-    * eq1: take moments of continuous Maxwellian and transform back to pdf space
-    * eq2: take cumulants of continuous Maxwellian, transform to moments then transform to pdf space
-    * eq3: take discrete equilibrium from LBM literature
-    * eq4: same as eq1 but with built-in function
-    """
-    for stencil in [get_stencil('D2Q9'), get_stencil('D3Q27')]:
-        dim = len(stencil[0])
-        u = sp.symbols("u_:{dim}".format(dim=dim))
-        indices = tuple(moments_up_to_component_order(2, dim=dim))
-        c_s_sq = sp.Rational(1, 3)
-        eq_moments1 = get_moments_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
-        eq_cumulants = get_cumulants_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
-        eq_cumulants = {idx: c for idx, c in zip(indices, eq_cumulants)}
-        eq_moments2 = [raw_moment_as_function_of_cumulants(idx, eq_cumulants) for idx in indices]
-        pdfs_to_moments = moment_matrix(indices, stencil)
-
-        def normalize(expressions):
-            return [remove_higher_order_terms(e.expand(), symbols=u, order=3) for e in expressions]
-
-        eq1 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments1))
-        eq2 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments2))
-        eq3 = normalize(discrete_maxwellian_equilibrium(stencil, order=3, c_s_sq=c_s_sq, compressible=True))
-        eq4 = normalize(generate_equilibrium_by_matching_moments(stencil, indices, c_s_sq=c_s_sq))
-
-        assert eq1 == eq2
-        assert eq2 == eq3
-        assert eq3 == eq4
-
-
-def test_moment_cumulant_continuous_equivalence_polynomial_formulation():
-    """Same as test above, but instead of index tuples, the polynomial formulation is used."""
-    for stencil in [get_stencil('D2Q9'), get_stencil('D3Q27')]:
-        dim = len(stencil[0])
-        u = sp.symbols(f"u_:{dim}")
-        index_tuples = tuple(moments_up_to_component_order(2, dim=dim))
-        indices = exponents_to_polynomial_representations(index_tuples)
-        c_s_sq = sp.Rational(1, 3)
-        eq_moments1 = get_moments_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
-        eq_cumulants = get_cumulants_of_continuous_maxwellian_equilibrium(indices, dim=dim, u=u, c_s_sq=c_s_sq)
-        eq_cumulants = {idx: c for idx, c in zip(index_tuples, eq_cumulants)}
-        eq_moments2 = [raw_moment_as_function_of_cumulants(idx, eq_cumulants) for idx in index_tuples]
-        pdfs_to_moments = moment_matrix(indices, stencil)
-
-        def normalize(expressions):
-            return [remove_higher_order_terms(e.expand(), symbols=u, order=3) for e in expressions]
-
-        eq1 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments1))
-        eq2 = normalize(pdfs_to_moments.inv() * sp.Matrix(eq_moments2))
-        eq3 = normalize(discrete_maxwellian_equilibrium(stencil, order=3, c_s_sq=c_s_sq, compressible=True))
-        eq4 = normalize(generate_equilibrium_by_matching_moments(stencil, indices, c_s_sq=c_s_sq))
-
-        assert eq1 == eq2
-        assert eq2 == eq3
-        assert eq3 == eq4
--- a/lbmpy_tests/test_n_phase_boyer_noncoupled.ipynb
+++ b/lbmpy_tests/test_n_phase_boyer_noncoupled.ipynb
-%% Cell type:code id: tags:
-
-``` python
-import pytest
-pytest.importorskip('pycuda')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-from lbmpy.phasefield.contact_angle_circle_fitting import *
-from scipy.ndimage.filters import gaussian_filter
-from pystencils.simp import sympy_cse_on_assignment_list
-one = sp.sympify(1)
-
-import pyximport
-pyximport.install(language_level=3)
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-```
-
-%% Cell type:markdown id: tags:
-
-# Simulation arbitrary surface tension case
-
-%% Cell type:code id: tags:
-
-``` python
-n = 4
-dx, dt = 1, 1
-mobility = 2e-3
-domain_size = (150, 150)
-ε = one * 4
-penalty_factor = 0
-stabilization_factor = 10
-
-κ = (one,  one/2, one/3, one/4)
-sigma_factor = one / 15
-σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-#σ
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_target='gpu')
-c = dh.add_array('c', values_per_cell=n)
-c_tmp = dh.add_array_like('c_tmp', 'c')
-
-μ = dh.add_array('mu', values_per_cell=n)
-
-cvec = c.center_vector
-μvec = μ.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-α, _ = diffusion_coefficients(σ)
-
-f = lambda c: c**2 * ( 1 - c ) **2
-a, b = compute_ab(f)
-
-capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
-
-f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
-f_if = free_energy_interfacial(cvec, σ, a, ε)
-f = f_bulk + f_if
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#f_bulk
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_assignments = mu_kernel(f, cvec, c, μ)
-μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
-μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-def lapl(e):
-    return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-rhs = α * μvec
-discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
-                   for i, rhs_i in enumerate(rhs)]
-c_assignments = [Assignment(lhs, rhs)
-                 for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#c_assignments
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_sync = dh.synchronization_function(μ.name)
-c_sync = dh.synchronization_function(c.name)
-optimization = {'cpu_openmp': 4, 'cpu_vectorize_info': None}
-μ_kernel = create_kernel(μ_assignments, target=dh.default_target, **optimization).compile()
-c_kernel = create_kernel(c_assignments, target=dh.default_target, **optimization).compile()
-
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c.name]
-        arr[..., : ] = values
-
-def smooth():
-    for block in dh.iterate(ghost_layers=True):
-        c_arr = block[c.name]
-        for i in range(n):
-            gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
-
-def time_loop(steps):
-    dh.all_to_gpu()
-    for t in range(steps):
-        c_sync()
-        dh.run_kernel(μ_kernel)
-        μ_sync()
-        dh.run_kernel(c_kernel)
-        dh.swap(c.name, c_tmp.name)
-        #simplex_projection_2d(dh.cpu_arrays[c.name])
-    dh.all_to_cpu()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-set_c(make_slice[:, :], [0, 0, 0, 0])
-set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-smooth()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#dh.load_all('n_phases_state_size200_stab10.npz')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot(dh.gather_array(c.name))
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 146.44269023807928, \  117.81813928465394, \  95.73917047726678\right]$
-    [146.44269023807928, 117.81813928465394, 95.73917047726678]
-
-%% Cell type:code id: tags:
-
-``` python
-import time
-for i in range(10):
-    start = time.perf_counter()
-    time_loop(1_000)
-    end = time.perf_counter()
-
-    try:
-        print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
-    except Exception:
-        print(i, end - start, "none found")
-```
-
-%% Output
-
-    0 0.30607624700132874 [83.97888710273061, 100.48794346625529, 175.5331694310141]
-    1 0.2600655169990205 [83.73094685534376, 100.65854574856168, 175.6105073960945]
-    2 0.2601136189987301 [83.49914818603683, 100.82173327601079, 175.67911853795226]
-    3 0.25987518599868054 [83.31519592224448, 100.94468140501989, 175.74012267273554]
-    4 0.2651959220020217 [83.14239972296966, 101.06100094405181, 175.79659933297853]
-    5 0.25910847799968906 [82.984481834461, 101.16731750637399, 175.8482006591651]
-    6 0.259863024999504 [82.84781128433397, 101.2570276449976, 175.89516107066834]
-    7 0.2606479199966998 [82.7456965110742, 101.31687551766585, 175.93742797125986]
-    8 0.25991897900166805 [82.67010885583116, 101.35099855297112, 175.97889259119805]
-    9 0.2590353729974595 [75.9000280154447, 108.9652166787719, 175.1347553057833]
-
-%% Cell type:code id: tags:
-
-``` python
-plt.subplot(1,3,1)
-t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
-plt.plot(t);
-
-plt.subplot(1,3,2)
-plt.phase_plot(dh.gather_array(c.name), linewidth=1)
-
-plt.subplot(1,3,3)
-plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert not np.isnan(dh.max(c.name))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
-plt.hlines(0.5, 0, 30)
-plt.plot(t);
-```
-
-%% Output
-
-
-%% Cell type:code id: tags:
-
-``` python
-import pytest
-pytest.importorskip('pycuda')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.n_phase_boyer import *
-from lbmpy.phasefield.kerneleqs import *
-from lbmpy.phasefield.contact_angle_circle_fitting import *
-from scipy.ndimage.filters import gaussian_filter
-from pystencils.simp import sympy_cse_on_assignment_list
-one = sp.sympify(1)
-
-import pyximport
-pyximport.install(language_level=3)
-from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
-```
-
-%% Cell type:markdown id: tags:
-
-# Simulation arbitrary surface tension case
-
-%% Cell type:code id: tags:
-
-``` python
-n = 4
-dx, dt = 1, 1
-mobility = 2e-3
-domain_size = (150, 150)
-ε = one * 4
-penalty_factor = 0
-stabilization_factor = 10
-
-κ = (one,  one/2, one/3, one/4)
-sigma_factor = one / 15
-σ = sp.ImmutableDenseMatrix(n, n, lambda i,j: sigma_factor* (κ[i] + κ[j]) if i != j else 0 )
-#σ
-```
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_target='gpu')
-c = dh.add_array('c', values_per_cell=n)
-c_tmp = dh.add_array_like('c_tmp', 'c')
-
-μ = dh.add_array('mu', values_per_cell=n)
-
-cvec = c.center_vector
-μvec = μ.center_vector
-```
-
-%% Cell type:code id: tags:
-
-``` python
-α, _ = diffusion_coefficients(σ)
-
-f = lambda c: c**2 * ( 1 - c ) **2
-a, b = compute_ab(f)
-
-capital_f = capital_f0(cvec, σ) + correction_g(cvec, σ) + stabilization_factor * stabilization_term(cvec, α)
-
-f_bulk = free_energy_bulk(capital_f, b, ε) + penalty_factor * (one - sum(cvec))
-f_if = free_energy_interfacial(cvec, σ, a, ε)
-f = f_bulk + f_if
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#f_bulk
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_assignments = mu_kernel(f, cvec, c, μ)
-μ_assignments = [Assignment(a.lhs, a.rhs.doit()) for a in μ_assignments]
-μ_assignments = sympy_cse_on_assignment_list(μ_assignments)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-discretize = fd.Discretization2ndOrder(dx=dx, dt=dt)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-def lapl(e):
-    return sum(ps.fd.diff(e, d, d) for d in range(dh.dim))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-rhs = α * μvec
-discretized_rhs = [discretize(fd.expand_diff_full( lapl(mobility * rhs_i) + fd.transient(cvec[i], idx=i), functions=μvec))
-                   for i, rhs_i in enumerate(rhs)]
-c_assignments = [Assignment(lhs, rhs)
-                 for lhs, rhs in zip(c_tmp.center_vector, discretized_rhs)]
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#c_assignments
-```
-
-%% Cell type:code id: tags:
-
-``` python
-μ_sync = dh.synchronization_function(μ.name)
-c_sync = dh.synchronization_function(c.name)
-optimization = {'cpu_openmp': 4, 'cpu_vectorize_info': None}
-μ_kernel = create_kernel(μ_assignments, target=dh.default_target, **optimization).compile()
-c_kernel = create_kernel(c_assignments, target=dh.default_target, **optimization).compile()
-
-def set_c(slice_obj, values):
-    for block in dh.iterate(slice_obj):
-        arr = block[c.name]
-        arr[..., : ] = values
-
-def smooth():
-    for block in dh.iterate(ghost_layers=True):
-        c_arr = block[c.name]
-        for i in range(n):
-            gaussian_filter(c_arr[..., i], sigma=2, output=c_arr[..., i])
-
-def time_loop(steps):
-    dh.all_to_gpu()
-    for t in range(steps):
-        c_sync()
-        dh.run_kernel(μ_kernel)
-        μ_sync()
-        dh.run_kernel(c_kernel)
-        dh.swap(c.name, c_tmp.name)
-        #simplex_projection_2d(dh.cpu_arrays[c.name])
-    dh.all_to_cpu()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-set_c(make_slice[:, :], [0, 0, 0, 0])
-set_c(make_slice[:, 0.5:], [1, 0, 0, 0])
-set_c(make_slice[:, :0.5], [0, 1, 0, 0])
-set_c(make_slice[0.3:0.7, 0.3:0.7], [0, 0, 1, 0])
-smooth()
-```
-
-%% Cell type:code id: tags:
-
-``` python
-#dh.load_all('n_phases_state_size200_stab10.npz')
-```
-
-%% Cell type:code id: tags:
-
-``` python
-plt.phase_plot(dh.gather_array(c.name))
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-neumann_angles_from_surface_tensions(lambda i, j: float(σ[i, j]))
-```
-
-%% Output
-
-
-    $\displaystyle \left[ 146.44269023807928, \  117.81813928465394, \  95.73917047726678\right]$
-    [146.44269023807928, 117.81813928465394, 95.73917047726678]
-
-%% Cell type:code id: tags:
-
-``` python
-import time
-for i in range(10):
-    start = time.perf_counter()
-    time_loop(1_000)
-    end = time.perf_counter()
-
-    try:
-        print(i, end - start, liquid_lens_neumann_angles(dh.gather_array(c.name)))
-    except Exception:
-        print(i, end - start, "none found")
-```
-
-%% Output
-
-    0 0.30607624700132874 [83.97888710273061, 100.48794346625529, 175.5331694310141]
-    1 0.2600655169990205 [83.73094685534376, 100.65854574856168, 175.6105073960945]
-    2 0.2601136189987301 [83.49914818603683, 100.82173327601079, 175.67911853795226]
-    3 0.25987518599868054 [83.31519592224448, 100.94468140501989, 175.74012267273554]
-    4 0.2651959220020217 [83.14239972296966, 101.06100094405181, 175.79659933297853]
-    5 0.25910847799968906 [82.984481834461, 101.16731750637399, 175.8482006591651]
-    6 0.259863024999504 [82.84781128433397, 101.2570276449976, 175.89516107066834]
-    7 0.2606479199966998 [82.7456965110742, 101.31687551766585, 175.93742797125986]
-    8 0.25991897900166805 [82.67010885583116, 101.35099855297112, 175.97889259119805]
-    9 0.2590353729974595 [75.9000280154447, 108.9652166787719, 175.1347553057833]
-
-%% Cell type:code id: tags:
-
-``` python
-plt.subplot(1,3,1)
-t = dh.gather_array(c.name, make_slice[25, :]).squeeze()
-plt.plot(t);
-
-plt.subplot(1,3,2)
-plt.phase_plot(dh.gather_array(c.name), linewidth=1)
-
-plt.subplot(1,3,3)
-plt.scalar_field(dh.gather_array(μ.name)[:, :, 2])
-plt.colorbar();
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-assert not np.isnan(dh.max(c.name))
-```
-
-%% Cell type:code id: tags:
-
-``` python
-t = dh.gather_array(c.name, make_slice[25, 55:90]).squeeze()
-plt.hlines(0.5, 0, 30)
-plt.plot(t);
-```
-
-%% Output
-
-
No results found