diff --git a/Python/actuator_line/actuator_line.py b/Python/actuator_line/actuator_line.py
index cd0e06149a52430cbdc811106b505bc6049e2221..6e3c8608617df1267535984d53307dea9184c6ab 100644
--- a/Python/actuator_line/actuator_line.py
+++ b/Python/actuator_line/actuator_line.py
@@ -16,8 +16,8 @@ sim_name = "ActuatorLine"
 config_file = Path(__file__).parent/Path("config.txt")
 output_path = Path(__file__).parent/Path("output")
 output_path.mkdir(exist_ok=True)
-timeStepOut = 500
-t_end = 50
+t_out = 100.
+t_end = 500.
 
 #%%
 logger.Logger.initialize_logger()
@@ -25,15 +25,18 @@ basics.logger.Logger.add_stdout()
 basics.logger.Logger.set_debug_level(basics.logger.Level.INFO_LOW)
 basics.logger.Logger.time_stamp(basics.logger.TimeStamp.ENABLE)
 basics.logger.Logger.enable_printed_rank_numbers(True)
+# %%
+comm = gpu.Communicator.get_instance()
+#%%
+grid_factory = gpu.grid_generator.GridFactory.make()
+grid_builder = gpu.grid_generator.MultipleGridBuilder.make_shared(grid_factory)
+
 #%%
-grid_builder = gpu.MultipleGridBuilder.make_shared()
 dx = reference_diameter/nodes_per_diameter
 
 grid_builder.add_coarse_grid(0.0, 0.0, 0.0, *length, dx)
 grid_builder.set_periodic_boundary_condition(False, False, False)
 grid_builder.build_grids(basics.LbmOrGks.LBM, False)
-# %%
-comm = gpu.Communicator.get_instance()
 #%%
 config = basics.ConfigurationFile()
 config.load(str(config_file))
@@ -55,6 +58,7 @@ para.set_max_level(1)
 para.set_velocity(velocity_lb)
 para.set_viscosity(viscosity_lb)    
 para.set_velocity_ratio(dx/dt)
+para.set_viscosity_ratio(dx*dx/dt)
 para.set_main_kernel("TurbulentViscosityCumulantK17CompChim")
 para.set_use_AMD(True)
 para.set_SGS_constant(0.083)
@@ -63,7 +67,7 @@ def init_func(coord_x, coord_y, coord_z):
     return [0.0, velocity_lb, 0.0, 0.0]
 
 para.set_initial_condition(init_func)
-para.set_t_out(timeStepOut)
+para.set_t_out(int(t_out/dt))
 para.set_t_end(int(t_end/dt))
 para.set_is_body_force(True)
 
