diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 7b08340a7fa7fae99f97e6f3bcc9360a588dc3db..52e3b9b9dd3b1aa2c238bf8566aa7634bbf94a4d 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -3,6 +3,7 @@ image: irmb/virtualfluids-python-deps:latest
 stages:
   - build
   - build_python
+  - container_upload
   - test
   - benchmark
   - analyze
@@ -141,26 +142,6 @@ msvc_16:
     paths:
       - $CI_PROJECT_DIR/$env:BUILD_FOLDER/
 
-###############################################################################
-build_singularity_image:
-  stage: build
-
-  tags:
-    - priviliged
-    - linux
-
-  rules:
-    - if: $CI_COMMIT_TAG
-
-  artifacts:
-    expire_in: 1 hrs
-    paths:
-      - Containers/VirtualFluidsOpenMPI.sif
-
-  script:
-    - singularity build --fakeroot Containers/VirtualFluidsOpenMPI.sif Containers/VirtualFluidsOpenMPI.def
-    - ls -sh Containers/VirtualFluidsOpenMPI.sif
-
 
 ###############################################################################
 ##                             Build Python                                  ##
@@ -189,6 +170,26 @@ gcc_9_python:
   script:
     - python3 setup.py bdist_wheel
 
+###############################################################################
+##                            Container Upload                               ##
+###############################################################################
+build_singularity_image:
+  stage: container_upload
+
+  needs:
+    - gcc_9_python
+
+  tags:
+    - linux
+    - privileged
+
+  rules:
+    - if: $CI_COMMIT_TAG
+
+  script:
+    - singularity build Containers/VirtualFluidsPython.sif Containers/VirtualFluidsPython.def
+    - singularity push --docker-username "${CI_REGISTRY_USER}" --docker-password "${CI_REGISTRY_PASSWORD}" Containers/VirtualFluidsPython.sif oras://"$CI_REGISTRY_IMAGE"/"$CI_PROJECT_NAME":"$CI_COMMIT_TAG"
+
 ###############################################################################
 ##                                Tests                                      ##
 ###############################################################################
@@ -236,6 +237,47 @@ gcc_9_python_bindings_test:
   script:
     - python3 -m unittest discover -s Python -v
 
+###############################################################################
+gcc_9_python_slurm_test:
+  stage: test
+
+  needs: ["gcc_9_python"]
+
+  rules:
+    - if: $PHOENIX_PRIVATE_KEY
+
+  tags:
+    - linux
+    - privileged
+
+  variables:
+    SSH_KEY: "$PHOENIX_PRIVATE_KEY"
+    HOST: "$PHOENIX_HOSTNAME"
+    USER: "$PHOENIX_USER"
+
+  before_script:
+    - 'command -v ssh-agent >/dev/null || ( apt-get update -y && apt-get install openssh-client -y )'
+    - apt-get install -y rsync
+    - mkdir -p ~/.ssh
+    - chmod 700 ~/.ssh
+    - eval $(ssh-agent -s)
+    - echo "$SSH_KEY" | tr -d '\r' | ssh-add -
+    - echo $SSH_KEY >> ansible/private_key
+    - ssh-keyscan -t rsa $HOST >> ~/.ssh/known_hosts
+    - ssh $USER@$HOST "rm -rf output-*"
+    - ssh $USER@$HOST "rm -f *.out"
+    - pip3 install git+git://github.com/SvenMarcus/ssh-slurm-runner
+
+  script:
+    - singularity build PoiseuilleTestContainer.sif Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def
+    - scp PoiseuilleTestContainer.sif $USER@$HOST:PoiseuilleTestContainer.sif
+    - scp Python/SlurmTests/poiseuille/slurm.job $USER@$HOST:slurm.job
+    - python3 -m ssh_slurm_runner slurm.job --host $HOST --user $USER --keyfile ansible/private_key
+    - ssh $USER@$HOST "rm -rf output-*"
+    - ssh $USER@$HOST "rm -f *.out"
+    - ssh $USER@$HOST "rm PoiseuilleTestContainer.sif"
+    - ssh $USER@$HOST "rm slurm.job"
+
 ###############################################################################
 ##                            Benchmark                                      ##
 ###############################################################################
@@ -549,10 +591,10 @@ vf_wheel_to_jupyterhub:
   needs: ["gcc_9_python", "gcc_9_unit_tests", "gcc_9_python_bindings_test"]
 
   variables:
-    HOST: "finrod.irmb.bau.tu-bs.de"
-    SSH_KEY: "$SSH_PRIVATE_KEY_JUPYTER_HOST_AT_FINROD"
-    REMOTE_USER: "jupyter_host"
-    jupyter_host: "jupyter_host"
+    HOST: "gitlab-runner01.irmb.bau.tu-bs.de"
+    SSH_KEY: "$SSH_PRIVATE_KEY"
+    REMOTE_USER: "runner"
+    jupyter_host: "runner"
 
   script:
     - ansible-playbook -i ansible/hosts.cfg -u $REMOTE_USER ansible/playbook_jupyter_update.yml
