diff --git a/apps/gpu/AtmosphericBoundaryLayer/AtmosphericBoundaryLayer.cpp b/apps/gpu/AtmosphericBoundaryLayer/AtmosphericBoundaryLayer.cpp
index 3940c1f51a67138c6eebd87d30e1948e059082ba..81f0b354b0fea4cccfeee4dfe2caca63f4469fc9 100644
--- a/apps/gpu/AtmosphericBoundaryLayer/AtmosphericBoundaryLayer.cpp
+++ b/apps/gpu/AtmosphericBoundaryLayer/AtmosphericBoundaryLayer.cpp
@@ -36,6 +36,7 @@
 #include <basics/DataTypes.h>
 #include <basics/config/ConfigurationFile.h>
 #include <basics/constants/NumericConstants.h>
+#include <basics/geometry3d/Axis.h>
 
 #include <logger/Logger.h>
 
@@ -336,7 +337,7 @@ void run(const vf::basics::ConfigurationFile& config)
     if (!usePrecursorInflow && (isFirstSubDomain || !isMultiGPU)) {
         const auto planarAverageProbe = std::make_shared<PlanarAverageProbe>(
             para, cudaMemoryManager, para->getOutputPath(), "planarAverageProbe", timeStepStartAveraging,
-            timeStepStartTemporalAveraging, timeStepAveraging, timeStepStartOutProbe, timeStepOutProbe, PlanarAverageProbe::PlaneNormal::z, true);
+            timeStepStartTemporalAveraging, timeStepAveraging, timeStepStartOutProbe, timeStepOutProbe, Axis::z, true, false);
         planarAverageProbe->addAllAvailableStatistics();
         planarAverageProbe->setFileNameToNOut();
         para->addSampler(planarAverageProbe);
diff --git a/pythonbindings/pyfluids-stubs/__init__.pyi b/pythonbindings/pyfluids-stubs/__init__.pyi
index d2ee6ab7c90fa6401897669f5b68ddd062d2c5c8..72d78df769a9475c6104aad91283b3ac7f9f9499 100644
--- a/pythonbindings/pyfluids-stubs/__init__.pyi
+++ b/pythonbindings/pyfluids-stubs/__init__.pyi
@@ -30,9 +30,9 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 from . import basics as basics
 from . import logger as logger
 from . import lbm as lbm
 from . import gpu as gpu
 from . import cpu as cpu
-
diff --git a/pythonbindings/pyfluids-stubs/basics/__init__.pyi b/pythonbindings/pyfluids-stubs/basics/__init__.pyi
index 30f26e31893006454655cf1cf286d48474aa18ed..23ef2ea2c07a9c668f660f1357644bf80329eac9 100644
--- a/pythonbindings/pyfluids-stubs/basics/__init__.pyi
+++ b/pythonbindings/pyfluids-stubs/basics/__init__.pyi
@@ -30,10 +30,13 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 from __future__ import annotations
 
 from typing import overload
 
+from . import geometry3d as geometry3d
+
 class ConfigurationFile:
     def __init__(self) -> None: ...
     def contains(self, key: str) -> bool: ...
diff --git a/pythonbindings/pyfluids-stubs/basics/geometry3d.pyi b/pythonbindings/pyfluids-stubs/basics/geometry3d.pyi
new file mode 100644
index 0000000000000000000000000000000000000000..051a7f3b0805974538e0e1f23d46b78363de20ae
--- /dev/null
+++ b/pythonbindings/pyfluids-stubs/basics/geometry3d.pyi
@@ -0,0 +1,40 @@
+r"""
+=======================================================================================
+ ____          ____    __    ______     __________   __      __       __        __
+ \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+      \    \  |    |   ________________________________________________________________
+       \    \ |    |  |  ______________________________________________________________|
+        \    \|    |  |  |         __          __     __     __     ______      _______
+         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+
+  This file is part of VirtualFluids. VirtualFluids is free software: you can
+  redistribute it and/or modify it under the terms of the GNU General Public
+  License as published by the Free Software Foundation, either version 3 of
+  the License, or (at your option) any later version.
+
+  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+  for more details.
+
+  SPDX-License-Identifier: GPL-3.0-or-later
+  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+
+! \author Henry Korb
+=======================================================================================
+"""
+
+from __future__ import annotations
+from enum import Enum
+
+class Axis(Enum):
+    x: int
+    y: int
+    z: int
diff --git a/pythonbindings/pyfluids-stubs/gpu/__init__.pyi b/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
index cfc38142b9ac6be063d0370a948f94cbeb9ad9b8..e0f643e392dc5c40b65b453da49a063c52e07a00 100644
--- a/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/__init__.pyi
@@ -30,13 +30,14 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 from __future__ import annotations
 from typing import Callable, ClassVar, List, Optional
 
 from typing import overload, Union
 import numpy as np
 import numpy.typing as npt
-import basics, parallel
+from .. import basics, parallel
 
 from . import grid_generator as grid_generator
 from . import probes as probes
@@ -48,12 +49,26 @@ class PreCollisionInteractor:
 class Sampler:
     def __init__(self, *args, **kwargs) -> None: ...
 
-
 class FileCollection:
     def __init__(self, *args, **kwargs) -> None: ...
 
 class ActuatorFarm(PreCollisionInteractor):
-    def __init__(self, para: Parameter, cuda_memory_manager: CudaMemoryManager,  diameter: float, blade_radii: npt.NDArray[np.float32],turbine_positions_x: npt.NDArray[np.float32], turbine_positions_y: npt.NDArray[np.float32], turbine_positions_z: npt.NDArray[np.float32], density: float, smearing_width: float, level: int, delta_t: float, delta_x: float, use_host_arrays: bool) -> None: ...
+    def __init__(
+        self,
+        para: Parameter,
+        cuda_memory_manager: CudaMemoryManager,
+        diameter: float,
+        blade_radii: npt.NDArray[np.float32],
+        turbine_positions_x: npt.NDArray[np.float32],
+        turbine_positions_y: npt.NDArray[np.float32],
+        turbine_positions_z: npt.NDArray[np.float32],
+        density: float,
+        smearing_width: float,
+        level: int,
+        delta_t: float,
+        delta_x: float,
+        use_host_arrays: bool,
+    ) -> None: ...
     def update_forces_and_coordinates(self) -> None: ...
     def get_all_blade_coords_x(self) -> npt.NDArray[np.float32]: ...
     def get_all_blade_coords_x_device(self) -> int: ...
@@ -94,12 +109,45 @@ class ActuatorFarm(PreCollisionInteractor):
     def get_turbine_blade_velocities_y_device(self, turbine: int) -> int: ...
     def get_turbine_blade_velocities_z(self, turbine: int) -> npt.NDArray[np.float32]: ...
     def get_turbine_blade_velocities_z_device(self, turbine: int) -> int: ...
-    def set_all_blade_coords(self, blade_coords_x: npt.NDArray[np.float32], blade_coords_y: npt.NDArray[np.float32], blade_coords_z: npt.NDArray[np.float32]) -> None: ...
-    def set_all_blade_forces(self, blade_forces_x: npt.NDArray[np.float32], blade_forces_y: npt.NDArray[np.float32], blade_forces_z: npt.NDArray[np.float32]) -> None: ...
-    def set_all_blade_velocities(self, blade_velocities_x: npt.NDArray[np.float32], blade_velocities_y: npt.NDArray[np.float32], blade_velocities_z: npt.NDArray[np.float32]) -> None: ...
-    def set_turbine_blade_coords(self, turbine: int, blade_coords_x: npt.NDArray[np.float32], blade_coords_y: npt.NDArray[np.float32], blade_coords_z: npt.NDArray[np.float32]) -> None: ...
-    def set_turbine_blade_forces(self, turbine: int, blade_forces_x: npt.NDArray[np.float32], blade_forces_y: npt.NDArray[np.float32], blade_forces_z: npt.NDArray[np.float32]) -> None: ...
-    def set_turbine_blade_velocities(self, turbine: int, blade_velocities_x: npt.NDArray[np.float32], blade_velocities_y: npt.NDArray[np.float32], blade_velocities_z: npt.NDArray[np.float32]) -> None: ...
+    def set_all_blade_coords(
+        self,
+        blade_coords_x: npt.NDArray[np.float32],
+        blade_coords_y: npt.NDArray[np.float32],
+        blade_coords_z: npt.NDArray[np.float32],
+    ) -> None: ...
+    def set_all_blade_forces(
+        self,
+        blade_forces_x: npt.NDArray[np.float32],
+        blade_forces_y: npt.NDArray[np.float32],
+        blade_forces_z: npt.NDArray[np.float32],
+    ) -> None: ...
+    def set_all_blade_velocities(
+        self,
+        blade_velocities_x: npt.NDArray[np.float32],
+        blade_velocities_y: npt.NDArray[np.float32],
+        blade_velocities_z: npt.NDArray[np.float32],
+    ) -> None: ...
+    def set_turbine_blade_coords(
+        self,
+        turbine: int,
+        blade_coords_x: npt.NDArray[np.float32],
+        blade_coords_y: npt.NDArray[np.float32],
+        blade_coords_z: npt.NDArray[np.float32],
+    ) -> None: ...
+    def set_turbine_blade_forces(
+        self,
+        turbine: int,
+        blade_forces_x: npt.NDArray[np.float32],
+        blade_forces_y: npt.NDArray[np.float32],
+        blade_forces_z: npt.NDArray[np.float32],
+    ) -> None: ...
+    def set_turbine_blade_velocities(
+        self,
+        turbine: int,
+        blade_velocities_x: npt.NDArray[np.float32],
+        blade_velocities_y: npt.NDArray[np.float32],
+        blade_velocities_z: npt.NDArray[np.float32],
+    ) -> None: ...
     def set_turbine_azimuth(self, turbine: int, azimuth: float) -> None: ...
     def enable_output(self, output_name: str, t_start_out: int, t_out: int) -> None: ...
     @property
@@ -120,7 +168,22 @@ class ActuatorFarm(PreCollisionInteractor):
     def number_of_turbines(self) -> int: ...
 
 class ActuatorFarmStandalone(ActuatorFarm):
-    def __init__(self, para: Parameter, cuda_memory_manager: CudaMemoryManager, diameter: float, number_of_nodes_per_blade: int, turbine_positions_x: npt.NDArray[np.float32], turbine_positions_y: npt.NDArray[np.float32], turbine_positions_z: npt.NDArray[np.float32], rotor_speeds: npt.NDArray[np.float32], density: float, smearing_width: float, level: int, delta_t: float, delta_x: float) -> None: ...
+    def __init__(
+        self,
+        para: Parameter,
+        cuda_memory_manager: CudaMemoryManager,
+        diameter: float,
+        number_of_nodes_per_blade: int,
+        turbine_positions_x: npt.NDArray[np.float32],
+        turbine_positions_y: npt.NDArray[np.float32],
+        turbine_positions_z: npt.NDArray[np.float32],
+        rotor_speeds: npt.NDArray[np.float32],
+        density: float,
+        smearing_width: float,
+        level: int,
+        delta_t: float,
+        delta_x: float,
+    ) -> None: ...
 
 class BoundaryConditionFactory:
     def __init__(self) -> None: ...
@@ -132,11 +195,9 @@ class BoundaryConditionFactory:
     def set_stress_boundary_condition(self, boundary_condition_type: StressBC) -> None: ...
     def set_velocity_boundary_condition(self, boundary_condition_type: VelocityBC) -> None: ...
 
-
 class CudaMemoryManager:
     def __init__(self, parameter: Parameter) -> None: ...
 
-
 class TransientBCFileType:
     __members__: ClassVar[dict] = ...  # read-only
     VTK: ClassVar[TransientBCFileType] = ...
@@ -152,14 +213,15 @@ class TransientBCFileType:
     @property
     def name(self) -> str: ...
 
-
 class GridProvider:
     def __init__(self, *args, **kwargs) -> None: ...
     @staticmethod
-    def make_grid_generator(builder: grid_generator.GridBuilder, para: Parameter, cuda_memory_manager: CudaMemoryManager, communicator: parallel.Communicator) -> GridProvider: ...
-
-class MultipleGridBuilder:
-    def __init__(self) -> None: ...
+    def make_grid_generator(
+        builder: grid_generator.GridBuilder,
+        para: Parameter,
+        cuda_memory_manager: CudaMemoryManager,
+        communicator: parallel.Communicator,
+    ) -> GridProvider: ...
 
 class GridScaling:
     __members__: ClassVar[dict] = ...  # read-only
@@ -178,12 +240,10 @@ class GridScaling:
     @property
     def name(self) -> str: ...
 
-
 class GridScalingFactory:
     def __init__(self) -> None: ...
     def set_scaling_factory(self, scaling_type) -> None: ...
 
-
 class NoSlipBC:
     __members__: ClassVar[dict] = ...  # read-only
     NoSlipDelayBounceBack: ClassVar[NoSlipBC] = ...
@@ -202,7 +262,6 @@ class NoSlipBC:
     @property
     def name(self) -> str: ...
 
-
 class OutputVariable:
     __members__: ClassVar[dict] = ...  # read-only
     Distributions: ClassVar[OutputVariable] = ...
@@ -219,10 +278,11 @@ class OutputVariable:
     @property
     def name(self) -> str: ...
 
-
 class Parameter:
     @overload
-    def __init__(self, number_of_processes: int, my_ID: int, config_data: Optional[basics.ConfigurationFile]) -> None: ...
+    def __init__(
+        self, number_of_processes: int, my_ID: int, config_data: Optional[basics.ConfigurationFile]
+    ) -> None: ...
     @overload
     def __init__(self, number_of_processes: int, my_ID: int) -> None: ...
     @overload
@@ -241,7 +301,6 @@ class Parameter:
     def get_viscosity_ratio(self) -> float: ...
     def set_AD_kernel(self, ad_kernel: str) -> None: ...
     def set_calc_turbulence_intensity(self, calc_velocity_and_fluctuations: bool) -> None: ...
-    def set_comp_on(self, is_comp: bool) -> None: ...
     def set_density_ratio(self, density_ratio: float) -> None: ...
     def set_devices(self, devices: List[int]) -> None: ...
     def set_diff_on(self, is_diff: bool) -> None: ...
@@ -249,7 +308,9 @@ class Parameter:
     def set_has_wall_model_monitor(self, has_wall_monitor: bool) -> None: ...
     def set_initial_condition(self, init_func: Callable[[float, float, float], List[float]]) -> None: ...
     def set_initial_condition_log_law(self, u_star: float, z0: float, velocity_ratio: float) -> None: ...