@@ -79,8 +83,8 @@ grid_builder.set_velocity_boundary_condition(gpu.SideType.PZ, velocity_lb, 0.0,
 grid_builder.set_pressure_boundary_condition(gpu.SideType.PX, 0.0)
 
 #%%
-cuda_memory_manager = gpu.CudaMemoryManager.make(para)
-grid_generator = gpu.GridProvider.make_grid_generator(grid_builder, para, cuda_memory_manager)
+cuda_memory_manager = gpu.CudaMemoryManager(para)
+grid_generator = gpu.GridProvider.make_grid_generator(grid_builder, para, cuda_memory_manager, comm)
 #%%
 turb_pos = np.array([3,3,3])*reference_diameter
 epsilon = 5
@@ -91,21 +95,17 @@ n_blade_nodes = 32
 alm = gpu.ActuatorLine(n_blades, density, n_blade_nodes, epsilon, *turb_pos, reference_diameter, level, dt, dx)
 para.add_actuator(alm)
 #%%
-point_probe = gpu.probes.PointProbe("pointProbe", str(output_path), 100, 500, 100)
+point_probe = gpu.probes.PointProbe("pointProbe", str(output_path), 100, 1, 500, 100)
 point_probe.add_probe_points_from_list(np.array([1,2,5])*reference_diameter, np.array([3,3,3])*reference_diameter, np.array([3,3,3])*reference_diameter)
-point_probe.add_post_processing_variable(gpu.probes.PostProcessingVariable.Means)
+point_probe.add_statistic(gpu.probes.Statistic.Means)
 
 para.add_probe(point_probe)
 
-plane_probe = gpu.probes.PlaneProbe("planeProbe", str(output_path), 100, 500, 100)
+plane_probe = gpu.probes.PlaneProbe("planeProbe", str(output_path), 100, 1, 500, 100)
 plane_probe.set_probe_plane(5*reference_diameter, 0, 0, dx, length[1], length[2])
 para.add_probe(plane_probe)
 #%%
-sim = gpu.Simulation(comm)
-kernel_factory = gpu.KernelFactory.get_instance()
-sim.set_factories(kernel_factory, gpu.PreProcessorFactory.get_instance())
-sim.init(para, grid_generator, gpu.FileWriter(), cuda_memory_manager)
+sim = gpu.Simulation(para, cuda_memory_manager, comm, grid_generator)
 #%%
 sim.run()
-sim.free()
 MPI.Finalize()
\ No newline at end of file
diff --git a/Python/boundary_layer/boundary_layer.py b/Python/boundary_layer/boundary_layer.py
index cf941a9418e5c3ec5d94864f119de20401601622..1c01f50946b49bc0ddab7e50065a24aab4ae869f 100644
--- a/Python/boundary_layer/boundary_layer.py
+++ b/Python/boundary_layer/boundary_layer.py
@@ -4,20 +4,34 @@ from pathlib import Path
 from mpi4py import MPI
 from pyfluids import basics, gpu, logger
 #%%
-reference_diameter = 126
+reference_height = 1000 # boundary layer height in m
 
-length = np.array([30,8,8])*reference_diameter
+length = np.array([6,4,1])*reference_height
 viscosity = 1.56e-5
-velocity = 9
 mach = 0.1
-nodes_per_diameter = 32
+nodes_per_height = 32
+
+z_0 = 0.1
+u_star = 0.4
+kappa = 0.4
+
+velocity = 0.5*u_star/kappa*np.log(length[2]/z_0+1)
+flow_through_time = length[0]/velocity
+use_AMD = True
+
 
 sim_name = "BoundaryLayer"
 config_file = Path(__file__).parent/Path("config.txt")
 output_path = Path(__file__).parent/Path("output")
 output_path.mkdir(exist_ok=True)
-timeStepOut = 500
-t_end = 50
+t_out = 1000.
+t_end = 5000.
+
+t_start_averaging = 0
+t_start_tmp_averaging =  100_000
+t_averaging = 200
+t_start_out_probe = 0
+t_out_probe = 1000
 
 #%%
 logger.Logger.initialize_logger()
@@ -25,24 +39,37 @@ basics.logger.Logger.add_stdout()
 basics.logger.Logger.set_debug_level(basics.logger.Level.INFO_LOW)
 basics.logger.Logger.time_stamp(basics.logger.TimeStamp.ENABLE)
 basics.logger.Logger.enable_printed_rank_numbers(True)
-#%%
-grid_builder = gpu.MultipleGridBuilder.make_shared()
-dx = reference_diameter/nodes_per_diameter
-
-grid_builder.add_coarse_grid(0.0, 0.0, 0.0, *length, dx)
-grid_builder.set_periodic_boundary_condition(False, False, False)
-grid_builder.build_grids(basics.LbmOrGks.LBM, False)
 # %%
 comm = gpu.Communicator.get_instance()
+#%%
+grid_factory = gpu.grid_generator.GridFactory.make()
+grid_builder = gpu.grid_generator.MultipleGridBuilder.make_shared(grid_factory)
+
+#%%
+dx = reference_height/nodes_per_height
+dt = dx * mach / (np.sqrt(3) * velocity)
+velocity_lb = velocity * dt / dx # LB units
+viscosity_lb = viscosity * dt / (dx * dx) # LB units
+
+pressure_gradient = u_star**2 / reference_height
+pressure_gradient_lb = pressure_gradient * dt**2 / dx
+
+logger.vf_log_info(f"velocity    = {velocity_lb:1.6} dx/dt")
+logger.vf_log_info(f"dt          = {dt:1.6}")
+logger.vf_log_info(f"dx          = {dx:1.6}")
+logger.vf_log_info(f"u*          = {u_star:1.6}")
+logger.vf_log_info(f"dpdx        = {pressure_gradient:1.6}")
+logger.vf_log_info(f"dpdx        = {pressure_gradient_lb:1.6} dx/dt^2")
+logger.vf_log_info(f"viscosity   = {viscosity_lb:1.6} dx^2/dt")
+
+
 #%%
 config = basics.ConfigurationFile()
 config.load(str(config_file))
 #%%
 para = gpu.Parameter(config, comm.get_number_of_process(), comm.get_pid())
 
-dt = dx * mach / (np.sqrt(3) * velocity)
-velocity_lb = velocity * dt / dx # LB units
-viscosity_lb = viscosity * dt / (dx * dx) # LB units
+
 
 #%%
 para.set_devices([0])
@@ -55,54 +82,52 @@ para.set_max_level(1)
 para.set_velocity(velocity_lb)
 para.set_viscosity(viscosity_lb)    
 para.set_velocity_ratio(dx/dt)
-para.set_main_kernel("CumulantK17CompChim")
+para.set_viscosity_ratio(dx*dx/dt)
+para.set_use_AMD(use_AMD)
+
+para.set_main_kernel("TurbulentViscosityCumulantK17CompChim" if para.get_use_AMD() else "CummulantK17CompChim")
+
+para.set_SGS_constant(0.083)
 
 def init_func(coord_x, coord_y, coord_z):
-    return [0.0, velocity_lb, 0.0, 0.0]
+    return [
+        0.0, 
+        (u_star/kappa*np.log(max(coord_z/z_0,0)+1) + 2*np.sin(np.pi*16*coord_x/length[0])*np.sin(np.pi*8*coord_z/length[2]))/((coord_z/reference_height)**2+0.1)*dt/dx, 
+        2*np.sin(np.pi*16*coord_x/length[0])*np.sin(np.pi*8*coord_z/length[2])/((coord_z/reference_height)**2+0.1)*dt/dx, 
+        8*u_star/kappa*(np.sin(np.pi*8*coord_y/reference_height)*np.sin(np.pi*8*coord_z/reference_height)+np.sin(np.pi*8*coord_x/length[0]))/((length[2]/2-coord_z)**2+0.1)*dt/dx
+        ]
 
 para.set_initial_condition(init_func)
-para.set_t_out(timeStepOut)
+para.set_t_out(int(t_out/dt))
 para.set_t_end(int(t_end/dt))
 para.set_is_body_force(True)
+para.set_has_wall_model_monitor(True)
 
-#%%
-grid_builder.set_velocity_boundary_condition(gpu.SideType.MX, velocity_lb, 0.0, 0.0)
-grid_builder.set_velocity_boundary_condition(gpu.SideType.PX, velocity_lb, 0.0, 0.0)
-
-grid_builder.set_velocity_boundary_condition(gpu.SideType.MY, velocity_lb, 0.0, 0.0)
-grid_builder.set_velocity_boundary_condition(gpu.SideType.PY, velocity_lb, 0.0, 0.0)
 
-grid_builder.set_velocity_boundary_condition(gpu.SideType.MZ, velocity_lb, 0.0, 0.0)
-grid_builder.set_velocity_boundary_condition(gpu.SideType.PZ, velocity_lb, 0.0, 0.0)
-
-#%%
-cuda_memory_manager = gpu.CudaMemoryManager.make(para)
-grid_generator = gpu.GridProvider.make_grid_generator(grid_builder, para, cuda_memory_manager)
+grid_builder.add_coarse_grid(0.0, 0.0, 0.0, *length, dx)
+grid_builder.set_periodic_boundary_condition(True, True, False)
+grid_builder.build_grids(basics.LbmOrGks.LBM, False)
 #%%
-turb_pos = np.array([3,3,3])*reference_diameter
-epsilon = 5
-density = 1.225
-level = 0
-n_blades = 3
-n_blade_nodes = 32
-alm = gpu.ActuatorLine(n_blades, density, n_blade_nodes, epsilon, *turb_pos, reference_diameter, level, dt, dx)
-para.add_actuator(alm)
+sampling_offset = 2
+grid_builder.set_stress_boundary_condition(gpu.SideType.MZ, 0.0, 0.0, 1.0, sampling_offset, z_0/dx)
+grid_builder.set_slip_boundary_condition(gpu.SideType.PZ, 0.0, 0.0, 0.0)
+
 #%%
-point_probe = gpu.probes.PointProbe("pointProbe", str(output_path), 100, 500, 100)
-point_probe.add_probe_points_from_list(np.array([1,2,5])*reference_diameter, np.array([3,3,3])*reference_diameter, np.array([3,3,3])*reference_diameter)
-point_probe.add_post_processing_variable(gpu.probes.PostProcessingVariable.Means)
+cuda_memory_manager = gpu.CudaMemoryManager(para)
+grid_generator = gpu.GridProvider.make_grid_generator(grid_builder, para, cuda_memory_manager, comm)
 
-para.add_probe(point_probe)
+#%%
+wall_probe = gpu.probes.WallModelProbe("wallModelProbe", str(output_path), int(t_start_averaging/dt), int(t_start_tmp_averaging/dt), int(t_averaging/dt/4), int(t_start_out_probe/dt), int(t_out_probe/dt))
+wall_probe.add_all_available_statistics()
+wall_probe.set_file_name_to_n_out()
+wall_probe.set_force_output_to_stress(True)
+if para.get_is_body_force():
+    wall_probe.set_evaluate_pressure_gradient(True)
+planar_probe = gpu.probes.PlanarAverageProbe("planarAverageProbe", str(output_path), int(t_start_averaging/dt), int(t_start_tmp_averaging/dt), int(t_averaging/dt), int(t_start_out_probe/dt), int(t_out_probe/dt), "z")
+para.add_probe(wall_probe)
 
-plane_probe = gpu.probes.PlaneProbe("planeProbe", str(output_path), 100, 500, 100)
-plane_probe.set_probe_plane(5*reference_diameter, 0, 0, dx, length[1], length[2])
-para.add_probe(plane_probe)
 #%%
-sim = gpu.Simulation(comm)
-kernel_factory = gpu.KernelFactory.get_instance()
-sim.set_factories(kernel_factory, gpu.PreProcessorFactory.get_instance())
-sim.init(para, grid_generator, gpu.FileWriter(), cuda_memory_manager)
+sim = gpu.Simulation(para, cuda_memory_manager, comm, grid_generator)
 #%%
 sim.run()
-sim.free()
 MPI.Finalize()
\ No newline at end of file
diff --git a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
index 20b2bd7786736a44b21550ce5a92dc0c077ee87a..9a44362d1b0753736e719aeefa100890174ad799 100644
--- a/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
+++ b/apps/gpu/LBM/ActuatorLine/ActuatorLine.cpp
@@ -194,13 +194,14 @@ void multipleLevel(const std::string& configPath)
     std::vector<real> probeCoordsZ = {3*reference_diameter,3*reference_diameter,3*reference_diameter};
     pointProbe->addProbePointsFromList(probeCoordsX, probeCoordsY, probeCoordsZ);
     // pointProbe->addProbePointsFromXNormalPlane(2*D, 0.0, 0.0, L_y, L_z, (uint)L_y/dx, (uint)L_z/dx);
-    pointProbe->addPostProcessingVariable(PostProcessingVariable::Means);
-    pointProbe->addPostProcessingVariable(PostProcessingVariable::Variances);
+    
+    pointProbe->addStatistic(Statistic::Means);
+    pointProbe->addStatistic(Statistic::Variances);
     para->addProbe( pointProbe );
 
-    SPtr<PlaneProbe> planeProbe = SPtr<PlaneProbe>( new PlaneProbe("planeProbe", para->getOutputPath(), 100, 500, 100) );
+    SPtr<PlaneProbe> planeProbe = SPtr<PlaneProbe>( new PlaneProbe("planeProbe", para->getOutputPath(), 100, 500, 100, 100) );
     planeProbe->setProbePlane(5*reference_diameter, 0, 0, dx, L_y, L_z);
-    planeProbe->addPostProcessingVariable(PostProcessingVariable::Means);
+    planeProbe->addStatistic(Statistic::Means);
     para->addProbe( planeProbe );
 
 
diff --git a/pythonbindings/src/gpu/gpu.cpp b/pythonbindings/src/gpu/gpu.cpp
index 1dd960dbb3ff9c9ef21cecb36e5df90e74360726..dc110cd5e19a9aad4937f9c2133ddf74c0ddf9bf 100644
--- a/pythonbindings/src/gpu/gpu.cpp
+++ b/pythonbindings/src/gpu/gpu.cpp
@@ -5,14 +5,10 @@
 #include "submodules/parameter.cpp"
 #include "submodules/boundary_conditions.cpp"
 #include "submodules/communicator.cpp"
-#include "submodules/grid_builder.cpp"
 #include "submodules/cuda_memory_manager.cpp"
 #include "submodules/grid_provider.cpp"
-#include "submodules/probes.cpp"
-#include "submodules/kernel_factory.cpp"
-#include "submodules/pre_processor_factory.cpp"
-#include "submodules/file_writer.cpp"
 #include "submodules/grid_generator.cpp"
+#include "submodules/probes.cpp"
 
 namespace gpu
 {
@@ -27,13 +23,9 @@ namespace gpu
         actuator_line::makeModule(gpuModule);
         boundary_conditions::makeModule(gpuModule);
         communicator::makeModule(gpuModule); 
-        grid_builder::makeModule(gpuModule);
         cuda_memory_manager::makeModule(gpuModule);
         grid_provider::makeModule(gpuModule);
         probes::makeModule(gpuModule);
-        kernel_factory::makeModule(gpuModule);
-        pre_processor_factory::makeModule(gpuModule);
-        file_writer::makeModule(gpuModule);
         grid_generator::makeModule(gpuModule);
         return gpuModule;
     }
diff --git a/pythonbindings/src/gpu/submodules/cuda_memory_manager.cpp b/pythonbindings/src/gpu/submodules/cuda_memory_manager.cpp
index f337aeb306fd1836e5bc2962e85f45d593185836..bf27080cb3cd050343ba42b0571827ed58870cfd 100644
--- a/pythonbindings/src/gpu/submodules/cuda_memory_manager.cpp
+++ b/pythonbindings/src/gpu/submodules/cuda_memory_manager.cpp
@@ -10,6 +10,6 @@ namespace cuda_memory_manager
     void makeModule(py::module_ &parentModule)
     {
         py::class_<CudaMemoryManager, std::shared_ptr<CudaMemoryManager>>(parentModule, "CudaMemoryManager")
-        .def("make", &CudaMemoryManager::make, py::return_value_policy::reference);
+        .def(py::init<std::shared_ptr<Parameter>>(), "parameter");
     }
 }
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/file_writer.cpp b/pythonbindings/src/gpu/submodules/file_writer.cpp
deleted file mode 100644
index 2ad90fe7381be215b2d257b5be99caf25db1e0ae..0000000000000000000000000000000000000000
--- a/pythonbindings/src/gpu/submodules/file_writer.cpp
+++ /dev/null
@@ -1,17 +0,0 @@
-#include <pybind11/pybind11.h>
-#include <gpu/VirtualFluids_GPU/Output/FileWriter.h>
-#include <gpu/VirtualFluids_GPU/Output/DataWriter.h>
-
-
-namespace file_writer
-{
-    namespace py = pybind11;
-
-    void makeModule(py::module_ &parentModule)
-    {
-        py::class_<DataWriter, std::shared_ptr<DataWriter>>(parentModule, "_DataWriter");
-
-        py::class_<FileWriter, DataWriter, std::shared_ptr<FileWriter>>(parentModule, "FileWriter")
-        .def(py::init<>());
-    }
-}
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/grid_builder.cpp b/pythonbindings/src/gpu/submodules/grid_builder.cpp
deleted file mode 100644
index 09241f71cbe9cd013600c34e4d0a07e53a1bc79f..0000000000000000000000000000000000000000
--- a/pythonbindings/src/gpu/submodules/grid_builder.cpp
+++ /dev/null
@@ -1,36 +0,0 @@
-#include <pybind11/pybind11.h>
-#include "gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
-#include "gpu/GridGenerator/grid/GridBuilder/GridBuilder.h"
-#include "gpu/GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
-#include "gpu/GridGenerator/geometries/Object.h"
-#include "gpu/GridGenerator/grid/BoundaryConditions/Side.h"
-
-namespace grid_builder
-{
-    namespace py = pybind11;
-
-    void makeModule(py::module_ &parentModule)
-    {        
-        py::class_<GridBuilder, std::shared_ptr<GridBuilder>>(parentModule, "GridBuilder")
-        .def("get_number_of_grid_levels", &GridBuilder::getNumberOfGridLevels)
-        .def("get_grid", &GridBuilder::getGrid);
-
-        py::class_<LevelGridBuilder, GridBuilder, std::shared_ptr<LevelGridBuilder>>(parentModule, "LevelGridBuilder")
-        .def("get_grid", py::overload_cast<int, int>(&LevelGridBuilder::getGrid))
-        .def("set_slip_boundary_condition", &LevelGridBuilder::setSlipBoundaryCondition)
-        .def("set_velocity_boundary_condition", &LevelGridBuilder::setVelocityBoundaryCondition)
-        .def("set_pressure_boundary_condition", &LevelGridBuilder::setPressureBoundaryCondition)
-        .def("set_periodic_boundary_condition", &LevelGridBuilder::setPeriodicBoundaryCondition)
-        .def("set_no_slip_boundary_condition", &LevelGridBuilder::setNoSlipBoundaryCondition);
-
-        py::class_<MultipleGridBuilder, LevelGridBuilder, std::shared_ptr<MultipleGridBuilder>>(parentModule, "MultipleGridBuilder")
-        .def("make_shared", &MultipleGridBuilder::makeShared, py::return_value_policy::reference)
-        .def("add_coarse_grid", &MultipleGridBuilder::addCoarseGrid)
-        .def("add_grid", py::overload_cast<Object*>(&MultipleGridBuilder::addGrid))
-        .def("add_grid", py::overload_cast<Object*, uint>(&MultipleGridBuilder::addGrid))
-        .def("add_geometry", py::overload_cast<Object*>(&MultipleGridBuilder::addGeometry))
-        .def("add_geometry", py::overload_cast<Object*, uint>(&MultipleGridBuilder::addGeometry))
-        .def("get_number_of_levels", &MultipleGridBuilder::getNumberOfLevels)
-        .def("build_grids", &MultipleGridBuilder::buildGrids);
-    }
-}
\ 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 b9babe1d6b45f69d516a4cd57377eaa6c1b37149..579c06c4e00cae9646ced8b554d71631eeb7e793 100644
--- a/pythonbindings/src/gpu/submodules/grid_generator.cpp
+++ b/pythonbindings/src/gpu/submodules/grid_generator.cpp
@@ -5,6 +5,10 @@
 #include "gpu/GridGenerator/geometries/Cuboid/Cuboid.h"
 #include "gpu/GridGenerator/geometries/Sphere/Sphere.h"
 #include "gpu/GridGenerator/geometries/TriangularMesh/TriangularMesh.h"
