890811d8 · 890811d8 · 890811d8 · 890811d8 · 890811d8 · 890811d8
--- a/lbmpy_tests/phasefield/test_analytical_3phase_model.ipynb
+++ b/lbmpy_tests/phasefield/test_analytical_3phase_model.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from pystencils.fd import evaluate_diffs
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for 3-Phase model
-
-Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
-
-%% Cell type:code id: tags:
-
-``` python
-order_parameters = sp.symbols("rho phi psi")
-rho, phi, psi = order_parameters
-F, _ = free_energy_functional_3_phases(include_bulk=True,
-                                       include_interface=True,
-                                       transformed=True,
-                                       expand_derivatives=False)
-F
-```
-
-%% Output
-
-    $$\frac{\alpha^{2} \kappa_{0}}{2} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{1}}{2} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{2}}{2} {\partial \psi}^{2} + \frac{\kappa_{0}}{2} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{1}}{2} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{2} \psi^{2}}{2} \left(- \psi + 1\right)^{2}$$
-    
-    
-     2                            2    2                             2    2
-    α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2)    α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2)    α ⋅κ₂⋅D(p
-    ─────────────────────────────── + ──────────────────────────────── + ─────────
-                   2                                 2                         2
-    
-                         2                  2                   2                2
-              ⎛φ   ψ   ρ⎞  ⎛  φ   ψ   ρ    ⎞       ⎛  φ   ψ   ρ⎞  ⎛φ   ψ   ρ    ⎞
-       2   κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟    κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
-    si)       ⎝2   2   2⎠  ⎝  2   2   2    ⎠       ⎝  2   2   2⎠  ⎝2   2   2    ⎠
-    ──── + ────────────────────────────────── + ──────────────────────────────────
-                           2                                    2
-    
-    
-    
-           2         2
-       κ₂⋅ψ ⋅(-ψ + 1)
-     + ───────────────
-              2
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-%% Cell type:markdown id: tags:
-
-Automatically deriving chemical potential as functional derivative of free energy
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$- \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \rho}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \rho}} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0}}{8} \phi^{2} \psi + \frac{3 \kappa_{0}}{8} \phi^{2} \rho - \frac{3 \kappa_{0}}{8} \phi^{2} + \frac{3 \kappa_{0}}{8} \phi \psi^{2} - \frac{3 \kappa_{0}}{4} \phi \psi \rho + \frac{3 \kappa_{0}}{4} \phi \psi + \frac{3 \kappa_{0}}{8} \phi \rho^{2} - \frac{3 \kappa_{0}}{4} \phi \rho + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0}}{8} \psi^{2} \rho - \frac{3 \kappa_{0}}{8} \psi^{2} - \frac{3 \kappa_{0}}{8} \psi \rho^{2} + \frac{3 \kappa_{0}}{4} \psi \rho - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0}}{8} \rho^{2} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1}}{8} \phi^{2} \psi + \frac{3 \kappa_{1}}{8} \phi^{2} \rho - \frac{3 \kappa_{1}}{8} \phi^{2} - \frac{3 \kappa_{1}}{8} \phi \psi^{2} + \frac{3 \kappa_{1}}{4} \phi \psi \rho - \frac{3 \kappa_{1}}{4} \phi \psi - \frac{3 \kappa_{1}}{8} \phi \rho^{2} + \frac{3 \kappa_{1}}{4} \phi \rho - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1}}{8} \psi^{2} \rho - \frac{3 \kappa_{1}}{8} \psi^{2} - \frac{3 \kappa_{1}}{8} \psi \rho^{2} + \frac{3 \kappa_{1}}{4} \psi \rho - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1}}{8} \rho^{2} + \frac{\kappa_{1} \rho}{4}$$
-       2                 2                 2                 2                 2
-      α ⋅κ₀⋅D(D(phi))   α ⋅κ₀⋅D(D(psi))   α ⋅κ₀⋅D(D(rho))   α ⋅κ₁⋅D(D(phi))   α ⋅κ
-    - ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
-             4                 4                 4                 4
-    
-                   2                    3         2           2           2
-    ₁⋅D(D(psi))   α ⋅κ₁⋅D(D(rho))   κ₀⋅φ    3⋅κ₀⋅φ ⋅ψ   3⋅κ₀⋅φ ⋅ρ   3⋅κ₀⋅φ    3⋅κ₀
-    ─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
-       4                 4            8         8           8          8
-    
-        2                                   2                         3         2
-    ⋅φ⋅ψ    3⋅κ₀⋅φ⋅ψ⋅ρ   3⋅κ₀⋅φ⋅ψ   3⋅κ₀⋅φ⋅ρ    3⋅κ₀⋅φ⋅ρ   κ₀⋅φ   κ₀⋅ψ    3⋅κ₀⋅ψ ⋅
-    ───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
-    8           4           4           8          4        4       8         8
-    
-              2           2                         3         2              3
-    ρ   3⋅κ₀⋅ψ    3⋅κ₀⋅ψ⋅ρ    3⋅κ₀⋅ψ⋅ρ   κ₀⋅ψ   κ₀⋅ρ    3⋅κ₀⋅ρ    κ₀⋅ρ   κ₁⋅φ    3
-    ─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
-           8          8          4        4       8        8       4       8
-    
-         2           2           2           2                                   2
-    ⋅κ₁⋅φ ⋅ψ   3⋅κ₁⋅φ ⋅ρ   3⋅κ₁⋅φ    3⋅κ₁⋅φ⋅ψ    3⋅κ₁⋅φ⋅ψ⋅ρ   3⋅κ₁⋅φ⋅ψ   3⋅κ₁⋅φ⋅ρ
-    ──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
-       8           8          8          8           4           4           8
-    
-                             3         2           2           2
-       3⋅κ₁⋅φ⋅ρ   κ₁⋅φ   κ₁⋅ψ    3⋅κ₁⋅ψ ⋅ρ   3⋅κ₁⋅ψ    3⋅κ₁⋅ψ⋅ρ    3⋅κ₁⋅ψ⋅ρ   κ₁⋅ψ
-     + ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
-          4        4       8         8          8          8          4        4
-    
-           3         2
-       κ₁⋅ρ    3⋅κ₁⋅ρ    κ₁⋅ρ
-     + ───── - ─────── + ────
-         8        8       4
-
-%% Cell type:markdown id: tags:
-
-Checking if expected profile is a solution of the differential equation:
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expectedProfile = analytic_interface_profile(x)
-expectedProfile
-```
-
-%% 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:
-
-Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
-
-%% Cell type:code id: tags:
-
-``` python
-for eq in mu_diff_eq:
-    eq = eq.subs({rho: 1,
-                  phi: 1 - expectedProfile,
-                  psi: expectedProfile})
-    eq = evaluate_diffs(eq, x).expand()
-    assert eq == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-F = expand_diff_linear(F, functions=order_parameters)  # expand derivatives using product rule
-two_phase_free_energy = F.subs({rho: 1,
-                                phi: 1 - expectedProfile,
-                                psi: expectedProfile})
-
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-two_phase_free_energy
-```
-
-%% Output
-
-
-    $$\frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}$$
-       κ₀ + κ₂
-    ─────────────
-           4⎛ x ⎞
-    16⋅cosh ⎜───⎟
-            ⎝2⋅α⎠
-
-%% Cell type:code id: tags:
-
-``` python
-gamma = cosh_integral(two_phase_free_energy, x)
-gamma
-```
-
-%% Output
-
-
-    $$\frac{\alpha}{6} \left(\kappa_{0} + \kappa_{2}\right)$$
-    α⋅(κ₀ + κ₂)
-    ───────────
-         6
-
-%% Cell type:code id: tags:
-
-``` python
-alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
-assert gamma == alpha/6 * (k0 + k2)
-```
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from pystencils.fd import evaluate_diffs
-```
-
-%% Cell type:markdown id: tags:
-
-# Analytical checks for 3-Phase model
-
-Here you can inspect the components of the free energy. View only bulk or interface contributions, before and after transformation from $U \rightarrow (\rho, \phi,\psi)$:
-
-%% Cell type:code id: tags:
-
-``` python
-order_parameters = sp.symbols("rho phi psi")
-rho, phi, psi = order_parameters
-F, _ = free_energy_functional_3_phases(include_bulk=True,
-                                       include_interface=True,
-                                       transformed=True,
-                                       expand_derivatives=False)
-F
-```
-
-%% Output
-
-    $$\frac{\alpha^{2} \kappa_{0}}{2} {\partial (\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{1}}{2} {\partial (- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}) }^{2} + \frac{\alpha^{2} \kappa_{2}}{2} {\partial \psi}^{2} + \frac{\kappa_{0}}{2} \left(\frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(- \frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{1}}{2} \left(- \frac{\phi}{2} - \frac{\psi}{2} + \frac{\rho}{2}\right)^{2} \left(\frac{\phi}{2} + \frac{\psi}{2} - \frac{\rho}{2} + 1\right)^{2} + \frac{\kappa_{2} \psi^{2}}{2} \left(- \psi + 1\right)^{2}$$
-    
-    
-     2                            2    2                             2    2
-    α ⋅κ₀⋅D(phi/2 - psi/2 + rho/2)    α ⋅κ₁⋅D(-phi/2 - psi/2 + rho/2)    α ⋅κ₂⋅D(p
-    ─────────────────────────────── + ──────────────────────────────── + ─────────
-                   2                                 2                         2
-    
-                         2                  2                   2                2
-              ⎛φ   ψ   ρ⎞  ⎛  φ   ψ   ρ    ⎞       ⎛  φ   ψ   ρ⎞  ⎛φ   ψ   ρ    ⎞
-       2   κ₀⋅⎜─ - ─ + ─⎟ ⋅⎜- ─ + ─ - ─ + 1⎟    κ₁⋅⎜- ─ - ─ + ─⎟ ⋅⎜─ + ─ - ─ + 1⎟
-    si)       ⎝2   2   2⎠  ⎝  2   2   2    ⎠       ⎝  2   2   2⎠  ⎝2   2   2    ⎠
-    ──── + ────────────────────────────────── + ──────────────────────────────────
-                           2                                    2
-    
-    
-    
-           2         2
-       κ₂⋅ψ ⋅(-ψ + 1)
-     + ───────────────
-              2
-
-%% Cell type:markdown id: tags:
-
-### Analytically checking the phase transition profile
-
-%% Cell type:markdown id: tags:
-
-Automatically deriving chemical potential as functional derivative of free energy
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-
-mu_diff_eq = chemical_potentials_from_free_energy(F, order_parameters)
-mu_diff_eq[0]
-```
-
-%% Output
-
-    $$- \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{0}}{4} {\partial {\partial \rho}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \phi}} + \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \psi}} - \frac{\alpha^{2} \kappa_{1}}{4} {\partial {\partial \rho}} + \frac{\kappa_{0} \phi^{3}}{8} - \frac{3 \kappa_{0}}{8} \phi^{2} \psi + \frac{3 \kappa_{0}}{8} \phi^{2} \rho - \frac{3 \kappa_{0}}{8} \phi^{2} + \frac{3 \kappa_{0}}{8} \phi \psi^{2} - \frac{3 \kappa_{0}}{4} \phi \psi \rho + \frac{3 \kappa_{0}}{4} \phi \psi + \frac{3 \kappa_{0}}{8} \phi \rho^{2} - \frac{3 \kappa_{0}}{4} \phi \rho + \frac{\kappa_{0} \phi}{4} - \frac{\kappa_{0} \psi^{3}}{8} + \frac{3 \kappa_{0}}{8} \psi^{2} \rho - \frac{3 \kappa_{0}}{8} \psi^{2} - \frac{3 \kappa_{0}}{8} \psi \rho^{2} + \frac{3 \kappa_{0}}{4} \psi \rho - \frac{\kappa_{0} \psi}{4} + \frac{\kappa_{0} \rho^{3}}{8} - \frac{3 \kappa_{0}}{8} \rho^{2} + \frac{\kappa_{0} \rho}{4} - \frac{\kappa_{1} \phi^{3}}{8} - \frac{3 \kappa_{1}}{8} \phi^{2} \psi + \frac{3 \kappa_{1}}{8} \phi^{2} \rho - \frac{3 \kappa_{1}}{8} \phi^{2} - \frac{3 \kappa_{1}}{8} \phi \psi^{2} + \frac{3 \kappa_{1}}{4} \phi \psi \rho - \frac{3 \kappa_{1}}{4} \phi \psi - \frac{3 \kappa_{1}}{8} \phi \rho^{2} + \frac{3 \kappa_{1}}{4} \phi \rho - \frac{\kappa_{1} \phi}{4} - \frac{\kappa_{1} \psi^{3}}{8} + \frac{3 \kappa_{1}}{8} \psi^{2} \rho - \frac{3 \kappa_{1}}{8} \psi^{2} - \frac{3 \kappa_{1}}{8} \psi \rho^{2} + \frac{3 \kappa_{1}}{4} \psi \rho - \frac{\kappa_{1} \psi}{4} + \frac{\kappa_{1} \rho^{3}}{8} - \frac{3 \kappa_{1}}{8} \rho^{2} + \frac{\kappa_{1} \rho}{4}$$
-       2                 2                 2                 2                 2
-      α ⋅κ₀⋅D(D(phi))   α ⋅κ₀⋅D(D(psi))   α ⋅κ₀⋅D(D(rho))   α ⋅κ₁⋅D(D(phi))   α ⋅κ
-    - ─────────────── + ─────────────── - ─────────────── + ─────────────── + ────
-             4                 4                 4                 4
-    
-                   2                    3         2           2           2
-    ₁⋅D(D(psi))   α ⋅κ₁⋅D(D(rho))   κ₀⋅φ    3⋅κ₀⋅φ ⋅ψ   3⋅κ₀⋅φ ⋅ρ   3⋅κ₀⋅φ    3⋅κ₀
-    ─────────── - ─────────────── + ───── - ───────── + ───────── - ─────── + ────
-       4                 4            8         8           8          8
-    
-        2                                   2                         3         2
-    ⋅φ⋅ψ    3⋅κ₀⋅φ⋅ψ⋅ρ   3⋅κ₀⋅φ⋅ψ   3⋅κ₀⋅φ⋅ρ    3⋅κ₀⋅φ⋅ρ   κ₀⋅φ   κ₀⋅ψ    3⋅κ₀⋅ψ ⋅
-    ───── - ────────── + ──────── + ───────── - ──────── + ──── - ───── + ────────
-    8           4           4           8          4        4       8         8
-    
-              2           2                         3         2              3
-    ρ   3⋅κ₀⋅ψ    3⋅κ₀⋅ψ⋅ρ    3⋅κ₀⋅ψ⋅ρ   κ₀⋅ψ   κ₀⋅ρ    3⋅κ₀⋅ρ    κ₀⋅ρ   κ₁⋅φ    3
-    ─ - ─────── - ───────── + ──────── - ──── + ───── - ─────── + ──── - ───── - ─
-           8          8          4        4       8        8       4       8
-    
-         2           2           2           2                                   2
-    ⋅κ₁⋅φ ⋅ψ   3⋅κ₁⋅φ ⋅ρ   3⋅κ₁⋅φ    3⋅κ₁⋅φ⋅ψ    3⋅κ₁⋅φ⋅ψ⋅ρ   3⋅κ₁⋅φ⋅ψ   3⋅κ₁⋅φ⋅ρ
-    ──────── + ───────── - ─────── - ───────── + ────────── - ──────── - ─────────
-       8           8          8          8           4           4           8
-    
-                             3         2           2           2
-       3⋅κ₁⋅φ⋅ρ   κ₁⋅φ   κ₁⋅ψ    3⋅κ₁⋅ψ ⋅ρ   3⋅κ₁⋅ψ    3⋅κ₁⋅ψ⋅ρ    3⋅κ₁⋅ψ⋅ρ   κ₁⋅ψ
-     + ──────── - ──── - ───── + ───────── - ─────── - ───────── + ──────── - ────
-          4        4       8         8          8          8          4        4
-    
-           3         2
-       κ₁⋅ρ    3⋅κ₁⋅ρ    κ₁⋅ρ
-     + ───── - ─────── + ────
-         8        8       4
-
-%% Cell type:markdown id: tags:
-
-Checking if expected profile is a solution of the differential equation:
-
-%% Cell type:code id: tags:
-
-``` python
-x = sp.symbols("x")
-expectedProfile = analytic_interface_profile(x)
-expectedProfile
-```
-
-%% 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:
-
-Checking a phase transition from $C_0$ to $C_2$. This means that $\rho=1$ while $phi$ and $psi$ are the analytical profile or 1-analytical profile
-
-%% Cell type:code id: tags:
-
-``` python
-for eq in mu_diff_eq:
-    eq = eq.subs({rho: 1,
-                  phi: 1 - expectedProfile,
-                  psi: expectedProfile})
-    eq = evaluate_diffs(eq, x).expand()
-    assert eq == 0
-```
-
-%% Cell type:markdown id: tags:
-
-### Checking the surface tensions parameters
-
-%% Cell type:code id: tags:
-
-``` python
-F, _ = free_energy_functional_3_phases(order_parameters)
-F = expand_diff_linear(F, functions=order_parameters)  # expand derivatives using product rule
-two_phase_free_energy = F.subs({rho: 1,
-                                phi: 1 - expectedProfile,
-                                psi: expectedProfile})
-
-two_phase_free_energy = sp.simplify(evaluate_diffs(two_phase_free_energy, x))
-two_phase_free_energy
-```
-
-%% Output
-
-
-    $$\frac{\kappa_{0} + \kappa_{2}}{16 \cosh^{4}{\left (\frac{x}{2 \alpha} \right )}}$$
-       κ₀ + κ₂
-    ─────────────
-           4⎛ x ⎞
-    16⋅cosh ⎜───⎟
-            ⎝2⋅α⎠
-
-%% Cell type:code id: tags:
-
-``` python
-gamma = cosh_integral(two_phase_free_energy, x)
-gamma
-```
-
-%% Output
-
-
-    $$\frac{\alpha}{6} \left(\kappa_{0} + \kappa_{2}\right)$$
-    α⋅(κ₀ + κ₂)
-    ───────────
-         6
-
-%% Cell type:code id: tags:
-
-``` python
-alpha, k0, k2 = sp.symbols("alpha, kappa_0, kappa_2")
-assert gamma == alpha/6 * (k0 + k2)
-```
--- 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_entropic_model.ipynb
+++ b/lbmpy_tests/phasefield/test_entropic_model.ipynb
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from lbmpy.phasefield.eos import *
-from lbmpy.phasefield.high_density_ratio_model import *
-from pystencils.fd.derivative import replace_generic_laplacian
-from pystencils.fd.spatial import discretize_spatial, \
-        fd_stencils_standard, fd_stencils_isotropic, fd_stencils_isotropic_high_density_code
-from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
-from lbmpy.forcemodels import EDM
-from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
-from scipy.ndimage.filters import gaussian_filter
-```
-
-%% Cell type:markdown id: tags:
-
-# Implementation of high density difference model
-
-According to *"Ternary free-energy entropic lattice Boltzmann model with high density ratio" by Wöhrwag, Semprebon, Mazloomi, Karlin and Kusumaatmaja*
-
-Up front we define all necessary parameters in one place:
-
-%% Cell type:code id: tags:
-
-``` python
-a = 0.037
-b = 0.2
-reduced_temperature = 0.61
-gas_constant = 1
-κ = (0.01, 1, 1)
-λ = (0.6, 1, 1)
-χ = 5
-φ_relaxation_rate = 1
-ρ_relaxation_rate = 1
-external_force = (0, 0)
-clipping = False
-
-domain_size = (100, 100)
-stencil = get_stencil('D2Q9', ordering='uk')
-
-fd_discretization = fd_stencils_isotropic_high_density_code
-target = 'cpu'
-threads = 4
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 1: Free Energy
-
-The free energy of this model contains a term that is derived from an equation of state. The equation of state and its parametrization determines the density of the liquid and gas phase.
-
-Here we use the Carnahan-Starling EOS, with the parametrization from the paper
-
-%% Cell type:code id: tags:
-
-``` python
-ρ, φ, c_l1, c_l2 = sp.symbols("rho, phi, c_l1, c_l2")
-
-critical_temperature = carnahan_starling_critical_temperature(a, b, gas_constant)
-temperature = reduced_temperature * critical_temperature
-
-eos = carnahan_starling_eos(ρ, gas_constant, temperature, a, b)
-eos
-```
-
-%% Output
-
-
-    $$- 0.037 \rho^{2} + \frac{0.042578305 \rho \left(- 0.000125 \rho^{3} + 0.0025 \rho^{2} + 0.05 \rho + 1\right)}{\left(- 0.05 \rho + 1\right)^{3}}$$
-                               ⎛            3           2             ⎞
-             2   0.042578305⋅ρ⋅⎝- 0.000125⋅ρ  + 0.0025⋅ρ  + 0.05⋅ρ + 1⎠
-    - 0.037⋅ρ  + ──────────────────────────────────────────────────────
-                                                  3
-                                     (-0.05⋅ρ + 1)
-
-%% Cell type:markdown id: tags:
-
-Next, we use a function that determines the gas and liquid density, using the Maxwell construction rule. This function uses an iterative procedure, that terminates when the equal-area condition is fulfilled with a certain tolerance.
-
-%% Cell type:code id: tags:
-
-``` python
-ρ_g, ρ_l = maxwell_construction(eos, tolerance=1e-3)
-(ρ_g, ρ_l, ρ_l / ρ_g)
-```
-
-%% Output
-
-
-    $$\left ( 0.0695273860315274, \quad 8.02904149705209, \quad 115.480272671423\right )$$
-    (0.0695273860315274, 8.02904149705209, 115.480272671423)
-
-%% Cell type:markdown id: tags:
-
-With this information we call a function that assembles the full free energy density:
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy = free_energy_high_density_ratio(eos, ρ, ρ_g, ρ_l, c_l1, c_l2, λ, κ)
-free_energy
-```
-
-%% Output
-
-    $$0.5 c_{l1}^{2} \left(- c_{l1} + 1\right)^{2} + 0.5 c_{l2}^{2} \left(- c_{l2} + 1\right)^{2} + 0.3 \rho \left(- 0.037 \rho - \frac{1.0 \left(1.7031322 \rho - 51.0939660000001\right)}{1.0 \rho^{2} - 40.0 \rho + 400.0} + 0.042578305 \log{\left (1.0 \rho \right )} - 0.0530922164415325\right) + 0.5 {\partial c_{l1}}^{2} + 0.5 {\partial c_{l2}}^{2} + 0.005 {\partial \rho}^{2} + 0.00085206629489322$$
-           2           2          2           2         ⎛           1.0⋅(1.7031322
-    0.5⋅cₗ₁ ⋅(-cₗ₁ + 1)  + 0.5⋅cₗ₂ ⋅(-cₗ₂ + 1)  + 0.3⋅ρ⋅⎜-0.037⋅ρ - ──────────────
-                                                        ⎜                      2
-                                                        ⎝                 1.0⋅ρ  -
-    
-    ⋅ρ - 51.0939660000001)                                              ⎞
-    ────────────────────── + 0.042578305⋅log(1.0⋅ρ) - 0.0530922164415325⎟ + 0.5⋅D(
-                                                                        ⎟
-     40.0⋅ρ + 400.0                                                     ⎠
-    
-         2              2               2
-    c_l1)  + 0.5⋅D(c_l2)  + 0.005⋅D(rho)  + 0.00085206629489322
-    
-    
-
-%% Cell type:markdown id: tags:
-
-This is the free energy expressed in the order parameters $\rho, c_{l1}, c_{l2}$. Next we have to transform it into coordinates $\rho, \phi$.
-
-%% Cell type:code id: tags:
-
-``` python
-transformation_eqs = [ c_l1 - (1 + φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2,
-                       c_l2 - (1 - φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2]
-transform_forward_substitutions = sp.solve(transformation_eqs, [c_l1, c_l2])
-transform_backward_substitutions = sp.solve(transformation_eqs, [ρ, φ])
-```
-
-%% Cell type:markdown id: tags:
-
-To do the transformation, we use the substitutions dict.
-After the substitutions the differentials have to be expanded again.
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy_transformed = free_energy.subs(transform_forward_substitutions)
-free_energy_transformed = expand_diff_full(free_energy_transformed, functions=(ρ, φ))
-free_energy_transformed.atoms(sp.Symbol)
-```
-
-%% Output
-
-
-    $$\left\{\phi, \rho\right\}$$
-    {φ, ρ}
-
-%% Cell type:markdown id: tags:
-
-Now the free energy depends only on ρ and φ. This transformed form is later used to derive expressions for the chemical potential, pressure tensor and force computations.
-
-%% Cell type:markdown id: tags:
-
-## Part 2: Data setup
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_target=target)
-
-# Fields for order parameters
-ρ_field = dh.add_array("rho")
-φ_field = dh.add_array("phi")
-c_field = dh.add_array("c", values_per_cell=2)
-
-# Chemical potential, pressure tensor, forces and velocities
-μ_phi_field = dh.add_array("mu_phi", latex_name=r"\mu_{\phi}")
-pbs_field = dh.add_array("pbs")
-pressure_tensor_field = dh.add_array("p", len(symmetric_tensor_linearization(dh.dim)))
-force_field = dh.add_array("force", values_per_cell=dh.dim, latex_name="F")
-vel_field = dh.add_array("velocity", values_per_cell=dh.dim)
-
-# PDF fields for lattice Boltzmann schemes
-pdf_src_rho = dh.add_array("pdf_src_rho", values_per_cell=len(stencil))
-pdf_dst_rho = dh.add_array_like("pdf_dst_rho", "pdf_src_rho")
-
-pdf_src_phi = dh.add_array("pdf_src_phi", values_per_cell=len(stencil))
-pdf_dst_phi = dh.add_array_like("pdf_dst_phi", "pdf_src_phi")
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 3a: Compute kernels and time loop
-
-We define one function that takes an expression with derivative objects in it, substitutes the spatial derivatives with finite differences using the strategy defined in the `fd_discretization` function and compiles a kernel from it.
-
-%% Cell type:code id: tags:
-
-``` python
-def make_kernel(assignments):
-    # assignments may be using the symbols ρ and φ
-    # these is substituted with the access to the corresponding fields here
-    field_substitutions = {
-        ρ: ρ_field.center,
-        φ: φ_field.center
-    }
-
-    processed_assignments = []
-    for a in assignments:
-        new_rhs = a.rhs.subs(field_substitutions)
-
-        # ∂∂f representing the laplacian of f is replaced by the explicit carteisan form
-        # ∂_0 ∂_0 f + ∂_1 ∂_1 f     (example for 2D)
-        # otherwise the discretization would not do the correct thing
-        new_rhs = replace_generic_laplacian(new_rhs)
-
-        # Next the "∂" objects are replaced using finite differences
-        new_rhs = discretize_spatial(new_rhs, dx=1, stencil=fd_discretization)
-        processed_assignments.append(Assignment(a.lhs, new_rhs))
-
-    return create_kernel(processed_assignments, target=target, cpu_openmp=threads).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-#### Chemical Potential
-
-In the next cell the kernel to compute the chemical potential is created. First an analytic expression for μ is obtained using the free energy, which is then passed to the discretization function above to create a kernel from it. We only have to store the chemical potential of the φ coordinate explicitly, which enters the Cahn-Hilliard lattice Boltzmann for φ.
-
-%% Cell type:code id: tags:
-
-``` python
-μ_ρ, μ_φ = chemical_potentials_from_free_energy(free_energy_transformed,
-                                                order_parameters=(ρ, φ))
-μ_phi_assignment = Assignment(μ_phi_field.center, μ_φ)
-μ_kernel = make_kernel([μ_phi_assignment])
-μ_phi_assignment
-```
-
-%% Output
-
-    $${{\mu_{\phi}}_{(0,0)}} \leftarrow 0.0004 \phi^{3} + 0.000473530700220886 \phi \rho^{2} - 0.00760399528440326 \phi \rho + 0.0205263968409311 \phi - 0.02 {\partial {\partial \phi}}$$
-                     3                           2
-    μ_φ_C := 0.0004⋅φ  + 0.000473530700220886⋅φ⋅ρ  - 0.00760399528440326⋅φ⋅ρ + 0.0
-    
-    
-    205263968409311⋅φ - 0.02⋅D(D(phi))
-
-%% Cell type:markdown id: tags:
-
-#### Pressure tensor and force computation
-
-For the pressure tensor a trick for enhancing numerical stability is used: the bulk component is not stored directly in the pressure tensor field, but the related quantity called `pbs` is stored in a separate field.
-
-$ pbs = \sqrt{|ρ  c_s^2 - p_{bulk} |} $
-
-The force is then calculated as $ \nabla \cdot  P_{if} + 2  (\nabla pbs) pbs$
-
-In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression.
-
-%% Cell type:code id: tags:
-
-``` python
-# Bulk part
-pressure_assignments = [
-    Assignment(pbs_field.center,
-               pressure_tensor_bulk_sqrt_term(free_energy_transformed, (ρ, φ), ρ)),
-]
-
-# Interface part
-P_if = pressure_tensor_from_free_energy(free_energy_transformed, (ρ, φ),
-                                        dim=dh.dim, include_bulk=False)
-index_map = symmetric_tensor_linearization(dh.dim)
-
-pressure_assignments += [
-    Assignment(pressure_tensor_field(index_1d), P_if[index_2d])
-    for index_2d, index_1d in index_map.items()
-]
-pressure_kernel = make_kernel(pressure_assignments)
-
-
-# Force kernel
-pressure_tensor_sym = sp.Matrix(dh.dim, dh.dim,
-                                lambda i, j: pressure_tensor_field(index_map[i, j]
-                                                                   if i < j else index_map[j, i]))
-force_term = force_from_pressure_tensor(pressure_tensor_sym,
-                                        functions=[ρ, φ],
-                                        pbs=pbs_field.center)
-force_assignments = [
-    Assignment(force_field(i),
-               force_term[i] + external_force[i] *  ρ_field.center / ρ_l)
-    for i in range(dh.dim)
-]
-force_kernel = make_kernel(force_assignments)
-```
-
-%% Cell type:markdown id: tags:
-
-#### Lattice Boltzmann schemes for time evolution of ρ and φ
-
- ρ is handled by a normal LB method (compressible, entropic equilibrium)
- stream and collide are splitted into separate kernels
- macroscopic values are computed after the stream, but inside the stream kernel
- velocity field stores the velocity which was not corrected for the forces yet
- the φ collision kernel corrects the velocity itself, because then the updated forces are used for the correction. When u is computed, the updated forces are not computed yet
- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system
- exact difference method is used to couple the force into the ρ-LBM
-
-%% Cell type:markdown id: tags:
-
-The following cell handles the clipping of the order parameters:
-
-%% Cell type:code id: tags:
-
-``` python
-if clipping:
-    def clip(ac, symbol, min_value, max_value):
-        """Function to clip the value of a symbol which is on one of lhs of the assignments
-        in an assignment collection"""
-        assert symbol in ac.bound_symbols
-        for i in range(len(ac.subexpressions)):
-            a = ac.subexpressions[i]
-            if a.lhs == symbol:
-                new_assignment = Assignment(symbol, sp.Piecewise((max_value, a.rhs > max_value),
-                                                                 (min_value, a.rhs < min_value),
-                                                                 (a.rhs, True)))
-                ac.subexpressions[i] = new_assignment
-                break
-
-    # TODO: how can this 'densgin' be derived automatically?
-    tred = reduced_temperature
-    densgin = -67.098 \
-              + 549.69 * tred \
-              - 1850.6 * tred * tred \
-              + 3281 * tred * tred * tred \
-              - 3237.3 * tred * tred * tred * tred \
-              + 1687.6 * tred * tred * tred * tred * tred \
-              - 361.51 * tred * tred * tred * tred * tred * tred
-    ρ_clip_min, ρ_clip_max = densgin * 0.5, 1.2 * ρ_l
-    φ_clip_min, φ_clip_max = -χ * 1.5, χ * 1.5
-```
-
-%% Cell type:markdown id: tags:
-
-Next, the collide and stream kernels for the ρ lattice Boltzmann are created
-
-%% Cell type:code id: tags:
-
-``` python
-force_model = EDM(force_field.center_vector)
-
-ρ_lbm_params = {
-    'stencil': stencil,
-    'method' : 'trt-kbc-n2',
-    'compressible': True,
-    'relaxation_rate': ρ_relaxation_rate,
-    'optimization': {
-        'symbolic_field': pdf_src_rho,
-        'symbolic_temporary_field': pdf_dst_rho,
-        'openmp': threads,
-        'target': target
-    }
-}
-
-# Standard collision step, that does not compute ρ and u from pdfs, but reads
-# them from fields - this is necessary because ρ may have been clipped before
-# the velocity field is not force corrected, which is the correct for the EDM model
-# but might be wrong for other force models
-ρ_collide = create_lb_function(kernel_type='collide_only',
-                               density_input=ρ_field,
-                               force_model=force_model,
-                               velocity_input=vel_field,
-                               **ρ_lbm_params)
-
-
-# First the assignments are created, then the density is clipped
-# then a kernel is created from the clipped assignments
-ρ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
-                                    force_model=None, #save uncorrected velocity
-                                    output={'density': ρ_field,
-                                            'velocity': vel_field},
-                                    **ρ_lbm_params)
-
-if clipping:
-    clip(ρ_stream_ur, sp.Symbol("rho"), ρ_clip_min, ρ_clip_max)
-ρ_stream = create_lb_function(update_rule=ρ_stream_ur, **ρ_lbm_params)
-```
-
-%% Cell type:markdown id: tags:
-
-The φ lattice Boltzmann solve the Cahn-Hilliard equation.
-
-%% Cell type:code id: tags:
-
-``` python
-φ_lb_method = cahn_hilliard_lb_method(stencil=stencil,
-                                      mu=μ_phi_field.center,
-                                      relaxation_rate=φ_relaxation_rate,
-                                      gamma=1)
-φ_lbm_params = {
-    'lb_method' : φ_lb_method,
-    'compressible': True,
-    'optimization': {
-        'symbolic_field': pdf_src_phi,
-        'symbolic_temporary_field': pdf_dst_phi,
-        'openmp': threads,
-        'target': target
-    }
-}
-
-ch_method = cahn_hilliard_lb_method(stencil=stencil,
-                                    mu=μ_phi_field.center,
-                                    relaxation_rate=φ_relaxation_rate,
-                                    gamma=1)
-
-corrected_vel = vel_field.center_vector + sp.Matrix(force_model.macroscopic_velocity_shift(ρ_field.center))
-φ_collide = create_lb_function(kernel_type='collide_only',
-                               density_input=φ_field,
-                               velocity_input=corrected_vel,
-                               **φ_lbm_params)
-
-φ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
-                                    output={'density': φ_field},
-                                    **φ_lbm_params)
-if clipping:
-    clip(φ_stream_ur, sp.Symbol("rho"), φ_clip_min, φ_clip_max)
-φ_stream = create_lb_function(update_rule=φ_stream_ur, **φ_lbm_params)
-```
-
-%% Cell type:markdown id: tags:
-
-#### Time loop
-
-Now we can put all kernels together into a time loop function
-
-%% Cell type:code id: tags:
-
-``` python
-op_sync = dh.synchronization_function([ρ_field.name, φ_field.name])
-p_sync = dh.synchronization_function([pbs_field.name, pressure_tensor_field.name])
-pdf_sync = dh.synchronization_function([pdf_src_phi.name, pdf_src_rho.name])
-
-def time_loop(steps):
-    for t in range(steps):
-        op_sync()
-        dh.run_kernel(μ_kernel)
-        dh.run_kernel(pressure_kernel)
-
-        p_sync()
-        dh.run_kernel(force_kernel)
-
-        dh.run_kernel(ρ_collide)
-        dh.run_kernel(φ_collide)
-
-        pdf_sync()
-        dh.run_kernel(ρ_stream)
-        dh.run_kernel(φ_stream)
-        dh.swap(pdf_dst_phi.name, pdf_src_phi.name)
-        dh.swap(pdf_dst_rho.name, pdf_src_rho.name)
-    return dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 3b: Compiling getter & setter kernels
-
-%% Cell type:markdown id: tags:
-
-The setter kernel computes ρ, φ from C and sets the pdfs to equilibrium using the values in the order parameter and velocity fields.
-
-%% Cell type:code id: tags:
-
-``` python
-init_assignments = [
-    Assignment( φ_field.center, transform_backward_substitutions[φ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
-    Assignment( ρ_field.center, transform_backward_substitutions[ρ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
-]
-
-init_rho = pdf_initialization_assignments(ρ_collide.method,
-                                          density=ρ_field.center,
-                                          velocity=vel_field.center_vector,
-                                          pdfs=pdf_src_rho.center_vector)
-init_rho = init_rho.new_without_subexpressions()
-init_assignments += init_rho.all_assignments
-
-init_phi = pdf_initialization_assignments(φ_collide.method,
-                                          density=φ_field.center,
-                                          velocity=(0,0,0),
-                                          pdfs=pdf_src_phi.center_vector)
-init_phi = init_phi.new_without_subexpressions().new_with_substitutions({μ_phi_field.center: 0})
-init_assignments += init_phi.all_assignments
-
-init_pdfs_assignments = init_rho.all_assignments + init_phi.all_assignments
-
-init_pdfs_kernel = create_kernel(init_pdfs_assignments).compile()
-init_kernel = create_kernel(init_assignments).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 4: Geometry initialization & plotting
-
-%% Cell type:code id: tags:
-
-``` python
-def plot_status():
-    plt.figure(figsize=(25, 6))
-    plt.subplot(1, 4, 1)
-    ρ_arr = dh.cpu_arrays[ρ_field.name][1:-1, 1:-1]
-    plt.scalar_field(ρ_arr)
-    plt.title("ρ ({:.2f}, {:.2f})".format(np.min(ρ_arr), np.max(ρ_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 2)
-    φ_arr = dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
-    plt.scalar_field(φ_arr)
-    plt.title("φ ({:.2f}, {:.2f})".format(np.min(φ_arr), np.max(φ_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 3)
-    f_arr = dh.cpu_arrays[force_field.name][1:-1, 1:-1]
-    plt.vector_field_magnitude(f_arr)
-    plt.title("F ({:.2f}, {:.2f})".format(np.min(f_arr), np.max(f_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 4)
-    μ_arr = dh.cpu_arrays[μ_phi_field.name][1:-1, 1:-1]
-    plt.scalar_field(μ_arr)
-    plt.title("μ_φ ({:.2f}, {:.2f})".format(np.min(μ_arr), np.max(μ_arr)))
-    plt.colorbar()
-
-def init_drop():
-    radius = dh.shape[0] // 5
-    mid1 = [dh.shape[0] // 2, dh.shape[1] // 2]
-
-    for block in dh.iterate(ghost_layers=True):
-        x, y = block.midpoint_arrays
-        mask1 = (x - mid1[0]) ** 2 + (y - mid1[1])**2 < radius ** 2
-
-        block[force_field.name].fill(0)
-        block[vel_field.name].fill(0)
-        block[μ_phi_field.name].fill(0)
-
-        c_arr = block[c_field.name]
-        c_arr[:, :].fill(0.0)
-        c_arr[mask1, 0] = 1.0
-
-        gaussian_filter(c_arr[..., 0], sigma=3, output=c_arr[..., 0])
-        gaussian_filter(c_arr[..., 1], sigma=3, output=c_arr[..., 1])
-
-    dh.run_kernel(init_kernel)
-```
-
-%% Cell type:markdown id: tags:
-
-# Part 5: Putting it all together
-
-%% Cell type:code id: tags:
-
-``` python
-init_drop()
-time_loop(1000)
-plot_status()
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-ρ_arr = dh.gather_array(ρ_field.name)
-assert not np.isnan(np.min(ρ_arr))
-assert not np.isnan(np.max(ρ_arr))
-```
-%% Cell type:code id: tags:
-
-``` python
-from lbmpy.session import *
-from lbmpy.phasefield.analytical import *
-from lbmpy.phasefield.eos import *
-from lbmpy.phasefield.high_density_ratio_model import *
-from pystencils.fd.derivative import replace_generic_laplacian
-from pystencils.fd.spatial import discretize_spatial, \
-        fd_stencils_standard, fd_stencils_isotropic, fd_stencils_isotropic_high_density_code
-from lbmpy.phasefield.cahn_hilliard_lbm import cahn_hilliard_lb_method
-from lbmpy.forcemodels import EDM
-from lbmpy.macroscopic_value_kernels import pdf_initialization_assignments
-from scipy.ndimage.filters import gaussian_filter
-```
-
-%% Cell type:markdown id: tags:
-
-# Implementation of high density difference model
-
-According to *"Ternary free-energy entropic lattice Boltzmann model with high density ratio" by Wöhrwag, Semprebon, Mazloomi, Karlin and Kusumaatmaja*
-
-Up front we define all necessary parameters in one place:
-
-%% Cell type:code id: tags:
-
-``` python
-a = 0.037
-b = 0.2
-reduced_temperature = 0.61
-gas_constant = 1
-κ = (0.01, 1, 1)
-λ = (0.6, 1, 1)
-χ = 5
-φ_relaxation_rate = 1
-ρ_relaxation_rate = 1
-external_force = (0, 0)
-clipping = False
-
-domain_size = (100, 100)
-stencil = get_stencil('D2Q9', ordering='uk')
-
-fd_discretization = fd_stencils_isotropic_high_density_code
-target = 'cpu'
-threads = 4
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 1: Free Energy
-
-The free energy of this model contains a term that is derived from an equation of state. The equation of state and its parametrization determines the density of the liquid and gas phase.
-
-Here we use the Carnahan-Starling EOS, with the parametrization from the paper
-
-%% Cell type:code id: tags:
-
-``` python
-ρ, φ, c_l1, c_l2 = sp.symbols("rho, phi, c_l1, c_l2")
-
-critical_temperature = carnahan_starling_critical_temperature(a, b, gas_constant)
-temperature = reduced_temperature * critical_temperature
-
-eos = carnahan_starling_eos(ρ, gas_constant, temperature, a, b)
-eos
-```
-
-%% Output
-
-
-    $$- 0.037 \rho^{2} + \frac{0.042578305 \rho \left(- 0.000125 \rho^{3} + 0.0025 \rho^{2} + 0.05 \rho + 1\right)}{\left(- 0.05 \rho + 1\right)^{3}}$$
-                               ⎛            3           2             ⎞
-             2   0.042578305⋅ρ⋅⎝- 0.000125⋅ρ  + 0.0025⋅ρ  + 0.05⋅ρ + 1⎠
-    - 0.037⋅ρ  + ──────────────────────────────────────────────────────
-                                                  3
-                                     (-0.05⋅ρ + 1)
-
-%% Cell type:markdown id: tags:
-
-Next, we use a function that determines the gas and liquid density, using the Maxwell construction rule. This function uses an iterative procedure, that terminates when the equal-area condition is fulfilled with a certain tolerance.
-
-%% Cell type:code id: tags:
-
-``` python
-ρ_g, ρ_l = maxwell_construction(eos, tolerance=1e-3)
-(ρ_g, ρ_l, ρ_l / ρ_g)
-```
-
-%% Output
-
-
-    $$\left ( 0.0695273860315274, \quad 8.02904149705209, \quad 115.480272671423\right )$$
-    (0.0695273860315274, 8.02904149705209, 115.480272671423)
-
-%% Cell type:markdown id: tags:
-
-With this information we call a function that assembles the full free energy density:
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy = free_energy_high_density_ratio(eos, ρ, ρ_g, ρ_l, c_l1, c_l2, λ, κ)
-free_energy
-```
-
-%% Output
-
-    $$0.5 c_{l1}^{2} \left(- c_{l1} + 1\right)^{2} + 0.5 c_{l2}^{2} \left(- c_{l2} + 1\right)^{2} + 0.3 \rho \left(- 0.037 \rho - \frac{1.0 \left(1.7031322 \rho - 51.0939660000001\right)}{1.0 \rho^{2} - 40.0 \rho + 400.0} + 0.042578305 \log{\left (1.0 \rho \right )} - 0.0530922164415325\right) + 0.5 {\partial c_{l1}}^{2} + 0.5 {\partial c_{l2}}^{2} + 0.005 {\partial \rho}^{2} + 0.00085206629489322$$
-           2           2          2           2         ⎛           1.0⋅(1.7031322
-    0.5⋅cₗ₁ ⋅(-cₗ₁ + 1)  + 0.5⋅cₗ₂ ⋅(-cₗ₂ + 1)  + 0.3⋅ρ⋅⎜-0.037⋅ρ - ──────────────
-                                                        ⎜                      2
-                                                        ⎝                 1.0⋅ρ  -
-    
-    ⋅ρ - 51.0939660000001)                                              ⎞
-    ────────────────────── + 0.042578305⋅log(1.0⋅ρ) - 0.0530922164415325⎟ + 0.5⋅D(
-                                                                        ⎟
-     40.0⋅ρ + 400.0                                                     ⎠
-    
-         2              2               2
-    c_l1)  + 0.5⋅D(c_l2)  + 0.005⋅D(rho)  + 0.00085206629489322
-    
-    
-
-%% Cell type:markdown id: tags:
-
-This is the free energy expressed in the order parameters $\rho, c_{l1}, c_{l2}$. Next we have to transform it into coordinates $\rho, \phi$.
-
-%% Cell type:code id: tags:
-
-``` python
-transformation_eqs = [ c_l1 - (1 + φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2,
-                       c_l2 - (1 - φ/χ - (ρ - ρ_l)/(ρ_g - ρ_l)) / 2]
-transform_forward_substitutions = sp.solve(transformation_eqs, [c_l1, c_l2])
-transform_backward_substitutions = sp.solve(transformation_eqs, [ρ, φ])
-```
-
-%% Cell type:markdown id: tags:
-
-To do the transformation, we use the substitutions dict.
-After the substitutions the differentials have to be expanded again.
-
-%% Cell type:code id: tags:
-
-``` python
-free_energy_transformed = free_energy.subs(transform_forward_substitutions)
-free_energy_transformed = expand_diff_full(free_energy_transformed, functions=(ρ, φ))
-free_energy_transformed.atoms(sp.Symbol)
-```
-
-%% Output
-
-
-    $$\left\{\phi, \rho\right\}$$
-    {φ, ρ}
-
-%% Cell type:markdown id: tags:
-
-Now the free energy depends only on ρ and φ. This transformed form is later used to derive expressions for the chemical potential, pressure tensor and force computations.
-
-%% Cell type:markdown id: tags:
-
-## Part 2: Data setup
-
-%% Cell type:code id: tags:
-
-``` python
-dh = create_data_handling(domain_size, periodicity=True, default_target=target)
-
-# Fields for order parameters
-ρ_field = dh.add_array("rho")
-φ_field = dh.add_array("phi")
-c_field = dh.add_array("c", values_per_cell=2)
-
-# Chemical potential, pressure tensor, forces and velocities
-μ_phi_field = dh.add_array("mu_phi", latex_name=r"\mu_{\phi}")
-pbs_field = dh.add_array("pbs")
-pressure_tensor_field = dh.add_array("p", len(symmetric_tensor_linearization(dh.dim)))
-force_field = dh.add_array("force", values_per_cell=dh.dim, latex_name="F")
-vel_field = dh.add_array("velocity", values_per_cell=dh.dim)
-
-# PDF fields for lattice Boltzmann schemes
-pdf_src_rho = dh.add_array("pdf_src_rho", values_per_cell=len(stencil))
-pdf_dst_rho = dh.add_array_like("pdf_dst_rho", "pdf_src_rho")
-
-pdf_src_phi = dh.add_array("pdf_src_phi", values_per_cell=len(stencil))
-pdf_dst_phi = dh.add_array_like("pdf_dst_phi", "pdf_src_phi")
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 3a: Compute kernels and time loop
-
-We define one function that takes an expression with derivative objects in it, substitutes the spatial derivatives with finite differences using the strategy defined in the `fd_discretization` function and compiles a kernel from it.
-
-%% Cell type:code id: tags:
-
-``` python
-def make_kernel(assignments):
-    # assignments may be using the symbols ρ and φ
-    # these is substituted with the access to the corresponding fields here
-    field_substitutions = {
-        ρ: ρ_field.center,
-        φ: φ_field.center
-    }
-
-    processed_assignments = []
-    for a in assignments:
-        new_rhs = a.rhs.subs(field_substitutions)
-
-        # ∂∂f representing the laplacian of f is replaced by the explicit carteisan form
-        # ∂_0 ∂_0 f + ∂_1 ∂_1 f     (example for 2D)
-        # otherwise the discretization would not do the correct thing
-        new_rhs = replace_generic_laplacian(new_rhs)
-
-        # Next the "∂" objects are replaced using finite differences
-        new_rhs = discretize_spatial(new_rhs, dx=1, stencil=fd_discretization)
-        processed_assignments.append(Assignment(a.lhs, new_rhs))
-
-    return create_kernel(processed_assignments, target=target, cpu_openmp=threads).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-#### Chemical Potential
-
-In the next cell the kernel to compute the chemical potential is created. First an analytic expression for μ is obtained using the free energy, which is then passed to the discretization function above to create a kernel from it. We only have to store the chemical potential of the φ coordinate explicitly, which enters the Cahn-Hilliard lattice Boltzmann for φ.
-
-%% Cell type:code id: tags:
-
-``` python
-μ_ρ, μ_φ = chemical_potentials_from_free_energy(free_energy_transformed,
-                                                order_parameters=(ρ, φ))
-μ_phi_assignment = Assignment(μ_phi_field.center, μ_φ)
-μ_kernel = make_kernel([μ_phi_assignment])
-μ_phi_assignment
-```
-
-%% Output
-
-    $${{\mu_{\phi}}_{(0,0)}} \leftarrow 0.0004 \phi^{3} + 0.000473530700220886 \phi \rho^{2} - 0.00760399528440326 \phi \rho + 0.0205263968409311 \phi - 0.02 {\partial {\partial \phi}}$$
-                     3                           2
-    μ_φ_C := 0.0004⋅φ  + 0.000473530700220886⋅φ⋅ρ  - 0.00760399528440326⋅φ⋅ρ + 0.0
-    
-    
-    205263968409311⋅φ - 0.02⋅D(D(phi))
-
-%% Cell type:markdown id: tags:
-
-#### Pressure tensor and force computation
-
-For the pressure tensor a trick for enhancing numerical stability is used: the bulk component is not stored directly in the pressure tensor field, but the related quantity called `pbs` is stored in a separate field.
-
-$ pbs = \sqrt{|ρ  c_s^2 - p_{bulk} |} $
-
-The force is then calculated as $ \nabla \cdot  P_{if} + 2  (\nabla pbs) pbs$
-
-In the following kernel the pressure tensor field is filled with $P_{if}$ and the pbs field with above expression.
-
-%% Cell type:code id: tags:
-
-``` python
-# Bulk part
-pressure_assignments = [
-    Assignment(pbs_field.center,
-               pressure_tensor_bulk_sqrt_term(free_energy_transformed, (ρ, φ), ρ)),
-]
-
-# Interface part
-P_if = pressure_tensor_from_free_energy(free_energy_transformed, (ρ, φ),
-                                        dim=dh.dim, include_bulk=False)
-index_map = symmetric_tensor_linearization(dh.dim)
-
-pressure_assignments += [
-    Assignment(pressure_tensor_field(index_1d), P_if[index_2d])
-    for index_2d, index_1d in index_map.items()
-]
-pressure_kernel = make_kernel(pressure_assignments)
-
-
-# Force kernel
-pressure_tensor_sym = sp.Matrix(dh.dim, dh.dim,
-                                lambda i, j: pressure_tensor_field(index_map[i, j]
-                                                                   if i < j else index_map[j, i]))
-force_term = force_from_pressure_tensor(pressure_tensor_sym,
-                                        functions=[ρ, φ],
-                                        pbs=pbs_field.center)
-force_assignments = [
-    Assignment(force_field(i),
-               force_term[i] + external_force[i] *  ρ_field.center / ρ_l)
-    for i in range(dh.dim)
-]
-force_kernel = make_kernel(force_assignments)
-```
-
-%% Cell type:markdown id: tags:
-
-#### Lattice Boltzmann schemes for time evolution of ρ and φ
-
- ρ is handled by a normal LB method (compressible, entropic equilibrium)
- stream and collide are splitted into separate kernels
- macroscopic values are computed after the stream, but inside the stream kernel
- velocity field stores the velocity which was not corrected for the forces yet
- the φ collision kernel corrects the velocity itself, because then the updated forces are used for the correction. When u is computed, the updated forces are not computed yet
- when ρ and φ are updated, they are clipped to a valid region, this clipping should be only necessary during equilibration of the system
- exact difference method is used to couple the force into the ρ-LBM
-
-%% Cell type:markdown id: tags:
-
-The following cell handles the clipping of the order parameters:
-
-%% Cell type:code id: tags:
-
-``` python
-if clipping:
-    def clip(ac, symbol, min_value, max_value):
-        """Function to clip the value of a symbol which is on one of lhs of the assignments
-        in an assignment collection"""
-        assert symbol in ac.bound_symbols
-        for i in range(len(ac.subexpressions)):
-            a = ac.subexpressions[i]
-            if a.lhs == symbol:
-                new_assignment = Assignment(symbol, sp.Piecewise((max_value, a.rhs > max_value),
-                                                                 (min_value, a.rhs < min_value),
-                                                                 (a.rhs, True)))
-                ac.subexpressions[i] = new_assignment
-                break
-
-    # TODO: how can this 'densgin' be derived automatically?
-    tred = reduced_temperature
-    densgin = -67.098 \
-              + 549.69 * tred \
-              - 1850.6 * tred * tred \
-              + 3281 * tred * tred * tred \
-              - 3237.3 * tred * tred * tred * tred \
-              + 1687.6 * tred * tred * tred * tred * tred \
-              - 361.51 * tred * tred * tred * tred * tred * tred
-    ρ_clip_min, ρ_clip_max = densgin * 0.5, 1.2 * ρ_l
-    φ_clip_min, φ_clip_max = -χ * 1.5, χ * 1.5
-```
-
-%% Cell type:markdown id: tags:
-
-Next, the collide and stream kernels for the ρ lattice Boltzmann are created
-
-%% Cell type:code id: tags:
-
-``` python
-force_model = EDM(force_field.center_vector)
-
-ρ_lbm_params = {
-    'stencil': stencil,
-    'method' : 'trt-kbc-n2',
-    'compressible': True,
-    'relaxation_rate': ρ_relaxation_rate,
-    'optimization': {
-        'symbolic_field': pdf_src_rho,
-        'symbolic_temporary_field': pdf_dst_rho,
-        'openmp': threads,
-        'target': target
-    }
-}
-
-# Standard collision step, that does not compute ρ and u from pdfs, but reads
-# them from fields - this is necessary because ρ may have been clipped before
-# the velocity field is not force corrected, which is the correct for the EDM model
-# but might be wrong for other force models
-ρ_collide = create_lb_function(kernel_type='collide_only',
-                               density_input=ρ_field,
-                               force_model=force_model,
-                               velocity_input=vel_field,
-                               **ρ_lbm_params)
-
-
-# First the assignments are created, then the density is clipped
-# then a kernel is created from the clipped assignments
-ρ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
-                                    force_model=None, #save uncorrected velocity
-                                    output={'density': ρ_field,
-                                            'velocity': vel_field},
-                                    **ρ_lbm_params)
-
-if clipping:
-    clip(ρ_stream_ur, sp.Symbol("rho"), ρ_clip_min, ρ_clip_max)
-ρ_stream = create_lb_function(update_rule=ρ_stream_ur, **ρ_lbm_params)
-```
-
-%% Cell type:markdown id: tags:
-
-The φ lattice Boltzmann solve the Cahn-Hilliard equation.
-
-%% Cell type:code id: tags:
-
-``` python
-φ_lb_method = cahn_hilliard_lb_method(stencil=stencil,
-                                      mu=μ_phi_field.center,
-                                      relaxation_rate=φ_relaxation_rate,
-                                      gamma=1)
-φ_lbm_params = {
-    'lb_method' : φ_lb_method,
-    'compressible': True,
-    'optimization': {
-        'symbolic_field': pdf_src_phi,
-        'symbolic_temporary_field': pdf_dst_phi,
-        'openmp': threads,
-        'target': target
-    }
-}
-
-ch_method = cahn_hilliard_lb_method(stencil=stencil,
-                                    mu=μ_phi_field.center,
-                                    relaxation_rate=φ_relaxation_rate,
-                                    gamma=1)
-
-corrected_vel = vel_field.center_vector + sp.Matrix(force_model.macroscopic_velocity_shift(ρ_field.center))
-φ_collide = create_lb_function(kernel_type='collide_only',
-                               density_input=φ_field,
-                               velocity_input=corrected_vel,
-                               **φ_lbm_params)
-
-φ_stream_ur = create_lb_update_rule(kernel_type='stream_pull_only',
-                                    output={'density': φ_field},
-                                    **φ_lbm_params)
-if clipping:
-    clip(φ_stream_ur, sp.Symbol("rho"), φ_clip_min, φ_clip_max)
-φ_stream = create_lb_function(update_rule=φ_stream_ur, **φ_lbm_params)
-```
-
-%% Cell type:markdown id: tags:
-
-#### Time loop
-
-Now we can put all kernels together into a time loop function
-
-%% Cell type:code id: tags:
-
-``` python
-op_sync = dh.synchronization_function([ρ_field.name, φ_field.name])
-p_sync = dh.synchronization_function([pbs_field.name, pressure_tensor_field.name])
-pdf_sync = dh.synchronization_function([pdf_src_phi.name, pdf_src_rho.name])
-
-def time_loop(steps):
-    for t in range(steps):
-        op_sync()
-        dh.run_kernel(μ_kernel)
-        dh.run_kernel(pressure_kernel)
-
-        p_sync()
-        dh.run_kernel(force_kernel)
-
-        dh.run_kernel(ρ_collide)
-        dh.run_kernel(φ_collide)
-
-        pdf_sync()
-        dh.run_kernel(ρ_stream)
-        dh.run_kernel(φ_stream)
-        dh.swap(pdf_dst_phi.name, pdf_src_phi.name)
-        dh.swap(pdf_dst_rho.name, pdf_src_rho.name)
-    return dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 3b: Compiling getter & setter kernels
-
-%% Cell type:markdown id: tags:
-
-The setter kernel computes ρ, φ from C and sets the pdfs to equilibrium using the values in the order parameter and velocity fields.
-
-%% Cell type:code id: tags:
-
-``` python
-init_assignments = [
-    Assignment( φ_field.center, transform_backward_substitutions[φ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
-    Assignment( ρ_field.center, transform_backward_substitutions[ρ].subs({ c_l1: c_field(0), c_l2: c_field(1)} )),
-]
-
-init_rho = pdf_initialization_assignments(ρ_collide.method,
-                                          density=ρ_field.center,
-                                          velocity=vel_field.center_vector,
-                                          pdfs=pdf_src_rho.center_vector)
-init_rho = init_rho.new_without_subexpressions()
-init_assignments += init_rho.all_assignments
-
-init_phi = pdf_initialization_assignments(φ_collide.method,
-                                          density=φ_field.center,
-                                          velocity=(0,0,0),
-                                          pdfs=pdf_src_phi.center_vector)
-init_phi = init_phi.new_without_subexpressions().new_with_substitutions({μ_phi_field.center: 0})
-init_assignments += init_phi.all_assignments
-
-init_pdfs_assignments = init_rho.all_assignments + init_phi.all_assignments
-
-init_pdfs_kernel = create_kernel(init_pdfs_assignments).compile()
-init_kernel = create_kernel(init_assignments).compile()
-```
-
-%% Cell type:markdown id: tags:
-
-## Part 4: Geometry initialization & plotting
-
-%% Cell type:code id: tags:
-
-``` python
-def plot_status():
-    plt.figure(figsize=(25, 6))
-    plt.subplot(1, 4, 1)
-    ρ_arr = dh.cpu_arrays[ρ_field.name][1:-1, 1:-1]
-    plt.scalar_field(ρ_arr)
-    plt.title("ρ ({:.2f}, {:.2f})".format(np.min(ρ_arr), np.max(ρ_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 2)
-    φ_arr = dh.cpu_arrays[φ_field.name][1:-1, 1:-1]
-    plt.scalar_field(φ_arr)
-    plt.title("φ ({:.2f}, {:.2f})".format(np.min(φ_arr), np.max(φ_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 3)
-    f_arr = dh.cpu_arrays[force_field.name][1:-1, 1:-1]
-    plt.vector_field_magnitude(f_arr)
-    plt.title("F ({:.2f}, {:.2f})".format(np.min(f_arr), np.max(f_arr)))
-    plt.colorbar()
-
-    plt.subplot(1, 4, 4)
-    μ_arr = dh.cpu_arrays[μ_phi_field.name][1:-1, 1:-1]
-    plt.scalar_field(μ_arr)
-    plt.title("μ_φ ({:.2f}, {:.2f})".format(np.min(μ_arr), np.max(μ_arr)))
-    plt.colorbar()
-
-def init_drop():
-    radius = dh.shape[0] // 5
-    mid1 = [dh.shape[0] // 2, dh.shape[1] // 2]
-
-    for block in dh.iterate(ghost_layers=True):
-        x, y = block.midpoint_arrays
-        mask1 = (x - mid1[0]) ** 2 + (y - mid1[1])**2 < radius ** 2
-
-        block[force_field.name].fill(0)
-        block[vel_field.name].fill(0)
-        block[μ_phi_field.name].fill(0)
-
-        c_arr = block[c_field.name]
-        c_arr[:, :].fill(0.0)
-        c_arr[mask1, 0] = 1.0
-
-        gaussian_filter(c_arr[..., 0], sigma=3, output=c_arr[..., 0])
-        gaussian_filter(c_arr[..., 1], sigma=3, output=c_arr[..., 1])
-
-    dh.run_kernel(init_kernel)
-```
-
-%% Cell type:markdown id: tags:
-
-# Part 5: Putting it all together
-
-%% Cell type:code id: tags:
-
-``` python
-init_drop()
-time_loop(1000)
-plot_status()
-```
-
-%% Output
-
-
-
-%% Cell type:code id: tags:
-
-``` python
-ρ_arr = dh.gather_array(ρ_field.name)
-assert not np.isnan(np.min(ρ_arr))
-assert not np.isnan(np.max(ρ_arr))
-```
--- 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
-
-    $$0.005 c_{0}^{2} \left(- c_{0} + 1\right)^{2} + 0.005 c_{1}^{2} \left(- c_{1} + 1\right)^{2} + 0.005 c_{2}^{2} \left(- c_{2} + 1\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₀ ⋅(-c₀ + 1)  + 0.005⋅c₁ ⋅(-c₁ + 1)  + 0.005⋅c₂ ⋅(-c₂ + 1)  + 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(1000):
-    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
-
-
-    $$\left [ 107.7400812901783, \quad 99.91508694115207, \quad 152.34483176866968\right ]$$
-    [107.7400812901783, 99.91508694115207, 152.34483176866968]
-%% 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
-
-    $$0.005 c_{0}^{2} \left(- c_{0} + 1\right)^{2} + 0.005 c_{1}^{2} \left(- c_{1} + 1\right)^{2} + 0.005 c_{2}^{2} \left(- c_{2} + 1\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₀ ⋅(-c₀ + 1)  + 0.005⋅c₁ ⋅(-c₁ + 1)  + 0.005⋅c₂ ⋅(-c₂ + 1)  + 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(1000):
-    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
-
-
-    $$\left [ 107.7400812901783, \quad 99.91508694115207, \quad 152.34483176866968\right ]$$
-    [107.7400812901783, 99.91508694115207, 152.34483176866968]
--- 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
-from lbmpy.phasefield_allen_cahn.analytical import analytic_rising_speed
-from lbmpy.phasefield_allen_cahn.parameter_calculation import calculate_dimensionless_rising_bubble, \
-    calculate_parameters_rti
-from pystencils import fields, AssignmentCollection
-from pystencils.simp import sympy_cse
-
-from lbmpy.creationfunctions import create_lb_method, create_lb_update_rule
-from lbmpy.stencils import get_stencil
-from lbmpy.methods.creationfunctions import create_with_discrete_maxwellian_eq_moments
-
-from lbmpy.phasefield_allen_cahn.kernel_equations import initializer_kernel_hydro_lb, \
-    initializer_kernel_phase_field_lb, get_collision_assignments_hydro, interface_tracking_force, hydrodynamic_force
-from lbmpy.phasefield_allen_cahn.force_model import MultiphaseForceModel
-
-from collections import OrderedDict
-
-import numpy as np
-
-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')
-force = fields("force(" + str(dimensions) + "): [" + 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)
-
-mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, relaxation_rate, 1, 1, 1, 1])
-rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
-
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=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 = create_lb_method(stencil=stencil_phase,
-                                method='srt',
-                                relaxation_rate=w_c,
-                                compressible=True,
-                                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)
-allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
-
-# ---------------------------------------------------------------------------------------------------------
-
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
-
-hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                       density=density,
-                                                       velocity_input=u,
-                                                       force=force_g,
-                                                       optimization={"symbolic_field": g,
-                                                                     "symbolic_temporary_field": g_tmp},
-                                                       kernel_type='collide_only')
-
-# streaming of the hydrodynamic distribution
-stream_hydro = create_lb_update_rule(stencil=stencil_hydro,
-                                     optimization={"symbolic_field": g,
-                                                   "symbolic_temporary_field": g_tmp},
-                                     kernel_type='stream_pull_only')
-
--- a/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.ipynb
+++ b/lbmpy_tests/phasefield_allen_cahn/test_phase_field_allen_cahn.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
-Nx = L0
-Ny = 4 * L0
-X0 = (Nx/2) + 1
-Y0 = (Ny/2) + 1
-
-# time step
-tf = 10001
-reference_time = 16000
-
-# force acting on the bubble
-body_force = [0, 0, 0]
-rho_H = 1.0
-
-
-parameters = calculate_parameters_rti(reference_length=L0,
-                                      reference_time=reference_time,
-                                      density_heavy=rho_H,
-                                      capillary_number=0.26,
-                                      reynolds_number=3000,
-                                      atwood_number=0.5,
-                                      peclet_number=1000,
-                                      density_ratio=3,
-                                      viscosity_ratio=1)
-
-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")
-
-
-# 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))
-
-
-# fields
-domain_size = (Nx, Ny)
-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)
-
-force = dh.add_array("force",values_per_cell=dh.dim)
-dh.fill("force", 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[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)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-if len(stencil_hydro) == 9:
-    # for D2Q9 a self defined method is used
-    x, y, z = MOMENT_SYMBOLS
-    moment_table = [sp.Integer(1),
-                    -4 + 3*(x**2 + y**2),
-                    4 - sp.Rational(21, 2)*(x**2 + y**2) + sp.Rational(9, 2)*(x**2 + y**2)**2,
-                    x,
-                    (-5 + 3*(x**2 + y**2))*x,
-                    y,
-                    (-5 + 3*(x**2 + y**2))*y ,
-                    x**2 - y**2,
-                    x*y]
-
-    relax_table = [1, 1, 1, 1, 1, 1, 1, s8, s8]
-    rr_dict = OrderedDict(zip(moment_table, relax_table))
-elif len(stencil_hydro) == 19:
-    mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, 1, 1, s8, 1, 1])
-    rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
-else:
-    mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, s8, 1, 1, 1, 1])
-    rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
-
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# initialize the domain
-def Initialize_distributions():
-
-    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
-        x, y = block.midpoint_arrays
-        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)
-force_model_g = MultiphaseForceModel(force=force_g, rho=rho)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# create the allen cahl lattice boltzmann step for the h field as well as the hydrodynamic
-# lattice boltzmann step for the g field
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
-
-method_phase = create_lb_method(stencil=stencil_phase,
-                                method='srt',
-                                relaxation_rate=w_c,
-                                compressible = True,
-                                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}})
-
-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()
-#------------------------------------------------------------------------------------------------------------------------
-
-
-## collision of g
-#force_model = forcemodels.Guo(force_g(0, MRT_collision_op))
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
-
-hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                       density=rho,
-                                                       velocity_input=u,
-                                                       force = force_g,
-                                                       optimization={"symbolic_field": g,
-                                                                     "symbolic_temporary_field": g_tmp},
-                                                       kernel_type='collide_only')
-
-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()
-
-pos_liquid_front = np.zeros((2, tf))
-pos_bubble_front = np.zeros((2, tf))
-for i in range(0, tf):
-    # write out the phasefield
-    if gpu:
-        dh.to_cpu("C")
-
-    pos_liquid_front[0, i] = i / reference_time
-    pos_liquid_front[1, i] = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
-
-    pos_bubble_front[0, i] = i / reference_time
-    pos_bubble_front[1, i] = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
-
-    # do one timestep
-    Improved_PhaseField_h_g()
-```
-
-%% Output
-
-    /home/markus/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: Numpy has detected that you (may be) writing to an array with
-    overlapping memory from np.broadcast_arrays. If this is intentional
-    set the WRITEABLE flag True or make a copy immediately before writing.
-    
-
-%% Cell type:code id: tags:
-
-``` python
-assert(pos_liquid_front[1, -1] == -0.1953125)
-assert(pos_bubble_front[1, -1] == 0.1875)
-```
-%% 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
-Nx = L0
-Ny = 4 * L0
-X0 = (Nx/2) + 1
-Y0 = (Ny/2) + 1
-
-# time step
-tf = 10001
-reference_time = 16000
-
-# force acting on the bubble
-body_force = [0, 0, 0]
-rho_H = 1.0
-
-
-parameters = calculate_parameters_rti(reference_length=L0,
-                                      reference_time=reference_time,
-                                      density_heavy=rho_H,
-                                      capillary_number=0.26,
-                                      reynolds_number=3000,
-                                      atwood_number=0.5,
-                                      peclet_number=1000,
-                                      density_ratio=3,
-                                      viscosity_ratio=1)
-
-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")
-
-
-# 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))
-
-
-# fields
-domain_size = (Nx, Ny)
-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)
-
-force = dh.add_array("force",values_per_cell=dh.dim)
-dh.fill("force", 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[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)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-if len(stencil_hydro) == 9:
-    # for D2Q9 a self defined method is used
-    x, y, z = MOMENT_SYMBOLS
-    moment_table = [sp.Integer(1),
-                    -4 + 3*(x**2 + y**2),
-                    4 - sp.Rational(21, 2)*(x**2 + y**2) + sp.Rational(9, 2)*(x**2 + y**2)**2,
-                    x,
-                    (-5 + 3*(x**2 + y**2))*x,
-                    y,
-                    (-5 + 3*(x**2 + y**2))*y ,
-                    x**2 - y**2,
-                    x*y]
-
-    relax_table = [1, 1, 1, 1, 1, 1, 1, s8, s8]
-    rr_dict = OrderedDict(zip(moment_table, relax_table))
-elif len(stencil_hydro) == 19:
-    mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, 1, 1, s8, 1, 1])
-    rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
-else:
-    mrt = create_lb_method(method="mrt", stencil=stencil_hydro, relaxation_rates=[1, 1, s8, 1, 1, 1, 1])
-    rr_dict = OrderedDict(zip(mrt.moments, mrt.relaxation_rates))
-
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, compressible=False)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# initialize the domain
-def Initialize_distributions():
-
-    for block in dh.iterate(ghost_layers=True, inner_ghost_layers=False):
-        x, y = block.midpoint_arrays
-        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)
-force_model_g = MultiphaseForceModel(force=force_g, rho=rho)
-```
-
-%% Cell type:code id: tags:
-
-``` python
-# create the allen cahl lattice boltzmann step for the h field as well as the hydrodynamic
-# lattice boltzmann step for the g field
-h_tmp_symbol_list = [h_tmp.center(i) for i, _ in enumerate(stencil_phase)]
-sum_h = np.sum(h_tmp_symbol_list[:])
-
-method_phase = create_lb_method(stencil=stencil_phase,
-                                method='srt',
-                                relaxation_rate=w_c,
-                                compressible = True,
-                                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}})
-
-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()
-#------------------------------------------------------------------------------------------------------------------------
-
-
-## collision of g
-#force_model = forcemodels.Guo(force_g(0, MRT_collision_op))
-method_hydro = create_with_discrete_maxwellian_eq_moments(stencil_hydro, rr_dict, force_model=force_model_g)
-
-hydro_lb_update_rule = get_collision_assignments_hydro(lb_method=method_hydro,
-                                                       density=rho,
-                                                       velocity_input=u,
-                                                       force = force_g,
-                                                       optimization={"symbolic_field": g,
-                                                                     "symbolic_temporary_field": g_tmp},
-                                                       kernel_type='collide_only')
-
-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()
-
-pos_liquid_front = np.zeros((2, tf))
-pos_bubble_front = np.zeros((2, tf))
-for i in range(0, tf):
-    # write out the phasefield
-    if gpu:
-        dh.to_cpu("C")
-
-    pos_liquid_front[0, i] = i / reference_time
-    pos_liquid_front[1, i] = (np.argmax(dh.cpu_arrays["C"][L0//2, :] > 0.5) - Ny//2)/L0
-
-    pos_bubble_front[0, i] = i / reference_time
-    pos_bubble_front[1, i] = (np.argmax(dh.cpu_arrays["C"][0, :] > 0.5) - Ny//2)/L0
-
-    # do one timestep
-    Improved_PhaseField_h_g()
-```
-
-%% Output
-
-    /home/markus/miniconda3/lib/python3.7/site-packages/ipykernel_launcher.py:6: DeprecationWarning: Numpy has detected that you (may be) writing to an array with
-    overlapping memory from np.broadcast_arrays. If this is intentional
-    set the WRITEABLE flag True or make a copy immediately before writing.
-    
-
-%% Cell type:code id: tags:
-
-``` python
-assert(pos_liquid_front[1, -1] == -0.1953125)
-assert(pos_bubble_front[1, -1] == 0.1875)
-```
--- a/lbmpy_tests/test_boundary_handling.py
+++ b/lbmpy_tests/test_boundary_handling.py
-import numpy as np
-
-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
-
-
-def test_simple():
-    dh = create_data_handling((10, 5), parallel=False)
-    dh.add_array('pdfs', values_per_cell=9, cpu=True, gpu=True)
-    lb_func = create_lb_function(stencil='D2Q9', compressible=False, relaxation_rate=1.8)
-
-    bh = LatticeBoltzmannBoundaryHandling(lb_func.method, dh, 'pdfs')
-
-    wall = NoSlip()
-    moving_wall = UBB((0.001, 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()
-
-
-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_boundary_indexlist_creation.py
+++ b/lbmpy_tests/test_boundary_indexlist_creation.py
-import numpy as np
-
-import pystencils.boundaries.createindexlist as cil
-from lbmpy.stencils import get_stencil
-
-
-def test_equivalence_cython_python_version():
-    if not cil.cython_funcs_available:
-        return
-
-    stencil_2d = get_stencil("D2Q9")
-    stencil_3d = get_stencil("D3Q19")
-
-    for dtype in [int, np.int16, np.uint32]:
-        fluid_mask = dtype(1)
-        mask = dtype(2)
-        flag_field_2d = np.ones([15, 16], dtype=dtype) * fluid_mask
-        flag_field_3d = np.ones([15, 16, 17], dtype=dtype) * fluid_mask
-
-        flag_field_2d[0, :] = mask
-        flag_field_2d[-1, :] = mask
-        flag_field_2d[7, 7] = mask
-
-        flag_field_3d[0, :, :] = mask
-        flag_field_3d[-1, :, :] = mask
-        flag_field_3d[7, 7, 7] = mask
-
-        result_python_2d = cil._create_boundary_neighbor_index_list_python(flag_field_2d, 1, mask, fluid_mask,
-                                                                           stencil_2d, False)
-
-        result_python_3d = cil._create_boundary_neighbor_index_list_python(flag_field_3d, 1, mask, fluid_mask,
-                                                                           stencil_3d, False)
-
-        result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask,
-                                                          fluid_mask, 1, True, False)
-        result_cython_3d = cil.create_boundary_index_list(flag_field_3d, stencil_3d, mask,
-                                                          fluid_mask, 1, True, False)
-
-        np.testing.assert_equal(result_python_2d, result_cython_2d)
-        np.testing.assert_equal(result_python_3d, result_cython_3d)
-
-
-def test_equivalence_cell_idx_list_cython_python_version():
-    if not cil.cython_funcs_available:
-        return
-
-    stencil_2d = get_stencil("D2Q9")
-    stencil_3d = get_stencil("D3Q19")
-
-    for dtype in [int, np.int16, np.uint32]:
-        fluid_mask = dtype(1)
-        mask = dtype(2)
-        flag_field_2d = np.ones([15, 16], dtype=dtype) * fluid_mask
-        flag_field_3d = np.ones([15, 16, 17], dtype=dtype) * fluid_mask
-
-        flag_field_2d[0, :] = mask
-        flag_field_2d[-1, :] = mask
-        flag_field_2d[7, 7] = mask
-
-        flag_field_3d[0, :, :] = mask
-        flag_field_3d[-1, :, :] = mask
-        flag_field_3d[7, 7, 7] = mask
-
-        result_python_2d = cil._create_boundary_cell_index_list_python(flag_field_2d, mask, fluid_mask,
-                                                                       stencil_2d, False)
-
-        result_python_3d = cil._create_boundary_cell_index_list_python(flag_field_3d, mask, fluid_mask,
-                                                                       stencil_3d, False)
-
-        result_cython_2d = cil.create_boundary_index_list(flag_field_2d, stencil_2d, mask, fluid_mask, None,
-                                                          False, False)
-        result_cython_3d = cil.create_boundary_index_list(flag_field_3d, stencil_3d, mask, fluid_mask, None,
-                                                          False, False)
-
-        np.testing.assert_equal(result_python_2d, result_cython_2d)
-        np.testing.assert_equal(result_python_3d, result_cython_3d)
--- 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
-from pystencils.llvm.llvmjit import generate_llvm
-
-
-def test_hash_equivalence_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_conserved_quantity_relaxation_invariance.py
+++ b/lbmpy_tests/test_conserved_quantity_relaxation_invariance.py
-"""
-The update equations should not change if a relaxation rate of a conserved quantity (density/velocity)
-changes. This test checks that for moment-based methods
-"""
-from copy import copy
-
-import pytest
-import sympy as sp
-
-from lbmpy.methods.creationfunctions import RelaxationInfo, create_srt, create_trt, create_trt_kbc
-from lbmpy.methods.cumulantbased import CumulantBasedLbMethod
-from lbmpy.methods.momentbased import MomentBasedLbMethod
-from lbmpy.moments import MOMENT_SYMBOLS
-from lbmpy.simplificationfactory import create_simplification_strategy
-from lbmpy.stencils import get_stencil
-
-
-def __change_relaxation_rate_of_conserved_moments(method, new_relaxation_rate=sp.Symbol("test_omega")):
-    conserved_moments = (sp.Rational(1, 1),) + MOMENT_SYMBOLS[:method.dim]
-
-    rr_dict = copy(method.relaxation_info_dict)
-    for conserved_moment in conserved_moments:
-        prev = rr_dict[conserved_moment]
-        rr_dict[conserved_moment] = RelaxationInfo(prev.equilibrium_value, new_relaxation_rate)
-
-    if isinstance(method, MomentBasedLbMethod):
-        changed_method = MomentBasedLbMethod(method.stencil, rr_dict, method.conserved_quantity_computation,
-                                             force_model=method.force_model)
-    elif isinstance(method, CumulantBasedLbMethod):
-        changed_method = CumulantBasedLbMethod(method.stencil, rr_dict, method.conserved_quantity_computation,
-                                               force_model=method.force_model)
-    else:
-        raise ValueError("Not a moment or cumulant-based method")
-
-    return changed_method
-
-
-def check_for_collision_rule_equivalence(collision_rule1, collision_rule2):
-    collision_rule1 = collision_rule1.new_without_subexpressions()
-    collision_rule2 = collision_rule2.new_without_subexpressions()
-    for eq1, eq2 in zip(collision_rule1.main_assignments, collision_rule2.main_assignments):
-        diff = sp.cancel(sp.expand(eq1.rhs - eq2.rhs))
-        assert diff == 0
-
-
-def check_method_equivalence(m1, m2, do_simplifications):
-    cr1 = m1.get_collision_rule()
-    cr2 = m2.get_collision_rule()
-    if do_simplifications:
-        cr1 = create_simplification_strategy(m1)(cr1)
-        cr2 = create_simplification_strategy(m2)(cr2)
-    check_for_collision_rule_equivalence(cr1, cr2)
-
-
-@pytest.mark.longrun
-def test_srt():
-    for stencil_name in ('D2Q9',):
-        for cumulant in (False, True):
-            stencil = get_stencil(stencil_name)
-            original_method = create_srt(stencil, sp.Symbol("omega"), cumulant=cumulant, compressible=True,
-                                         maxwellian_moments=True)
-            changed_method = __change_relaxation_rate_of_conserved_moments(original_method)
-
-            check_method_equivalence(original_method, changed_method, True)
-            check_method_equivalence(original_method, changed_method, False)
-
-
-def test_srt_short():
-    stencil = get_stencil("D2Q9")
-    original_method = create_srt(stencil, sp.Symbol("omega"), compressible=True,
-                                 maxwellian_moments=True)
-    changed_method = __change_relaxation_rate_of_conserved_moments(original_method)
-
-    check_method_equivalence(original_method, changed_method, True)
-    check_method_equivalence(original_method, changed_method, False)
-
-
-@pytest.mark.longrun
-def test_trt():
-    for stencil_name in ("D2Q9", "D3Q19", "D3Q27"):
-        for continuous_moments in (False, True):
-            stencil = get_stencil(stencil_name)
-            original_method = create_trt(stencil, sp.Symbol("omega1"), sp.Symbol("omega2"),
-                                         maxwellian_moments=continuous_moments)
-            changed_method = __change_relaxation_rate_of_conserved_moments(original_method)
-
-            check_method_equivalence(original_method, changed_method, True)
-            check_method_equivalence(original_method, changed_method, False)
-
-
-@pytest.mark.longrun
-def test_trt_kbc_long():
-    for dim in (2, 3):
-        for method_name in ("KBC-N1", "KBC-N2", "KBC-N3", "KBC-N4"):
-            original_method = create_trt_kbc(dim, sp.Symbol("omega1"), sp.Symbol("omega2"), method_name=method_name,
-                                             maxwellian_moments=False)
-            changed_method = __change_relaxation_rate_of_conserved_moments(original_method)
-            check_method_equivalence(original_method, changed_method, True)
-            check_method_equivalence(original_method, changed_method, False)
-
-
-def test_trt_kbc_short():
-    for dim, method_name in [(2, "KBC-N2")]:
-        original_method = create_trt_kbc(dim, sp.Symbol("omega1"), sp.Symbol("omega2"), method_name=method_name,
-                                         maxwellian_moments=False)
-        changed_method = __change_relaxation_rate_of_conserved_moments(original_method)
-        check_method_equivalence(original_method, changed_method, True)
-        check_method_equivalence(original_method, changed_method, False)
--- 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],
-        '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), 'mrt3', 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],
-        'method': 'mrt3',
-        'pressure_difference': 0.001,
-        'optimization': {},
-    }
-
-    scenarios = []
-
-    # Different methods
-    for method in ('srt', 'mrt3'):
-        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'] = 'mrt3'
-        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_cumulant_methods.py
+++ b/lbmpy_tests/test_cumulant_methods.py
-from lbmpy.methods import create_srt
-from lbmpy.stencils import get_stencil
-from pystencils.sympyextensions import remove_higher_order_terms
-
-
-def test_weights():
-    for stencil_name in ('D2Q9', 'D3Q19', 'D3Q27'):
-        stencil = get_stencil(stencil_name)
-        cumulant_method = create_srt(stencil, 1, cumulant=True, compressible=True,
-                                     maxwellian_moments=True)
-        moment_method = create_srt(stencil, 1, cumulant=False, compressible=True,
-                                   maxwellian_moments=True)
-        assert cumulant_method.weights == moment_method.weights
-
-
-def test_equilibrium_equivalence():
-    for stencil_name in ('D2Q9', 'D3Q19', 'D3Q27'):
-        stencil = get_stencil(stencil_name)
-        cumulant_method = create_srt(stencil, 1, cumulant=True, compressible=True,
-                                     maxwellian_moments=True)
-        moment_method = create_srt(stencil, 1, cumulant=False, compressible=True,
-                                   maxwellian_moments=True)
-        moment_eq = moment_method.get_equilibrium()
-        cumulant_eq = cumulant_method.get_equilibrium()
-        u = moment_method.first_order_equilibrium_moment_symbols
-        for mom_eq, cum_eq in zip(moment_eq.main_assignments, cumulant_eq.main_assignments):
-            diff = cum_eq.rhs - mom_eq.rhs
-            assert remove_higher_order_terms(diff.expand(), order=2, symbols=u) == 0
--- 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
No results found