diff --git a/Containers/VirtualFluidsPython.def b/Containers/VirtualFluidsPython.def
new file mode 100644
index 0000000000000000000000000000000000000000..d54066bc634cf25f4340b1e659eae72515467fa8
--- /dev/null
+++ b/Containers/VirtualFluidsPython.def
@@ -0,0 +1,33 @@
+BootStrap: docker
+From: ubuntu:20.04
+
+%files
+    Python Python
+    dist dist
+
+
+%post
+    export DEBIAN_FRONTEND=noninteractive
+    apt-get update &&          \
+    apt-get install -y         \
+    build-essential            \
+    cmake=3.16.3-1ubuntu1      \
+    python3                    \
+    python3-dev                \
+    python3-pip                \
+    mpich                      \
+    libomp-dev
+
+    pip3 install setuptools wheel $(find dist/*.whl)
+
+%environment
+    export PYTHONPATH=/Python
+
+%runscript
+    python3 /Python/liddrivencavity/simulation.py
+
+%appenv poiseuille
+    export PYTHONPATH=Python
+
+%apprun poisueille
+    python3 /Python/poiseuille/poiseuille_hpc.py
diff --git a/Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def b/Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def
new file mode 100644
index 0000000000000000000000000000000000000000..a3836e7906b9be66ec79f68bf53ccc079db9d9ef
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/PoiseuilleTestContainer.def
@@ -0,0 +1,42 @@
+BootStrap: docker
+From: ubuntu:20.04
+
+%files
+    3rdParty 3rdParty
+    apps apps
+    CMake CMake
+    Python Python
+    src src
+    CMakeLists.txt CMakeLists.txt
+    cpu.cmake cpu.cmake
+    gpu.cmake gpu.cmake
+    setup.py setup.py
+    pyproject.toml pyproject.toml
+
+
+%post
+    export DEBIAN_FRONTEND=noninteractive
+    apt-get update &&          \
+    apt-get install -y         \
+    build-essential            \
+    cmake=3.16.3-1ubuntu1      \
+    python3                    \
+    python3-dev                \
+    python3-pip                \
+    mpich                      \
+    libomp-dev                 \
+    libgl1
+
+    pip3 install setuptools wheel numpy scipy pyvista
+
+    export PYTHONPATH=Python
+    python3 /setup.py install
+
+%environment
+    export PYTHONPATH=/Python
+
+%apprun testcase
+    python3 /Python/SlurmTests/poiseuille/simulation_runner.py
+
+%apprun evaluate
+    python3 /Python/SlurmTests/poiseuille/evaluator.py
\ No newline at end of file
diff --git a/Python/SlurmTests/poiseuille/evaluator.py b/Python/SlurmTests/poiseuille/evaluator.py
new file mode 100644
index 0000000000000000000000000000000000000000..74602846b67bb82f7e3c3d3ca3015fbe00a63041
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/evaluator.py
@@ -0,0 +1,20 @@
+import numpy as np
+import scipy.stats as stats
+import errors
+from SlurmTests.poiseuille.result_collector import collect_results
+from SlurmTests.poiseuille.settings import Scaling
+
+analytical_results, numerical_results = collect_results()
+normalized_l2_errors = [errors.normalized_l2_error(analytical, numerical)
+                        for analytical, numerical in zip(analytical_results, numerical_results)]
+
+nodes_in_x3_per_run = []
+for simulation_run in range(0, 3):
+    grid_params, _, _, _ = Scaling.configuration_for_scale_level(simulation_run)
+    nodes_in_x3_per_run.append(grid_params.number_of_nodes_per_direction[2])
+
+nodes_as_log = [np.log10(node) for node in nodes_in_x3_per_run]
+l2_norms_as_log = [np.log10(l2) for l2 in normalized_l2_errors]
+res = stats.linregress(nodes_as_log, l2_norms_as_log)
+
+assert res.slope <= -2, f"Expected slope of l2 error to be <= -2, but was {res.slope}"
diff --git a/Python/SlurmTests/poiseuille/result_collector.py b/Python/SlurmTests/poiseuille/result_collector.py
new file mode 100644
index 0000000000000000000000000000000000000000..06efa481c8c010647531426f2af2bec2c2d7eaee
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/result_collector.py
@@ -0,0 +1,73 @@
+from typing import Collection, List
+
+import pyvista as pv
+from poiseuille.analytical import PoiseuilleSettings, poiseuille_at_heights
+from vtk_utilities import vertical_column_from_mesh, get_values_from_indices
+from SlurmTests.poiseuille.settings import Scaling
+
+
+def get_output_file_name(output_folder, runtime_params):
+    timesteps = runtime_params.number_of_timesteps
+    file_name = f"{output_folder}/mq/mq{timesteps}/mq0_{timesteps}.bin.vtu"
+
+    return file_name
+
+
+def get_mesh_for_last_timestep(output_folder, runtime_params):
+    file_name_of_last_timestep = get_output_file_name(output_folder, runtime_params)
+    mesh_of_last_timestep = pv.read(file_name_of_last_timestep)
+    return mesh_of_last_timestep
+
+
+def get_heights_from_indices(mesh, indices):
+    return [mesh.points[index][2] for index in indices]
+
+
+def get_heights(output_folder, runtime_params):
+    mesh_of_last_timestep = get_mesh_for_last_timestep(output_folder, runtime_params)
+    column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
+    heights = get_heights_from_indices(mesh_of_last_timestep, column_indices)
+    return heights
+
+
+def get_numerical_results(runtime_params, output_folder):
+    mesh_of_last_timestep = get_mesh_for_last_timestep(output_folder, runtime_params)
+    velocities_in_x_direction = mesh_of_last_timestep.get_array("Vx")
+    column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
+    numerical_results = get_values_from_indices(velocities_in_x_direction, column_indices)
+
+    return numerical_results
+
+
+def get_analytical_results(grid_params, physical_params, kernel, height_values):
+    channel_height = grid_params.number_of_nodes_per_direction[2]
+    settings = get_analytical_poiseuille_settings(channel_height, physical_params, kernel)
+    max_grid_height = channel_height * grid_params.node_distance
+    adjusted_height_values = [value / max_grid_height * channel_height for value in height_values]
+    analytical_results = poiseuille_at_heights(settings, adjusted_height_values)
+    return analytical_results
+
+
+def get_analytical_poiseuille_settings(height, physical_params, kernel):
+    settings = PoiseuilleSettings()
+    settings.height = height
+    settings.viscosity = physical_params.lattice_viscosity
+    settings.density = 1
+    settings.force = kernel.forcing_in_x1
+
+    return settings
+
+
+def collect_results() -> (List[List[float]], List[List[float]]):
+    analytical_results = []
+    numerical_results = []
+
+    for simulation_run in range(0, 3):
+        output_folder = f"output-{simulation_run}"
+        grid_params, physical_params, runtime_params, kernel = Scaling.configuration_for_scale_level(simulation_run)
+        heights = get_heights(output_folder, runtime_params)
+        analytical_results.append(
+            get_analytical_results(grid_params, physical_params, kernel, heights))
+        numerical_results.append(get_numerical_results(runtime_params, output_folder))
+
+    return analytical_results, numerical_results
diff --git a/Python/SlurmTests/poiseuille/settings.py b/Python/SlurmTests/poiseuille/settings.py
new file mode 100644
index 0000000000000000000000000000000000000000..f75c2b1d7133323880dd5520de0a96cb8fa87860
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/settings.py
@@ -0,0 +1,26 @@
+import os
+from acousticscaling import OneDirectionalAcousticScaling
+from pyfluids.kernel import LBMKernel, KernelType
+from pyfluids.parameters import RuntimeParameters, GridParameters, PhysicalParameters
+
+
+grid_params = GridParameters()
+grid_params.node_distance = 1
+grid_params.number_of_nodes_per_direction = [1, 1, 16]
+grid_params.blocks_per_direction = [1, 1, 4]
+grid_params.periodic_boundary_in_x1 = True
+grid_params.periodic_boundary_in_x2 = True
+
+physical_params = PhysicalParameters()
+physical_params.lattice_viscosity = 1e-4
+
+runtime_params = RuntimeParameters()
+runtime_params.number_of_threads = int(os.environ["PYFLUIDS_NUM_THREADS"])
+runtime_params.number_of_timesteps = 4_000_000
+runtime_params.timestep_log_interval = 1_000_000
+
+kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+kernel.use_forcing = True
+kernel.forcing_in_x1 = 5e-10
+
+Scaling = OneDirectionalAcousticScaling(grid_params, physical_params, runtime_params, kernel)
diff --git a/Python/SlurmTests/poiseuille/simulation_runner.py b/Python/SlurmTests/poiseuille/simulation_runner.py
new file mode 100644
index 0000000000000000000000000000000000000000..0b75de40b6a8f11ccd76f97f2ed9d709dc5362dd
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/simulation_runner.py
@@ -0,0 +1,19 @@
+import os
+
+from SlurmTests.poiseuille.settings import Scaling
+from poiseuille.simulation import run_simulation
+from pyfluids.writer import Writer, OutputFormat
+
+
+scale_level = int(os.environ["PYFLUIDS_SCALE_LEVEL"])
+grid_params, physical_params, runtime_params, kernel = Scaling.configuration_for_scale_level(scale_level)
+
+writer = Writer()
+writer.output_format = OutputFormat.BINARY
+writer.output_path = "./output-" + str(scale_level)
+
+run_simulation(grid_params=grid_params,
+               physical_params=physical_params,
+               runtime_params=runtime_params,
+               kernel=kernel,
+               writer=writer)
diff --git a/Python/SlurmTests/poiseuille/slurm.job b/Python/SlurmTests/poiseuille/slurm.job
new file mode 100644
index 0000000000000000000000000000000000000000..488fc9a42f261d69a8212cff389721fdfb9cbf6e
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/slurm.job
@@ -0,0 +1,26 @@
+#!/bin/bash
+#SBATCH -J PyFluidsTest
+#SBATCH --nodes=1
+#SBATCH --ntasks-per-node=1
+#SBATCH --cpus-per-task=20
+
+#SBATCH --mem-per-cpu=3000
+#SBATCH --time=07:00:00
+#SBATCH --partition=standard
+
+source $HOME/.bashrc
+
+echo "PyFluids Poiseuille Test Case"
+echo "Number of tasks: ${SLURM_NTASKS}"
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=0
+export SINGULARITYENV_PYFLUIDS_NUM_THREADS=4
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=1
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+export SINGULARITYENV_PYFLUIDS_SCALE_LEVEL=2
+srun singularity run --app testcase PoiseuilleTestContainer.sif
+
+srun singularity run --app evaluate PoiseuilleTestContainer.sif
diff --git a/Python/acousticscaling.py b/Python/acousticscaling.py
new file mode 100644
index 0000000000000000000000000000000000000000..50b81db064251fa269f29bf72a561567ddedafbc
--- /dev/null
+++ b/Python/acousticscaling.py
@@ -0,0 +1,85 @@
+from pyfluids.kernel import LBMKernel
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+
+class OneDirectionalAcousticScaling:
+
+    def __init__(self, grid_parameters: GridParameters,
+                 physical_parameters: PhysicalParameters,
+                 runtime_parameters: RuntimeParameters,
+                 kernel: LBMKernel):
+        self._grid_params = grid_parameters
+        self._physical_params = physical_parameters
+        self._runtime_params = runtime_parameters
+        self._kernel = kernel
+
+    def configuration_for_scale_level(self, level: int = 1) -> (GridParameters,
+                                                                PhysicalParameters,
+                                                                RuntimeParameters,
+                                                                LBMKernel):
+        if level < 0:
+            raise ValueError("level must be >= 0")
+
+        grid_params = self.clone_grid_params_for_level(level)
+        physical_params = self.clone_physical_parameters(level)
+        runtime_params = self.clone_runtime_params_for_level(level)
+        kernel = self.clone_kernel_for_level(level)
+
+        return grid_params, physical_params, runtime_params, kernel
+
+    def clone_grid_params_for_level(self, level) -> GridParameters:
+        grid_params = GridParameters()
+        grid_params.reference_direction_index = self._grid_params.reference_direction_index
+        grid_params.periodic_boundary_in_x1 = self._grid_params.periodic_boundary_in_x1
+        grid_params.periodic_boundary_in_x2 = self._grid_params.periodic_boundary_in_x2
+        grid_params.periodic_boundary_in_x3 = self._grid_params.periodic_boundary_in_x3
+
+        grid_params.number_of_nodes_per_direction = list(self._grid_params.number_of_nodes_per_direction)
+        grid_params.blocks_per_direction = list(self._grid_params.blocks_per_direction)
+        grid_params.node_distance = self._grid_params.node_distance
+
+        if level > 0:
+            grid_params.node_distance /= (level * 2)
+            grid_params.number_of_nodes_per_direction = [grid_params.number_of_nodes_per_direction[0],
+                                                         grid_params.number_of_nodes_per_direction[1],
+                                                         grid_params.number_of_nodes_per_direction[2] * (level * 2)]
+
+            grid_params.blocks_per_direction = [grid_params.blocks_per_direction[0],
+                                                grid_params.blocks_per_direction[1],
+                                                grid_params.blocks_per_direction[2] * (level * 2)]
+
+        return grid_params
+
+    def clone_physical_parameters(self, level):
+        physical_params = PhysicalParameters()
+        physical_params.lattice_viscosity = self._physical_params.lattice_viscosity
+
+        if level > 0:
+            physical_params.lattice_viscosity *= (level * 2)
+
+        return physical_params
+
+    def clone_runtime_params_for_level(self, level):
+        runtime_params = RuntimeParameters()
+        runtime_params.number_of_timesteps = self._runtime_params.number_of_timesteps
+        runtime_params.number_of_threads = self._runtime_params.number_of_threads
+        runtime_params.timestep_log_interval = self._runtime_params.timestep_log_interval
+
+        if level > 0:
+            runtime_params.number_of_timesteps *= (level * 2)
+
+        return runtime_params
+
+    def clone_kernel_for_level(self, level):
+        kernel = LBMKernel(self._kernel.type)
+        kernel.use_forcing = self._kernel.use_forcing
+        kernel.forcing_in_x1 = self._kernel.forcing_in_x1
+        kernel.forcing_in_x2 = self._kernel.forcing_in_x2
+        kernel.forcing_in_x3 = self._kernel.forcing_in_x3
+
+        if level > 0:
+            kernel.forcing_in_x1 /= (level * 2)
+            kernel.forcing_in_x2 /= (level * 2)
+            kernel.forcing_in_x3 /= (level * 2)
+
+        return kernel
diff --git a/Python/norms.py b/Python/errors.py
similarity index 96%
rename from Python/norms.py
rename to Python/errors.py
index 78ae344591e4d91f28b9a98cf1e28ef447e2f62a..16e8c48ab9f0b7a46ed1372ef0b4d45738cccb1b 100644
--- a/Python/norms.py
+++ b/Python/errors.py
@@ -42,8 +42,8 @@ def mean_squared_error(real_values, numerical_values):
     return sum_of_squared_distances / num_values
 
 
-def l2_norm_error(real_values, numerical_values):
+def normalized_l2_error(real_values, numerical_values):
     sum_of_squared_distances = get_sum_of_squared_distances(real_values, numerical_values)
     sum_of_squared_real_values = sum(real_value ** 2 for real_value in real_values)
 
-    return math.sqrt(sum_of_squared_distances / sum_of_squared_real_values)
\ No newline at end of file
+    return math.sqrt(sum_of_squared_distances / sum_of_squared_real_values)
diff --git a/Python/poiseuille/analytical.py b/Python/poiseuille/analytical.py
index bca1a3ff95c1f28bca68ebb0c14efee48a1a5984..33e67595d94c50a1eb98751b7f10df9a031800e8 100644
--- a/Python/poiseuille/analytical.py
+++ b/Python/poiseuille/analytical.py
@@ -1,39 +1,52 @@
 from dataclasses import dataclass
 
 
+@dataclass
 class PoiseuilleSettings:
-
-    def __init__(self):
-        self.density = 1
-        self.viscosity = 0.005
-        self.height = 10
-        self.length = 1
-        self.pressure_in = 0
-        self.pressure_out = 0
-        self.force = 0
+    density = 1
+    viscosity = 0.005
+    height = 10
+    length = 1
+    pressure_in = 0
+    pressure_out = 0
+    force = 0
 
 
 def poiseuille_at_z(settings: PoiseuilleSettings, z: float):
     pressure_grad = ((settings.pressure_out - settings.pressure_in) / settings.length)
 
-    return (1 / settings.viscosity
+    return ((1 / settings.viscosity)
             * (- pressure_grad + settings.density * settings.force)
-            * z / 2 * (settings.height - z))
+            * (z / 2)
+            * (settings.height - z))
 
 
 def poiseuille_at_heights(settings: PoiseuilleSettings, heights):
     return [poiseuille_at_z(settings, z) for z in heights]
 
 
-if __name__ == '__main__':
-    # h1 = [1, 2, 3, 4, 5, 6, 7, 8, 9]
-    # h2 = [0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.5]
-    settings = PoiseuilleSettings()
-    settings.force = 1e-8
-    settings.height = 32
+def reynolds_number(settings: PoiseuilleSettings):
+    max_v = poiseuille_at_z(settings, settings.height / 2)
+    return max_v * settings.height / settings.viscosity
 
-    # print(max(poiseuille_at_heights(settings, h1)))
-    # print(max(poiseuille_at_heights(settings, h2)))
 
-    v = poiseuille_at_z(settings, 16)
-    print(v)
\ No newline at end of file
+if __name__ == '__main__':
+    sim_settings = PoiseuilleSettings()
+
+    sim_settings.force = 2e-7
+    sim_settings.viscosity = 1e-3
+    sim_settings.height = 16
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))
+
+    sim_settings.viscosity *= 2
+    sim_settings.height *= 2
+    sim_settings.force /= 2
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))
+
+    sim_settings.viscosity *= 2
+    sim_settings.height *= 2
+    sim_settings.force /= 2
+    print(f"v_max = ", poiseuille_at_z(sim_settings, sim_settings.height / 2))
+    print(f"Re =", reynolds_number(sim_settings))
diff --git a/Python/poiseuille/simulation.py b/Python/poiseuille/simulation.py
index f7f6d468f19820993c61b458b3b0138b8c886139..31ceb1ab9ef90fa4fd606bde4f47c45b8f7d7567 100644
--- a/Python/poiseuille/simulation.py
+++ b/Python/poiseuille/simulation.py
@@ -17,23 +17,30 @@ default_physical_params.lattice_viscosity = 0.005
 
 default_runtime_params = RuntimeParameters()
 default_runtime_params.number_of_threads = 4
-default_runtime_params.number_of_timesteps = 100000
-default_runtime_params.timestep_log_interval = 10000
+default_runtime_params.number_of_timesteps = 10000
+default_runtime_params.timestep_log_interval = 1000
+
+default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel.use_forcing = True
+default_kernel.forcing_in_x1 = 1e-8
+
+default_writer = Writer()
+default_writer.output_path = "./output"
+default_writer.output_format = OutputFormat.BINARY
+
+
+default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel.use_forcing = True
+default_kernel.forcing_in_x1 = 1e-8
 
 
 def run_simulation(physical_params=default_physical_params,
                    grid_params=default_grid_params,
-                   runtime_params=default_runtime_params):
+                   runtime_params=default_runtime_params,
+                   kernel=default_kernel,
+                   writer=default_writer):
     simulation = Simulation()
 
-    kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
-    kernel.use_forcing = True
-    kernel.forcing_in_x1 = 1e-8
-
-    writer = Writer()
-    writer.output_path = "./output"
-    writer.output_format = OutputFormat.BINARY
-
     simulation.set_kernel_config(kernel)
     simulation.set_physical_parameters(physical_params)
     simulation.set_grid_parameters(grid_params)
@@ -42,26 +49,26 @@ def run_simulation(physical_params=default_physical_params,
 
     no_slip_bc = NoSlipBoundaryCondition()
 
-    block_width = 3 * grid_params.node_distance
+    block_thickness = 3 * grid_params.node_distance
     simulation.add_object(
         GbCuboid3D(
-            grid_params.bounding_box.min_x1 - block_width,
-            grid_params.bounding_box.min_x2 - block_width,
-            grid_params.bounding_box.min_x3 - block_width,
-            grid_params.bounding_box.max_x1 + block_width,
-            grid_params.bounding_box.max_x2 + block_width,
+            grid_params.bounding_box.min_x1 - block_thickness,
+            grid_params.bounding_box.min_x2 - block_thickness,
+            grid_params.bounding_box.min_x3 - block_thickness,
+            grid_params.bounding_box.max_x1 + block_thickness,
+            grid_params.bounding_box.max_x2 + block_thickness,
             grid_params.bounding_box.min_x3),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMin")
 
     simulation.add_object(
         GbCuboid3D(
-            grid_params.bounding_box.min_x1 - block_width,
-            grid_params.bounding_box.min_x2 - block_width,
+            grid_params.bounding_box.min_x1 - block_thickness,
+            grid_params.bounding_box.min_x2 - block_thickness,
             grid_params.bounding_box.max_x3,
-            grid_params.bounding_box.max_x1 + block_width,
-            grid_params.bounding_box.max_x2 + block_width,
-            grid_params.bounding_box.max_x3 + block_width),
+            grid_params.bounding_box.max_x1 + block_thickness,
+            grid_params.bounding_box.max_x2 + block_thickness,
+            grid_params.bounding_box.max_x3 + block_thickness),
         no_slip_bc,
         State.SOLID, "/geo/addWallZMax")
 
diff --git a/Python/poiseuille/test_poiseuille_l2.py b/Python/poiseuille/test_poiseuille_l2.py
index 508fda5d6f551d428d314ffbfa657b4b6b2a230a..39c8b6dffe05e3c352e7fd340857e43d8d5a3dc8 100644
--- a/Python/poiseuille/test_poiseuille_l2.py
+++ b/Python/poiseuille/test_poiseuille_l2.py
@@ -1,68 +1,90 @@
+import os
 import shutil
 import unittest
 
 import matplotlib.pyplot as plt
+import numpy as np
 import pyvista as pv
+from pyfluids.kernel import LBMKernel, KernelType
 from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+from scipy import stats
 
-from norms import l2_norm_error
+from errors import normalized_l2_error
 from poiseuille.analytical import poiseuille_at_heights, PoiseuilleSettings
 from poiseuille.simulation import run_simulation
 from vtk_utilities import vertical_column_from_mesh, get_values_from_indices
 
 
 class TestPoiseuilleFlow(unittest.TestCase):
+    node_distances = [1, 0.5, 0.25]
+    number_of_nodes = [16, 32, 64]
+    number_of_timesteps = [2_500_000, 5_000_000, 10_000_000]
+    forcings = [1e-9, 5e-10, 2.5e-10]
+    viscosities = [1e-3, 2e-3, 4e-3]
+
+    def zipped_settings(self):
+        return zip(self.node_distances,
+                   self.number_of_nodes,
+                   self.number_of_timesteps,
+                   self.forcings,
+                   self.viscosities)
 
     def test_poiseuille_flow(self):
-        self.skipTest("Skipping test! This test has not been implemented correctly yet")
+        self.skipTest("This test is not implemented correctly yet")
         plt.ion()
 
-        channel_height = 10
-        number_of_nodes = [8, 16, 32]
-        number_of_timesteps = [10_000, 20_000, 40_000]
-        viscosities = [5e-3, 1e-2, 2e-2]
-        l2_norm_results = []
-
         physical_params = PhysicalParameters()
 
         runtime_params = RuntimeParameters()
-        runtime_params.number_of_threads = 4
-        runtime_params.timestep_log_interval = 1000
-
-        for test_number, nodes_in_column in enumerate(number_of_nodes):
-            runtime_params.number_of_timesteps = number_of_timesteps[test_number]
-            physical_params.lattice_viscosity = viscosities[test_number]
-            delta_x = channel_height / nodes_in_column
-            grid_params = create_grid_params_with_nodes_in_column(nodes_in_column, delta_x)
-            l2_norm_result = get_l2_norm_for_simulation(grid_params, physical_params, runtime_params)
-            l2_norm_results.append(l2_norm_result)
-
-        plt.plot(number_of_nodes, l2_norm_results)
+        runtime_params.number_of_threads = os.cpu_count()
+        runtime_params.timestep_log_interval = 10000
+
+        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel.use_forcing = True
+
+        normalized_l2_errors = []
+        for delta_x, nodes, timesteps, forcing, viscosity in self.zipped_settings():
+            physical_params.lattice_viscosity = viscosity
+            runtime_params.number_of_timesteps = timesteps
+            kernel.forcing_in_x1 = forcing
+
+            grid_params = create_grid_params_with_nodes_in_column(nodes, delta_x)
+            l2_error = get_l2_error_for_simulation(grid_params, physical_params, runtime_params, kernel)
+            normalized_l2_errors.append(l2_error)
+
+        nodes_as_log = [np.log10(node) for node in self.number_of_nodes]
+        l2_norms_as_log = [np.log10(l2) for l2 in normalized_l2_errors]
+        res = stats.linregress(nodes_as_log, l2_norms_as_log)
+
+        plt.xscale("log")
+        plt.yscale("log")
+        plt.plot(self.number_of_nodes, [np.power(10, res.intercept + res.slope * node) for node in nodes_as_log], 'r-')
+        plt.plot(self.number_of_nodes, normalized_l2_errors, "x:")
         plt.show()
 
-        self.assertTrue(l2_norm_results[1] <= l2_norm_results[0])
-        self.assertTrue(l2_norm_results[2] <= l2_norm_results[1])
-
-
-def run_simulation_with_settings(grid_params, physical_params, runtime_params, output_folder):
-    remove_existing_output_directory(output_folder)
-    run_simulation(physical_params, grid_params, runtime_params)
+        print(normalized_l2_errors)
+        self.assertAlmostEqual(res.slope, -2, places=2)
 
 
-def get_l2_norm_for_simulation(grid_params, physical_params, runtime_params):
+def get_l2_error_for_simulation(grid_params, physical_params, runtime_params, kernel):
     output_folder = "./output"
-    run_simulation_with_settings(grid_params, physical_params, runtime_params, output_folder)
+    run_simulation_with_settings(grid_params, physical_params, runtime_params, kernel, output_folder)
     heights = get_heights(output_folder, runtime_params)
 
-    numerical_results = get_numerical_results(runtime_params, output_folder, heights)
-    analytical_results = get_analytical_results(physical_params, heights, grid_params.number_of_nodes_per_direction[2])
+    numerical_results = get_numerical_results(runtime_params, output_folder)
+    analytical_results = get_analytical_results(grid_params, physical_params, kernel, heights)
 
     plt.plot(heights, numerical_results)
     plt.plot(heights, analytical_results)
     plt.legend(["numerical", "analytical"])
     plt.show()
 
-    return l2_norm_error(analytical_results, numerical_results)
+    return normalized_l2_error(analytical_results, numerical_results)
+
+
+def run_simulation_with_settings(grid_params, physical_params, runtime_params, kernel, output_folder):
+    shutil.rmtree(output_folder, ignore_errors=True)
+    run_simulation(physical_params, grid_params, runtime_params, kernel)
 
 
 def get_heights(output_folder, runtime_params):
@@ -72,7 +94,7 @@ def get_heights(output_folder, runtime_params):
     return heights
 
 
-def get_numerical_results(runtime_params, output_folder, heights):
+def get_numerical_results(runtime_params, output_folder):
     mesh_of_last_timestep = get_mesh_for_last_timestep(output_folder, runtime_params)
     velocities_in_x_direction = mesh_of_last_timestep.get_array("Vx")
     column_indices = vertical_column_from_mesh(mesh_of_last_timestep)
@@ -81,16 +103,12 @@ def get_numerical_results(runtime_params, output_folder, heights):
     return numerical_results
 
 
-def calculate_analytical_results(physical_params, height_values, channel_height):
-    settings = get_analytical_poiseuille_settings(channel_height, physical_params)
-    max_height = max(height_values)
-    height_values = [value / max_height * channel_height for value in height_values]
-    analytical_results = poiseuille_at_heights(settings, height_values)
-    return analytical_results
-
-
-def get_analytical_results(physical_params, heights, channel_height):
-    analytical_results = calculate_analytical_results(physical_params, heights, channel_height)
+def get_analytical_results(grid_params, physical_params, kernel, height_values):
+    channel_height = grid_params.number_of_nodes_per_direction[2]
+    settings = get_analytical_poiseuille_settings(channel_height, physical_params, kernel)
+    max_grid_height = channel_height * grid_params.node_distance
+    adjusted_height_values = [value / max_grid_height * channel_height for value in height_values]
+    analytical_results = poiseuille_at_heights(settings, adjusted_height_values)
     return analytical_results
 
 
@@ -100,18 +118,12 @@ def get_mesh_for_last_timestep(output_folder, runtime_params):
     return mesh_of_last_timestep
 
 
-def remove_existing_output_directory(output_dir):
-    shutil.rmtree(output_dir, ignore_errors=True)
-
-
-def get_analytical_poiseuille_settings(height, physical_params):
+def get_analytical_poiseuille_settings(height, physical_params, kernel):
     settings = PoiseuilleSettings()
     settings.height = height
     settings.viscosity = physical_params.lattice_viscosity
     settings.density = 1
-    settings.force = 1e-8
-
-    # print(settings)
+    settings.force = kernel.forcing_in_x1
 
     return settings
 
@@ -130,13 +142,10 @@ def get_heights_from_indices(mesh, indices):
 def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x):
     grid_params = GridParameters()
     grid_params.node_distance = delta_x
-    grid_params.number_of_nodes_per_direction = [2, 2, nodes_in_column]
-    grid_params.blocks_per_direction = [1, 1, 6]
+    grid_params.number_of_nodes_per_direction = [1, 1, nodes_in_column]
+    grid_params.blocks_per_direction = [1, 1, 8]
     grid_params.periodic_boundary_in_x1 = True
     grid_params.periodic_boundary_in_x2 = True
     grid_params.periodic_boundary_in_x3 = False
 
-    print(f"GridParameters.node_distance = {grid_params.node_distance}")
-    print(f"GridParameters.number_of_nodes_per_direction = {grid_params.number_of_nodes_per_direction}")
-
     return grid_params
diff --git a/Python/requirements.txt b/Python/requirements.txt
index ded26051abb63145842f304ba865dc8487a8be73..8628634d1b85ebc0b07328d563d479f35641be97 100644
--- a/Python/requirements.txt
+++ b/Python/requirements.txt
@@ -14,8 +14,10 @@ py==1.10.0
 pyparsing==2.4.7
 pytest==6.2.1
 python-dateutil==2.8.1
-pyvista==0.27.4
+pyvista==0.28.1
+scipy==1.6.1
 scooby==0.5.6
 six==1.15.0
 toml==0.10.2
+transforms3d==0.3.1
 vtk==9.0.1
diff --git a/Python/tests/test_acousticscaling.py b/Python/tests/test_acousticscaling.py
new file mode 100644
index 0000000000000000000000000000000000000000..2da5314529f9559f9ac316f2d1bb3f1a9d0e1211
--- /dev/null
+++ b/Python/tests/test_acousticscaling.py
@@ -0,0 +1,115 @@
+import unittest
+from typing import List
+
+from pyfluids.kernel import LBMKernel, KernelType
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+from acousticscaling import OneDirectionalAcousticScaling
+
+
+class OneDirectionalAcousticScalingTest(unittest.TestCase):
+
+    def setUp(self) -> None:
+        self.grid_params = self.make_grid_params()
+        self.physical_params = self.make_physical_params()
+        self.runtime_params = self.make_runtime_params()
+        self.kernel = self.make_kernel()
+
+        self.sut = OneDirectionalAcousticScaling(self.grid_params, self.physical_params, self.runtime_params,
+                                                 self.kernel)
+
+    def test_given_sim_parameters__when_scaling_level_zero__should_return_equal_sim_parameters(self):
+        factor = 1
+        actual_params = self.sut.configuration_for_scale_level(0)
+        actual_grid_params = actual_params[0]
+        actual_physical_params = actual_params[factor]
+        actual_runtime_params = actual_params[2]
+        actual_kernel = actual_params[3]
+
+        self.assert_parameters_scaled_by_factor(actual_grid_params, actual_kernel,
+                                                actual_physical_params, actual_runtime_params, factor)
+
+    def test_given_sim_parameters__when_scaling_level_one__should_return_sim_parameters_scaled_by_two(self):
+        actual_params = self.sut.configuration_for_scale_level(1)
+        actual_grid_params = actual_params[0]
+        actual_physical_params = actual_params[1]
+        actual_runtime_params = actual_params[2]
+        actual_kernel = actual_params[3]
+
+        self.assert_parameters_scaled_by_factor(actual_grid_params, actual_kernel,
+                                                actual_physical_params, actual_runtime_params, 2)
+
+    def assert_parameters_scaled_by_factor(self, actual_grid_params, actual_kernel,
+                                           actual_physical_params, actual_runtime_params, factor):
+        self.assert_grid_params_scaled_by_factor(actual_grid_params, factor=factor)
+        self.assert_physical_params_scaled_by_factor(actual_physical_params, factor=factor)
+        self.assert_runtime_params_scaled_by_factor(actual_runtime_params, factor=factor)
+        self.assert_kernel_forcing_scaled_by_factor(actual_kernel, factor=factor)
+
+    def assert_grid_params_scaled_by_factor(self, actual_grid_params: GridParameters, factor: int):
+        expected_nodes_per_direction = self.scaled_list(self.grid_params.number_of_nodes_per_direction, factor)
+        expected_blocks_per_direction = self.scaled_list(self.grid_params.blocks_per_direction, factor)
+        expected_node_distance = self.grid_params.node_distance / factor
+        self.assertEqual(expected_node_distance, actual_grid_params.node_distance)
+        self.assertEqual(expected_nodes_per_direction, actual_grid_params.number_of_nodes_per_direction)
+        self.assertEqual(expected_blocks_per_direction, actual_grid_params.blocks_per_direction)
+        self.assertEqual(self.grid_params.reference_direction_index, actual_grid_params.reference_direction_index)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x1, actual_grid_params.periodic_boundary_in_x1)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x2, actual_grid_params.periodic_boundary_in_x2)
+        self.assertEqual(self.grid_params.periodic_boundary_in_x3, actual_grid_params.periodic_boundary_in_x3)
+
+    def assert_physical_params_scaled_by_factor(self, actual_params: PhysicalParameters, factor: int):
+        self.assertEqual(self.physical_params.lattice_viscosity * factor, actual_params.lattice_viscosity)
+        self.assertEqual(self.physical_params.bulk_viscosity_factor, actual_params.bulk_viscosity_factor)
+
+    def assert_runtime_params_scaled_by_factor(self, actual_params: RuntimeParameters, factor: int):
+        self.assertEqual(self.runtime_params.number_of_timesteps * factor, actual_params.number_of_timesteps)
+        self.assertEqual(self.runtime_params.number_of_threads, actual_params.number_of_threads)
+        self.assertEqual(self.runtime_params.timestep_log_interval, actual_params.timestep_log_interval)
+
+    def assert_kernel_forcing_scaled_by_factor(self, actual_kernel: LBMKernel, factor: int):
+        self.assertEqual(self.kernel.type, actual_kernel.type)
+        self.assertEqual(self.kernel.use_forcing, actual_kernel.use_forcing)
+        self.assertAlmostEqual(self.kernel.forcing_in_x1 / factor, actual_kernel.forcing_in_x1)
+        self.assertAlmostEqual(self.kernel.forcing_in_x2, actual_kernel.forcing_in_x2)
+        self.assertAlmostEqual(self.kernel.forcing_in_x3, actual_kernel.forcing_in_x3)
+
+    @staticmethod
+    def scaled_list(list_to_scale: List[int], factor: int) -> List[int]:
+        return [list_to_scale[0], list_to_scale[1], list_to_scale[2] * factor]
+
+    @staticmethod
+    def make_kernel():
+        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel.use_forcing = True
+        kernel.forcing_in_x1 = 5e-10
+        return kernel
+
+    @staticmethod
+    def make_runtime_params():
+        runtime_params = RuntimeParameters()
+        runtime_params.number_of_threads = 4
+        runtime_params.number_of_timesteps = 4_000_000
+        runtime_params.timestep_log_interval = 1_000_000
+        return runtime_params
+
+    @staticmethod
+    def make_physical_params():
+        physical_params = PhysicalParameters()
+        physical_params.lattice_viscosity = 1e-4
+        return physical_params
+
+    @staticmethod
+    def make_grid_params():
+        grid_params = GridParameters()
+        grid_params.node_distance = 1
+        grid_params.number_of_nodes_per_direction = [1, 1, 16]
+        grid_params.blocks_per_direction = [1, 1, 16]
+        grid_params.periodic_boundary_in_x1 = True
+        grid_params.periodic_boundary_in_x2 = True
+
+        return grid_params
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/Python/tests/test_boundaryconditions.py b/Python/tests/test_boundaryconditions.py
new file mode 100644
index 0000000000000000000000000000000000000000..5a7d61f36337398fc5621540951f15b72262b17b
--- /dev/null
+++ b/Python/tests/test_boundaryconditions.py
@@ -0,0 +1,61 @@
+import unittest
+from pyfluids.boundaryconditions import *
+
+
+class BoundaryConditionsTest(unittest.TestCase):
+
+    def test__can_create_no_slip_bc(self):
+        """
+        Should be able to create NoSlipBoundaryCondition
+        """
+        sut = NoSlipBoundaryCondition()
+
+    def test__can_create_velocity_bc(self):
+        """
+        Should be able to create VelocityBoundaryCondition
+        """
+        sut = VelocityBoundaryCondition()
+
+    def test__can_create_velocity_bc_with_directions_function_and_time(self):
+        """
+        Should be able to create VelocityBoundaryCondition with directions, function and start/end time
+        """
+        from pymuparser import Parser
+
+        parser = Parser()
+        parser.expression = "1"
+        sut = VelocityBoundaryCondition(True, True, True, parser, 0, 1)
+
+    def test__can_create_velocity_bc_with_directions__function_per_direction__and__time(self):
+        """
+        Should be able to create VelocityBoundaryCondition with directions, function per direction and start/end time
+        """
+        from pymuparser import Parser
+
+        f1 = Parser()
+        f1.expression = "1"
+
+        f2 = Parser()
+        f2.expression = "1"
+
+        f3 = Parser()
+        f3.expression = "1"
+        sut = VelocityBoundaryCondition(True, True, True, f1, f2, f3, 0, 1)
+
+    def test__can_create_velocity_bc_with_speeds_and_times_per_direction(self):
+        """
+        Should be able to create VelocityBoundaryCondition with speeds and start/end times per direction
+        """
+        vx1, vx2, vx3 = 1, 2, 3
+        start1, end1 = 0, 1
+        start2, end2 = 1, 2
+        start3, end3 = 2, 3
+
+        sut = VelocityBoundaryCondition(vx1, start1, end1, vx2, start2, end2, vx3, start3, end3)
+
+    def test__can_create_non_reflecting_outflow(self):
+        """
+        Should be able to create NonReflectingOutflow
+        """
+
+        sut = NonReflectingOutflow()
diff --git a/Python/test_geometry.py b/Python/tests/test_geometry.py
similarity index 100%
rename from Python/test_geometry.py
rename to Python/tests/test_geometry.py
diff --git a/Python/test_kernel.py b/Python/tests/test_kernel.py
similarity index 85%
rename from Python/test_kernel.py
rename to Python/tests/test_kernel.py
index 5e5487903a279501b0acdba1ab192ac8dd381383..b1016f15308a77c9025787e061355819cbca3874 100644
--- a/Python/test_kernel.py
+++ b/Python/tests/test_kernel.py
@@ -51,3 +51,10 @@ class TestLBMKernel(unittest.TestCase):
         self.assertEqual(self.sut.forcing_in_x2, 8)
         self.assertEqual(self.sut.forcing_in_x3, 5)
 
+    def test_lbm_kernel__when_getting_type__should_equal_kernel_type_enum_value(self) -> None:
+        """
+        WHEN getting the kernel type IT should equal the corresponding KernelType enum value
+        """
+
+        actual = self.sut.type
+        self.assertEqual(KernelType.BGK, actual)
diff --git a/Python/vtk_utilities.py b/Python/vtk_utilities.py
index e8c16ce02c85518edacde024a012c4ce261fbc14..14d13a40d85e669a463c0b5b31ca2a21e771af75 100644
--- a/Python/vtk_utilities.py
+++ b/Python/vtk_utilities.py
@@ -1,20 +1,29 @@
 import math
+from typing import List
+
 import pyvista as pv
 
 
 def vertical_column_from_mesh(mesh):
     last_seen = math.inf
     relevant_indices = []
+    first_x = 0
+    first_y = 0
     for index, point in enumerate(mesh.points):
-        if point[2] == last_seen:
+        if index == 0:
+            first_x = point[0]
+            first_y = point[1]
+
+        if (point[0] != first_x or point[1] != first_y) and point[2] == last_seen:
             continue
+
         relevant_indices.append(index)
         last_seen = point[2]
 
     return relevant_indices
 
 
-def get_values_from_indices(array, indices):
+def get_values_from_indices(array, indices) -> List[float]:
     return [array[index] for index in indices]
 
 
diff --git a/setup.py b/setup.py
index 6f4ff71b014dbafa8800af234cf715f5207fb2a8..11676929cf55ad17742d37a7679a4d3e91ebb436 100644
--- a/setup.py
+++ b/setup.py
@@ -11,6 +11,7 @@ from distutils.version import LooseVersion
 vf_cmake_args = [
     "-DBUILD_VF_PYTHON_BINDINGS=ON",
     "-DBUILD_VF_CPU:BOOL=ON",
+    "-DBUILD_VF_GPU:BOOL=OFF",
     "-DUSE_METIS=ON",
     "-DUSE_MPI=ON",
     "-DBUILD_SHARED_LIBS=OFF",
diff --git a/src/cpu/pythonbindings/src/VirtualfluidsModule.cpp b/src/cpu/pythonbindings/src/VirtualfluidsModule.cpp
index 61c6005b9fbdd1108d9ffff374763cad95a5c40e..564dc1838d48a92340fa5491779177b299bcb270 100644
--- a/src/cpu/pythonbindings/src/VirtualfluidsModule.cpp
+++ b/src/cpu/pythonbindings/src/VirtualfluidsModule.cpp
@@ -6,14 +6,17 @@
 #include "submodules/simulationparameters.cpp"
 #include "submodules/writer.cpp"
 
-namespace py = pybind11;
-
-PYBIND11_MODULE(pyfluids, m)
+namespace py_bindings
 {
-    makeBoundaryConditionsModule(m);
-    makeSimulationModule(m);
-    makeGeometryModule(m);
-    makeKernelModule(m);
-    makeParametersModule(m);
-    makeWriterModule(m);
+    namespace py = pybind11;
+
+    PYBIND11_MODULE(pyfluids, m)
+    {
+        boundaryconditions::makeModule(m);
+        simulation::makeModule(m);
+        geometry::makeModule(m);
+        kernel::makeModule(m);
+        parameters::makeModule(m);
+        writer::makeModule(m);
+    }
 }
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/boundaryconditions.cpp b/src/cpu/pythonbindings/src/submodules/boundaryconditions.cpp
index ef8b923fd9aa2feea9c81da686df38cb891a4f84..3bff7bc069ca20fe1c0cf3d1847b9714e0381505 100644
--- a/src/cpu/pythonbindings/src/submodules/boundaryconditions.cpp
+++ b/src/cpu/pythonbindings/src/submodules/boundaryconditions.cpp
@@ -7,56 +7,58 @@
 #include <BoundaryConditions/VelocityBCAdapter.h>
 #include <BoundaryConditions/NoSlipBCAlgorithm.h>
 #include <BoundaryConditions/VelocityBCAlgorithm.h>
+#include <BoundaryConditions/HighViscosityNoSlipBCAlgorithm.h>
 
-namespace py = pybind11;
-using namespace py::literals;
+namespace boundaryconditions
+{
+    namespace py = pybind11;
+    using namespace py::literals;
 
-template<class Adapter>
-using py_bc_class = py::class_<Adapter, BCAdapter, std::shared_ptr<Adapter>>;
+    template<class adapter, class algorithm,
+            class = std::enable_if_t<std::is_base_of<BCAdapter, adapter>::value>,
+            class = std::enable_if_t<std::is_base_of<BCAlgorithm, algorithm>::value>>
+    class PyBoundaryCondition : public adapter
+    {
+    public:
+        template<typename ...Args>
+        PyBoundaryCondition(Args &&... args) : adapter(std::forward<Args>(args)...)
+        {
+            this->setBcAlgorithm(std::make_shared<algorithm>());
+        }
+    };
 
-template<class Adapter, class Algorithm, typename ...Args>
-std::shared_ptr<Adapter> create(Args... args)
-{
-    auto adapter = std::make_shared<Adapter>(args...);
-    adapter->setBcAlgorithm(std::make_shared<Algorithm>());
-    return adapter;
-}
+    template<class adapter, class algorithm>
+    using bc_class = py::class_<PyBoundaryCondition<adapter, algorithm>, BCAdapter,
+            std::shared_ptr<PyBoundaryCondition<adapter, algorithm>>>;
 
-template<class Algorithm>
-void add_constructors_to_velocity_bc(py_bc_class<VelocityBCAdapter> &cls)
-{
-    auto firstConstructor = &create<VelocityBCAdapter, Algorithm, bool &, bool &, bool &, mu::Parser &, double &, double &>;
-    auto secondConstructor = &create<VelocityBCAdapter, Algorithm,
-            bool &, bool &, bool &, mu::Parser &, mu::Parser &, mu::Parser &, double &, double &>;
-    auto thirdConstructor = &create<VelocityBCAdapter, Algorithm,
-            double &, double &, double &, double &, double &, double &, double &, double &, double &>;
-
-    cls.def(py::init(&create<VelocityBCAdapter, Algorithm>))
-            .def(py::init(firstConstructor),
-                 "vx1"_a, "vx2"_a, "vx3"_a, "function"_a, "start_time"_a, "end_time"_a)
-            .def(py::init(secondConstructor),
-                 "vx1"_a, "vx2"_a, "vx3"_a,
-                 "function_vx1"_a, "function_vx2"_a, "function_vx2"_a,
-                 "start_time"_a, "end_time"_a)
-            .def(py::init(thirdConstructor),
-                 "vx1"_a, "vx1_start_time"_a, "vx1_end_time"_a,
-                 "vx2"_a, "vx2_start_time"_a, "vx2_end_time"_a,
-                 "vx3"_a, "vx3_start_time"_a, "vx3_end_time"_a);
-}
-
-void makeBoundaryConditionsModule(py::module_ &parentModule)
-{
-    py::module_ bcModule = parentModule.def_submodule("boundaryconditions");
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module_ bcModule = parentModule.def_submodule("boundaryconditions");
+
+        auto _ = py::class_<BCAdapter, std::shared_ptr<BCAdapter>>(bcModule, "BCAdapter");
 
-    py::class_<BCAdapter, std::shared_ptr<BCAdapter>>(bcModule, "BCAdapter");
+        bc_class<NoSlipBCAdapter, NoSlipBCAlgorithm>(bcModule, "NoSlipBoundaryCondition")
+                .def(py::init());
 
-    py_bc_class<NoSlipBCAdapter>(bcModule, "NoSlipBoundaryCondition")
-            .def(py::init(&create<NoSlipBCAdapter, NoSlipBCAlgorithm>));
+        bc_class<NoSlipBCAdapter, HighViscosityNoSlipBCAlgorithm>(bcModule, "HighViscosityNoSlipBoundaryCondition")
+                .def(py::init());
 
-    auto velocityBoundaryCondition = py_bc_class<VelocityBCAdapter>(bcModule, "VelocityBoundaryCondition");
-    add_constructors_to_velocity_bc<VelocityBCAlgorithm>(velocityBoundaryCondition);
+        bc_class<VelocityBCAdapter, VelocityBCAlgorithm>(bcModule, "VelocityBoundaryCondition")
+                .def(py::init())
+                .def(py::init<bool &, bool &, bool &, mu::Parser &, double &, double &>(),
+                     "vx1"_a, "vx2"_a, "vx3"_a,
+                     "function"_a, "start_time"_a, "end_time"_a)
+                .def(py::init<bool &, bool &, bool &, mu::Parser &, mu::Parser &, mu::Parser &, double &, double &>(),
+                     "vx1"_a, "vx2"_a, "vx3"_a,
+                     "function_vx1"_a, "function_vx2"_a, "function_vx2"_a,
+                     "start_time"_a, "end_time"_a)
+                .def(py::init<double &, double &, double &, double &, double &, double &, double &, double &, double &>(),
+                     "vx1"_a, "vx1_start_time"_a, "vx1_end_time"_a,
+                     "vx2"_a, "vx2_start_time"_a, "vx2_end_time"_a,
+                     "vx3"_a, "vx3_start_time"_a, "vx3_end_time"_a);
 
-    py_bc_class<DensityBCAdapter>(bcModule, "DensityBoundaryCondition")
-            .def(py::init(&create<DensityBCAdapter, NonReflectingOutflowBCAlgorithm>));
-}
+        bc_class<DensityBCAdapter, NonReflectingOutflowBCAlgorithm>(bcModule, "NonReflectingOutflow")
+                .def(py::init());
+    }
 
+}
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/geometry.cpp b/src/cpu/pythonbindings/src/submodules/geometry.cpp
index 884ced7b92ddae8a30a8f482e4b22dcbfa37beec..b7ff4dd761258d41687589d2dd89c3479093753e 100644
--- a/src/cpu/pythonbindings/src/submodules/geometry.cpp
+++ b/src/cpu/pythonbindings/src/submodules/geometry.cpp
@@ -5,76 +5,80 @@
 #include <geometry3d/GbLine3D.h>
 #include <Interactors/Interactor3D.h>
 
-namespace py = pybind11;
 
-template<class GeoObject>
-using py_geometry = py::class_<GeoObject, GbObject3D, std::shared_ptr<GeoObject>>;
-
-std::string GbPoint3D_repr_(const GbPoint3D &instance)
+namespace geometry
 {
-    std::ostringstream stream;
-    stream << "<GbPoint3D"
-           << " x1: " << instance.getX1Coordinate()
-           << " x2: " << instance.getX2Coordinate()
-           << " x3: " << instance.getX3Coordinate() << ">";
+    namespace py = pybind11;
 
-    return stream.str();
-}
+    template<class GeoObject>
+    using py_geometry = py::class_<GeoObject, GbObject3D, std::shared_ptr<GeoObject>>;
 
-void makeGeometryModule(py::module_ &parentModule)
-{
-    py::module geometry = parentModule.def_submodule("geometry");
+    std::string GbPoint3D_repr_(const GbPoint3D &instance)
+    {
+        std::ostringstream stream;
+        stream << "<GbPoint3D"
+               << " x1: " << instance.getX1Coordinate()
+               << " x2: " << instance.getX2Coordinate()
+               << " x3: " << instance.getX3Coordinate() << ">";
+
+        return stream.str();
+    }
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module geometry = parentModule.def_submodule("geometry");
+
+        py::class_<GbObject3D, std::shared_ptr<GbObject3D>>(geometry, "GbObject3D");
 
-    py::class_<GbObject3D, std::shared_ptr<GbObject3D>>(geometry, "GbObject3D");
+        py_geometry<GbPoint3D>(geometry, "GbPoint3D")
+                .def(py::init())
+                .def(py::init<double &, double &, double &>())
+                .def(py::init<GbPoint3D *>())
+                .def_property("x1", &GbPoint3D::getX1Coordinate, &GbPoint3D::setX1)
+                .def_property("x2", &GbPoint3D::getX2Coordinate, &GbPoint3D::setX2)
+                .def_property("x3", &GbPoint3D::getX3Coordinate, &GbPoint3D::setX3)
+                .def("get_distance", &GbPoint3D::getDistance)
+                .def("__repr__", &GbPoint3D_repr_);
 
-    py_geometry<GbPoint3D>(geometry, "GbPoint3D")
-            .def(py::init())
-            .def(py::init<double &, double &, double &>())
-            .def(py::init<GbPoint3D *>())
-            .def_property("x1", &GbPoint3D::getX1Coordinate, &GbPoint3D::setX1)
-            .def_property("x2", &GbPoint3D::getX2Coordinate, &GbPoint3D::setX2)
-            .def_property("x3", &GbPoint3D::getX3Coordinate, &GbPoint3D::setX3)
-            .def("get_distance", &GbPoint3D::getDistance)
-            .def("__repr__", [&](const GbPoint3D &instance)
-            { return GbPoint3D_repr_(instance); });
+        py_geometry<GbCuboid3D>(geometry, "GbCuboid3D")
+                .def(py::init())
+                .def(py::init<double &, double &, double &, double &, double &, double &>())
+                .def(py::init<GbPoint3D *, GbPoint3D *>())
+                .def(py::init<GbCuboid3D *>())
+                .def_property("point1", &GbCuboid3D::getPoint1, &GbCuboid3D::setPoint1)
+                .def_property("point2", &GbCuboid3D::getPoint2, &GbCuboid3D::setPoint2)
+                .def("__repr__", [&](GbCuboid3D &instance)
+                {
+                    std::ostringstream stream;
+                    stream << "<GbCuboid3D" << std::endl
+                           << "point1: " << GbPoint3D_repr_(instance.getPoint1()) << std::endl
+                           << "point2: " << GbPoint3D_repr_(instance.getPoint2()) << ">";
+                    return stream.str();
+                });
 
-    py_geometry<GbCuboid3D>(geometry, "GbCuboid3D")
-            .def(py::init())
-            .def(py::init<double &, double &, double &, double &, double &, double &>())
-            .def(py::init<GbPoint3D *, GbPoint3D *>())
-            .def(py::init<GbCuboid3D *>())
-            .def_property("point1", &GbCuboid3D::getPoint1, &GbCuboid3D::setPoint1)
-            .def_property("point2", &GbCuboid3D::getPoint2, &GbCuboid3D::setPoint2)
-            .def("__repr__", [&](GbCuboid3D instance)
-            {
-                std::ostringstream stream;
-                stream << "<GbCuboid3D" << std::endl
-                       << "point1: " << GbPoint3D_repr_(instance.getPoint1()) << std::endl
-                       << "point2: " << GbPoint3D_repr_(instance.getPoint2()) << ">";
-                return stream.str();
-            });
+        py_geometry<GbLine3D>(geometry, "GbLine3D")
+                .def(py::init())
+                .def(py::init<GbPoint3D *, GbPoint3D *>())
+                .def(py::init<GbLine3D>())
+                .def_property("point1", &GbLine3D::getPoint1, &GbLine3D::setPoint1)
+                .def_property("point2", &GbLine3D::getPoint2, &GbLine3D::setPoint2)
+                .def("__repr__", [&](GbLine3D &instance)
+                {
+                    std::ostringstream stream;
+                    stream << "<GbLine3D" << std::endl
+                           << "point1: " << GbPoint3D_repr_(instance.getPoint1()) << std::endl
+                           << "point2: " << GbPoint3D_repr_(instance.getPoint2()) << ">";
+                    return stream.str();
+                });
 
-    py_geometry<GbLine3D>(geometry, "GbLine3D")
-            .def(py::init())
-            .def(py::init<GbPoint3D *, GbPoint3D *>())
-            .def(py::init<GbLine3D>())
-            .def_property("point1", &GbLine3D::getPoint1, &GbLine3D::setPoint1)
-            .def_property("point2", &GbLine3D::getPoint2, &GbLine3D::setPoint2)
-            .def("__repr__", [&](GbLine3D instance)
-            {
-                std::ostringstream stream;
-                stream << "<GbLine3D" << std::endl
-                       << "point1: " << GbPoint3D_repr_(instance.getPoint1()) << std::endl
-                       << "point2: " << GbPoint3D_repr_(instance.getPoint2()) << ">";
-                return stream.str();
-            });
 
+        py::class_<Interactor3D, std::shared_ptr<Interactor3D>>(geometry, "State")
+                .def_readonly_static("SOLID", &Interactor3D::SOLID)
+                .def_readonly_static("INVERSESOLID", &Interactor3D::INVERSESOLID)
+                .def_readonly_static("TIMEDEPENDENT", &Interactor3D::TIMEDEPENDENT)
+                .def_readonly_static("FLUID", &Interactor3D::FLUID)
+                .def_readonly_static("MOVEABLE", &Interactor3D::MOVEABLE)
+                .def_readonly_static("CHANGENOTNECESSARY", &Interactor3D::CHANGENOTNECESSARY);
+    }
 
-    py::class_<Interactor3D, std::shared_ptr<Interactor3D>>(geometry, "State")
-            .def_readonly_static("SOLID", &Interactor3D::SOLID)
-            .def_readonly_static("INVERSESOLID", &Interactor3D::INVERSESOLID)
-            .def_readonly_static("TIMEDEPENDENT", &Interactor3D::TIMEDEPENDENT)
-            .def_readonly_static("FLUID", &Interactor3D::FLUID)
-            .def_readonly_static("MOVEABLE", &Interactor3D::MOVEABLE)
-            .def_readonly_static("CHANGENOTNECESSARY", &Interactor3D::CHANGENOTNECESSARY);
 }
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/kernel.cpp b/src/cpu/pythonbindings/src/submodules/kernel.cpp
index 0e2e23c8c6a8ff56ea01c277e9825e7eb78c5c02..fb291790632cc2041410f60a14fca8d966283343 100644
--- a/src/cpu/pythonbindings/src/submodules/kernel.cpp
+++ b/src/cpu/pythonbindings/src/submodules/kernel.cpp
@@ -3,39 +3,43 @@
 #include <simulationconfig/KernelFactory.h>
 #include <simulationconfig/KernelConfigStructs.h>
 
-
-namespace py = pybind11;
-
-
-void makeKernelModule(py::module_ &parentModule)
+namespace kernel
 {
-    py::module kernelModule = parentModule.def_submodule("kernel");
-
-    py::enum_<KernelFactory::KernelType>(kernelModule, "KernelType")
-            .value("BGK", KernelFactory::BGK)
-            .value("CompressibleCumulantFourthOrderViscosity",
-                   KernelFactory::COMPRESSIBLE_CUMULANT_4TH_ORDER_VISCOSITY);
-
-
-    py::class_<LBMKernelConfiguration, std::shared_ptr<LBMKernelConfiguration>>(kernelModule, "LBMKernel")
-            .def(py::init<KernelFactory::KernelType>())
-            .def_readwrite("use_forcing", &LBMKernelConfiguration::useForcing)
-            .def_readwrite("forcing_in_x1", &LBMKernelConfiguration::forcingX1)
-            .def_readwrite("forcing_in_x2", &LBMKernelConfiguration::forcingX2)
-            .def_readwrite("forcing_in_x3", &LBMKernelConfiguration::forcingX3)
-            .def("set_forcing", [](LBMKernelConfiguration &kernelConfig, double x1, double x2, double x3) {
-                kernelConfig.forcingX1 = x1;
-                kernelConfig.forcingX2 = x2;
-                kernelConfig.forcingX3 = x3;
-            })
-            .def("__repr__", [](LBMKernelConfiguration &kernelConfig) {
-                std::ostringstream stream;
-                stream << "<" << kernelConfig.kernelType << std::endl
-                       << "Use forcing: " << kernelConfig.useForcing << std::endl
-                       << "Forcing in x1: " << kernelConfig.forcingX1 << std::endl
-                       << "Forcing in x2: " << kernelConfig.forcingX2 << std::endl
-                       << "Forcing in x3: " << kernelConfig.forcingX3 << ">" << std::endl;
+    namespace py = pybind11;
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module kernelModule = parentModule.def_submodule("kernel");
+
+        py::enum_<KernelFactory::KernelType>(kernelModule, "KernelType")
+                .value("BGK", KernelFactory::BGK)
+                .value("CompressibleCumulantFourthOrderViscosity",
+                       KernelFactory::COMPRESSIBLE_CUMULANT_4TH_ORDER_VISCOSITY);
+
+        py::class_<LBMKernelConfiguration, std::shared_ptr<LBMKernelConfiguration>>(kernelModule, "LBMKernel")
+                .def(py::init<KernelFactory::KernelType>())
+                .def_readwrite("type", &LBMKernelConfiguration::kernelType)
+                .def_readwrite("use_forcing", &LBMKernelConfiguration::useForcing)
+                .def_readwrite("forcing_in_x1", &LBMKernelConfiguration::forcingX1)
+                .def_readwrite("forcing_in_x2", &LBMKernelConfiguration::forcingX2)
+                .def_readwrite("forcing_in_x3", &LBMKernelConfiguration::forcingX3)
+                .def("set_forcing", [](LBMKernelConfiguration &kernelConfig, double x1, double x2, double x3)
+                {
+                    kernelConfig.forcingX1 = x1;
+                    kernelConfig.forcingX2 = x2;
+                    kernelConfig.forcingX3 = x3;
+                })
+                .def("__repr__", [](LBMKernelConfiguration &kernelConfig)
+                {
+                    std::ostringstream stream;
+                    stream << "<" << kernelConfig.kernelType << std::endl
+                           << "Use forcing: " << kernelConfig.useForcing << std::endl
+                           << "Forcing in x1: " << kernelConfig.forcingX1 << std::endl
+                           << "Forcing in x2: " << kernelConfig.forcingX2 << std::endl
+                           << "Forcing in x3: " << kernelConfig.forcingX3 << ">" << std::endl;
+
+                    return stream.str();
+                });
+    }
 
-                return stream.str();
-            });
 }
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/simulationconfig.cpp b/src/cpu/pythonbindings/src/submodules/simulationconfig.cpp
index a525691bdf8dd9a4a0a3994b607816fba509c7a5..60af4e36af4dca67e9262dd9f5ee1f46d5b7bb58 100644
--- a/src/cpu/pythonbindings/src/submodules/simulationconfig.cpp
+++ b/src/cpu/pythonbindings/src/submodules/simulationconfig.cpp
@@ -1,18 +1,22 @@
 #include <pybind11/pybind11.h>
 #include <simulationconfig/Simulation.h>
 
-namespace py = pybind11;
-
-void makeSimulationModule(py::module_ &parentModule)
+namespace simulation
 {
-    py::class_<Simulation, std::shared_ptr<Simulation>>(parentModule, "Simulation")
-            .def(py::init())
-            .def("set_writer", &Simulation::setWriterConfiguration)
-            .def("set_grid_parameters", &Simulation::setGridParameters)
-            .def("set_physical_parameters", &Simulation::setPhysicalParameters)
-            .def("set_runtime_parameters", &Simulation::setRuntimeParameters)
-            .def("set_kernel_config", &Simulation::setKernelConfiguration)
-            .def("add_object", &Simulation::addObject)
-            .def("add_bc_adapter", &Simulation::addBCAdapter)
-            .def("run_simulation", &Simulation::run);
+    namespace py = pybind11;
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::class_<Simulation, std::shared_ptr<Simulation>>(parentModule, "Simulation")
+                .def(py::init())
+                .def("set_writer", &Simulation::setWriterConfiguration)
+                .def("set_grid_parameters", &Simulation::setGridParameters)
+                .def("set_physical_parameters", &Simulation::setPhysicalParameters)
+                .def("set_runtime_parameters", &Simulation::setRuntimeParameters)
+                .def("set_kernel_config", &Simulation::setKernelConfiguration)
+                .def("add_object", &Simulation::addObject)
+                .def("add_bc_adapter", &Simulation::addBCAdapter)
+                .def("run_simulation", &Simulation::run);
+    }
+
 }
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/simulationparameters.cpp b/src/cpu/pythonbindings/src/submodules/simulationparameters.cpp
index f59d1c0ec0473c537f0cc334044bcd113f822687..acc272f2ee412cfbafd9007b4b18610cfd0a1e9b 100644
--- a/src/cpu/pythonbindings/src/submodules/simulationparameters.cpp
+++ b/src/cpu/pythonbindings/src/submodules/simulationparameters.cpp
@@ -3,54 +3,57 @@
 #include <complex>
 #include <simulationconfig/SimulationParameters.h>
 
-namespace py = pybind11;
-
-void makeParametersModule(py::module_ &parentModule)
+namespace parameters
 {
-    py::module parametersModule = parentModule.def_submodule("parameters");
-
-    py::class_<PhysicalParameters, std::shared_ptr<PhysicalParameters>>(parametersModule, "PhysicalParameters")
-            .def(py::init())
-            .def_readwrite("bulk_viscosity_factor", &PhysicalParameters::bulkViscosityFactor,
-                           "The viscosity of the fluid will be multiplied with this factor to calculate its bulk viscosity. Default is 1.0")
-            .def_readwrite("lattice_viscosity", &PhysicalParameters::latticeViscosity, "Lattice viscosity");
-
-    py::class_<GridParameters, std::shared_ptr<GridParameters>>(parametersModule, "GridParameters")
-            .def(py::init())
-            .def_readwrite("node_distance", &GridParameters::nodeDistance)
-            .def_readwrite("reference_direction_index", &GridParameters::referenceDirectionIndex)
-            .def_readwrite("number_of_nodes_per_direction", &GridParameters::numberOfNodesPerDirection)
-            .def_readwrite("blocks_per_direction", &GridParameters::blocksPerDirection)
-            .def_readwrite("periodic_boundary_in_x1", &GridParameters::periodicBoundaryInX1)
-            .def_readwrite("periodic_boundary_in_x2", &GridParameters::periodicBoundaryInX2)
-            .def_readwrite("periodic_boundary_in_x3", &GridParameters::periodicBoundaryInX3)
-            .def_property_readonly("bounding_box", &GridParameters::boundingBox);
-
-    py::class_<BoundingBox, std::shared_ptr<BoundingBox>>(parametersModule, "BoundingBox")
-            .def_readonly("min_x1", &BoundingBox::minX1)
-            .def_readonly("min_x2", &BoundingBox::minX2)
-            .def_readonly("min_x3", &BoundingBox::minX3)
-            .def_readonly("max_x1", &BoundingBox::maxX1)
-            .def_readonly("max_x2", &BoundingBox::maxX2)
-            .def_readonly("max_x3", &BoundingBox::maxX3)
-            .def("__repr__", [](BoundingBox &self)
-            {
-                std::ostringstream stream;
-                stream << "<BoundingBox" << std::endl
-                       << "min x1: " << self.minX1 << std::endl
-                       << "min x2: " << self.minX2 << std::endl
-                       << "min x3: " << self.minX3 << std::endl
-                       << "max x1: " << self.maxX1 << std::endl
-                       << "max x2: " << self.maxX2 << std::endl
-                       << "max x3: " << self.maxX3 << std::endl << ">";
-
-                return stream.str();
-            });
-
-    py::class_<RuntimeParameters, std::shared_ptr<RuntimeParameters>>(parametersModule, "RuntimeParameters")
-            .def(py::init())
-            .def_readwrite("number_of_timesteps", &RuntimeParameters::numberOfTimeSteps)
-            .def_readwrite("timestep_log_interval", &RuntimeParameters::timeStepLogInterval)
-            .def_readwrite("number_of_threads", &RuntimeParameters::numberOfThreads);
-
+    namespace py = pybind11;
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module parametersModule = parentModule.def_submodule("parameters");
+
+        py::class_<PhysicalParameters, std::shared_ptr<PhysicalParameters>>(parametersModule, "PhysicalParameters")
+                .def(py::init())
+                .def_readwrite("bulk_viscosity_factor", &PhysicalParameters::bulkViscosityFactor,
+                               "The viscosity of the fluid will be multiplied with this factor to calculate its bulk viscosity. Default is 1.0")
+                .def_readwrite("lattice_viscosity", &PhysicalParameters::latticeViscosity, "Lattice viscosity");
+
+        py::class_<GridParameters, std::shared_ptr<GridParameters>>(parametersModule, "GridParameters")
+                .def(py::init())
+                .def_readwrite("node_distance", &GridParameters::nodeDistance)
+                .def_readwrite("reference_direction_index", &GridParameters::referenceDirectionIndex)
+                .def_readwrite("number_of_nodes_per_direction", &GridParameters::numberOfNodesPerDirection)
+                .def_readwrite("blocks_per_direction", &GridParameters::blocksPerDirection)
+                .def_readwrite("periodic_boundary_in_x1", &GridParameters::periodicBoundaryInX1)
+                .def_readwrite("periodic_boundary_in_x2", &GridParameters::periodicBoundaryInX2)
+                .def_readwrite("periodic_boundary_in_x3", &GridParameters::periodicBoundaryInX3)
+                .def_property_readonly("bounding_box", &GridParameters::boundingBox);
+
+        py::class_<BoundingBox, std::shared_ptr<BoundingBox>>(parametersModule, "BoundingBox")
+                .def_readonly("min_x1", &BoundingBox::minX1)
+                .def_readonly("min_x2", &BoundingBox::minX2)
+                .def_readonly("min_x3", &BoundingBox::minX3)
+                .def_readonly("max_x1", &BoundingBox::maxX1)
+                .def_readonly("max_x2", &BoundingBox::maxX2)
+                .def_readonly("max_x3", &BoundingBox::maxX3)
+                .def("__repr__", [](BoundingBox &self)
+                {
+                    std::ostringstream stream;
+                    stream << "<BoundingBox" << std::endl
+                           << "min x1: " << self.minX1 << std::endl
+                           << "min x2: " << self.minX2 << std::endl
+                           << "min x3: " << self.minX3 << std::endl
+                           << "max x1: " << self.maxX1 << std::endl
+                           << "max x2: " << self.maxX2 << std::endl
+                           << "max x3: " << self.maxX3 << std::endl << ">";
+
+                    return stream.str();
+                });
+
+        py::class_<RuntimeParameters, std::shared_ptr<RuntimeParameters>>(parametersModule, "RuntimeParameters")
+                .def(py::init())
+                .def_readwrite("number_of_timesteps", &RuntimeParameters::numberOfTimeSteps)
+                .def_readwrite("timestep_log_interval", &RuntimeParameters::timeStepLogInterval)
+                .def_readwrite("number_of_threads", &RuntimeParameters::numberOfThreads);
+
+    }
 }
\ No newline at end of file
diff --git a/src/cpu/pythonbindings/src/submodules/writer.cpp b/src/cpu/pythonbindings/src/submodules/writer.cpp
index 40819e4766eb30a442967067f954fef5508a4707..d5ec527a27caf63d9a3066c51e1f675b307fe0b2 100644
--- a/src/cpu/pythonbindings/src/submodules/writer.cpp
+++ b/src/cpu/pythonbindings/src/submodules/writer.cpp
@@ -1,18 +1,21 @@
 #include <pybind11/pybind11.h>
 #include <simulationconfig/WriterConfiguration.h>
 
-namespace py = pybind11;
-
-void makeWriterModule(py::module_ &parentModule)
+namespace writer
 {
-    py::module writerModule = parentModule.def_submodule("writer");
+    namespace py = pybind11;
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module writerModule = parentModule.def_submodule("writer");
 
-    py::enum_<OutputFormat>(writerModule, "OutputFormat")
-            .value("ASCII", OutputFormat::ASCII)
-            .value("BINARY", OutputFormat::BINARY);
+        py::enum_<OutputFormat>(writerModule, "OutputFormat")
+                .value("ASCII", OutputFormat::ASCII)
+                .value("BINARY", OutputFormat::BINARY);
 
-    py::class_<WriterConfiguration>(writerModule, "Writer")
-            .def(py::init())
-            .def_readwrite("output_path", &WriterConfiguration::outputPath)
-            .def_readwrite("output_format", &WriterConfiguration::outputFormat);
+        py::class_<WriterConfiguration>(writerModule, "Writer")
+                .def(py::init())
+                .def_readwrite("output_path", &WriterConfiguration::outputPath)
+                .def_readwrite("output_format", &WriterConfiguration::outputFormat);
+    }
 }
\ No newline at end of file
diff --git a/src/cpu/simulationconfig/include/simulationconfig/Simulation.h b/src/cpu/simulationconfig/include/simulationconfig/Simulation.h
index 47f8b0e6ef1c844a70efcf0faedd5cbcdb7dbc05..4bf800c2375347ab7040424bdb2e4c53d5cc2bf8 100644
--- a/src/cpu/simulationconfig/include/simulationconfig/Simulation.h
+++ b/src/cpu/simulationconfig/include/simulationconfig/Simulation.h
@@ -63,7 +63,9 @@ public:
     void run();
 
 private:
-    std::shared_ptr<GbObject3D> makeSimulationBoundingBox() const;
+    bool isMainProcess();
+
+    std::shared_ptr<GbObject3D> makeSimulationBoundingBox();
 
     void writeBlocksToFile() const;
 
@@ -84,6 +86,12 @@ private:
 
 
     void setKernelForcing(const std::shared_ptr<LBMKernel> &kernel, std::shared_ptr<LBMUnitConverter> &converter) const;
+
+    void setConnectors();
+
+    void initializeDistributions();
+
+    std::shared_ptr<CoProcessor> makeNupsCoProcessor() const;
 };
 
 #endif
\ No newline at end of file
diff --git a/src/cpu/simulationconfig/src/Simulation.cpp b/src/cpu/simulationconfig/src/Simulation.cpp
index 6d80a1cfa2cc62d2575d1e8f67193b879f217b90..2ac83701cbd17d0c790022552ace5c7c804a7acd 100644
--- a/src/cpu/simulationconfig/src/Simulation.cpp
+++ b/src/cpu/simulationconfig/src/Simulation.cpp
@@ -65,6 +65,8 @@ Simulation::addObject(const std::shared_ptr<GbObject3D> &object, const std::shar
     const bool is_in = registeredAdapters.find(bcAdapter) != registeredAdapters.end();
     if (!is_in) addBCAdapter(bcAdapter);
     this->interactors.push_back(lbmSystem->makeInteractor(object, this->grid, bcAdapter, state));
+    if (communicator->getProcessID() != 0) return;
+
     GbSystem3D::writeGeoObject(object, writerConfig.outputPath + folderPath, writerConfig.getWriter());
 }
 
@@ -105,13 +107,14 @@ void Simulation::run()
     int &nodesInX1 = gridParameters->numberOfNodesPerDirection[0];
     int &nodesInX2 = gridParameters->numberOfNodesPerDirection[1];
     int &nodesInX3 = gridParameters->numberOfNodesPerDirection[2];
-    logSimulationData(nodesInX1, nodesInX2, nodesInX3);
+
+    if (isMainProcess())
+        logSimulationData(nodesInX1, nodesInX2, nodesInX3);
 
     setBlockSize(nodesInX1, nodesInX2, nodesInX3);
     auto gridCube = makeSimulationBoundingBox();
 
     generateBlockGrid(gridCube);
-
     setKernelForcing(lbmKernel, converter);
     setBoundaryConditionProcessor(lbmKernel);
 
@@ -125,55 +128,44 @@ void Simulation::run()
 
     intHelper.selectBlocks();
 
-
     int numberOfProcesses = communicator->getNumberOfProcesses();
     SetKernelBlockVisitor kernelVisitor(lbmKernel, physicalParameters->latticeViscosity,
                                         numberOfProcesses);
     grid->accept(kernelVisitor);
     intHelper.setBC();
 
-    double bulkViscosity = physicalParameters->latticeViscosity * physicalParameters->bulkViscosityFactor;
-    //auto iProcessor = std::make_shared<CompressibleOffsetMomentsInterpolationProcessor>();
-    //iProcessor->setBulkViscosity(physicalParameters->latticeViscosity, bulkViscosity);
-
-    //SetConnectorsBlockVisitor setConnsVisitor(communicator, true,
-    //                                          lbmSystem->getNumberOfDirections(),
-    //                                          physicalParameters->latticeViscosity, iProcessor);
 
-    OneDistributionSetConnectorsBlockVisitor setConnsVisitor(communicator);
-    grid->accept(setConnsVisitor);
-
-    InitDistributionsBlockVisitor initVisitor;
-    grid->accept(initVisitor);
-    grid->accept(setConnsVisitor);
+    writeBlocksToFile(); // important: run this after metis & intHelper.selectBlocks()
+    setConnectors();
+    initializeDistributions();
     grid->accept(bcVisitor);
-
     writeBoundaryConditions();
-    // important: run this after metis & intHelper.selectBlocks()
-    writeBlocksToFile();
+
+#ifdef _OPENMP
+    omp_set_num_threads(simulationParameters->numberOfThreads);
+    if (isMainProcess())
+        UBLOG(logINFO, "OpenMP is set to run with " << simulationParameters->numberOfThreads << " threads")
+#endif
 
     auto visualizationScheduler = std::make_shared<UbScheduler>(simulationParameters->timeStepLogInterval);
     auto mqCoProcessor = makeMacroscopicQuantitiesCoProcessor(converter,
                                                               visualizationScheduler);
 
-    std::shared_ptr<UbScheduler> nupsScheduler(new UbScheduler(100, 100));
-    std::shared_ptr<CoProcessor> nupsCoProcessor(
-            new NUPSCounterCoProcessor(grid, nupsScheduler, simulationParameters->numberOfThreads, communicator));
-
-
-#ifdef _OPENMP
-    omp_set_num_threads(simulationParameters->numberOfThreads);
-    UBLOG(logINFO, "OpenMP is set to run with " << omp_get_num_threads() << " threads")
-#endif
+    auto nupsCoProcessor = makeNupsCoProcessor();
 
     auto calculator = std::make_shared<BasicCalculator>(grid, visualizationScheduler,
                                                         simulationParameters->numberOfTimeSteps);
     calculator->addCoProcessor(nupsCoProcessor);
     calculator->addCoProcessor(mqCoProcessor);
 
-    UBLOG(logINFO, "Simulation-start")
+    if (isMainProcess()) UBLOG(logINFO, "Simulation-start")
     calculator->calculate();
-    UBLOG(logINFO, "Simulation-end")
+    if (isMainProcess()) UBLOG(logINFO, "Simulation-end")
+}
+
+bool Simulation::isMainProcess()
+{
+    return communicator->getProcessID() == 0;
 }
 
 void
@@ -224,18 +216,7 @@ Simulation::makeLBMUnitConverter()
     return std::make_shared<LBMUnitConverter>();
 }
 
-std::shared_ptr<CoProcessor>
-Simulation::makeMacroscopicQuantitiesCoProcessor(const std::shared_ptr<LBMUnitConverter> &converter,
-                                                 const std::shared_ptr<UbScheduler> &visualizationScheduler) const
-{
-    auto mqCoProcessor = std::make_shared<WriteMacroscopicQuantitiesCoProcessor>(grid, visualizationScheduler,
-                                                                                 writerConfig.outputPath,
-                                                                                 writerConfig.getWriter(),
-                                                                                 converter,
-                                                                                 communicator);
-    mqCoProcessor->process(0);
-    return mqCoProcessor;
-}
+
 
 void Simulation::writeBoundaryConditions() const
 {
@@ -248,8 +229,9 @@ void Simulation::writeBoundaryConditions() const
 void Simulation::writeBlocksToFile() const
 {
     UBLOG(logINFO, "Write block grid to VTK-file")
+    auto scheduler = std::make_shared<UbScheduler>(1);
     auto ppblocks = std::make_shared<WriteBlocksCoProcessor>(grid,
-                                                             std::make_shared<UbScheduler>(1),
+                                                             scheduler,
                                                              writerConfig.outputPath,
                                                              writerConfig.getWriter(),
                                                              communicator);
@@ -258,17 +240,56 @@ void Simulation::writeBlocksToFile() const
 }
 
 std::shared_ptr<GbObject3D>
-Simulation::makeSimulationBoundingBox() const
+Simulation::makeSimulationBoundingBox()
 {
     auto box = gridParameters->boundingBox();
+    auto gridCube = std::make_shared<GbCuboid3D>(box->minX1, box->minX2, box->minX3, box->maxX1, box->maxX2,
+                                                 box->maxX3);
 
-    UBLOG(logINFO, "Bounding box dimensions = [("
-            << box->minX1 << ", " << box->minX2 << ", " << box->minX3 << "); ("
-            << box->maxX1 << ", " << box->maxX2 << ", " << box->maxX3 << ")]")
+    if (isMainProcess())
+    {
+        UBLOG(logINFO, "Bounding box dimensions = [("
+                << box->minX1 << ", " << box->minX2 << ", " << box->minX3 << "); ("
+                << box->maxX1 << ", " << box->maxX2 << ", " << box->maxX3 << ")]")
+
+        GbSystem3D::writeGeoObject(gridCube.get(), writerConfig.outputPath + "/geo/gridCube", writerConfig.getWriter());
+    }
 
-    auto gridCube = std::make_shared<GbCuboid3D>(box->minX1, box->minX2, box->minX3, box->maxX1, box->maxX2, box->maxX3);
-    GbSystem3D::writeGeoObject(gridCube.get(), writerConfig.outputPath + "/geo/gridCube", writerConfig.getWriter());
     return gridCube;
 }
 
+void Simulation::setConnectors()
+{
+    OneDistributionSetConnectorsBlockVisitor setConnsVisitor(communicator);
+    grid->accept(setConnsVisitor);
+}
+
+void Simulation::initializeDistributions()
+{
+    InitDistributionsBlockVisitor initVisitor;
+    grid->accept(initVisitor);
+}
+
+std::shared_ptr<CoProcessor>
+Simulation::makeMacroscopicQuantitiesCoProcessor(const std::shared_ptr<LBMUnitConverter> &converter,
+                                                 const std::shared_ptr<UbScheduler> &visualizationScheduler) const
+{
+    auto mqCoProcessor = std::make_shared<WriteMacroscopicQuantitiesCoProcessor>(grid, visualizationScheduler,
+                                                                                 writerConfig.outputPath,
+                                                                                 writerConfig.getWriter(),
+                                                                                 converter,
+                                                                                 communicator);
+    mqCoProcessor->process(0);
+    return mqCoProcessor;
+}
+
+std::shared_ptr<CoProcessor> Simulation::makeNupsCoProcessor() const
+{
+    auto scheduler = std::make_shared<UbScheduler>(100, 100);
+    return std::make_shared<NUPSCounterCoProcessor>(grid, scheduler,
+                                                    simulationParameters->numberOfThreads,
+                                                    communicator);
+}
+
+
 Simulation::~Simulation() = default;