+#include "gpu/GridGenerator/grid/GridFactory.h"
+#include "gpu/GridGenerator/grid/GridBuilder/GridBuilder.h"
+#include "gpu/GridGenerator/grid/GridBuilder/LevelGridBuilder.h"
+#include "gpu/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h"
 
 namespace grid_generator
 {
@@ -13,6 +17,9 @@ namespace grid_generator
     {  
         py::module gridGeneratorModule = parentModule.def_submodule("grid_generator");
 
+        py::class_<GridFactory, std::shared_ptr<GridFactory>>(gridGeneratorModule, "GridFactory")
+        .def("make", &GridFactory::make, py::return_value_policy::reference);
+
         py::class_<BoundingBox>(gridGeneratorModule, "BoundingBox")
         .def(py::init<real, real, real, real, real, real>(),"min_x","max_x","min_y","max_y","min_z","max_z");
 
@@ -33,6 +40,29 @@ namespace grid_generator
         py::class_<TriangularMesh, Object, std::shared_ptr<TriangularMesh>>(gridGeneratorModule, "TriangularMesh")
         .def("make", &TriangularMesh::make, py::return_value_policy::reference);
 
+        py::class_<GridBuilder, std::shared_ptr<GridBuilder>>(gridGeneratorModule, "GridBuilder")
+        .def("get_number_of_grid_levels", &GridBuilder::getNumberOfGridLevels)
+        .def("get_grid", &GridBuilder::getGrid);
+
+        py::class_<LevelGridBuilder, GridBuilder, std::shared_ptr<LevelGridBuilder>>(gridGeneratorModule, "LevelGridBuilder")
+        .def("get_grid", py::overload_cast<int, int>(&LevelGridBuilder::getGrid))
+        .def("set_slip_boundary_condition", &LevelGridBuilder::setSlipBoundaryCondition)
+        .def("set_velocity_boundary_condition", &LevelGridBuilder::setVelocityBoundaryCondition)
+        .def("set_pressure_boundary_condition", &LevelGridBuilder::setPressureBoundaryCondition)
+        .def("set_periodic_boundary_condition", &LevelGridBuilder::setPeriodicBoundaryCondition)
+        .def("set_no_slip_boundary_condition", &LevelGridBuilder::setNoSlipBoundaryCondition)
+        .def("set_stress_boundary_condition", &LevelGridBuilder::setStressBoundaryCondition);
+
+        py::class_<MultipleGridBuilder, LevelGridBuilder, std::shared_ptr<MultipleGridBuilder>>(gridGeneratorModule, "MultipleGridBuilder")
+        .def("make_shared", &MultipleGridBuilder::makeShared, py::return_value_policy::reference)
+        .def("add_coarse_grid", &MultipleGridBuilder::addCoarseGrid)
+        .def("add_grid", py::overload_cast<Object*>(&MultipleGridBuilder::addGrid))
+        .def("add_grid", py::overload_cast<Object*, uint>(&MultipleGridBuilder::addGrid))
+        .def("add_geometry", py::overload_cast<Object*>(&MultipleGridBuilder::addGeometry))
+        .def("add_geometry", py::overload_cast<Object*, uint>(&MultipleGridBuilder::addGeometry))
+        .def("get_number_of_levels", &MultipleGridBuilder::getNumberOfLevels)
+        .def("build_grids", &MultipleGridBuilder::buildGrids);
+
         return gridGeneratorModule;
     }
 }
