.gitignore

@@ -16,4 +16,10 @@ _local_tmp
 doc/bibtex.json
 /db
 /lbmpy/phasefield/simplex_projection.*.so
-/lbmpy/phasefield/simplex_projection.c
\ No newline at end of file
+/lbmpy/phasefield/simplex_projection.c
+
+# macOS
+**/.DS_Store
+
+# benchmark database
+/lbmpy_tests/benchmark/db
\ No newline at end of file

.gitlab-ci.yml

 stages:
+  - pretest
   - test
   - deploy
 
 
-# --------------------------  Tests ------------------------------------------------------------------------------------
+# --------------------------  Pre Tests --------------------------------------------------------------------------------
 
 # Normal test - runs on every commit all but "long run" tests
 tests-and-coverage:
-  stage: test
+  stage: pretest
   except:
     variables:
       - $ENABLE_NIGHTLY_BUILDS
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   script:
-    - env
-    - pip list
+    # - pip install sympy --upgrade
     - export NUM_CORES=$(nproc --all)
     - mkdir -p ~/.config/matplotlib
     - echo "backend:template" > ~/.config/matplotlib/matplotlibrc
@@ -36,6 +36,33 @@ tests-and-coverage:
       cobertura: coverage.xml
       junit: report.xml
 
+minimal-conda:
+  stage: pretest
+  except:
+    variables:
+      - $ENABLE_NIGHTLY_BUILDS
+  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
+  script:
+    - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
+    - python setup.py quicktest
+  tags:
+    - docker
+
+# Linter for code formatting
+flake8-lint:
+  stage: pretest
+  except:
+    variables:
+      - $ENABLE_NIGHTLY_BUILDS
+  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
+  script:
+    - flake8 lbmpy
+  tags:
+    - docker
+    - cuda11
+
+# -------------------------- Tests -------------------------------------------------------------------------------------
+
 # pipeline with latest python version
 latest-python:
   stage: test
@@ -106,6 +133,22 @@ minimal-windows:
     - python -c "import numpy"
     - py.test -v -n $NUM_CORES -m "not (notebook or longrun)"
 
+minimal-sympy-master:
+  stage: test
+  except:
+    variables:
+      - $ENABLE_NIGHTLY_BUILDS
+  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
+  script:
+    - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
+    - python -m pip install --upgrade git+https://github.com/sympy/sympy.git
+    - pip list
+    - python setup.py quicktest
+  allow_failure: true
+  tags:
+    - docker
+    - cuda
+
 ubuntu:
   stage: test
   except:
@@ -113,9 +156,9 @@ ubuntu:
       - $ENABLE_NIGHTLY_BUILDS
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/ubuntu
   before_script:
-    - apt-get -y remove python3-sympy
+    # - apt-get -y remove python3-sympy
     - ln -s /usr/include/locale.h /usr/include/xlocale.h
-    - pip3 install `grep -Eo 'sympy[>=]+[0-9\.]+' setup.py | sed 's/>/=/g'`
+    # - pip3 install `grep -Eo 'sympy[>=]+[0-9\.]+' setup.py | sed 's/>/=/g'`
     - pip3 install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
   script:
     - export NUM_CORES=$(nproc --all)
@@ -132,18 +175,6 @@ ubuntu:
     reports:
       junit: report.xml
 
-minimal-conda:
-  stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
-  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/minimal_conda
-  script:
-    - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils
-    - python setup.py quicktest
-  tags:
-    - docker
-
 pycodegen-integration:
   image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
   stage: test
@@ -179,25 +210,11 @@ pycodegen-integration:
     - cuda11
     - AVX
 
-# -------------------- Linter & Documentation --------------------------------------------------------------------------
-
-
-flake8-lint:
-  stage: test
-  except:
-    variables:
-      - $ENABLE_NIGHTLY_BUILDS
-  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
-  script:
-    - flake8 lbmpy
-  tags:
-    - docker
-    - cuda11
-
+# -------------------- Documentation and deploy ------------------------------------------------------------------------
 
 build-documentation:
   stage: test