-    def set_initial_condition_perturbed_log_law(self, u_star: float, z0: float, length_x: float, length_z: float, height: float, velocity_ratio: float) -> None: ...
+    def set_initial_condition_perturbed_log_law(
+        self, u_star: float, z0: float, length_x: float, length_z: float, height: float, velocity_ratio: float
+    ) -> None: ...
     def set_initial_condition_uniform(self, velocity_x: float, velocity_y: float, velocity_z: float) -> None: ...
     def set_is_body_force(self, is_body_force: bool) -> None: ...
     def configure_main_kernel(self, kernel: str) -> None: ...
@@ -259,9 +320,9 @@ class Parameter:
     def set_output_path(self, o_path: str) -> None: ...
     def set_output_prefix(self, o_prefix: str) -> None: ...
     def set_print_files(self, print_files: bool) -> None: ...
-    def set_quadric_limiters(self, quadric_limiter_p: float, quadric_limiter_m: float, quadric_limiter_d: float) -> None: ...
-    def set_temperature_BC(self, temp_bc: float) -> None: ...
-    def set_temperature_init(self, temp: float) -> None: ...
+    def set_quadric_limiters(
+        self, quadric_limiter_p: float, quadric_limiter_m: float, quadric_limiter_d: float
+    ) -> None: ...
     def set_timestep_end(self, tend: int) -> None: ...
     def set_timestep_of_coarse_level(self, timestep: int) -> None: ...
     def set_timestep_out(self, tout: int) -> None: ...
@@ -272,7 +333,6 @@ class Parameter:
     def set_viscosity_LB(self, viscosity: float) -> None: ...
     def set_viscosity_ratio(self, viscosity_ratio: float) -> None: ...
 
-
 class PrecursorBC:
     __members__: ClassVar[dict] = ...  # read-only
     DistributionsPrecursor: ClassVar[PrecursorBC] = ...
@@ -290,10 +350,23 @@ class PrecursorBC:
     @property
     def name(self) -> str: ...
 
-
 class PrecursorWriter(Sampler):
-    def __init__(self, para: Parameter, cuda_memory_manager: CudaMemoryManager, output_path: str, filename: str, x_pos: float, y_min: float, y_max: float, z_min: float, z_max: float, t_start_out: int, t_save: int, output_variable: OutputVariable, max_timesteps_per_file: int) -> None: ...
-
+    def __init__(
+        self,
+        para: Parameter,
+        cuda_memory_manager: CudaMemoryManager,
+        output_path: str,
+        filename: str,
+        x_pos: float,
+        y_min: float,
+        y_max: float,
+        z_min: float,
+        z_max: float,
+        t_start_out: int,
+        t_save: int,
+        output_variable: OutputVariable,
+        max_timesteps_per_file: int,
+    ) -> None: ...
 
 class PressureBC:
     __members__: ClassVar[dict] = ...  # read-only
@@ -314,7 +387,6 @@ class PressureBC:
     @property
     def name(self) -> str: ...
 
-
 class SideType:
     __members__: ClassVar[dict] = ...  # read-only
     GEOMETRY: ClassVar[SideType] = ...
@@ -336,19 +408,41 @@ class SideType:
     @property
     def name(self) -> str: ...
 
-
 class Simulation:
     @overload
-    def __init__(self, parameter: Parameter, memoryManager: CudaMemoryManager, communicator, gridProvider: GridProvider, bcFactory: BoundaryConditionFactory, gridScalingFactory: GridScalingFactory) -> None: ...
+    def __init__(
+        self,
+        parameter: Parameter,
+        memoryManager: CudaMemoryManager,
+        communicator,
+        gridProvider: GridProvider,
+        bcFactory: BoundaryConditionFactory,
+        gridScalingFactory: GridScalingFactory,
+    ) -> None: ...
     @overload
-    def __init__(self, parameter: Parameter, memoryManager: CudaMemoryManager, communicator, gridProvider: GridProvider, bcFactory: BoundaryConditionFactory) -> None: ...
+    def __init__(
+        self,
+        parameter: Parameter,
+        memoryManager: CudaMemoryManager,
+        communicator,
+        gridProvider: GridProvider,
+        bcFactory: BoundaryConditionFactory,
+    ) -> None: ...
     @overload
-    def __init__(self, parameter: Parameter, memoryManager: CudaMemoryManager, communicator, gridProvider: GridProvider, bcFactory: BoundaryConditionFactory, tmFactory: TurbulenceModelFactory, gridScalingFactory: GridScalingFactory) -> None: ...
+    def __init__(
+        self,
+        parameter: Parameter,
+        memoryManager: CudaMemoryManager,
+        communicator,
+        gridProvider: GridProvider,
+        bcFactory: BoundaryConditionFactory,
+        tmFactory: TurbulenceModelFactory,
+        gridScalingFactory: GridScalingFactory,
+    ) -> None: ...
     def addEnstrophyAnalyzer(self, t_analyse: int) -> None: ...
     def addKineticEnergyAnalyzer(self, t_analyse: int) -> None: ...
     def run(self) -> None: ...
 
-
 class SlipBC:
     __members__: ClassVar[dict] = ...  # read-only
     NotSpecified: ClassVar[SlipBC] = ...
@@ -366,7 +460,6 @@ class SlipBC:
     @property
     def name(self) -> str: ...
 
-
 class StressBC:
     __members__: ClassVar[dict] = ...  # read-only
     NotSpecified: ClassVar[StressBC] = ...
@@ -385,7 +478,6 @@ class StressBC:
     @property
     def name(self) -> str: ...
 
-
 class TurbulenceModel:
     __members__: ClassVar[dict] = ...  # read-only
     AMD: ClassVar[TurbulenceModel] = ...
@@ -404,18 +496,15 @@ class TurbulenceModel:
     @property
     def name(self) -> str: ...
 
-
 class TurbulenceModelFactory:
     def __init__(self, para: Parameter) -> None: ...
     def read_config_file(self, config_data: basics.ConfigurationFile) -> None: ...
     def set_model_constant(self, model_constant: float) -> None: ...
     def set_turbulence_model(self, turbulence_model: TurbulenceModel) -> None: ...
 
-
 class VTKFileCollection(FileCollection):
     def __init__(self, prefix: str) -> None: ...
 
-
 class VelocityBC:
     __members__: ClassVar[dict] = ...  # read-only
     NotSpecified: ClassVar[VelocityBC] = ...
@@ -435,5 +524,4 @@ class VelocityBC:
     @property
     def name(self) -> str: ...
 
-
 def create_file_collection(prefix: str, type: TransientBCFileType) -> FileCollection: ...
diff --git a/pythonbindings/pyfluids-stubs/gpu/grid_generator.pyi b/pythonbindings/pyfluids-stubs/gpu/grid_generator.pyi
index 413fec385f706df58e948980976093356f51997b..8edb9c9e8fd47e38dfde7f7ce7343f1273f5027c 100644
--- a/pythonbindings/pyfluids-stubs/gpu/grid_generator.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/grid_generator.pyi
@@ -30,22 +30,20 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 from __future__ import annotations
 
 from typing import List
 
 from typing import overload
-import gpu
-
+from .. import gpu
 
 class BoundingBox:
     def __init__(self, min_x: float, max_x: float, min_y: float, max_y: float, min_z: float, max_z: float) -> None: ...
 
-
 class Object:
     def __init__(self, *args, **kwargs) -> None: ...
 
-
 class Conglomerate(Object):
     def __init__(self, *args, **kwargs) -> None: ...
     def add(self, object: Object) -> None: ...
@@ -53,32 +51,49 @@ class Conglomerate(Object):
     def make_shared() -> Conglomerate: ...
     def subtract(self, object: Object) -> None: ...
 
-
 class Cuboid(Object):
-    def __init__(self, min_x1: float, min_x2: float, min_x3: float, max_x1: float, max_x2: float, max_x3: float) -> None: ...
-
+    def __init__(
+        self, min_x1: float, min_x2: float, min_x3: float, max_x1: float, max_x2: float, max_x3: float
+    ) -> None: ...
 
 class GridBuilder:
     def __init__(self, *args, **kwargs) -> None: ...
     def get_number_of_grid_levels(self) -> int: ...
 
-
-class GridFactory:
-    def __init__(self, *args, **kwargs) -> None: ...
-    @staticmethod
-    def make() -> GridFactory: ...
-
 class GridDimensions:
-    def __init__(self, min_x: float, max_x: float, min_y: float, max_y: float, min_z: float, max_z: float, delta: float) -> None: ...
+    def __init__(
+        self, min_x: float, max_x: float, min_y: float, max_y: float, min_z: float, max_z: float, delta: float
+    ) -> None: ...
 
 class LevelGridBuilder(GridBuilder):
     def __init__(self, *args, **kwargs) -> None: ...
     def set_no_slip_boundary_condition(self, side_type: gpu.SideType) -> None: ...
     def set_periodic_boundary_condition(self, periodic_x: bool, periodic_y: bool, periodic_z: bool) -> None: ...
-    def set_precursor_boundary_condition(self, side_type: gpu.SideType, file_collection: gpu.FileCollection, n_t_read: int, velocity_x: float = ..., velocity_y: float = ..., velocity_z: float = ..., file_level_to_grid_level_map: List[int] = ...) -> None: ...
+    def set_precursor_boundary_condition(
+        self,
+        side_type: gpu.SideType,
+        file_collection: gpu.FileCollection,
+        n_t_read: int,
+        velocity_x: float = ...,
+        velocity_y: float = ...,
+        velocity_z: float = ...,
+        file_level_to_grid_level_map: List[int] = ...,
+    ) -> None: ...
     def set_pressure_boundary_condition(self, side_type: gpu.SideType, rho: float) -> None: ...
-    def set_slip_boundary_condition(self, side_type: gpu.SideType, normal_x: float, normal_y: float, normal_z: float) -> None: ...
-    def set_stress_boundary_condition(self, side_type: gpu.SideType, normal_x: float, normal_y: float, normal_z: float, sampling_offset: int, z0: float, dx: float, q: float = ...) -> None: ...
+    def set_slip_boundary_condition(
+        self, side_type: gpu.SideType, normal_x: float, normal_y: float, normal_z: float
+    ) -> None: ...
+    def set_stress_boundary_condition(
+        self,
+        side_type: gpu.SideType,
+        normal_x: float,
+        normal_y: float,
+        normal_z: float,
+        sampling_offset: int,
+        z0: float,
+        dx: float,
+        q: float = ...,
+    ) -> None: ...
     def set_velocity_boundary_condition(self, side_type: gpu.SideType, vx: float, vy: float, vz: float) -> None: ...
     def set_periodic_shift_on_x_boundary_in_y_direction(self, shift: float) -> None: ...
     def set_periodic_shift_on_x_boundary_in_z_direction(self, shift: float) -> None: ...
@@ -88,14 +103,12 @@ class LevelGridBuilder(GridBuilder):
     def set_periodic_shift_on_z_boundary_in_y_direction(self, shift: float) -> None: ...
     def set_communication_process(self, direction: int, process: int) -> None: ...
 
-
-
 class MultipleGridBuilder(LevelGridBuilder):
     def __init__(self, *args, **kwargs) -> None: ...
-    @staticmethod
-    def make_shared(grid_factory: GridFactory) -> MultipleGridBuilder: ...
     @overload
-    def add_coarse_grid(self, start_x: float, start_y: float, start_z: float, end_x: float, end_y: float, end_z: float, delta: float) -> None: ...
+    def add_coarse_grid(
+        self, start_x: float, start_y: float, start_z: float, end_x: float, end_y: float, end_z: float, delta: float
+    ) -> None: ...
     @overload
     def add_coarse_grid(self, gridDimensions: GridDimensions) -> None: ...
     @overload
@@ -111,13 +124,11 @@ class MultipleGridBuilder(LevelGridBuilder):
     def find_communication_indices(self, direction: int, do_shift: bool = False) -> None: ...
     def get_number_of_levels(self) -> int: ...
 
-
 class Sphere(Object):
     def __init__(self, *args, **kwargs) -> None: ...
     @staticmethod
     def make_shared() -> Sphere: ...
 
-
 class TriangularMesh(Object):
     def __init__(self, *args, **kwargs) -> None: ...
     @staticmethod
diff --git a/pythonbindings/pyfluids-stubs/gpu/kernel/__init__.pyi b/pythonbindings/pyfluids-stubs/gpu/kernel/__init__.pyi
index 4f2a6f283019d6eed3d92be65ee959a10e7bb064..ab3699018255e75e4b3b90613338d37c1a5269eb 100644
--- a/pythonbindings/pyfluids-stubs/gpu/kernel/__init__.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/kernel/__init__.pyi
@@ -30,5 +30,6 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 from . import compressible as compressible
-from . import incompressible as incompressible
\ No newline at end of file
+from . import incompressible as incompressible
diff --git a/pythonbindings/pyfluids-stubs/gpu/kernel/compressible.pyi b/pythonbindings/pyfluids-stubs/gpu/kernel/compressible.pyi
index eda5289d766d53ef82303ae6f3ee22db78cdb108..fb04e03a79caae76bb7fa5833f8fde0556a56907 100644
--- a/pythonbindings/pyfluids-stubs/gpu/kernel/compressible.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/kernel/compressible.pyi
@@ -30,6 +30,7 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 BGK: str
 BGKPlus: str
 K17CompressibleNavierStokes: str
diff --git a/pythonbindings/pyfluids-stubs/gpu/kernel/incompressible.pyi b/pythonbindings/pyfluids-stubs/gpu/kernel/incompressible.pyi
index c7ed80a1b8d8c0a0df28a478442c2d6b49478980..910326c8140e7df0ded7ec33fe8cd1b7c0e044b9 100644
--- a/pythonbindings/pyfluids-stubs/gpu/kernel/incompressible.pyi
+++ b/pythonbindings/pyfluids-stubs/gpu/kernel/incompressible.pyi
@@ -30,6 +30,7 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
+
 BGK: str
 BGKPlus: str