diff --git a/pythonbindings/src/gpu/submodules/grid_provider.cpp b/pythonbindings/src/gpu/submodules/grid_provider.cpp
index 5a5514dd8b14d24f4818740e115c8504bc973726..02ff273e2cd1a2022943e19c9a48a447d9dfe54b 100644
--- a/pythonbindings/src/gpu/submodules/grid_provider.cpp
+++ b/pythonbindings/src/gpu/submodules/grid_provider.cpp
@@ -1,8 +1,8 @@
 #include <pybind11/pybind11.h>
 #include "gpu/VirtualFluids_GPU/DataStructureInitializer/GridProvider.h"
-#include <gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h>
-#include <gpu/VirtualFluids_GPU/Parameter/Parameter.h>
-#include "gpu/GridGenerator/grid/GridBuilder/GridBuilder.h"
+// #include <gpu/VirtualFluids_GPU/GPU/CudaMemoryManager.h>
+// #include <gpu/VirtualFluids_GPU/Parameter/Parameter.h>
+// #include "gpu/GridGenerator/grid/GridBuilder/GridBuilder.h"
 
 namespace grid_provider
 {
diff --git a/pythonbindings/src/gpu/submodules/kernel_factory.cpp b/pythonbindings/src/gpu/submodules/kernel_factory.cpp
deleted file mode 100644
index af710948aa75900ba9f4051360ec26dd8510118a..0000000000000000000000000000000000000000
--- a/pythonbindings/src/gpu/submodules/kernel_factory.cpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#include <pybind11/pybind11.h>
-#include <gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.h>
-#include <gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactory.h>
-
-namespace kernel_factory
-{
-    namespace py = pybind11;
-
-    void makeModule(py::module_ &parentModule)
-    {
-        py::class_<KernelFactory, std::shared_ptr<KernelFactory>>(parentModule, "_KernelFactory");
-        
-        py::class_<KernelFactoryImp, KernelFactory, std::shared_ptr<KernelFactoryImp>>(parentModule, "KernelFactory")
-        .def("get_instance", &KernelFactoryImp::getInstance, py::return_value_policy::reference);
-    }
-}
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/parameter.cpp b/pythonbindings/src/gpu/submodules/parameter.cpp
index 6e65c9d4c7382298739b1254c4294e65855567da..7b4e67f101e3928abbd4262557864ea1d0f45b02 100644
--- a/pythonbindings/src/gpu/submodules/parameter.cpp
+++ b/pythonbindings/src/gpu/submodules/parameter.cpp
@@ -70,6 +70,11 @@ namespace parameter
         .def("get_viscosity_ratio", &Parameter::getViscosityRatio)
         .def("get_density_ratio", &Parameter::getDensityRatio)
         .def("get_force_ratio", &Parameter::getForceRatio)
+        .def("get_use_AMD", &Parameter::getUseAMD)
+        .def("get_use_Wale", &Parameter::getUseWale)
+        .def("get_SGS_constant", &Parameter::getSGSConstant)
+        .def("get_is_body_force", &Parameter::getIsBodyForce)
+        .def("set_has_wall_model_monitor", &Parameter::setHasWallModelMonitor)
         ;
     }
 }
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/pre_processor_factory.cpp b/pythonbindings/src/gpu/submodules/pre_processor_factory.cpp
deleted file mode 100644
index b76ae285b413594847f6672e7f3fc0e656da3cec..0000000000000000000000000000000000000000
--- a/pythonbindings/src/gpu/submodules/pre_processor_factory.cpp
+++ /dev/null
@@ -1,16 +0,0 @@
-#include <pybind11/pybind11.h>
-#include <gpu/VirtualFluids_GPU/PreProcessor/PreProcessorFactory/PreProcessorFactory.h>
-#include <gpu/VirtualFluids_GPU/PreProcessor/PreProcessorFactory/PreProcessorFactoryImp.h>
-
-namespace pre_processor_factory
-{
-    namespace py = pybind11;
-
-    void makeModule(py::module_ &parentModule)
-    {
-        py::class_<PreProcessorFactory, std::shared_ptr<PreProcessorFactory>>(parentModule, "_PreProcessorFactory");
-
-        py::class_<PreProcessorFactoryImp, PreProcessorFactory, std::shared_ptr<PreProcessorFactoryImp>>(parentModule, "PreProcessorFactory")
-        .def("get_instance", &PreProcessorFactoryImp::getInstance, py::return_value_policy::reference);
-    }
-}
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/probes.cpp b/pythonbindings/src/gpu/submodules/probes.cpp
index 9b751a513d23a8a17e4e9981ababc720fc7346b0..6993d9617d870922d7ed90ed9ecbebb8a797be25 100644
--- a/pythonbindings/src/gpu/submodules/probes.cpp
+++ b/pythonbindings/src/gpu/submodules/probes.cpp
@@ -3,6 +3,8 @@
 #include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/Probe.h>
 #include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PointProbe.h>
 #include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlaneProbe.h>
