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..7b2f13b835f0f5cd37bae2a2fafc6af206e504ca
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/settings.py
@@ -0,0 +1,26 @@
+import os
+from acousticscaling import AcousticScaling
+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, 8]
+grid_params.periodic_boundary_in_x1 = True
+grid_params.periodic_boundary_in_x2 = True
+
+physical_params = PhysicalParameters()
+physical_params.lattice_viscosity = 1e-3
+
+runtime_params = RuntimeParameters()
+runtime_params.number_of_threads = int(os.environ["SLURM_NTASKS"])
+runtime_params.number_of_timesteps = 2_500_000
+runtime_params.timestep_log_interval = 500_000
+
+kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+kernel.use_forcing = True
+kernel.forcing_in_x1 = 1e-9
+
+Scaling = AcousticScaling(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..1d7072123060f94981cb5d975fbbb5575bd71fda
--- /dev/null
+++ b/Python/SlurmTests/poiseuille/slurm.job
@@ -0,0 +1,25 @@
+#!/bin/bash
+#SBATCH -J PyFluidsTest
+#SBATCH --nodes=2
+#SBATCH --ntasks-per-node=1
+#SBATCH --cpus-per-task=8
+
+#SBATCH --mem-per-cpu=30000
+#SBATCH --time=02:00:00
+#SBATCH --partition=shortrun_small
+
+source .bashrc
+
+echo "PyFluids Poiseuille Test Case"
+echo "Number of tasks: ${SLURM_NTASKS}"
+
+export PYFLUIDS_SCALE_LEVEL=0
+mpiexec -np "$SLURM_NTASKS" singularity exec SingularityImages/VirtualFluidsMPICH.sif python3 /VirtualFluids_dev/Python/poiseuille/simulation_runner.py
+
+export PYFLUIDS_SCALE_LEVEL=1
+mpiexec -np "$SLURM_NTASKS" singularity exec SingularityImages/VirtualFluidsMPICH.sif python3 /VirtualFluids_dev/Python/poiseuille/simulation_runner.py
+
+export PYFLUIDS_SCALE_LEVEL=2
+mpiexec -np "$SLURM_NTASKS" singularity exec SingularityImages/VirtualFluidsMPICH.sif python3 /VirtualFluids_dev/Python/poiseuille/simulation_runner.py
+
+singularity exec SingularityImages/VirtualFluidsMPICH.sif python3 /VirtualFluids_dev/Python/poiseuille/SlurmTests/poiseuille/evaluator.py
diff --git a/Python/acousticscaling.py b/Python/acousticscaling.py
new file mode 100644
index 0000000000000000000000000000000000000000..ff10cbd49d344f87e6a7b32a991bf0b2d2c5e41a
--- /dev/null
+++ b/Python/acousticscaling.py
@@ -0,0 +1,85 @@
+from pyfluids.kernel import LBMKernel
+from pyfluids.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+
+
+class AcousticScaling:
+
+    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/poiseuille/simulation.py b/Python/poiseuille/simulation.py
index 898809fbb3e39e2303d47fc8cca4e3fe163560cd..31ceb1ab9ef90fa4fd606bde4f47c45b8f7d7567 100644
--- a/Python/poiseuille/simulation.py
+++ b/Python/poiseuille/simulation.py
@@ -20,11 +20,14 @@ default_runtime_params.number_of_threads = 4
 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
@@ -34,13 +37,10 @@ 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,
-                   kernel=default_kernel):
+                   kernel=default_kernel,
+                   writer=default_writer):
     simulation = Simulation()
 
-    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)
diff --git a/Python/test_boundaryconditions.py b/Python/test_boundaryconditions.py
index 4fd51680eefe565bef79657a3f020229ae906407..5a7d61f36337398fc5621540951f15b72262b17b 100644
--- a/Python/test_boundaryconditions.py
+++ b/Python/test_boundaryconditions.py
@@ -58,4 +58,4 @@ class BoundaryConditionsTest(unittest.TestCase):
         Should be able to create NonReflectingOutflow
         """
 
-        sut = NonReflectingOutflow()
\ No newline at end of file
+        sut = NonReflectingOutflow()
diff --git a/Python/test_kernel.py b/Python/test_kernel.py
index 5e5487903a279501b0acdba1ab192ac8dd381383..b1016f15308a77c9025787e061355819cbca3874 100644
--- a/Python/test_kernel.py
+++ b/Python/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 55b83aed1cf7dcd8bd430aa17546b26a8f9155e5..14d13a40d85e669a463c0b5b31ca2a21e771af75 100644
--- a/Python/vtk_utilities.py
+++ b/Python/vtk_utilities.py
@@ -1,4 +1,6 @@
 import math
+from typing import List
+
 import pyvista as pv
 
 
@@ -21,7 +23,7 @@ def vertical_column_from_mesh(mesh):
     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/src/cpu/pythonbindings/src/submodules/kernel.cpp b/src/cpu/pythonbindings/src/submodules/kernel.cpp
index 2418aa1d089bc9e0deff7ae1e79007ca424a7522..fb291790632cc2041410f60a14fca8d966283343 100644
--- a/src/cpu/pythonbindings/src/submodules/kernel.cpp
+++ b/src/cpu/pythonbindings/src/submodules/kernel.cpp
@@ -16,9 +16,9 @@ namespace kernel
                 .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)