-CumulantK15: str
\ No newline at end of file
+CumulantK15: str
diff --git a/pythonbindings/pyfluids-stubs/gpu/probes.pyi b/pythonbindings/pyfluids-stubs/gpu/probes.pyi
deleted file mode 100644
index 0073d2812024be26dc3eefb1baf4233aadaef30f..0000000000000000000000000000000000000000
--- a/pythonbindings/pyfluids-stubs/gpu/probes.pyi
+++ /dev/null
@@ -1,97 +0,0 @@
-r"""
-=======================================================================================
- ____          ____    __    ______     __________   __      __       __        __
- \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
-  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
-   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
-    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
-     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
-      \    \  |    |   ________________________________________________________________
-       \    \ |    |  |  ______________________________________________________________|
-        \    \|    |  |  |         __          __     __     __     ______      _______
-         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
-          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
-           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
-            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
-
-  This file is part of VirtualFluids. VirtualFluids is free software: you can
-  redistribute it and/or modify it under the terms of the GNU General Public
-  License as published by the Free Software Foundation, either version 3 of
-  the License, or (at your option) any later version.
-
-  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
-  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
-  for more details.
-
-  SPDX-License-Identifier: GPL-3.0-or-later
-  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
-
-! \author Henry Korb
-=======================================================================================
-"""
-from __future__ import annotations
-from typing import ClassVar, List
-from enum import Enum
-
-import gpu
-
-
-class Statistic:
-    __members__: ClassVar[dict] = ...  # read-only
-    Instantaneous: ClassVar[Statistic] = ...
-    Means: ClassVar[Statistic] = ...
-    SpatialCovariances: ClassVar[Statistic] = ...
-    SpatialFlatness: ClassVar[Statistic] = ...
-    SpatialMeans: ClassVar[Statistic] = ...
-    SpatialSkewness: ClassVar[Statistic] = ...
-    SpatioTemporalCovariances: ClassVar[Statistic] = ...
-    SpatioTemporalFlatness: ClassVar[Statistic] = ...
-    SpatioTemporalMeans: ClassVar[Statistic] = ...
-    SpatioTemporalSkewness: ClassVar[Statistic] = ...
-    Variances: ClassVar[Statistic] = ...
-    __entries: ClassVar[dict] = ...
-    def __init__(self, arg0: int) -> None: ...
-    def __eq__(self, arg0: object) -> bool: ...
-    def __getstate__(self) -> int: ...
-    def __hash__(self) -> int: ...
-    def __index__(self) -> int: ...
-    def __int__(self) -> int: ...
-    def __ne__(self, arg0: object) -> bool: ...
-    def __setstate__(self, arg0: int) -> None: ...
-    @property
-    def name(self) -> str: ...
-
-
-class Probe(gpu.Sampler):
-    def __init__(self, *args, **kwargs) -> None: ...
-    def add_all_available_statistics(self) -> None: ...
-    def add_statistic(self, variable: Statistic) -> None: ...
-    def set_file_name_to_n_out(self) -> None: ...
-
-class PlaneNormal(Enum):
-    X = ...
-    Y = ...
-    Z = ...
-
-
-class PlanarAverageProbe(Probe):
-    def __init__(self, para: gpu.Parameter, cuda_memory_manager: gpu.CudaMemoryManager, probe_name: str, output_path: str, t_start_avg: int, t_start_tmp_avg: int, t_avg: int, t_start_out: int, t_out: int, plane_normal: PlaneNormal, average_every_timestep: bool) -> None: ...
-
-
-class PlaneProbe(Probe):
-    def __init__(self, para: gpu.Parameter, cuda_memory_manager: gpu.CudaMemoryManager, probe_name: str, output_path: str, t_start_avg: int, t_avg: int, t_start_out: int, t_out: int) -> None: ...
-    def set_probe_plane(self, pos_x: float, pos_y: float, pos_z: float, delta_x: float, delta_y: float, delta_z: float) -> None: ...
-
-
-class PointProbe(Probe):
-    def __init__(self, para: gpu.Parameter, cuda_memory_manager: gpu.CudaMemoryManager, probe_name: str, output_path: str, t_start_avg: int, t_avg: int, t_start_out: int, t_out: int, output_timeseries: bool) -> None: ...
-    def add_probe_point(self, point_coord_x: float, point_coord_y: float, point_coord_z: float) -> None: ...
-    def add_probe_points_from_list(self, point_coords_x: List[float], point_coords_y: List[float], point_coords_z: List[float]) -> None: ...
-    def add_probe_points_from_x_normal_plane(self, pos_x: float, pos0_y: float, pos0_z: float, pos1_y: float, pos1_z: float, n_y: int, n_z: int) -> None: ...
-
-
-class WallModelProbe(gpu.Sampler):
-    def __init__(self, para: gpu.Parameter, cuda_memory_manager: gpu.CudaMemoryManager, probe_name: str, output_path: str, t_start_avg: int, t_start_tmp_avg: int, t_avg: int, t_start_out: int, t_out: int, average_every_timestep: bool, compute_temporal_averages: bool, output_stress: bool, evaluate_pressure_gradient: bool) -> None: ...
-    def set_evaluate_pressure_gradient(self, eval_press_grad: bool) -> None: ...
-    def set_force_output_to_stress(self, output_stress: bool) -> None: ...
diff --git a/pythonbindings/pyfluids-stubs/gpu/probes/__init__.pyi b/pythonbindings/pyfluids-stubs/gpu/probes/__init__.pyi
new file mode 100644
index 0000000000000000000000000000000000000000..d272dab3f358f77c20743a3fd4dfc27bb7108112
--- /dev/null
+++ b/pythonbindings/pyfluids-stubs/gpu/probes/__init__.pyi
@@ -0,0 +1,36 @@
+r"""
+=======================================================================================
+ ____          ____    __    ______     __________   __      __       __        __
+ \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+      \    \  |    |   ________________________________________________________________
+       \    \ |    |  |  ______________________________________________________________|
+        \    \|    |  |  |         __          __     __     __     ______      _______
+         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+
+  This file is part of VirtualFluids. VirtualFluids is free software: you can
+  redistribute it and/or modify it under the terms of the GNU General Public
+  License as published by the Free Software Foundation, either version 3 of
+  the License, or (at your option) any later version.
+
+  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+  for more details.
+
+  SPDX-License-Identifier: GPL-3.0-or-later
+  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+
+! \author Henry Korb
+=======================================================================================
+"""
+
+from __future__ import annotations
+from . import probe as probe
+from . import planar_average_probe as planar_average_probe
diff --git a/pythonbindings/pyfluids-stubs/gpu/probes/planar_average_probe.pyi b/pythonbindings/pyfluids-stubs/gpu/probes/planar_average_probe.pyi
new file mode 100644
index 0000000000000000000000000000000000000000..62930434aafb38f5078d03fba48413da395e7f04
--- /dev/null
+++ b/pythonbindings/pyfluids-stubs/gpu/probes/planar_average_probe.pyi
@@ -0,0 +1,64 @@
+r"""
+=======================================================================================
+ ____          ____    __    ______     __________   __      __       __        __
+ \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+      \    \  |    |   ________________________________________________________________
+       \    \ |    |  |  ______________________________________________________________|
+        \    \|    |  |  |         __          __     __     __     ______      _______
+         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+
+  This file is part of VirtualFluids. VirtualFluids is free software: you can
+  redistribute it and/or modify it under the terms of the GNU General Public
+  License as published by the Free Software Foundation, either version 3 of
+  the License, or (at your option) any later version.
+
+  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+  for more details.
+
+  SPDX-License-Identifier: GPL-3.0-or-later
+  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+
+! \author Henry Korb
+=======================================================================================
+"""
+
+from __future__ import annotations
+
+import enum
+from ... import gpu
+from ...basics.geometry3d import Axis
+
+class Statistic(enum.Enum):
+    Means: int
+    Covariances: int
+    Skewness: int
+    Flatness: int
+
+class PlanarAverageProbe(gpu.Sampler):
+    def __init__(
+        self,
+        para: gpu.Parameter,
+        cuda_memory_manager: gpu.CudaMemoryManager,
+        probe_name: str,
+        output_path: str,
+        t_start_avg: int,
+        t_start_tmp_avg: int,
+        t_avg: int,
+        t_start_out: int,
+        t_out: int,
+        plane_normal: Axis,
+        compute_time_averages: bool,
+        compute_statistics_of_concentration: bool,
+    ) -> None: ...
+    def add_statistic(self, variable: Statistic) -> None: ...
+    def add_all_available_statistics(self) -> None: ...
+    def set_file_name_to_n_out(self) -> None: ...
diff --git a/pythonbindings/pyfluids-stubs/gpu/probes/probe.pyi b/pythonbindings/pyfluids-stubs/gpu/probes/probe.pyi
new file mode 100644
index 0000000000000000000000000000000000000000..39b3a34cd713a00a79470568de58a1cc018984b4
--- /dev/null
+++ b/pythonbindings/pyfluids-stubs/gpu/probes/probe.pyi
@@ -0,0 +1,67 @@
+r"""
+=======================================================================================
+ ____          ____    __    ______     __________   __      __       __        __
+ \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+      \    \  |    |   ________________________________________________________________
+       \    \ |    |  |  ______________________________________________________________|
+        \    \|    |  |  |         __          __     __     __     ______      _______
+         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+
+  This file is part of VirtualFluids. VirtualFluids is free software: you can
+  redistribute it and/or modify it under the terms of the GNU General Public
+  License as published by the Free Software Foundation, either version 3 of
+  the License, or (at your option) any later version.
+
+  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+  for more details.
+
+  SPDX-License-Identifier: GPL-3.0-or-later
+  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+
+! \author Henry Korb
+=======================================================================================
+"""
+
+from __future__ import annotations
+from enum import Enum
+
+from ... import gpu
+
+class Statistic(Enum):
+    Instantaneous: int
+    Means: int
+    Variances: int
+
+class Probe(gpu.Sampler):
+    def __init__(
+        self,
+        para: gpu.Parameter,
+        cuda_memory_manager: gpu.CudaMemoryManager,
+        probe_name: str,
+        output_path: str,
+        t_start_avg: int,
+        t_start_tmp_avg: int,
+        t_avg: int,
+        t_start_out: int,
+        t_out: int,
+        output_timeseries: bool,
+        average_every_timestep: bool,
+    ) -> None: ...
+    def add_statistic(self, variable: Statistic) -> None: ...
+    def add_all_available_statistics(self) -> None: ...
+    def add_probe_point(self, point_coord_x: float, point_coord_y: float, point_coord_z: float) -> None: ...
+    def add_probe_points_from_list(
+        self, point_coords_x: list[float], point_coords_y: list[float], point_coords_z: list[float]
+    ) -> None: ...
+    def set_probe_plane(
+        self, pos_x: float, pos_y: float, pos_z: float, delta_x: float, delta_y: float, delta_z: float
+    ) -> None: ...
diff --git a/pythonbindings/pyfluids-stubs/gpu/probes/wall_model_probe.pyi b/pythonbindings/pyfluids-stubs/gpu/probes/wall_model_probe.pyi
new file mode 100644
index 0000000000000000000000000000000000000000..20dd87fe480bacafab62f32cebfb62f6cb81af22
--- /dev/null
+++ b/pythonbindings/pyfluids-stubs/gpu/probes/wall_model_probe.pyi
@@ -0,0 +1,54 @@
+r"""
+=======================================================================================
+ ____          ____    __    ______     __________   __      __       __        __
+ \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+      \    \  |    |   ________________________________________________________________
+       \    \ |    |  |  ______________________________________________________________|
+        \    \|    |  |  |         __          __     __     __     ______      _______
+         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+
+  This file is part of VirtualFluids. VirtualFluids is free software: you can
+  redistribute it and/or modify it under the terms of the GNU General Public
+  License as published by the Free Software Foundation, either version 3 of
+  the License, or (at your option) any later version.
+
+  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+  for more details.
+
+  SPDX-License-Identifier: GPL-3.0-or-later
+  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+
+! \author Henry Korb
+=======================================================================================
+"""
+
+from __future__ import annotations
+
+from ... import gpu
+
+class WallModelProbe(gpu.Sampler):
+    def __init__(
+        self,
+        para: gpu.Parameter,
+        cuda_memory_manager: gpu.CudaMemoryManager,
+        probe_name: str,
+        output_path: str,
+        t_start_avg: int,
+        t_start_tmp_avg: int,
+        t_avg: int,
+        t_start_out: int,
+        t_out: int,
+        average_every_timestep: bool,
+        compute_temporal_averages: bool,
+        output_stress: bool,
+        evaluate_pressure_gradient: bool,
+    ) -> None: ...
diff --git a/pythonbindings/pyfluids-stubs/lbm.pyi b/pythonbindings/pyfluids-stubs/lbm.pyi
index 1c576a0ecca269b8dc7db0a506dbb1136cf6e4f5..62f00efbc491ca3ecca9cab3f34fbf53ba69ce87 100644
--- a/pythonbindings/pyfluids-stubs/lbm.pyi
+++ b/pythonbindings/pyfluids-stubs/lbm.pyi
@@ -29,4 +29,4 @@ r"""
 
 ! \author Henry Korb
 =======================================================================================
-"""
\ No newline at end of file
+"""
diff --git a/pythonbindings/pyfluids-stubs/logger.pyi b/pythonbindings/pyfluids-stubs/logger.pyi
index 3d736abedf93e6200d9230fd76f0754cc23b76ad..17e64d52a517fb4e0cc7eb60512c385676cdf8db 100644
--- a/pythonbindings/pyfluids-stubs/logger.pyi
+++ b/pythonbindings/pyfluids-stubs/logger.pyi
@@ -30,8 +30,8 @@ r"""
 ! \author Henry Korb
 =======================================================================================
 """
-from __future__ import annotations
 
+from __future__ import annotations
 
 class Logger:
     @staticmethod
@@ -39,7 +39,6 @@ class Logger:
     @staticmethod
     def initialize_logger() -> None: ...
 
-
 def vf_log_critical(message: str) -> None: ...
 def vf_log_debug(message: str) -> None: ...
 def vf_log_info(message: str) -> None: ...
diff --git a/pythonbindings/pyfluids-stubs/parallel.pyi b/pythonbindings/pyfluids-stubs/parallel.pyi
index f5d20a55d982dbf80791c180d9c3a5f320782f71..5430788277815b2abd42337f7249023fc7fafff6 100644
--- a/pythonbindings/pyfluids-stubs/parallel.pyi
+++ b/pythonbindings/pyfluids-stubs/parallel.pyi
@@ -43,4 +43,4 @@ class Communicator:
 
 class MPICommunicator(Communicator):
     @staticmethod
-    def get_instance() -> MPICommunicator: ...
\ No newline at end of file
+    def get_instance() -> MPICommunicator: ...
diff --git a/pythonbindings/src/basics/basics.cpp b/pythonbindings/src/basics/basics.cpp
index 99bf8cc8d378dc21ab0064935e20a63aba890a3b..610c6789d874394c94bbee92f1b96ae73a550a9b 100644
--- a/pythonbindings/src/basics/basics.cpp
+++ b/pythonbindings/src/basics/basics.cpp
@@ -30,6 +30,7 @@
 //=======================================================================================
 #include <pybind11/pybind11.h>
 #include "submodules/configuration_file.cpp"
+#include "submodules/geometry3d.cpp"
 
 namespace basics
 {
@@ -38,5 +39,6 @@ namespace basics
     PYBIND11_MODULE(basics, m)
     {
         configuration::makeModule(m);
+        geometry3d::makeModule(m);
     }
 }
\ No newline at end of file
diff --git a/pythonbindings/src/basics/submodules/geometry3d.cpp b/pythonbindings/src/basics/submodules/geometry3d.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..96686485eacd634b10c23978e904beb216bf554a
--- /dev/null
+++ b/pythonbindings/src/basics/submodules/geometry3d.cpp
@@ -0,0 +1,45 @@
+//=======================================================================================
+// ____          ____    __    ______     __________   __      __       __        __
+// \    \       |    |  |  |  |   _   \  |___    ___| |  |    |  |     /  \      |  |
+//  \    \      |    |  |  |  |  |_)   |     |  |     |  |    |  |    /    \     |  |
+//   \    \     |    |  |  |  |   _   /      |  |     |  |    |  |   /  /\  \    |  |
+//    \    \    |    |  |  |  |  | \  \      |  |     |   \__/   |  /  ____  \   |  |____
+//     \    \   |    |  |__|  |__|  \__\     |__|      \________/  /__/    \__\  |_______|
+//      \    \  |    |   ________________________________________________________________
+//       \    \ |    |  |  ______________________________________________________________|
+//        \    \|    |  |  |         __          __     __     __     ______      _______
+//         \         |  |  |_____   |  |        |  |   |  |   |  |   |   _  \    /  _____)
+//          \        |  |   _____|  |  |        |  |   |  |   |  |   |  | \  \   \_______
+//           \       |  |  |        |  |_____   |   \_/   |   |  |   |  |_/  /    _____  |
+//            \ _____|  |__|        |________|   \_______/    |__|   |______/    (_______/
+//
+//  This file is part of VirtualFluids. VirtualFluids is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+//  for more details.
+//
+//  SPDX-License-Identifier: GPL-3.0-or-later
+//  SPDX-FileCopyrightText: Copyright © VirtualFluids Project contributors, see AUTHORS.md in root folder
+//
+//! \author Henry Korb
+//=======================================================================================
+#include <pybind11/pybind11.h>
+#include <basics/geometry3d/Axis.h>
+
+namespace geometry3d {
+    namespace py = pybind11;
+
+    void makeModule(py::module_ &parentModule)
+    {
+        py::module geometry3dModule = parentModule.def_submodule("geometry3d");
+        py::enum_<Axis>(geometry3dModule, "Axis")
+        .value("x", Axis::x)
+        .value("y", Axis::y)
+        .value("z", Axis::z);
+    }
+}
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/grid_generator.cpp b/pythonbindings/src/gpu/submodules/grid_generator.cpp
index 0bc9e80f39412a6e8e30bcaa94d68da55c1729d1..9ac3731ba91ceb97b634a83b496949b03bd284f1 100644
--- a/pythonbindings/src/gpu/submodules/grid_generator.cpp
+++ b/pythonbindings/src/gpu/submodules/grid_generator.cpp
@@ -81,7 +81,7 @@ namespace grid_generator
         py::class_<TriangularMesh, Object, std::shared_ptr<TriangularMesh>>(gridGeneratorModule, "TriangularMesh")
         .def_static("make", &TriangularMesh::make, py::return_value_policy::reference);
 
-        py::class_<GridDimensions, std::shared_ptr<GridDimensions>>(gridGeneratorModule, "GridDimension")
+        py::class_<GridDimensions, std::shared_ptr<GridDimensions>>(gridGeneratorModule, "GridDimensions")
         .def(py::init<real, real, real, real, real, real, real>(), py::arg("min_x"), py::arg("max_x"), py::arg("min_y"), py::arg("max_y"), py::arg("min_z"), py::arg("max_z"), py::arg("delta"));
 
         py::class_<GridBuilder, std::shared_ptr<GridBuilder>>(gridGeneratorModule, "GridBuilder")
diff --git a/pythonbindings/src/gpu/submodules/probes.cpp b/pythonbindings/src/gpu/submodules/probes.cpp
index 93f23cea40db16cd46d21b4480c4754f6a2b3a35..fe4562b45f7d18eaa54da3988f40477114bbb1c2 100644
--- a/pythonbindings/src/gpu/submodules/probes.cpp
+++ b/pythonbindings/src/gpu/submodules/probes.cpp
@@ -30,6 +30,7 @@
 //=======================================================================================
 #include <pybind11/pybind11.h>
 #include <pybind11/stl.h>
+#include <basics/geometry3d/Axis.h>
 #include <gpu/core/Samplers/Probe.h>
 #include <gpu/core/Samplers/WallModelProbe.h>
 #include <gpu/core/Samplers/PlanarAverageProbe.h>
@@ -43,7 +44,7 @@ namespace probes
     {
         py::module probeModule = parentModule.def_submodule("probes");
 
-        py::module probeProbeModule = probeModule.def_submodule("Probe");
+        py::module probeProbeModule = probeModule.def_submodule("probe");
 
         py::enum_<Probe::Statistic>(probeProbeModule, "Statistic")
         .value("Instantaneous", Probe::Statistic::Instantaneous)
@@ -78,12 +79,7 @@ namespace probes
         .def("add_probe_points_from_list", &Probe::addProbePointsFromList, py::arg("point_coords_x"), py::arg("point_coords_y"), py::arg("point_coords_z"))
         .def("set_probe_plane", &Probe::addProbePlane, py::arg("pos_x"), py::arg("pos_y"), py::arg("pos_z"), py::arg("delta_x"), py::arg("delta_y"), py::arg("delta_z"));
 
-        py::module planarAverageProbeModule = probeModule.def_submodule("PlanarAverageProbe");
-
-        py::enum_<PlanarAverageProbe::PlaneNormal>(planarAverageProbeModule, "PlaneNormal")
-        .value("x", PlanarAverageProbe::PlaneNormal::x)
-        .value("y", PlanarAverageProbe::PlaneNormal::y)
-        .value("z", PlanarAverageProbe::PlaneNormal::z);
+        py::module planarAverageProbeModule = probeModule.def_submodule("planar_average_probe");
 
         py::enum_<PlanarAverageProbe::Statistic>(planarAverageProbeModule, "Statistic")
         .value("Means", PlanarAverageProbe::Statistic::Means)
@@ -101,7 +97,8 @@ namespace probes
                         uint,
                         uint,
                         uint,
-                        PlanarAverageProbe::PlaneNormal,
+                        Axis,
+                        bool,
                         bool>(),
                         py::arg("para"),
                         py::arg("cuda_memory_manager"),
@@ -113,10 +110,16 @@ namespace probes
                         py::arg("t_start_out"),
                         py::arg("t_out"),
                         py::arg("plane_normal"),
-                        py::arg("compute_time_averages"));
+                        py::arg("compute_time_averages"),
+                        py::arg("compute_statistics_of_concentration"))
+        .def("add_statistic", &PlanarAverageProbe::addStatistic, py::arg("statistic"))
+        .def("add_all_available_statistics", &PlanarAverageProbe::addAllAvailableStatistics)
+        .def("set_file_name_to_n_out", &PlanarAverageProbe::setFileNameToNOut);
+
+        py::module wallModelProbeModule = probeModule.def_submodule("wall_model_probe");
 
 
-        py::class_<WallModelProbe, Sampler, std::shared_ptr<WallModelProbe>>(probeModule, "WallModelProbe")
+        py::class_<WallModelProbe, Sampler, std::shared_ptr<WallModelProbe>>(wallModelProbeModule, "WallModelProbe")
         .def(py::init<  SPtr<Parameter>,
                         SPtr<CudaMemoryManager>,
                         const std::string,
diff --git a/src/gpu/core/Samplers/PlanarAverageProbe.cu b/src/gpu/core/Samplers/PlanarAverageProbe.cu
index 7bb3003d0d905202fa8e8b70e5c0d4d1c25b1efe..9d32499f16ea30e48df1fa724e4a99f21a52d672 100644
--- a/src/gpu/core/Samplers/PlanarAverageProbe.cu
+++ b/src/gpu/core/Samplers/PlanarAverageProbe.cu
@@ -35,7 +35,7 @@
 #include <cmath>
 #include <stdexcept>
 #include <string>
-#include <tuple>
+#include <cstddef>
 
 #include <helper_cuda.h>
 
@@ -47,6 +47,7 @@
 
 #include <basics/DataTypes.h>
 #include <basics/constants/NumericConstants.h>
+#include <basics/geometry3d/Axis.h>
 #include <basics/utilities/UbTuple.h>
 #include <basics/writer/WbWriterVtkXmlBinary.h>
 
@@ -79,7 +80,7 @@ iterPair getPermutationIterators(real* values, unsigned long long* indices, uint
 struct shiftIndex
 {
     const uint* neighborNormal;
-    __host__ __device__ unsigned long long operator()(unsigned long long& pointIndices)
+    constexpr unsigned long long operator()(unsigned long long& pointIndices) const
     {
         return neighborNormal[pointIndices];
     }
@@ -93,7 +94,7 @@ struct covariance
     }
 
     template <typename Tuple>
-    __host__ __device__ real operator()(const Tuple& t) const
+    constexpr real operator()(const Tuple& t) const
     {
         return (thrust::get<0>(t) - mean_x) * (thrust::get<1>(t) - mean_y);
     }
@@ -105,7 +106,7 @@ struct skewness
     skewness(real mean) : mean(mean)
     {
     }
-    __host__ __device__ real operator()(const real& x) const
+    constexpr real operator()(const real& x) const
     {
         return (x - mean) * (x - mean) * (x - mean);
     }
@@ -117,7 +118,7 @@ struct flatness
     flatness(real mean) : mean(mean)
     {
     }
-    __host__ __device__ real operator()(const real& x) const
+    constexpr real operator()(const real& x) const
     {
         return (x - mean) * (x - mean) * (x - mean) * (x - mean);
     }
@@ -125,22 +126,23 @@ struct flatness
 
 struct Means
 {
-    real vx, vy, vz;
+    // phi is the scalar
+    real vx, vy, vz, phi;
 };
 
 struct Covariances
 {
-    real vxvx, vyvy, vzvz, vxvy, vxvz, vyvz;
+    real vxvx, vyvy, vzvz, vxvy, vxvz, vyvz, phiphi, vxphi, vyphi, vzphi;
 };
 
 struct Skewnesses
 {
-    real Sx, Sy, Sz;
+    real Sx, Sy, Sz, Sphi;
 };
 
 struct Flatnesses
 {
-    real Fx, Fy, Fz;
+    real Fx, Fy, Fz, Fphi;
 };
 
 ///////////////////////////////////////////////////////////////////////////////////
@@ -164,8 +166,7 @@ std::string getName(std::string quantity, bool timeAverage)
 }
 
 ///////////////////////////////////////////////////////////////////////////////////
-std::vector<std::string> PlanarAverageProbe::getVariableNames(PlanarAverageProbe::Statistic statistic,
-                                                              bool namesForTimeAverages)
+std::vector<std::string> PlanarAverageProbe::getVariableNames(Statistic statistic, bool namesForTimeAverages) const
 {
 
     std::vector<std::string> variableNames;
@@ -174,9 +175,10 @@ std::vector<std::string> PlanarAverageProbe::getVariableNames(PlanarAverageProbe
             variableNames.emplace_back(getName("vx", namesForTimeAverages));
             variableNames.emplace_back(getName("vy", namesForTimeAverages));
             variableNames.emplace_back(getName("vz", namesForTimeAverages));
-            if (para->getUseTurbulentViscosity()) {
+            if (para->getUseTurbulentViscosity())
                 variableNames.emplace_back(getName("EddyViscosity", namesForTimeAverages));
-            }
+            if (sampleScalar)
+                variableNames.emplace_back(getName("phi", namesForTimeAverages));
             break;
         case Statistic::Covariances:
             variableNames.emplace_back(getName("vxvx", namesForTimeAverages));
@@ -185,16 +187,26 @@ std::vector<std::string> PlanarAverageProbe::getVariableNames(PlanarAverageProbe
             variableNames.emplace_back(getName("vxvy", namesForTimeAverages));
             variableNames.emplace_back(getName("vxvz", namesForTimeAverages));
             variableNames.emplace_back(getName("vyvz", namesForTimeAverages));
+            if (sampleScalar) {
+                variableNames.emplace_back(getName("phiphi", namesForTimeAverages));
+                variableNames.emplace_back(getName("vxphi", namesForTimeAverages));
+                variableNames.emplace_back(getName("vyphi", namesForTimeAverages));
+                variableNames.emplace_back(getName("vzphi", namesForTimeAverages));
+            }
             break;
         case Statistic::Skewness:
             variableNames.emplace_back(getName("SkewnessX", namesForTimeAverages));
             variableNames.emplace_back(getName("SkewnessY", namesForTimeAverages));
             variableNames.emplace_back(getName("SkewnessZ", namesForTimeAverages));
+            if (sampleScalar)
+                variableNames.emplace_back(getName("SkewnessPhi", namesForTimeAverages));
             break;
         case Statistic::Flatness:
             variableNames.emplace_back(getName("FlatnessX", namesForTimeAverages));
             variableNames.emplace_back(getName("FlatnessY", namesForTimeAverages));
             variableNames.emplace_back(getName("FlatnessZ", namesForTimeAverages));
+            if (sampleScalar)
+                variableNames.emplace_back(getName("FlatnessPhi", namesForTimeAverages));
             break;
 
         default:
@@ -205,12 +217,29 @@ std::vector<std::string> PlanarAverageProbe::getVariableNames(PlanarAverageProbe
     return variableNames;
 }
 
-void PlanarAverageProbe::addStatistic(PlanarAverageProbe::Statistic statistic)
+void PlanarAverageProbe::addStatistic(Statistic statistic)
 {
     if (!isStatisticIn(statistic, statistics))
         statistics.push_back(statistic);
 }
 
+PlanarAverageProbe::PlanarAverageProbe(SPtr<Parameter> para, SPtr<CudaMemoryManager> cudaMemoryManager,
+                                       std::string outputPath, std::string probeName, uint tStartSampling,
+                                       uint tStartTemporalAveraging, uint tBetweenSamples, uint tStartWritingOutput,
+                                       uint tBetweenWriting, Axis planeNormal, bool computeTimeAverages, bool sampleScalar)
+    : para(para), cudaMemoryManager(cudaMemoryManager), tStartSampling(tStartSampling),
+      tStartTemporalAveraging(tStartTemporalAveraging), tBetweenSamples(tBetweenSamples),
+      tStartWritingOutput(tStartWritingOutput), tBetweenWriting(tBetweenWriting), computeTimeAverages(computeTimeAverages),
+      planeNormal(planeNormal), sampleScalar(sampleScalar), Sampler(outputPath, probeName)
+{
+    if (tStartTemporalAveraging < tStartSampling && computeTimeAverages)
+        throw std::runtime_error("PlaneAverageProbe: tStartTemporalAveraging must be larger than tStartSampling!");
+    if (tBetweenWriting == 0)
+        throw std::runtime_error("PlaneAverageProbe: tBetweenWriting must be larger than 0!");
+    if (sampleScalar && !para->getDiffOn())
+        throw std::runtime_error("PlaneAverageProbe: Scalar can only be sampled if diff is on!");
+}
+
 ///////////////////////////////////////////////////////////////////////////////////
 void PlanarAverageProbe::init()
 {
@@ -237,16 +266,16 @@ void PlanarAverageProbe::init()
 
 void PlanarAverageProbe::sample(int level, uint t)
 {
-    if (t < tStartAveraging)
+    if (t < tStartSampling)
         return;
 
     const uint tLevel = para->getTimeStep(level, t, true);
     const uint levelFactor = exp2(level);
-    const uint tAfterStartAveraging = tLevel - tStartAveraging * levelFactor;
+    const uint tAfterStartAveraging = tLevel - tStartSampling * levelFactor;
     const uint tAfterStartWriting = tLevel - tStartWritingOutput * levelFactor;
     const bool doTimeAverages = computeTimeAverages && tLevel > (tStartTemporalAveraging * levelFactor);
 
-    if (tAfterStartAveraging % (tBetweenAverages * levelFactor) == 0) {
+    if (tAfterStartAveraging % (tBetweenSamples * levelFactor) == 0) {
         calculateQuantities(level, doTimeAverages);
     }
 
@@ -273,11 +302,11 @@ std::vector<unsigned long long> PlanarAverageProbe::findIndicesInPlane(int level
 
     const real* coordinatesPlaneNormal = [&] {
         switch (planeNormal) {
-            case PlaneNormal::x:
+            case Axis::x:
                 return param->coordinateX;
-            case PlaneNormal::y:
+            case Axis::y:
                 return param->coordinateY;
-            case PlaneNormal::z:
+            case Axis::z:
                 return param->coordinateZ;
             default:
                 throw std::runtime_error("PlaneNormal not defined!");
@@ -304,19 +333,19 @@ void PlanarAverageProbe::findCoordinatesForPlanes(int level, std::vector<real>&
     unsigned long long nextIndex = startIndex;
     do {
         switch (planeNormal) {
-            case PlaneNormal::x:
+            case Axis::x:
                 coordinateX.push_back(para->getParH(level)->coordinateX[nextIndex]);
                 coordinateY.push_back(999999);
                 coordinateZ.push_back(999999);
                 nextIndex = para->getParH(level)->neighborX[nextIndex];
                 break;
-            case PlaneNormal::y:
+            case Axis::y:
                 coordinateX.push_back(999999);
                 coordinateY.push_back(para->getParH(level)->coordinateY[nextIndex]);
                 coordinateZ.push_back(999999);
                 nextIndex = para->getParH(level)->neighborY[nextIndex];
                 break;
-            case PlaneNormal::z:
+            case Axis::z:
                 coordinateX.push_back(999999);
                 coordinateY.push_back(999999);
                 coordinateZ.push_back(para->getParH(level)->coordinateZ[nextIndex]);
@@ -334,12 +363,14 @@ real computeMean(iterPair x, real invNPointsPerPlane)
     return sum * invNPointsPerPlane;
 }
 
-Means computeMeans(iterPair vx, iterPair vy, iterPair vz, real invNPointsPerPlane)
+Means computeMeans(iterPair vx, iterPair vy, iterPair vz, iterPair phi, real invNPointsPerPlane, bool sampleScalar)
 {
     Means means;
     means.vx = computeMean(vx, invNPointsPerPlane);
     means.vy = computeMean(vy, invNPointsPerPlane);
     means.vz = computeMean(vz, invNPointsPerPlane);
+    if (sampleScalar)
+        means.phi = computeMean(phi, invNPointsPerPlane);
     return means;
 }
 
@@ -350,7 +381,8 @@ real computeCovariance(iterPair x, iterPair y, real mean_x, real mean_y, real in
     return thrust::transform_reduce(begin, end, covariance(mean_x, mean_y), c0o1, thrust::plus<real>()) * invNPointsPerPlane;
 }
 
-Covariances computeCovariances(iterPair vx, iterPair vy, iterPair vz, Means means, real invNPointsPerPlane)
+Covariances computeCovariances(iterPair vx, iterPair vy, iterPair vz, iterPair phi, Means means, real invNPointsPerPlane,
+                               bool sampleScalar)
 {
     Covariances covariances;
 
@@ -360,6 +392,12 @@ Covariances computeCovariances(iterPair vx, iterPair vy, iterPair vz, Means mean
     covariances.vxvy = computeCovariance(vx, vy, means.vx, means.vy, invNPointsPerPlane);
     covariances.vxvz = computeCovariance(vx, vz, means.vx, means.vz, invNPointsPerPlane);
     covariances.vyvz = computeCovariance(vy, vz, means.vy, means.vz, invNPointsPerPlane);
+    if (sampleScalar) {
+        covariances.phiphi = computeCovariance(phi, phi, means.phi, means.phi, invNPointsPerPlane);
+        covariances.vxphi = computeCovariance(vx, phi, means.vx, means.phi, invNPointsPerPlane);
+        covariances.vyphi = computeCovariance(vy, phi, means.vy, means.phi, invNPointsPerPlane);
+        covariances.vzphi = computeCovariance(vz, phi, means.vz, means.phi, invNPointsPerPlane);
+    }
 
     return covariances;
 }
@@ -367,17 +405,19 @@ Covariances computeCovariances(iterPair vx, iterPair vy, iterPair vz, Means mean
 real computeSkewness(iterPair x, real mean, real covariance, real invNPointsPerPlane)
 {
     return thrust::transform_reduce(x.first, x.second, skewness(mean), c0o1, thrust::plus<real>()) * invNPointsPerPlane *
-           pow(covariance, -1.5f);
+           std::pow(covariance, -c3o2);
 }
 
-Skewnesses computeSkewnesses(Means means, Covariances covariances, iterPair vx, iterPair vy, iterPair vz,
-                             real invNPointsPerPlane)
+Skewnesses computeSkewnesses(Means means, Covariances covariances, iterPair vx, iterPair vy, iterPair vz, iterPair phi,
+                             real invNPointsPerPlane, bool sampleScalar)
 {
     Skewnesses skewnesses;
 
     skewnesses.Sx = computeSkewness(vx, means.vx, covariances.vxvx, invNPointsPerPlane);
     skewnesses.Sy = computeSkewness(vy, means.vy, covariances.vyvy, invNPointsPerPlane);
     skewnesses.Sz = computeSkewness(vz, means.vz, covariances.vzvz, invNPointsPerPlane);
+    if (sampleScalar)
+        skewnesses.Sphi = computeSkewness(phi, means.phi, covariances.phiphi, invNPointsPerPlane);
 
     return skewnesses;
 }
@@ -388,15 +428,16 @@ real computeFlatness(iterPair x, real mean, real covariance, real invNPointsPerP
            pow(covariance, -c2o1);
 }
 
-Flatnesses computeFlatnesses(iterPair vx, iterPair vy, iterPair vz, Means means, Covariances covariances,
-                             real invNPointsPerPlane)
+Flatnesses computeFlatnesses(iterPair vx, iterPair vy, iterPair vz, iterPair phi, Means means, Covariances covariances,
+                             real invNPointsPerPlane, bool sampleScalar)
 {
     Flatnesses flatnesses;
 
     flatnesses.Fx = computeFlatness(vx, means.vx, covariances.vxvx, invNPointsPerPlane);
     flatnesses.Fy = computeFlatness(vy, means.vy, covariances.vyvy, invNPointsPerPlane);
     flatnesses.Fz = computeFlatness(vz, means.vz, covariances.vzvz, invNPointsPerPlane);
-
+    if (sampleScalar)
+        flatnesses.Fphi = computeFlatness(phi, means.phi, covariances.phiphi, invNPointsPerPlane);
     return flatnesses;
 }
 
@@ -414,27 +455,32 @@ std::vector<real> PlanarAverageProbe::computePlaneStatistics(int level)
         getPermutationIterators(parameter->velocityZ, data->indicesOfFirstPlaneD, data->numberOfPointsPerPlane);
     const auto turbulentViscosity =
         getPermutationIterators(parameter->turbViscosity, data->indicesOfFirstPlaneD, data->numberOfPointsPerPlane);
+    const auto phi =
+        getPermutationIterators(parameter->concentration, data->indicesOfFirstPlaneD, data->numberOfPointsPerPlane);
 
     std::vector<real> averages;
 
-    if (!isStatisticIn(PlanarAverageProbe::Statistic::Means, statistics))
+    if (!isStatisticIn(Statistic::Means, statistics))
         return averages;
 
     const real velocityRatio = para->getScaledVelocityRatio(level);
     const real viscosityRatio = para->getScaledViscosityRatio(level);
     const real stressRatio = para->getScaledStressRatio(level);
 
-    const auto means = computeMeans(velocityX, velocityY, velocityZ, invNPointsPerPlane);
+    const auto means = computeMeans(velocityX, velocityY, velocityZ, phi, invNPointsPerPlane, sampleScalar);
     averages.push_back(means.vx * velocityRatio);
     averages.push_back(means.vy * velocityRatio);
     averages.push_back(means.vz * velocityRatio);
     if (para->getUseTurbulentViscosity())
         averages.push_back(computeMean(turbulentViscosity, invNPointsPerPlane) * viscosityRatio);
+    if (sampleScalar)
+        averages.push_back(means.phi);
 
-    if (!isStatisticIn(PlanarAverageProbe::Statistic::Covariances, statistics))
+    if (!isStatisticIn(Statistic::Covariances, statistics))
         return averages;
 
-    const auto covariances = computeCovariances(velocityX, velocityY, velocityZ, means, invNPointsPerPlane);
+    const auto covariances =
+        computeCovariances(velocityX, velocityY, velocityZ, phi, means, invNPointsPerPlane, sampleScalar);
     averages.push_back(covariances.vxvx * stressRatio);
     averages.push_back(covariances.vyvy * stressRatio);
     averages.push_back(covariances.vzvz * stressRatio);
@@ -442,21 +488,34 @@ std::vector<real> PlanarAverageProbe::computePlaneStatistics(int level)
     averages.push_back(covariances.vxvz * stressRatio);
     averages.push_back(covariances.vyvz * stressRatio);
 
-    if (!isStatisticIn(PlanarAverageProbe::Statistic::Skewness, statistics))
+    if (sampleScalar) {
+        averages.push_back(covariances.phiphi);
+        averages.push_back(covariances.vxphi);
+        averages.push_back(covariances.vyphi);
+        averages.push_back(covariances.vzphi);
+    }
+
+    if (!isStatisticIn(Statistic::Skewness, statistics))
         return averages;
 
-    const auto skewnesses = computeSkewnesses(means, covariances, velocityX, velocityY, velocityZ, invNPointsPerPlane);
+    const auto skewnesses =
+        computeSkewnesses(means, covariances, velocityX, velocityY, velocityZ, phi, invNPointsPerPlane, sampleScalar);
     averages.push_back(skewnesses.Sx);
     averages.push_back(skewnesses.Sy);
     averages.push_back(skewnesses.Sz);
+    if (sampleScalar)
+        averages.push_back(skewnesses.Sphi);
 
-    if (!isStatisticIn(PlanarAverageProbe::Statistic::Flatness, statistics))
+    if (!isStatisticIn(Statistic::Flatness, statistics))
         return averages;
 
-    const auto flatnesses = computeFlatnesses(velocityX, velocityY, velocityZ, means, covariances, invNPointsPerPlane);
+    const auto flatnesses =
+        computeFlatnesses(velocityX, velocityY, velocityZ, phi, means, covariances, invNPointsPerPlane, sampleScalar);
     averages.push_back(flatnesses.Fx);
     averages.push_back(flatnesses.Fy);
     averages.push_back(flatnesses.Fz);
+    if (sampleScalar)
+        averages.push_back(flatnesses.Fphi);
 
     return averages;
 }
@@ -495,6 +554,8 @@ void PlanarAverageProbe::calculateQuantities(int level, bool doTimeAverages)
 
         thrust::transform(indices, indices + data->numberOfPointsPerPlane, indices, shiftOp);
     }
+    if (doTimeAverages)
+        data->numberOfTimestepsInTimeAverage++;
     cudaMemoryManager->cudaCopyPlanarAverageProbeIndicesHtoD(this, level);
     getLastCudaError("PlanarAverageProbe::calculateQuantities execution failed");
 }
@@ -502,11 +563,11 @@ void PlanarAverageProbe::calculateQuantities(int level, bool doTimeAverages)
 const uint* PlanarAverageProbe::getNeighborIndicesInPlaneNormal(int level)
 {
     switch (planeNormal) {
-        case PlaneNormal::x:
+        case Axis::x:
             return para->getParD(level)->neighborX;
-        case PlaneNormal::y:
+        case Axis::y:
             return para->getParD(level)->neighborY;
-        case PlaneNormal::z:
+        case Axis::z:
             return para->getParD(level)->neighborZ;
     };
 
diff --git a/src/gpu/core/Samplers/PlanarAverageProbe.h b/src/gpu/core/Samplers/PlanarAverageProbe.h
index e3dfd14b37b8bc75ee6a865d6a67b74363e7cd4a..06a2c9ed3a6275c2b537d253ba29871559449031 100644
--- a/src/gpu/core/Samplers/PlanarAverageProbe.h
+++ b/src/gpu/core/Samplers/PlanarAverageProbe.h
@@ -29,12 +29,8 @@
 //! \addtogroup gpu_Samplers Samplers
 //! \ingroup gpu_core core
 //! \{
-//! \author Henrik Asmuth
+//! \author Henrik Asmuth, Henry Korb
 //! \date 13/05/2022
-//! \brief Probe computing statistics across planes spanning the entire domain
-//!
-
-//!
 //=======================================================================================
 
 #ifndef PlanarAverageProbe_H
@@ -42,77 +38,67 @@
 
 #include "Sampler.h"
 
-#include <stdexcept>
 #include <string>
 #include <vector>
 
 #include <basics/DataTypes.h>
+#include <basics/geometry3d/Axis.h>
 
 class Parameter;
 class CudaMemoryManager;
-struct PlanarAverageProbeLevelData
-{
-    unsigned long long *indicesOfFirstPlaneH, *indicesOfFirstPlaneD;
-    uint numberOfPlanes, numberOfPointsPerPlane, numberOfTimestepsInTimeAverage;
-    std::vector<real> coordinateX, coordinateY, coordinateZ;
-    std::vector<std::vector<real>> instantaneous;
-    std::vector<std::vector<real>> timeAverages;
-};
-
 
 //! \brief Computes spatial statistics across x, y or z-normal planes defined by planeNormal.
 //! The planes include all points of the domain at each respective position along that normal direction.
 //! The spatial statistics can additionally be averaged in time.
+//! The name phi is used to denote the scalar field.
 class PlanarAverageProbe : public Sampler
 {
 public:
-    enum class PlaneNormal { x, y, z };
     enum class Statistic {
         Means,
         Covariances,
         Skewness,
         Flatness,
     };
+    struct LevelData;
 
 public:
-    PlanarAverageProbe(SPtr<Parameter> para, SPtr<CudaMemoryManager> cudaMemoryManager, const std::string outputPath,
-                       const std::string probeName, uint tStartAveraging, uint tStartTemporalAveraging,
-                       uint tBetweenAverages, uint tStartWritingOutput, uint tBetweenWriting,
-                       PlanarAverageProbe::PlaneNormal planeNormal, bool computeTimeAverages)
-        : para(para), cudaMemoryManager(cudaMemoryManager), tStartAveraging(tStartAveraging), tStartTemporalAveraging(tStartTemporalAveraging),
-          tBetweenAverages(tBetweenAverages), tStartWritingOutput(tStartWritingOutput), tBetweenWriting(tBetweenWriting),
-          computeTimeAverages(computeTimeAverages), planeNormal(planeNormal),
-          Sampler(outputPath, probeName)
-    {
-        if (tStartTemporalAveraging < tStartAveraging && computeTimeAverages)
-            throw std::runtime_error("PlaneAverageProbe: tStartTemporalAveraging must be larger than tStartAveraging!");
-        if (tBetweenWriting == 0)
-            throw std::runtime_error("PlaneAverageProbe: tBetweenWriting must be larger than 0!");
-    }
+//! \param tStartSampling The first timestep at which the probe samples planar averages.
+//! \param tStartTemporalAveraging The first timestep at which temporal averaging starts.
+//! \param tBetweenSamples The number of timesteps between samples.
+//! \param tStartWritingOutput The first timestep at which the probe writes output.
+//! \param tBetweenWriting The number of timesteps between writing output.
+//! \param planeNormal The normal direction of the planes along which the probe samples.
+//! \param computeTimeAverages If true, the probe computes time averages.
+//! \param sampleScalar If true, the probe samples statistics related to the scalar.
+    PlanarAverageProbe(SPtr<Parameter> para, SPtr<CudaMemoryManager> cudaMemoryManager, std::string outputPath,
+                       std::string probeName, uint tStartSampling, uint tStartTemporalAveraging, uint tBetweenSamples,
+                       uint tStartWritingOutput, uint tBetweenWriting, Axis planeNormal, bool computeTimeAverages,
+                       bool sampleScalar);
     ~PlanarAverageProbe();
 
     void init() override;
     void sample(int level, uint t) override;
     void getTaggedFluidNodes(GridProvider* gridProvider) override {};
-    PlanarAverageProbeLevelData* getLevelData(int level)
+    LevelData* getLevelData(int level)
     {
         return &levelData[level];
     }
     void addAllAvailableStatistics()
     {
-        addStatistic(PlanarAverageProbe::Statistic::Means);
-        addStatistic(PlanarAverageProbe::Statistic::Covariances);
-        addStatistic(PlanarAverageProbe::Statistic::Skewness);
-        addStatistic(PlanarAverageProbe::Statistic::Flatness);
+        addStatistic(Statistic::Means);
+        addStatistic(Statistic::Covariances);
+        addStatistic(Statistic::Skewness);
+        addStatistic(Statistic::Flatness);
     }
-    void addStatistic(PlanarAverageProbe::Statistic statistic);
+    void addStatistic(Statistic statistic);
     void setFileNameToNOut()
     {
         nameFilesWithFileCount = true;
     }
 
 private:
-    std::vector<std::string> getVariableNames(PlanarAverageProbe::Statistic statistic, bool namesForTimeAverages);
+    std::vector<std::string> getVariableNames(Statistic statistic, bool namesForTimeAverages) const;
     void copyDataToNodedata(std::vector<std::vector<real>>& data, std::vector<std::vector<double>>& nodeData);
     void calculateQuantities(int level, bool doTimeAverages);
     std::vector<unsigned long long> findIndicesInPlane(int level);
@@ -128,15 +114,25 @@ private:
 private:
     SPtr<Parameter> para;
     SPtr<CudaMemoryManager> cudaMemoryManager;
-    uint tStartAveraging, tStartTemporalAveraging, tBetweenAverages, tStartWritingOutput, tBetweenWriting;
-    bool computeTimeAverages, nameFilesWithFileCount = false;
-    PlaneNormal planeNormal;
-    std::vector<PlanarAverageProbe::Statistic> statistics;
-    std::vector<PlanarAverageProbeLevelData> levelData;
+    const uint tStartSampling, tStartTemporalAveraging, tBetweenSamples, tStartWritingOutput, tBetweenWriting;
+    const bool computeTimeAverages, sampleScalar;
+    bool nameFilesWithFileCount = false;
+    const Axis planeNormal;
+    std::vector<Statistic> statistics;
+    std::vector<LevelData> levelData;
     std::vector<std::string> fileNamesForCollectionFile;
 };
 
 bool isStatisticIn(PlanarAverageProbe::Statistic statistic, std::vector<PlanarAverageProbe::Statistic> statistics);
 
+struct PlanarAverageProbe::LevelData
+{
+    unsigned long long *indicesOfFirstPlaneH, *indicesOfFirstPlaneD;
+    uint numberOfPlanes {}, numberOfPointsPerPlane {}, numberOfTimestepsInTimeAverage {};
+    std::vector<real> coordinateX, coordinateY, coordinateZ;
+    std::vector<std::vector<real>> instantaneous;
+    std::vector<std::vector<real>> timeAverages;
+};
+
 #endif
 //! \}
diff --git a/src/gpu/core/Samplers/Probe.cu b/src/gpu/core/Samplers/Probe.cu
index 202549e170315a34d59fb440f4bb8965f6c10a77..6f78a68679d6686d00b453e54e457e575b72243e 100644
--- a/src/gpu/core/Samplers/Probe.cu
+++ b/src/gpu/core/Samplers/Probe.cu
@@ -35,11 +35,17 @@
 #include "Probe.h"
 
 #include <cmath>
+#include <memory>
 #include <stdexcept>
 #include <string>
 
 #include <basics/DataTypes.h>
 #include <basics/constants/NumericConstants.h>
+#include <basics/geometry3d/GbCuboid3D.h>
+#include <basics/geometry3d/GbObject3D.h>
+#include <basics/geometry3d/GbPoint3D.h>
+
+#include <logger/Logger.h>
 
 #include "gpu/core/Calculation/Calculation.h"
 #include "gpu/core/Cuda/CudaMemoryManager.h"
@@ -50,11 +56,51 @@
 #include "gpu/cuda_helper/CudaGrid.h"
 #include "gpu/cuda_helper/CudaIndexCalculation.h"
 
+#include "Sampler.h"
 #include "Utilities.h"
 
 using namespace vf::basics::constant;
 
-__host__ __device__ int calcArrayIndex(int node, int nNodes, int timestep, int nTimesteps, int quantity)
+Probe::Probe(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> cudaMemoryManager, std::string outputPath,
+             std::string probeName, uint tStartSampling, uint tBetweenSamples, uint tStartWritingOutput,
+             uint tBetweenWriting, bool outputTimeSeries, bool sampleEveryTimestep, bool sampleScalar)
+    : para(para), cudaMemoryManager(cudaMemoryManager), tStartSampling(tStartSampling), tBetweenSamples(tBetweenSamples),
+      tStartWritingOutput(tStartWritingOutput), tBetweenWriting(tBetweenWriting), outputTimeSeries(outputTimeSeries),
+      sampleEveryTimestep(sampleEveryTimestep), Sampler(outputPath, probeName)
+{
+    if (tStartWritingOutput < tStartSampling)
+        throw std::runtime_error("Probe: tStartWritingOutput must be larger than tStartSampling!");
+    if (sampleEveryTimestep)
+        VF_LOG_INFO("Probe: sampleEveryTimestep is true, ignoring tBetweenSamples");
+    if (sampleScalar && !para->getDiffOn())
+        throw std::runtime_error("Probe: can only sample scalar if diffusion is enabled in parameter");
+}
+
+void Probe::addProbePlane(real startX, real startY, real startZ, real length, real width, real height)
+{
+    probeObjects.emplace_back(
+        std::make_shared<GbCuboid3D>(startX, startY, startZ, startX + length, startY + width, startZ + height));
+}
+
+void Probe::addProbePoint(real x, real y, real z)
+{
+    probeObjects.emplace_back(std::make_shared<GbPoint3D>(x, y, z));
+}
+
+void Probe::addProbeVolume(std::shared_ptr<GbObject3D> probeVolume)
+{
+    probeObjects.push_back(probeVolume);
+}
+
+void Probe::addProbePointsFromList(std::vector<real> coordsX, std::vector<real> coordY, std::vector<real> coordZ)
+{
+    if (coordsX.size() != coordY.size() || coordsX.size() != coordZ.size())
+        throw std::runtime_error("Probe: Point coordinates must have the same size!");
+    for (uint i = 0; i < coordsX.size(); i++)
+        addProbePoint(coordsX[i], coordY[i], coordZ[i]);
+}
+
+constexpr int calcArrayIndex(int node, int nNodes, int timestep, int nTimesteps, int quantity)
 {
     return node + nNodes * (timestep + nTimesteps * quantity);
 }
@@ -78,26 +124,26 @@ real* getStatisticArray(Probe::ProbeData probeData, Probe::Statistic statistic)
     }
 }
 
-__host__ __device__ real computeMean(real oldMean, real newValue, real inverseCount)
+constexpr real computeMean(real oldMean, real newValue, real inverseCount)
 {
     return oldMean + (newValue - oldMean) * inverseCount;
 }
 
-__host__ __device__ real computeAndSaveMean(real* quantityArray, real oldValue, uint index, real currentValue, real invCount)
+constexpr real computeAndSaveMean(real* quantityArray, real oldValue, uint index, real currentValue, real invCount)
 {
     const real newValue = computeMean(oldValue, currentValue, invCount);
     quantityArray[index] = newValue;
     return newValue;
 }
 
-__host__ __device__ real computeVariance(real oldVariance, real oldMean, real newMean, real currentValue,
-                                         uint numberOfAveragedValues, real inverseCount)
+constexpr real computeVariance(real oldVariance, real oldMean, real newMean, real currentValue, uint numberOfAveragedValues,
+                               real inverseCount)
 {
     return (numberOfAveragedValues * oldVariance + (currentValue - oldMean) * (currentValue - newMean)) * inverseCount;
 }
 
-__host__ __device__ real computeAndSaveVariance(real* quantityArray, real oldVariance, uint indexNew, real currentValue,
-                                                real oldMean, real newMean, uint numberOfAveragedValues, real inverseCount)
+constexpr real computeAndSaveVariance(real* quantityArray, real oldVariance, uint indexNew, real currentValue, real oldMean,
+                                      real newMean, uint numberOfAveragedValues, real inverseCount)
 {
     const real newVariance =
         computeVariance(oldVariance, oldMean, newMean, currentValue, numberOfAveragedValues, inverseCount);
@@ -105,10 +151,9 @@ __host__ __device__ real computeAndSaveVariance(real* quantityArray, real oldVar
     return newVariance;
 }
 
-__forceinline__ __device__ void computeStatistics(uint nAveragedValues, uint currentTimestep, uint lastTimestep,
-                                                  uint nPoints, uint nodeIndex, bool computeInstant, bool computeMean,
-                                                  bool computeVariance, uint iQuantity, real currentValue,
-                                                  const Probe::ProbeData& probeData, real invCount)
+constexpr void computeStatistics(uint nAveragedValues, uint currentTimestep, uint lastTimestep, uint nPoints, uint nodeIndex,
+                                 bool computeInstant, bool computeMean, bool computeVariance, uint iQuantity,
+                                 real currentValue, const Probe::ProbeData& probeData, real invCount)
 {
     const uint indexCurrent = calcArrayIndex(nodeIndex, nPoints, currentTimestep, probeData.numberOfTimesteps, iQuantity);
     const uint indexLast = calcArrayIndex(nodeIndex, nPoints, lastTimestep, probeData.numberOfTimesteps, iQuantity);
@@ -152,6 +197,10 @@ __global__ void calculateQuantitiesKernel(uint numberOfAveragedValues, Probe::Gr
     computeStatistics(numberOfAveragedValues, currentTimestep, lastTimestep, nPoints, nodeIndex,
                       probeData.computeInstantaneous, probeData.computeMeans, probeData.computeVariances, 3,
                       gridParams.density[gridNodeIndex], probeData, invCount);
+    if (probeData.sampleScalar)
+        computeStatistics(numberOfAveragedValues, currentTimestep, lastTimestep, nPoints, nodeIndex,
+                          probeData.computeInstantaneous, probeData.computeMeans, probeData.computeVariances, 4,
+                          gridParams.scalar[gridNodeIndex], probeData, invCount);
 }
 
 std::vector<Probe::PostProcessingVariable> Probe::getPostProcessingVariables(Statistic statistic, int level) const
@@ -166,18 +215,24 @@ std::vector<Probe::PostProcessingVariable> Probe::getPostProcessingVariables(Sta
             postProcessingVariables.emplace_back("vy", velocityRatio);
             postProcessingVariables.emplace_back("vz", velocityRatio);
             postProcessingVariables.emplace_back("rho", densityRatio);
+            if (sampleScalar)
+                postProcessingVariables.emplace_back("phi", c1o1);
             break;
         case Statistic::Means:
             postProcessingVariables.emplace_back("vx_mean", velocityRatio);
             postProcessingVariables.emplace_back("vy_mean", velocityRatio);
             postProcessingVariables.emplace_back("vz_mean", velocityRatio);
             postProcessingVariables.emplace_back("rho_mean", densityRatio);
+            if (sampleScalar)
+                postProcessingVariables.emplace_back("phi_mean", c1o1);
             break;
         case Statistic::Variances:
             postProcessingVariables.emplace_back("vx_var", stressRatio);
             postProcessingVariables.emplace_back("vy_var", stressRatio);
             postProcessingVariables.emplace_back("vz_var", stressRatio);
             postProcessingVariables.emplace_back("rho_var", densityRatio);
+            if (sampleScalar)
+                postProcessingVariables.emplace_back("phi_var", c1o1);
             break;
 
         default:
@@ -232,33 +287,20 @@ void Probe::addLevelData(int level)
     const real* nodeCoordinatesZ = para->getParH(level)->coordinateZ;
     const real deltaX = nodeCoordinatesX[para->getParH(level)->neighborX[1]] - nodeCoordinatesX[1];
     for (unsigned long long pos = 1; pos < para->getParH(level)->numberOfNodes; pos++) {
-        const real pointCoordX = nodeCoordinatesX[pos];
-        const real pointCoordY = nodeCoordinatesY[pos];
-        const real pointCoordZ = nodeCoordinatesZ[pos];
-        for (auto point : points) {
-            const real distX = point.x - pointCoordX;
-            const real distY = point.y - pointCoordY;
-            const real distZ = point.z - pointCoordZ;
-            if (distX <= deltaX && distY <= deltaX && distZ <= deltaX && distX > c0o1 && distY > c0o1 && distZ > c0o1 &&
+        const real nodeCoordX = nodeCoordinatesX[pos];
+        const real nodeCoordY = nodeCoordinatesY[pos];
+        const real nodeCoordZ = nodeCoordinatesZ[pos];
+        const real maxX = nodeCoordX + deltaX;
+        const real maxY = nodeCoordY + deltaX;
+        const real maxZ = nodeCoordZ + deltaX;
+        for (auto object : probeObjects) {
+            if ((object->isInsideCell(nodeCoordX, nodeCoordY, nodeCoordZ, maxX, maxY, maxZ) ||
+                 object->isPointInGbObject3D(nodeCoordX, nodeCoordY, nodeCoordZ)) &&
                 isValidProbePoint(pos, para.get(), level)) {
                 indices.push_back(static_cast<uint>(pos));
-                coordinatesX.push_back(pointCoordX);
-                coordinatesY.push_back(pointCoordY);
-                coordinatesZ.push_back(pointCoordZ);
-                continue;
-            }
-        }
-        for (auto plane : planes) {
-            const real distanceX = pointCoordX - plane.startX;
-            const real distanceY = pointCoordY - plane.startY;
-            const real distanceZ = pointCoordZ - plane.startZ;
-
-            if (distanceX <= plane.length && distanceY <= plane.width && distanceZ <= plane.height && distanceX >= c0o1 &&
-                distanceY >= c0o1 && distanceZ >= c0o1 && isValidProbePoint(pos, para.get(), level)) {
-                indices.push_back(static_cast<uint>(pos));
-                coordinatesX.push_back(pointCoordX);
-                coordinatesY.push_back(pointCoordY);
-                coordinatesZ.push_back(pointCoordZ);
+                coordinatesX.push_back(nodeCoordX);
+                coordinatesY.push_back(nodeCoordY);
+                coordinatesZ.push_back(nodeCoordZ);
                 continue;
             }
         }
@@ -266,7 +308,7 @@ void Probe::addLevelData(int level)
 
     const uint numberOfQuantities = static_cast<uint>(getPostProcessingVariables(Statistic::Instantaneous, 0).size());
 
-    ProbeData probeData(enableComputationInstantaneous, enableComputationMeans, enableComputationVariances,
+    ProbeData probeData(enableComputationInstantaneous, enableComputationMeans, enableComputationVariances, sampleScalar,
                         static_cast<uint>(indices.size()), numberOfQuantities, getNumberOfTimestepsInTimeseries(level));
 
     const uint sizeData = probeData.numberOfPoints * probeData.numberOfTimesteps * numberOfQuantities;
@@ -294,15 +336,15 @@ void Probe::sample(int level, uint t)
         return;
     const uint tLevel = para->getTimeStep(level, t, false);
 
-    //! if averageEveryTimestep the probe will be evaluated in every sub-timestep of each respective level
-    //! else, the probe will only be evaluated in each synchronous time step tBetweenAverages
+    //! if sampleEveryTimestep the probe will be evaluated in every sub-timestep of each respective level
+    //! else, the probe will only be evaluated in each synchronous time step tBetweenSamples
 
     const uint levelFactor = exp2(level);
 
-    const uint tStartAvgLevel = this->tStartAveraging * levelFactor;
+    const uint tStartAvgLevel = this->tStartSampling * levelFactor;
     const uint tAfterStartAvg = tLevel - tStartAvgLevel;
-    const uint tAvgLevel = this->tBetweenAverages * levelFactor;
-    const bool averageThisTimestep = this->averageEveryTimestep || (tAfterStartAvg % tAvgLevel == 0);
+    const uint tAvgLevel = this->tBetweenSamples * levelFactor;
+    const bool averageThisTimestep = this->sampleEveryTimestep || (tAfterStartAvg % tAvgLevel == 0);
 
     const uint tStartOutLevel = this->tStartWritingOutput * levelFactor;
     const uint tAfterStartOut = tLevel - tStartOutLevel;
@@ -313,7 +355,7 @@ void Probe::sample(int level, uint t)
 
     const vf::cuda::CudaGrid grid(para->getParD(level)->numberofthreads, levelData->probeDataD.numberOfPoints);
 
-    if ((t > this->tStartAveraging) && averageThisTimestep) {
+    if ((t > this->tStartSampling) && averageThisTimestep) {
         if (outputTimeSeries) {
             const uint lastTimestep = calcOldTimestep(levelData->timeseriesParams.currentTimestep,
                                                       levelData->timeseriesParams.lastTimestepInOldTimeseries);
@@ -480,7 +522,7 @@ std::vector<real> Probe::getTimestepData(real time, int timestep, int level)
 
 void Probe::appendTimeseriesFile(int level, int t)
 {
-    const uint tAvg_level = this->tBetweenAverages == 1 ? this->tBetweenAverages : this->tBetweenAverages * exp2(-level);
+    const uint tAvg_level = this->tBetweenSamples == 1 ? this->tBetweenSamples : this->tBetweenSamples * exp2(-level);
     const real deltaT = para->getTimeRatio() * tAvg_level;
     auto levelData = levelDatas[level];
 
@@ -516,12 +558,7 @@ void Probe::getTaggedFluidNodes(GridProvider* gridProvider)
 
 Probe::GridParams Probe::getGridParams(LBMSimulationParameter* para)
 {
-    return {
-        para->velocityX,
-        para->velocityY,
-        para->velocityZ,
-        para->rho,
-    };
+    return { para->velocityX, para->velocityY, para->velocityZ, para->rho, para->concentration };
 }
 
 bool isCoarseInterpolationCell(unsigned long long pointIndex, Parameter* para, int level)
diff --git a/src/gpu/core/Samplers/Probe.h b/src/gpu/core/Samplers/Probe.h
index 67c59af8d8f9c9a722ed548e48b70e44469d4ff4..ac6544377e387ba18924bc53b5dd421f2f77402d 100644
--- a/src/gpu/core/Samplers/Probe.h
+++ b/src/gpu/core/Samplers/Probe.h
@@ -31,7 +31,6 @@
 //! \{
 //! \author Henry Korb, Henrik Asmuth
 //! \date 13/05/2022
-
 //=======================================================================================
 
 #ifndef Probe_H
@@ -39,15 +38,14 @@
 
 #include "Samplers/Sampler.h"
 
-#include <cuda.h>
 #include <functional>
 #include <stdexcept>
 #include <string>
 
 #include <basics/DataTypes.h>
 #include <basics/PointerDefinitions.h>
+#include <basics/geometry3d/GbObject3D.h>
 #include <basics/writer/WbWriterVtkXmlBinary.h>
-#include <logger/Logger.h>
 
 struct LBMSimulationParameter;
 class Parameter;
@@ -55,128 +53,49 @@ class CudaMemoryManager;
 
 //! \brief Computes statistics of pointwise data. Data can be written to vtk-file or timeseries file.
 //! All points and planes are written to the same file. Use different probes to write to separate files.
-//! Data is sampled in synchronous timestep, unless averageEveryTimestep is set to true.
+//! Data is sampled in synchronous timestep, unless sampleEveryTimestep is set to true.
 class Probe : public Sampler
 {
 public:
     enum class Statistic { Instantaneous, Means, Variances };
-
-    Probe(SPtr<Parameter> para, SPtr<CudaMemoryManager> cudaMemoryManager, std::string outputPath, std::string probeName,
-          uint tStartAveraging, uint tBetweenAverages, uint tStartWritingOutput, uint tBetweenWriting, bool outputTimeSeries,
-          bool averageEveryTimestep)
-        : para(para), cudaMemoryManager(cudaMemoryManager), tStartAveraging(tStartAveraging),
-          tBetweenAverages(tBetweenAverages), tStartWritingOutput(tStartWritingOutput), tBetweenWriting(tBetweenWriting),
-          outputTimeSeries(outputTimeSeries), averageEveryTimestep(averageEveryTimestep), Sampler(outputPath, probeName)
-    {
-        if (tStartWritingOutput < tStartAveraging)
-            throw std::runtime_error("Probe: tStartWritingOutput must be larger than tStartAveraging!");
-        if (averageEveryTimestep)
-            VF_LOG_INFO("Probe: averageEveryTimestep is true, ignoring tBetweenAverages");
-    }
-
+    Probe(std::shared_ptr<Parameter> para, std::shared_ptr<CudaMemoryManager> cudaMemoryManager, std::string outputPath,
+          std::string probeName, uint tStartSampling, uint tBetweenSamples, uint tStartWritingOutput, uint tBetweenWriting,
+          bool outputTimeSeries, bool sampleEveryTimestep, bool sampleScalar = false);
     ~Probe();
 
-    void init() override;
-    void sample(int level, uint t) override;
-    void getTaggedFluidNodes(GridProvider* gridProvider) override;
-
-    void addProbePlane(real startX, real startY, real startZ, real length, real width, real height)
-    {
-        planes.emplace_back(Plane { startX, startY, startZ, length, width, height });
-    }
-
-    void addProbePoint(real x, real y, real z)
-    {
-        points.emplace_back(Point { x, y, z });
-    }
-
-    void addProbePointsFromList(std::vector<real> coordsX, std::vector<real> coordY, std::vector<real> coordZ)
-    {
-        if (coordsX.size() != coordY.size() || coordsX.size() != coordZ.size())
-            throw std::runtime_error("Probe: Point coordinates must have the same size!");
-        for (uint i = 0; i < coordsX.size(); i++)
-            addProbePoint(coordsX[i], coordY[i], coordZ[i]);
-    }
-
-    void addStatistic(Probe::Statistic variable);
+    void addProbePlane(real startX, real startY, real startZ, real length, real width, real height);
+    void addProbePoint(real x, real y, real z);
+    void addProbePointsFromList(std::vector<real> coordsX, std::vector<real> coordY, std::vector<real> coordZ);
+    void addProbeVolume(std::shared_ptr<GbObject3D> probeVolume);
+    void addStatistic(Statistic variable);
     void addAllAvailableStatistics()
     {
-        addStatistic(Probe::Statistic::Instantaneous);
-        addStatistic(Probe::Statistic::Means);
-        addStatistic(Probe::Statistic::Variances);
+        addStatistic(Statistic::Instantaneous);
+        addStatistic(Statistic::Means);
+        addStatistic(Statistic::Variances);
     }
 
     void setFileNameToNOut()
     {
         this->fileNameLU = false;
     }
+    void init() override;
+    void sample(int level, uint t) override;
+    void getTaggedFluidNodes(GridProvider* gridProvider) override;
 
-    struct PostProcessingVariable
-    {
-        std::string name;
-        real conversionFactor;
-        PostProcessingVariable(std::string name, real conversionFactor) : name(name), conversionFactor(conversionFactor)
-        {
-        }
-    };
-
-    struct ProbeData
-    {
-        real *instantaneous, *means, *variances;
-        bool computeInstantaneous, computeMeans, computeVariances;
-        uint numberOfPoints, numberOfQuantities, numberOfTimesteps;
-        uint* indices;
-        __device__ __host__ ProbeData(bool computeInstantaneous, bool computeMeans, bool computeVariance,
-                                      uint numberOfPoints, uint numberOfQuantities, uint numberOfTimesteps)
-            : computeInstantaneous(computeInstantaneous), computeMeans(computeMeans), computeVariances(computeVariance),
-              numberOfPoints(numberOfPoints), numberOfQuantities(numberOfQuantities), numberOfTimesteps(numberOfTimesteps)
-        {
-        }
-    };
-
-    struct TimeseriesParams
-    {
-        uint currentTimestep {}, numberOfTimesteps {}, lastTimestepInOldTimeseries {};
-    };
-
-    struct LevelData
-    {
-        ProbeData probeDataH, probeDataD;
-        TimeseriesParams timeseriesParams;
-        std::vector<real> coordinatesX, coordinatesY, coordinatesZ;
-        uint numberOfAveragedValues {};
-        LevelData(ProbeData probeDataH, ProbeData probeDataD, std::vector<real> coordinatesX, std::vector<real> coordinatesY,
-                  std::vector<real> coordinatesZ)
-            : probeDataH(probeDataH), probeDataD(probeDataD), coordinatesX(coordinatesX), coordinatesY(coordinatesY),
-              coordinatesZ(coordinatesZ)
-        {
-        }
-    };
+    struct PostProcessingVariable;
+    struct ProbeData;
+    struct TimeseriesParams;
+    struct LevelData;
 
     LevelData* getLevelData(int level)
     {
         return &levelDatas[level];
     }
 
-    struct GridParams
-    {
-        real *velocityX, *velocityY, *velocityZ, *density;
-    };
-
+    struct GridParams;
     static GridParams getGridParams(LBMSimulationParameter* para);
 
-    struct Plane
-    {
-        real startX, startY, startZ;
-        real length, width, height;
-    };
-
-    struct Point
-    {
-        real x, y, z;
-    };
-
-
 private:
     void addLevelData(int level);
     void addPointsToLevelData(LevelData& levelData, std::vector<uint>& indices, int level);
@@ -211,21 +130,68 @@ protected:
     SPtr<Parameter> para;
     SPtr<CudaMemoryManager> cudaMemoryManager;
     std::vector<LevelData> levelDatas;
-    bool enableComputationInstantaneous {}, enableComputationMeans {}, enableComputationVariances {};
+    bool enableComputationInstantaneous {}, enableComputationMeans {}, enableComputationVariances {}, sampleScalar {};
     //! flag initiating time series output.
-    const bool outputTimeSeries, averageEveryTimestep;
+    const bool outputTimeSeries, sampleEveryTimestep;
     std::vector<std::string> fileNamesForCollectionFile;
     std::vector<std::string> timeseriesFileNames;
 
     //! if true, written file name contains time step in LU, else is the number of the written probe files
     bool fileNameLU = true;
 
-    const uint tStartAveraging;
-    const uint tBetweenAverages;
+    const uint tStartSampling;
+    const uint tBetweenSamples;
     const uint tStartWritingOutput;
     const uint tBetweenWriting;
-    std::vector<Plane> planes;
-    std::vector<Point> points;
+    std::vector<std::shared_ptr<GbObject3D>> probeObjects;
+};
+
+struct Probe::PostProcessingVariable
+{
+    std::string name;
+    real conversionFactor;
+    PostProcessingVariable(std::string name, real conversionFactor) : name(name), conversionFactor(conversionFactor)
+    {
+    }
+};
+
+struct Probe::ProbeData
+{
+    real *instantaneous, *means, *variances;
+    bool computeInstantaneous, computeMeans, computeVariances, sampleScalar;
+    uint numberOfPoints, numberOfQuantities, numberOfTimesteps;
+    uint* indices;
+    __device__ __host__ ProbeData(bool computeInstantaneous, bool computeMeans, bool computeVariance, bool sampleScalar,
+                                  uint numberOfPoints, uint numberOfQuantities, uint numberOfTimesteps)
+        : computeInstantaneous(computeInstantaneous), computeMeans(computeMeans), computeVariances(computeVariance),
+          sampleScalar(sampleScalar), numberOfPoints(numberOfPoints), numberOfQuantities(numberOfQuantities),
+          numberOfTimesteps(numberOfTimesteps)
+    {
+    }
+};
+
+struct Probe::TimeseriesParams
+{
+    uint currentTimestep {}, numberOfTimesteps {}, lastTimestepInOldTimeseries {};
+};
+
+struct Probe::LevelData
+{
+    ProbeData probeDataH, probeDataD;
+    TimeseriesParams timeseriesParams;
+    std::vector<real> coordinatesX, coordinatesY, coordinatesZ;
+    uint numberOfAveragedValues {};
+    LevelData(ProbeData probeDataH, ProbeData probeDataD, std::vector<real> coordinatesX, std::vector<real> coordinatesY,
+              std::vector<real> coordinatesZ)
+        : probeDataH(probeDataH), probeDataD(probeDataD), coordinatesX(coordinatesX), coordinatesY(coordinatesY),
+          coordinatesZ(coordinatesZ)
+    {
+    }
+};
+
+struct Probe::GridParams
+{
+    real *velocityX, *velocityY, *velocityZ, *density, *scalar;
 };
 
 bool isValidProbePoint(unsigned long long pointIndex, Parameter* para, int level);
diff --git a/src/gpu/core/Samplers/Utilities.h b/src/gpu/core/Samplers/Utilities.h
index 0ac991ce9c8e75ffe4c0458c59cf91392d6e2ed3..76b1b5599a1f9431902c1240ce89d2c7ba3915b5 100644
--- a/src/gpu/core/Samplers/Utilities.h
+++ b/src/gpu/core/Samplers/Utilities.h
@@ -63,7 +63,7 @@ inline std::string makeTimeseriesFileName(const std::string& probeName, int leve
 }
 
 template <typename T>
-__host__ __device__ inline T computeNewTimeAverage(T oldAverage, T newValue, real inverseNumberOfTimesteps)
+constexpr inline T computeNewTimeAverage(T oldAverage, T newValue, real inverseNumberOfTimesteps)
 {
     return oldAverage + (newValue - oldAverage) * inverseNumberOfTimesteps;
 }
diff --git a/src/gpu/core/Samplers/WallModelProbe.cu b/src/gpu/core/Samplers/WallModelProbe.cu
index 8d8c770d5bed5114e10b3c886a092041ea23db48..c698d54868cb6abdd23bc527bd52221a099595d6 100644
--- a/src/gpu/core/Samplers/WallModelProbe.cu
+++ b/src/gpu/core/Samplers/WallModelProbe.cu
@@ -31,7 +31,9 @@
 //! \{
 #include "WallModelProbe.h"
 
+#include <cstddef>
 #include <functional>
+#include <string>
 #include <vector>
 
 #include <thrust/device_ptr.h>
@@ -95,7 +97,7 @@ std::vector<std::string> WallModelProbe::getVariableNames()
 void WallModelProbe::init()
 {
 
-    std::vector<std::string> variableNames  = getVariableNames();
+    std::vector<std::string> variableNames = getVariableNames();
 
     const int numberOfQuantities = getNumberOfInstantaneousQuantities();
 
@@ -137,18 +139,20 @@ template <typename T>
 T computeMean(T* devicePointer, uint numberOfPoints, T conversionFactor)
 {
     thrust::device_ptr<T> thrustPointer = thrust::device_pointer_cast(devicePointer);
-    return  thrust::reduce(thrustPointer, thrustPointer + numberOfPoints) / real(numberOfPoints) * conversionFactor;
+    return thrust::reduce(thrustPointer, thrustPointer + numberOfPoints) / real(numberOfPoints) * conversionFactor;
 }
 
-struct isValidNode {
-    __host__ __device__ real operator()(thrust::tuple<real, uint> x)
+struct isValidNode
+{
+    constexpr real operator()(thrust::tuple<real, uint> x)
     {
         return thrust::get<1>(x) == GEO_FLUID ? thrust::get<0>(x) : c0o1;
     }
 };
 
 template <typename T>
-T computeIndexBasedMean(T* devicePointer, uint* typeOfGridNode, uint numberOfNodes, uint numberOfFluidNodes, T conversionFactor)
+T computeIndexBasedMean(T* devicePointer, uint* typeOfGridNode, uint numberOfNodes, uint numberOfFluidNodes,
+                        T conversionFactor)
 {
     thrust::device_ptr<T> thrustPointer = thrust::device_pointer_cast(devicePointer);
     thrust::device_ptr<uint> typePointer = thrust::device_pointer_cast(typeOfGridNode);
@@ -170,7 +174,8 @@ template <typename T>
 void computeAndSaveIndexBasedMean(T* devicePointer, uint* typeOfGridNode, uint numberOfNodes, uint numberOfFluidNodes,
                                   std::vector<real>& quantitiesArray, T conversionFactor)
 {
-    quantitiesArray.push_back(computeIndexBasedMean(devicePointer, typeOfGridNode, numberOfNodes, numberOfFluidNodes, conversionFactor));
+    quantitiesArray.push_back(
+        computeIndexBasedMean(devicePointer, typeOfGridNode, numberOfNodes, numberOfFluidNodes, conversionFactor));
 }
 
 template <typename T>
@@ -185,7 +190,7 @@ uint WallModelProbe::countFluidNodes(int level)
     return std::count(typePointer, typePointer + para->getParH(level)->numberOfNodes, GEO_FLUID);
 }
 
-void WallModelProbe::calculateQuantities(WallModelProbeLevelData* data, uint t, int level)
+void WallModelProbe::calculateQuantities(WallModelProbe::LevelData* data, uint t, int level)
 {
     const uint nPoints = para->getParD(level)->stressBC.numberOfBCnodes;
     if (nPoints < 1)
@@ -227,12 +232,12 @@ void WallModelProbe::calculateQuantities(WallModelProbeLevelData* data, uint t,
     }
 
     if (computeTemporalAverages) {
-        if(t > tStartTemporalAveraging){
+        if (t > tStartTemporalAveraging) {
             const real inverseNumberOfAveragedValues = c1o1 / real(data->numberOfAveragedValues + 1);
             std::vector<real>& oldAverages = data->averagedData.back();
             std::vector<real> newAverages;
             newAverages.reserve(numberOfQuantities);
-            for(int i=0; i<numberOfQuantities; i++)
+            for (int i = 0; i < numberOfQuantities; i++)
                 computeTemporalAverage(newAverages, oldAverages[i], newInstantaneous[i], inverseNumberOfAveragedValues);
             data->averagedData.push_back(newAverages);
             data->numberOfAveragedValues++;
@@ -247,18 +252,18 @@ void WallModelProbe::write(int level)
 {
     auto data = &levelData[level];
 
-    if(data->timestepTime.empty())
+    if (data->timestepTime.empty())
         return;
 
     std::vector<std::vector<real>> dataToWrite;
-    for(size_t i=0; i<data->timestepTime.size(); i++)
-    {
+    for (size_t i = 0; i < data->timestepTime.size(); i++) {
         std::vector<real> row;
-        row.reserve(data->instantaneousData[i].size() + (computeTemporalAverages ? data->averagedData[i+1].size() : 0) + 1);
+        row.reserve(data->instantaneousData[i].size() + (computeTemporalAverages ? data->averagedData[i + 1].size() : 0) +
+                    1);
         row.push_back(data->timestepTime[i]);
         std::copy(data->instantaneousData[i].begin(), data->instantaneousData[i].end(), std::back_inserter(row));
-        if(computeTemporalAverages)
-            std::copy(data->averagedData[i+1].begin(), data->averagedData[i+1].end(), std::back_inserter(row));
+        if (computeTemporalAverages)
+            std::copy(data->averagedData[i + 1].begin(), data->averagedData[i + 1].end(), std::back_inserter(row));
         dataToWrite.push_back(row);
     }
 
@@ -266,12 +271,11 @@ void WallModelProbe::write(int level)
 
     data->timestepTime.clear();
     data->instantaneousData.clear();
-    if(computeTemporalAverages){
+    if (computeTemporalAverages) {
         auto lastTimestep = data->averagedData.back();
         data->averagedData.clear();
         data->averagedData.push_back(lastTimestep);
     }
 }
 
-
 //! \}
diff --git a/src/gpu/core/Samplers/WallModelProbe.h b/src/gpu/core/Samplers/WallModelProbe.h
index f5a7b5873b8940e24a68729dfa5c3b082def8de3..5bafda67226a63c2c286feecafd9fa9ea734c32c 100644
--- a/src/gpu/core/Samplers/WallModelProbe.h
+++ b/src/gpu/core/Samplers/WallModelProbe.h
@@ -29,10 +29,8 @@
 //! \addtogroup gpu_Samplers Samplers
 //! \ingroup gpu_core core
 //! \{
-//! \author Henrik Asmuth
+//! \author Henrik Asmuth, Henry Korb
 //! \date 13/05/2022
-//!
-//!
 //=======================================================================================
 
 #ifndef WallModelProbe_H
@@ -51,19 +49,6 @@
 class Parameter;
 class CudaMemoryManager;
 
-struct WallModelProbeLevelData
-{
-    uint numberOfAveragedValues, numberOfFluidNodes;
-    std::string timeseriesFileName;
-    std::vector<std::vector<real>> instantaneousData, averagedData;
-    std::vector<real> timestepTime;
-    bool firstWrite = true;
-    WallModelProbeLevelData(std::string fileName, uint numberOfFluidNodes)
-        : timeseriesFileName(fileName), numberOfFluidNodes(numberOfFluidNodes)
-    {
-    }
-};
-
 //! \brief Probe computing statistics of all relevant wall model quantities used in the StressBC kernels
 //! Computes spatial statistics for all grid points of the StressBC
 //! The spatial statistics can additionally be averaged in time.
@@ -94,13 +79,15 @@ public:
     void sample(int level, uint t) override;
     void getTaggedFluidNodes(GridProvider* gridProvider) override {};
 
+    struct LevelData;
+
 private:
     std::vector<std::string> getVariableNames();
-    int getNumberOfInstantaneousQuantities()
+    int getNumberOfInstantaneousQuantities() const
     {
         return evaluatePressureGradient ? 13 : 10;
     }
-    void calculateQuantities(WallModelProbeLevelData* levelData, uint t, int level);
+    void calculateQuantities(LevelData* levelData, uint t, int level);
     void write(int level);
     uint countFluidNodes(int level);
 
@@ -112,7 +99,20 @@ private:
     bool evaluatePressureGradient = false;
     bool computeTemporalAverages = false;
     bool averageEveryTimestep = false;
-    std::vector<WallModelProbeLevelData> levelData;
+    std::vector<LevelData> levelData;
+};
+
+struct WallModelProbe::LevelData
+{
+    uint numberOfAveragedValues {}, numberOfFluidNodes {};
+    std::string timeseriesFileName;
+    std::vector<std::vector<real>> instantaneousData, averagedData;
+    std::vector<real> timestepTime;
+    bool firstWrite = true;
+    LevelData(std::string fileName, uint numberOfFluidNodes)
+        : timeseriesFileName(fileName), numberOfFluidNodes(numberOfFluidNodes)
+    {
+    }
 };
 
 #endif