+#include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/WallModelProbe.h>
+#include <gpu/VirtualFluids_GPU/PreCollisionInteractor/Probes/PlanarAverageProbe.h>
 #include <gpu/VirtualFluids_GPU/PreCollisionInteractor/PreCollisionInteractor.h>
 
 namespace probes
@@ -13,13 +15,23 @@ namespace probes
     {
         py::module probeModule = parentModule.def_submodule("probes");
 
-        py::enum_<PostProcessingVariable>(probeModule, "PostProcessingVariable")
-        .value("Instantaneous", PostProcessingVariable::Instantaneous)
-        .value("Means", PostProcessingVariable::Means)
-        .value("Variances", PostProcessingVariable::Variances);
+        py::enum_<Statistic>(probeModule, "Statistic")
+        .value("Instantaneous", Statistic::Instantaneous)
+        .value("Means", Statistic::Means)
+        .value("Variances", Statistic::Variances)
+        .value("SpatialMeans", Statistic::SpatialMeans)
+        .value("SpatioTemporalMeans", Statistic::SpatioTemporalMeans)
+        .value("SpatialCovariances", Statistic::SpatialCovariances)
+        .value("SpatioTemporalCovariances", Statistic::SpatioTemporalCovariances)
+        .value("SpatialSkewness", Statistic::SpatialSkewness)
+        .value("SpatioTemporalSkewness", Statistic::SpatioTemporalSkewness)
+        .value("SpatialFlatness", Statistic::SpatialFlatness)
+        .value("SpatioTemporalFlatness", Statistic::SpatioTemporalFlatness);
 
         py::class_<Probe, PreCollisionInteractor, std::shared_ptr<Probe>>(probeModule, "Probe")
-        .def("add_post_processing_variable", &Probe::addPostProcessingVariable);
+        .def("add_statistic", &Probe::addStatistic)
+        .def("set_file_name_to_n_out", &Probe::setFileNameToNOut)
+        .def("add_all_available_statistics", &Probe::addAllAvailableStatistics);
 
         py::class_<PointProbe, Probe, std::shared_ptr<PointProbe>>(probeModule, "PointProbe")
         .def(py::init<
@@ -27,10 +39,12 @@ namespace probes
                         const std::string,
                         uint,
                         uint, 
+                        uint,
                         uint>(), 
                         "probe_name",
                         "output_path"
                         "t_start_avg",
+                        "t_avg",
                         "t_start_out",
                         "t_out")
         .def("add_probe_points_from_list", &PointProbe::addProbePointsFromList)
@@ -42,14 +56,55 @@ namespace probes
                         const std::string,
                         uint,
                         uint, 
+                        uint,
                         uint>(), 
                         "probe_name",
                         "output_path"
                         "t_start_avg",
+                        "t_avg",
                         "t_start_out",
                         "t_out")
         .def("set_probe_plane", &PlaneProbe::setProbePlane);
 