-  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/full
+  image: i10git.cs.fau.de:5005/pycodegen/pycodegen/documentation
   script:
     - export PYTHONPATH=`pwd`
     - pip install git+https://gitlab-ci-token:${CI_JOB_TOKEN}@i10git.cs.fau.de/pycodegen/pystencils.git@master#egg=pystencils

AUTHORS.txt

+
+Contributors:
+-------------
+
+  - Martin Bauer <martin.bauer@fau.de>
+  - Markus Holzer <markus.holzer@fau.de>
+  - Michael Kuron <mkuron@icp.uni-stuttgart.de>
+  - Stephan Seitz <stephan.seitz@fau.de>
+  - Frederik Hennig <frederik.hennig@fau.de>
+  - Helen Schottenhamml <helen.schottenhamml@fau.de>
+  - Rudolf Weeber <weeber@icp.uni-stuttgart.de>
+  - Christian Godenschwager <christian.godenschwager@fau.de>
+  - Jan Hönig <jan.hoenig@fau.de>

CHANGELOG.md

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

CONTRIBUTING.md

+# Contributing
+
+lbmpy is built on the open-source python framework [pystencils](https://pypi.org/project/pystencils/). Please consider the [contribution guideline](https://i10git.cs.fau.de/pycodegen/pystencils/-/blob/master/CONTRIBUTING.md) of pystencils for contributing to lbmpy.
\ No newline at end of file

MANIFEST.in

 include README.md
 include COPYING.txt
-include RELEASE-VERSION
+include AUTHORS.txt
+include CONTRIBUTING.md
 global-include *.pyx
 include versioneer.py
 include lbmpy/_version.py

README.md

@@ -16,11 +16,12 @@ It even comes with an integrated Chapman Enskog analysis based on sympy!
 
 Common test scenarios can be set up quickly:
 ```python
-from lbmpy.scenarios import create_channel
+from pystencils import Target
+from lbmpy.session import *
 
-ch = create_channel(domain_size=(300,100, 100), force=1e-7, method="trt",
+ch = create_channel(domain_size=(300, 100, 100), force=1e-7, method=Method.TRT,
                     equilibrium_order=2, compressible=True,
-                    relaxation_rates=[1.97, 1.6], optimization={'target': 'gpu'})
+                    relaxation_rates=[1.97, 1.6], optimization={'target': Target.GPU})
 ```
 
 To find out more, check out the interactive [tutorial notebooks online with binder](https://mybinder.org/v2/gh/mabau/lbmpy/master?filepath=doc%2Fnotebooks).
@@ -55,3 +56,27 @@ Documentation
 
 Read the docs [here](http://pycodegen.pages.i10git.cs.fau.de/lbmpy) and
 check out the Jupyter notebooks in `doc/notebooks`. 
+
+Contributing
+-------
+To see how to open issues, submit bug reports, create feature requests or submit your additions to lbmpy please refer to
+[contribution documentation](https://i10git.cs.fau.de/pycodegen/pystencils/-/blob/master/CONTRIBUTING.md) of pystencils since lbmpy is heavily build on pystencils.
+
+
+Many thanks go to the [contributors](https://i10git.cs.fau.de/pycodegen/lbmpy/-/blob/master/AUTHORS.txt) of lbmpy.
+
+
+### Please cite us
+
+If you use lbmpy in a publication, please cite the following articles:
+
+Overview:
+  - M. Bauer et al, lbmpy: Automatic code generation for efficient parallel lattice Boltzmann methods. Journal of Computational Science, 2021. https://doi.org/10.1016/j.jocs.2020.101269 ([Preprint](https://arxiv.org/abs/2001.11806))
+
+Multiphase:
+   - M. Holzer et al, Highly efficient lattice Boltzmann multiphase simulations of immiscible fluids at high-density ratios on CPUs and GPUs through code generation. The International Journal of High Performance Computing Applications, 2021. https://doi.org/10.1177/10943420211016525
+
+
+### Further Reading
+
+- F. Hennig et al, Automatic Code Generation for the Cumulant Lattice Boltzmann Method. ICMMES, 2021. [Poster Link](https://www.researchgate.net/publication/353224406_Automatic_Code_Generation_for_the_Cumulant_Lattice_Boltzmann_Method)

conftest.py

@@ -4,6 +4,9 @@ import tempfile
 import runpy
 import sys
 import warnings
+import nbformat
+from nbconvert import PythonExporter
+import sympy
 # Trigger config file reading / creation once - to avoid race conditions when multiple instances are creating it
 # at the same time
 from pystencils.cpu import cpujit
@@ -12,12 +15,12 @@ from pystencils.cpu import cpujit
 # at the same time
 try:
     import pyximport
+
     pyximport.install(language_level=3)
 except ImportError:
     pass
 from lbmpy.phasefield.simplex_projection import simplex_projection_2d  # NOQA
 
-
 SCRIPT_FOLDER = os.path.dirname(os.path.realpath(__file__))
 sys.path.insert(0, os.path.abspath('lbmpy'))
 
@@ -40,7 +43,6 @@ collect_ignore = [os.path.join(SCRIPT_FOLDER, "doc", "conf.py"),
 add_path_to_ignore('pystencils_tests/benchmark')
 add_path_to_ignore('_local_tmp')
 
-
 try:
     import pycuda
 except ImportError:
@@ -49,18 +51,17 @@ except ImportError:
 try:
     import waLBerla
 except ImportError:
-    collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_serial_scenarios.py")]
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_datahandling_parallel.py")]
 
 try:
     import blitzdb
 except ImportError:
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/benchmark"),
-                       os.path.join(SCRIPT_FOLDER, "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py")]
+                       os.path.join(SCRIPT_FOLDER,
+                                    "lbmpy_tests/full_scenarios/kida_vortex_flow/scenario_kida_vortex_flow.py")]
 
-import sympy
 sver = sympy.__version__.split(".")
-if (int(sver[0]) == 1 and int(sver[1]) < 2):
+if int(sver[0]) == 1 and int(sver[1]) < 2:
     add_path_to_ignore('lbmpy_tests/phasefield')
     collect_ignore += [os.path.join(SCRIPT_FOLDER, "lbmpy_tests/test_n_phase_boyer_noncoupled.ipynb")]
 
@@ -72,10 +73,6 @@ for root, sub_dirs, files in os.walk('.'):
             collect_ignore.append(f)
 
 
-import nbformat
-from nbconvert import PythonExporter
-
-
 class IPythonMockup:
     def run_line_magic(self, *args, **kwargs):
         pass

doc/conf.py

@@ -26,8 +26,8 @@ templates_path = ['_templates']
 source_suffix = '.rst'
 master_doc = 'index'
 
-copyright = f'{datetime.datetime.now().year}, Martin Bauer'
-author = 'Martin Bauer'
+copyright = f'{datetime.datetime.now().year}, Martin Bauer, Markus Holzer'
+author = 'Martin Bauer, Markus Holzer'
 # The short X.Y version (including .devXXXX, rcX, b1 suffixes if present)
 version = re.sub(r'(\d+\.\d+)\.\d+(.*)', r'\1\2', lbmpy.__version__)
 version = re.sub(r'(\.dev\d+).*?$', r'\1', version)

doc/notebooks/08_tutorial_shanchen_twophase.ipynb

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

doc/notebooks/10_tutorial_conservative_allen_cahn_two_phase.ipynb

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

doc/notebooks/demo_automatic_chapman_enskog_analysis.ipynb

 %% Cell type:code id: tags:
 
 ``` python
-import sympy as sp
+from lbmpy.session import *
 from lbmpy.chapman_enskog import ChapmanEnskogAnalysis
-from lbmpy.creationfunctions import create_lb_method
 sp.init_printing()
 ```
 
 %% Cell type:markdown id: tags:
 
 # Demo: Automatic Chapman Enskog Analysis
 
 
 First, we create a SRT lattice Boltzmann method. It is defined as the set of moments, together with one relaxation rate per moment.
 
 %% Cell type:code id: tags:
 
 ``` python
-method = create_lb_method(method='trt', stencil='D3Q19', compressible=False)
+method = create_lb_method(LBMConfig(method=Method.TRT, stencil=Stencil.D3Q19, compressible=False))
 method
 ```
 
 %% Output
 
-    <lbmpy.methods.momentbased.MomentBasedLbMethod at 0x7f9bbc455a90>
+    <lbmpy.methods.momentbased.momentbasedmethod.MomentBasedLbMethod at 0x7f5ea617b190>
 
 %% Cell type:markdown id: tags:
 
 Next, the Chapman Enskog analysis object is created. This may take a while...
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis = ChapmanEnskogAnalysis(method)
 ```
 
 %% Cell type:markdown id: tags:
 
 This object now information about the method, e.g. the relation of relaxation rate to viscosities, if the method approximates the compressible or incompressible continuity equation ...
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis.compressible
 ```
 
 %% Output
 
     False
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis.pressure_equation
 ```
 
 %% Output
 
-
-    $$\left [ \frac{\rho}{3}\right ]$$
+
+    $\displaystyle \left[ \frac{\rho}{3}\right]$
     ⎡ρ⎤
     ⎢─⎥
     ⎣3⎦
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis.get_kinematic_viscosity()
 ```
 
 %% Output
 
-
-    $$- \frac{\omega_{0} - 2}{6 \omega_{0}}$$
-    -(ω₀ - 2)
-    ──────────
-       6⋅ω₀
+
+    $\displaystyle - \frac{\omega - 2}{6 \omega}$
+    -(ω - 2)
+    ─────────
+       6⋅ω
 
 %% Cell type:code id: tags:
 
 ``` python
 analysis.get_bulk_viscosity()
 ```
 
 %% Output
 
-
-    $$- \frac{1}{9} + \frac{2}{9 \omega_{0}}$$
+
+    $\displaystyle - \frac{1}{9} + \frac{2}{9 \omega}$
       1    2
-    - ─ + ────
-      9   9⋅ω₀
+    - ─ + ───
+      9   9⋅ω
 
 %% Cell type:markdown id: tags:
 
 But also details of the analysis are available:
 
 %% Cell type:code id: tags:
 
 ``` python
 sp.Matrix(analysis.get_macroscopic_equations())
 ```
 
 %% Output
 
-    $$\left[\begin{matrix}{\partial_{t} \rho} + {\partial_{0} u_{0}} + {\partial_{1} u_{1}} + {\partial_{2} u_{2}}\\- \frac{\epsilon \omega_{0}}{2} {\partial_{2} {\Pi_{02}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{1} {\Pi_{01}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{0} {\Pi_{00}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{02}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{01}^{(1)}}} + \epsilon {\partial_{0} {\Pi_{00}^{(1)}}} + \frac{{\partial_{0} \rho}}{3} + {\partial_{t} u_{0}} + {\partial_{0} (u_{0}^{2}) } + {\partial_{1} (u_{0} u_{1}) } + {\partial_{2} (u_{0} u_{2}) }\\- \frac{\epsilon \omega_{0}}{2} {\partial_{2} {\Pi_{12}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{1} {\Pi_{11}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{0} {\Pi_{01}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{12}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{11}^{(1)}}} + \epsilon {\partial_{0} {\Pi_{01}^{(1)}}} + \frac{{\partial_{1} \rho}}{3} + {\partial_{t} u_{1}} + {\partial_{1} (u_{1}^{2}) } + {\partial_{0} (u_{0} u_{1}) } + {\partial_{2} (u_{1} u_{2}) }\\- \frac{\epsilon \omega_{0}}{2} {\partial_{2} {\Pi_{22}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{1} {\Pi_{12}^{(1)}}} - \frac{\epsilon \omega_{0}}{2} {\partial_{0} {\Pi_{02}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{22}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{12}^{(1)}}} + \epsilon {\partial_{0} {\Pi_{02}^{(1)}}} + \frac{{\partial_{2} \rho}}{3} + {\partial_{t} u_{2}} + {\partial_{2} (u_{2}^{2}) } + {\partial_{0} (u_{0} u_{2}) } + {\partial_{1} (u_{1} u_{2}) }\end{matrix}\right]$$
+
+    $\displaystyle \left[\begin{matrix}{\partial_{t} \rho} + {\partial_{0} u_{0}} + {\partial_{1} u_{1}} + {\partial_{2} u_{2}}\\- \frac{\epsilon \omega {\partial_{0} {\Pi_{00}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{1} {\Pi_{01}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{2} {\Pi_{02}^{(1)}}}}{2} + \epsilon {\partial_{0} {\Pi_{00}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{01}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{02}^{(1)}}} + \frac{{\partial_{0} \rho}}{3} + {\partial_{t} u_{0}} + {\partial_{0} (u_{0}^{2}) } + {\partial_{1} (u_{0} u_{1}) } + {\partial_{2} (u_{0} u_{2}) }\\- \frac{\epsilon \omega {\partial_{0} {\Pi_{01}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{1} {\Pi_{11}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{2} {\Pi_{12}^{(1)}}}}{2} + \epsilon {\partial_{0} {\Pi_{01}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{11}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{12}^{(1)}}} + \frac{{\partial_{1} \rho}}{3} + {\partial_{t} u_{1}} + {\partial_{1} (u_{1}^{2}) } + {\partial_{0} (u_{0} u_{1}) } + {\partial_{2} (u_{1} u_{2}) }\\- \frac{\epsilon \omega {\partial_{0} {\Pi_{02}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{1} {\Pi_{12}^{(1)}}}}{2} - \frac{\epsilon \omega {\partial_{2} {\Pi_{22}^{(1)}}}}{2} + \epsilon {\partial_{0} {\Pi_{02}^{(1)}}} + \epsilon {\partial_{1} {\Pi_{12}^{(1)}}} + \epsilon {\partial_{2} {\Pi_{22}^{(1)}}} + \frac{{\partial_{2} \rho}}{3} + {\partial_{t} u_{2}} + {\partial_{2} (u_{2}^{2}) } + {\partial_{0} (u_{0} u_{2}) } + {\partial_{1} (u_{1} u_{2}) }\end{matrix}\right]$
     ⎡
     ⎢
-    ⎢  ε⋅ω₀⋅D(\Pi_(1)_(1, 0, 1))   ε⋅ω₀⋅D(\Pi_(1)_(1, 1, 0))   ε⋅ω₀⋅D(\Pi_(1)_(2,
-    ⎢- ───────────────────────── - ───────────────────────── - ───────────────────
-    ⎢              2                           2                           2
+    ⎢  ε⋅ω⋅D(\Pi_(1)_(2, 0, 0))   ε⋅ω⋅D(\Pi_(1)_(1, 1, 0))   ε⋅ω⋅D(\Pi_(1)_(1, 0,
+    ⎢- ──────────────────────── - ──────────────────────── - ─────────────────────
+    ⎢             2                          2                          2
     ⎢
-    ⎢  ε⋅ω₀⋅D(\Pi_(1)_(0, 1, 1))   ε⋅ω₀⋅D(\Pi_(1)_(0, 2, 0))   ε⋅ω₀⋅D(\Pi_(1)_(1,
-    ⎢- ───────────────────────── - ───────────────────────── - ───────────────────
-    ⎢              2                           2                           2
+    ⎢  ε⋅ω⋅D(\Pi_(1)_(1, 1, 0))   ε⋅ω⋅D(\Pi_(1)_(0, 2, 0))   ε⋅ω⋅D(\Pi_(1)_(0, 1,
+    ⎢- ──────────────────────── - ──────────────────────── - ─────────────────────
+    ⎢             2                          2                          2
     ⎢
-    ⎢  ε⋅ω₀⋅D(\Pi_(1)_(0, 0, 2))   ε⋅ω₀⋅D(\Pi_(1)_(0, 1, 1))   ε⋅ω₀⋅D(\Pi_(1)_(1,
-    ⎢- ───────────────────────── - ───────────────────────── - ───────────────────
-    ⎣              2                           2                           2
+    ⎢  ε⋅ω⋅D(\Pi_(1)_(1, 0, 1))   ε⋅ω⋅D(\Pi_(1)_(0, 1, 1))   ε⋅ω⋅D(\Pi_(1)_(0, 0,
+    ⎢- ──────────────────────── - ──────────────────────── - ─────────────────────
+    ⎣             2                          2                          2
     
-                 D(rho) + D(u_0) + D(u_1) + D(u_2)
+                D(rho) + D(u_0) + D(u_1) + D(u_2)
     
-    0, 0))
-    ────── + ε⋅D(\Pi_(1)_(1, 0, 1)) + ε⋅D(\Pi_(1)_(1, 1, 0)) + ε⋅D(\Pi_(1)_(2, 0,
+    1))
+    ─── + ε⋅D(\Pi_(1)_(2, 0, 0)) + ε⋅D(\Pi_(1)_(1, 1, 0)) + ε⋅D(\Pi_(1)_(1, 0, 1))
     
     
-    1, 0))
-    ────── + ε⋅D(\Pi_(1)_(0, 1, 1)) + ε⋅D(\Pi_(1)_(0, 2, 0)) + ε⋅D(\Pi_(1)_(1, 1,
+    1))
+    ─── + ε⋅D(\Pi_(1)_(1, 1, 0)) + ε⋅D(\Pi_(1)_(0, 2, 0)) + ε⋅D(\Pi_(1)_(0, 1, 1))
     
     
-    0, 1))
-    ────── + ε⋅D(\Pi_(1)_(0, 0, 2)) + ε⋅D(\Pi_(1)_(0, 1, 1)) + ε⋅D(\Pi_(1)_(1, 0,
+    2))
+    ─── + ε⋅D(\Pi_(1)_(1, 0, 1)) + ε⋅D(\Pi_(1)_(0, 1, 1)) + ε⋅D(\Pi_(1)_(0, 0, 2))
     
     
-                                                               ⎤
-                                                               ⎥
-          D(rho)                                               ⎥
-    0)) + ────── + D(u_0) + D(u_0**2) + D(u_0*u_1) + D(u_0*u_2)⎥
-            3                                                  ⎥
-                                                               ⎥
-          D(rho)                                               ⎥
-    0)) + ────── + D(u_1) + D(u_1**2) + D(u_0*u_1) + D(u_1*u_2)⎥
-            3                                                  ⎥
-                                                               ⎥
-          D(rho)                                               ⎥
-    1)) + ────── + D(u_2) + D(u_2**2) + D(u_0*u_2) + D(u_1*u_2)⎥
-            3                                                  ⎦
+                                                            ⎤
+                                                            ⎥
+       D(rho)                                               ⎥
+     + ────── + D(u_0) + D(u_0**2) + D(u_0*u_1) + D(u_0*u_2)⎥
+         3                                                  ⎥
+                                                            ⎥
+       D(rho)                                               ⎥
+     + ────── + D(u_1) + D(u_1**2) + D(u_0*u_1) + D(u_1*u_2)⎥
+         3                                                  ⎥
+                                                            ⎥
+       D(rho)                                               ⎥
+     + ────── + D(u_2) + D(u_2**2) + D(u_0*u_2) + D(u_1*u_2)⎥
+         3                                                  ⎦