+        py::class_<PlanarAverageProbe, Probe, std::shared_ptr<PlanarAverageProbe>>(probeModule, "PlanarAverageProbe")
+        .def(py::init<
+                        const std::string,
+                        const std::string,
+                        uint,
+                        uint,
+                        uint,
+                        uint,
+                        uint,
+                        char>(),
+                        "probe_name",
+                        "output_path",
+                        "t_start_avg",
+                        "t_start_tmp_avg",
+                        "t_avg",
+                        "t_start_out",
+                        "t_out",
+                        "plane_normal");
+
+
+        py::class_<WallModelProbe, Probe, std::shared_ptr<WallModelProbe>>(probeModule, "WallModelProbe")
+        .def(py::init<
+                        const std::string,
+                        const std::string,
+                        uint,
+                        uint, 
+                        uint,
+                        uint,
+                        uint>(), 
+                        "probe_name",
+                        "output_path"
+                        "t_start_avg",
+                        "t_start_tmp_avg",
+                        "t_avg",
+                        "t_start_out",
+                        "t_out")
+        .def("set_force_output_to_stress", &WallModelProbe::setForceOutputToStress)
+        .def("set_evaluate_pressure_gradient", &WallModelProbe::setEvaluatePressureGradient);
+
         return probeModule;
     }
 }
\ No newline at end of file
diff --git a/pythonbindings/src/gpu/submodules/simulation.cpp b/pythonbindings/src/gpu/submodules/simulation.cpp
index 9683e8768417ad2422fb1a53ef6a428543bceffc..b775d604ba41530223f22738c72785b2c15348b3 100644
--- a/pythonbindings/src/gpu/submodules/simulation.cpp
+++ b/pythonbindings/src/gpu/submodules/simulation.cpp
@@ -15,11 +15,18 @@ namespace simulation
 
     void makeModule(py::module_ &parentModule)
     {
+        // missing setFactories and setDataWriter, not possible to wrap these functions as long as they take unique ptr arguments
         py::class_<Simulation>(parentModule, "Simulation")
-        .def(py::init<vf::gpu::Communicator&>(), "communicator")
-        .def("set_factories", &Simulation::setFactories)
-        .def("init", &Simulation::init)
+        .def(py::init<  std::shared_ptr<Parameter>,
+                        std::shared_ptr<CudaMemoryManager>,
+                        vf::gpu::Communicator &,
+                        GridProvider &>(), 
+                        "parameter",
+                        "memoryManager",
+                        "communicator",
+                        "gridProvider")
         .def("run", &Simulation::run)
-        .def("free", &Simulation::free);
+        .def("addKineticEnergyAnalyzer", &Simulation::addKineticEnergyAnalyzer)
+        .def("addEnstrophyAnalyzer", &Simulation::addEnstrophyAnalyzer);
     }
 }
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cfg b/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cfg
deleted file mode 100644
index e414d4f3173e555b8944fa9637ebbd2023ce393c..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cfg
+++ /dev/null
@@ -1,3 +0,0 @@
-# these two parameters need to be defined in each config file
-Path = /output/path
-GridPath = /path/to/grid
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cpp b/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cpp
index e445f59663157ccb67b6ebb9dc5a08d4af4c2679..473e8c1a5424cf7ddd05f6ed0a534814a3971dc6 100644
--- a/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cpp
+++ b/src/gpu/VirtualFluids_GPU/Communication/ExchangeData27Test.cpp
@@ -28,7 +28,7 @@ protected:
 
     void SetUp() override
     {
-        para = std::make_shared<Parameter>(1, 0);
+        para = std::make_shared<Parameter>();
         para->setMaxLevel(level + 1);       // setMaxLevel resizes parH
         para->initLBMSimulationParameter(); // init parH
 
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cfg b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cfg
deleted file mode 100644
index e414d4f3173e555b8944fa9637ebbd2023ce393c..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cfg
+++ /dev/null
@@ -1,3 +0,0 @@
-# these two parameters need to be defined in each config file
-Path = /output/path
-GridPath = /path/to/grid
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cpp b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cpp
index 6e777e2c13212d33251f563dcad525bf92f9d566..405370c905adc9937badde2f6e54f2d54942056b 100644
--- a/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cpp
+++ b/src/gpu/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cpp
@@ -214,7 +214,7 @@ private:
 
     void SetUp() override
     {
-        para        = std::make_shared<Parameter>(1, 0);
+        para        = std::make_shared<Parameter>();
         testSubject = createTestSubjectFCBorderBulk();
     }
 };
@@ -313,7 +313,7 @@ private:
 
     void SetUp() override
     {
-        para        = std::make_shared<Parameter>(1, 0);
+        para        = std::make_shared<Parameter>();
         testSubject = createTestSubjectReorderSendIndices();
     };
 };
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 2a6ed54f7691b2f356f18e3e6dac57b670bcdfa5..c6adcac6a9ae83308d843d3bdc819f68f0e4ca03 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -45,14 +45,11 @@
 #include "VirtualFluids_GPU_export.h"
 
 struct curandStateXORWOW;
-typedef struct curandStateXORWOW curandState;
-namespace vf
-{
-namespace basics
+using curandState = struct curandStateXORWOW;
+namespace vf:: basics
 {
 class ConfigurationFile;
 }
-} // namespace vf
 class CudaStreamManager;
 
 //! \struct LBMSimulationParameter
@@ -377,8 +374,8 @@ struct LBMSimulationParameter {
 class VIRTUALFLUIDS_GPU_EXPORT Parameter
 {
 public:
-    Parameter(const vf::basics::ConfigurationFile &configData, int numberOfProcesses, int myId);
-    Parameter(int numberOfProcesses, int myId);
+    Parameter(const vf::basics::ConfigurationFile &configData, const int numberOfProcesses = 1, const int myId = 0);
+    Parameter(const int numberOfProcesses = 1, const int myId = 0);
     ~Parameter();
     void initLBMSimulationParameter();
 
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp b/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
index 5d9485b5aa0715f7141d8db411b3e1ce5399bd41..9e05ed1332b34420656e6c0c81f07501da7c7aac 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/ParameterTest.cpp
@@ -39,7 +39,7 @@ TEST(ParameterTest, check_all_Parameter_CanBePassedToConstructor)
     vf::basics::ConfigurationFile config;
     config.load(filePath.string());
 
-    Parameter para(config, 1, 0);
+    Parameter para(config);
 
     // test optional parameter
     EXPECT_THAT(para.getOutputPath(), testing::Eq("/output/path/"));
@@ -151,7 +151,7 @@ TEST(ParameterTest, check_all_Parameter_CanBePassedToConstructor)
 
 TEST(ParameterTest, defaultGridPath)
 {
-    Parameter para(1, 0);
+    Parameter para;
     EXPECT_THAT(para.getGridPath(), testing::Eq("grid/"));
     EXPECT_THAT(para.getConcentration(), testing::Eq("grid/conc.dat"));
 }
@@ -190,7 +190,7 @@ TEST(ParameterTest, setGridPathOverridesConfigFile)
 
 TEST(ParameterTest, userMissedSlash)
 {
-    Parameter para(1, 0);
+    Parameter para;
     para.setGridPath("gridPathTest");
 
     EXPECT_THAT(para.getGridPath(), testing::Eq("gridPathTest/"));