diff --git a/Python/SlurmTests/poiseuille/settings.py b/Python/SlurmTests/poiseuille/settings.py
index 4b4a1e4e9cc7f6118a60c22a40c70b027e3ac4e2..d70a9574b3d7a3ed617069dd9aaa090c9f293ab0 100644
--- a/Python/SlurmTests/poiseuille/settings.py
+++ b/Python/SlurmTests/poiseuille/settings.py
@@ -1,25 +1,24 @@
 import os
 from acousticscaling import OneDirectionalAcousticScaling
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import RuntimeParameters, GridParameters, PhysicalParameters
+from pyfluids import cpu
 
 
-grid_params = GridParameters()
+grid_params = cpu.parameters.GridParameters()
 grid_params.node_distance = 1
 grid_params.number_of_nodes_per_direction = [1, 1, 16]
 grid_params.blocks_per_direction = [1, 1, 4]
 grid_params.periodic_boundary_in_x1 = True
 grid_params.periodic_boundary_in_x2 = True
 
-physical_params = PhysicalParameters()
+physical_params = cpu.parameters.PhysicalParameters()
 physical_params.lattice_viscosity = 1e-4
 
-runtime_params = RuntimeParameters()
+runtime_params = cpu.parameters.RuntimeParameters()
 runtime_params.number_of_threads = int(os.environ["PYFLUIDS_NUM_THREADS"])
 runtime_params.number_of_timesteps = 4_000_000
 runtime_params.timestep_log_interval = 1_000_000
 
-kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
 kernel.use_forcing = True
 kernel.forcing_in_x1 = 5e-10
 
diff --git a/Python/acousticscaling.py b/Python/acousticscaling.py
index a664b8e924d648b680562b9aef11bee87b3562b1..9ba6f87653777830b7e6d860ae569246ff9505b7 100644
--- a/Python/acousticscaling.py
+++ b/Python/acousticscaling.py
@@ -1,22 +1,21 @@
-from pyfluids.cpu.kernel import LBMKernel
-from pyfluids.cpu.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+from pyfluids import cpu
 
 
 class OneDirectionalAcousticScaling:
 
-    def __init__(self, grid_parameters: GridParameters,
-                 physical_parameters: PhysicalParameters,
-                 runtime_parameters: RuntimeParameters,
-                 kernel: LBMKernel):
+    def __init__(self, grid_parameters: cpu.parameters.GridParameters,
+                 physical_parameters: cpu.parameters.PhysicalParameters,
+                 runtime_parameters: cpu.parameters.RuntimeParameters,
+                 kernel: cpu.kernel.LBMKernel):
         self._grid_params = grid_parameters
         self._physical_params = physical_parameters
         self._runtime_params = runtime_parameters
         self._kernel = kernel
 
-    def configuration_for_scale_level(self, level: int = 1) -> tuple[GridParameters,
-                                                                PhysicalParameters,
-                                                                RuntimeParameters,
-                                                                LBMKernel]:
+    def configuration_for_scale_level(self, level: int = 1) -> tuple[cpu.parameters.GridParameters,
+                                                                cpu.parameters.PhysicalParameters,
+                                                                cpu.parameters.RuntimeParameters,
+                                                                cpu.kernel.LBMKernel]:
         if level < 0:
             raise ValueError("level must be >= 0")
 
@@ -27,8 +26,8 @@ class OneDirectionalAcousticScaling:
 
         return grid_params, physical_params, runtime_params, kernel
 
-    def clone_grid_params_for_level(self, level) -> GridParameters:
-        grid_params = GridParameters()
+    def clone_grid_params_for_level(self, level) -> cpu.parameters.GridParameters:
+        grid_params = cpu.parameters.GridParameters()
         grid_params.reference_direction_index = self._grid_params.reference_direction_index
         grid_params.periodic_boundary_in_x1 = self._grid_params.periodic_boundary_in_x1
         grid_params.periodic_boundary_in_x2 = self._grid_params.periodic_boundary_in_x2
@@ -51,7 +50,7 @@ class OneDirectionalAcousticScaling:
         return grid_params
 
     def clone_physical_parameters(self, level):
-        physical_params = PhysicalParameters()
+        physical_params = cpu.parameters.PhysicalParameters()
         physical_params.lattice_viscosity = self._physical_params.lattice_viscosity
 
         if level > 0:
@@ -60,7 +59,7 @@ class OneDirectionalAcousticScaling:
         return physical_params
 
     def clone_runtime_params_for_level(self, level):
-        runtime_params = RuntimeParameters()
+        runtime_params = cpu.parameters.RuntimeParameters()
         runtime_params.number_of_timesteps = self._runtime_params.number_of_timesteps
         runtime_params.number_of_threads = self._runtime_params.number_of_threads
         runtime_params.timestep_log_interval = self._runtime_params.timestep_log_interval
@@ -71,7 +70,7 @@ class OneDirectionalAcousticScaling:
         return runtime_params
 
     def clone_kernel_for_level(self, level):
-        kernel = LBMKernel(self._kernel.type)
+        kernel = cpu.kernel.LBMKernel(self._kernel.type)
         kernel.use_forcing = self._kernel.use_forcing
         kernel.forcing_in_x1 = self._kernel.forcing_in_x1
         kernel.forcing_in_x2 = self._kernel.forcing_in_x2
diff --git a/Python/cubeflow/simulation.py b/Python/cubeflow/simulation.py
index 9e77e8d747c072188d8d81150afa8e2ccb76a792..da8d4f3df9a07c086c643ca40bde1d24873b42b6 100644
--- a/Python/cubeflow/simulation.py
+++ b/Python/cubeflow/simulation.py
@@ -1,13 +1,8 @@
-from pyfluids.cpu import Simulation
-from pyfluids.cpu.boundaryconditions import NoSlipBoundaryCondition, VelocityBoundaryCondition, DensityBoundaryCondition
-from pyfluids.cpu.geometry import GbCuboid3D
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import PhysicalParameters, RuntimeParameters, GridParameters
-from pyfluids.cpu.writer import Writer, OutputFormat
-from pymuparser import Parser
-
 import os
 
+from pyfluids import cpu
+from pymuparser import Parser
+
 
 def get_max_length(number_of_nodes_per_direction, delta_x):
     return (number_of_nodes_per_direction[0] * delta_x,
@@ -15,10 +10,10 @@ def get_max_length(number_of_nodes_per_direction, delta_x):
             number_of_nodes_per_direction[2] * delta_x)
 
 
-physical_params = PhysicalParameters()
+physical_params = cpu.parameters.PhysicalParameters()
 physical_params.lattice_viscosity = 0.005
 
-grid_params = GridParameters()
+grid_params = cpu.parameters.GridParameters()
 grid_params.number_of_nodes_per_direction = [200, 120, 120]
 grid_params.blocks_per_direction = [2, 2, 2]
 grid_params.node_distance = 0.125
@@ -26,7 +21,7 @@ grid_params.periodic_boundary_in_x1 = False
 grid_params.periodic_boundary_in_x2 = True
 grid_params.periodic_boundary_in_x3 = True
 
-runtime_params = RuntimeParameters()
+runtime_params = cpu.parameters.RuntimeParameters()
 runtime_params.timestep_log_interval = 1000
 runtime_params.number_of_timesteps = 50000
 runtime_params.number_of_threads = int(os.environ.get("OMP_NUM_THREADS", 4))
@@ -39,46 +34,46 @@ def run_simulation(physical_parameters=physical_params, grid_parameters=grid_par
     min_x, min_y, min_z = 0, 0, 0
     max_x, max_y, max_z = get_max_length(grid_parameters.number_of_nodes_per_direction, grid_parameters.node_distance)
 
-    bottom_wall = GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, min_z, max_x + wall_thickness,
+    bottom_wall = cpu.geometry.GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, min_z, max_x + wall_thickness,
                              max_y + wall_thickness, min_z - wall_thickness)
 
-    top_wall = GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, max_z, max_x + wall_thickness,
+    top_wall = cpu.geometry.GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, max_z, max_x + wall_thickness,
                           max_y + wall_thickness,
                           max_z + wall_thickness)
 
-    left_wall = GbCuboid3D(min_x - wall_thickness, min_y, min_z - wall_thickness, max_x + wall_thickness,
+    left_wall = cpu.geometry.GbCuboid3D(min_x - wall_thickness, min_y, min_z - wall_thickness, max_x + wall_thickness,
                            min_y - wall_thickness,
                            max_z + wall_thickness)
 
-    right_wall = GbCuboid3D(min_x - wall_thickness, max_y, min_z - wall_thickness, max_x + wall_thickness,
+    right_wall = cpu.geometry.GbCuboid3D(min_x - wall_thickness, max_y, min_z - wall_thickness, max_x + wall_thickness,
                             max_y + wall_thickness, max_z + wall_thickness)
 
-    obstacle = GbCuboid3D(7, 7, 7, 8, 8, 8)
+    obstacle = cpu.geometry.GbCuboid3D(7, 7, 7, 8, 8, 8)
 
-    velocity_boundary = GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, min_z - wall_thickness, min_x,
+    velocity_boundary = cpu.geometry.GbCuboid3D(min_x - wall_thickness, min_y - wall_thickness, min_z - wall_thickness, min_x,
                                    max_y + wall_thickness, max_z + wall_thickness)
 
-    outflow_boundary = GbCuboid3D(max_x, min_y - wall_thickness, min_z - wall_thickness, max_x + wall_thickness,
+    outflow_boundary = cpu.geometry.GbCuboid3D(max_x, min_y - wall_thickness, min_z - wall_thickness, max_x + wall_thickness,
                                   max_y + wall_thickness, max_z + wall_thickness)
 
-    no_slip_bc = NoSlipBoundaryCondition()
+    no_slip_bc = cpu.boundaryconditions.NoSlipBoundaryCondition()
 
-    outflow_bc = DensityBoundaryCondition()
+    outflow_bc = cpu.boundaryconditions.DensityBoundaryCondition()
 
     velocity_function = Parser()
     velocity_function.define_constant("u", 0.07)
     velocity_function.expression = "u"
-    velocity_bc = VelocityBoundaryCondition(True, False, False, velocity_function, 0, -10)
+    velocity_bc = cpu.boundaryconditions.VelocityBoundaryCondition(True, False, False, velocity_function, 0, -10)
 
-    kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+    kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
     # kernel.use_forcing = True
     # kernel.forcing_in_x1 = 3e-6
 
-    writer = Writer()
+    writer = cpu.writer.Writer()
     writer.output_path = "./output"
-    writer.output_format = OutputFormat.BINARY
+    writer.output_format = cpu.writer.OutputFormat.BINARY
 
-    simulation = Simulation()
+    simulation = cpu.Simulation()
     simulation.set_writer(writer)
 
     simulation.set_physical_parameters(physical_parameters)
diff --git a/Python/liddrivencavity/simulation.py b/Python/liddrivencavity/simulation.py
index 155fad2f6f8aade0368c8a7006b88f7985f8822c..468cce6691bc29b78ccf3bd7b4dfe8821c1e546f 100644
--- a/Python/liddrivencavity/simulation.py
+++ b/Python/liddrivencavity/simulation.py
@@ -1,32 +1,27 @@
-from pyfluids.cpu import Simulation
-from pyfluids.cpu.boundaryconditions import NoSlipBoundaryCondition, VelocityBoundaryCondition
-from pyfluids.cpu.geometry import GbCuboid3D
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import GridParameters, PhysicalParameters, RuntimeParameters
-from pyfluids.cpu.writer import Writer, OutputFormat
+from pyfluids import cpu
 from pymuparser import Parser
 
-runtime_params = RuntimeParameters()
+runtime_params = cpu.parameters.RuntimeParameters()
 runtime_params.number_of_threads = 4
 runtime_params.number_of_timesteps = 10000
 runtime_params.timestep_log_interval = 1000
 
-physical_params = PhysicalParameters()
+physical_params = cpu.parameters.PhysicalParameters()
 physical_params.lattice_viscosity = 0.005
 
-grid_params = GridParameters()
+grid_params = cpu.parameters.GridParameters()
 grid_params.number_of_nodes_per_direction = [64, 64, 64]
 grid_params.blocks_per_direction = [2, 2, 2]
 grid_params.node_distance = 1 / 10
 
 
 def run_simulation(physical_params=physical_params, grid_params=grid_params, runtime_params=runtime_params):
-    simulation = Simulation()
-    kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+    simulation = cpu.Simulation()
+    kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
 
-    writer = Writer()
+    writer = cpu.writer.Writer()
     writer.output_path = "./output"
-    writer.output_format = OutputFormat.BINARY
+    writer.output_format = cpu.writer.OutputFormat.BINARY
 
     simulation.set_grid_parameters(grid_params)
     simulation.set_physical_parameters(physical_params)
@@ -34,12 +29,12 @@ def run_simulation(physical_params=physical_params, grid_params=grid_params, run
     simulation.set_kernel_config(kernel)
     simulation.set_writer(writer)
 
-    no_slip_bc_adapter = NoSlipBoundaryCondition()
+    no_slip_bc_adapter = cpu.boundaryconditions.NoSlipBoundaryCondition()
 
     fct = Parser()
     fct.expression = "u"
     fct.define_constant("u", 0.005)
-    velocity_bc_adapter = VelocityBoundaryCondition(True, True, False, fct, 0, -10.0)
+    velocity_bc_adapter = cpu.boundaryconditions.VelocityBoundaryCondition(True, True, False, fct, 0, -10.0)
 
     g_min_x1, g_min_x2, g_min_x3 = 0, 0, 0
     g_max_x1 = grid_params.number_of_nodes_per_direction[0] * grid_params.node_distance
@@ -48,12 +43,12 @@ def run_simulation(physical_params=physical_params, grid_params=grid_params, run
 
     dx = grid_params.node_distance
 
-    wall_x_min = GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_min_x1, g_max_x2 + dx, g_max_x3)
-    wall_x_max = GbCuboid3D(g_max_x1, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_max_x3)
-    wall_y_min = GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_min_x2, g_max_x3)
-    wall_y_max = GbCuboid3D(g_min_x1 - dx, g_max_x2, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_max_x3)
-    wall_z_min = GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_min_x3)
-    wall_z_max = GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_max_x3, g_max_x1 + dx, g_max_x2 + dx, g_max_x3 + dx)
+    wall_x_min = cpu.geometry.GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_min_x1, g_max_x2 + dx, g_max_x3)
+    wall_x_max = cpu.geometry.GbCuboid3D(g_max_x1, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_max_x3)
+    wall_y_min = cpu.geometry.GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_min_x2, g_max_x3)
+    wall_y_max = cpu.geometry.GbCuboid3D(g_min_x1 - dx, g_max_x2, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_max_x3)
+    wall_z_min = cpu.geometry.GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_min_x3 - dx, g_max_x1 + dx, g_max_x2 + dx, g_min_x3)
+    wall_z_max = cpu.geometry.GbCuboid3D(g_min_x1 - dx, g_min_x2 - dx, g_max_x3, g_max_x1 + dx, g_max_x2 + dx, g_max_x3 + dx)
 
     simulation.add_object(wall_x_min, no_slip_bc_adapter, 1, "/geo/wallXmin")
     simulation.add_object(wall_x_max, no_slip_bc_adapter, 1, "/geo/wallXmax")
diff --git a/Python/poiseuille/poiseuille_hpc.py b/Python/poiseuille/poiseuille_hpc.py
index f5f5a1387c9fe234abae0c6f979cc7d5b283d1a4..e5a89250f5d6df8c490dca6d9641f3b9a3f4075c 100644
--- a/Python/poiseuille/poiseuille_hpc.py
+++ b/Python/poiseuille/poiseuille_hpc.py
@@ -1,15 +1,15 @@
 from poiseuille.simulation import run_simulation
-from pyfluids.cpu.parameters import *
+from pyfluids import cpu
 
-grid_parameters = GridParameters()
+grid_parameters = cpu.prameters.GridParameters()
 grid_parameters.number_of_nodes_per_direction = [64, 64, 512]
 grid_parameters.node_distance = 1
 grid_parameters.blocks_per_direction = [1, 2, 2]
 
-physical_parameters = PhysicalParameters()
+physical_parameters = cpu.prameters.PhysicalParameters()
 physical_parameters.lattice_viscosity = 0.0005
 
-runtime_parameters = RuntimeParameters()
+runtime_parameters = cpu.prameters.RuntimeParameters()
 runtime_parameters.number_of_threads = 4
 runtime_parameters.number_of_timesteps = 1000
 runtime_parameters.timestep_log_interval = 100
diff --git a/Python/poiseuille/simulation.py b/Python/poiseuille/simulation.py
index d107801fa84cfe16d1d7e91d31dc3ff4b8671f02..0f793d0c0bcf81016d6f9d0c65ba3637f6317ffa 100644
--- a/Python/poiseuille/simulation.py
+++ b/Python/poiseuille/simulation.py
@@ -1,35 +1,31 @@
-from pyfluids.cpu import Simulation
-from pyfluids.cpu.boundaryconditions import NoSlipBoundaryCondition
-from pyfluids.cpu.geometry import GbCuboid3D, State
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import RuntimeParameters, GridParameters, PhysicalParameters
-from pyfluids.cpu.writer import Writer, OutputFormat
-
-default_grid_params = GridParameters()
+from pyfluids import cpu
+
+
+default_grid_params = cpu.parameters.GridParameters()
 default_grid_params.node_distance = 10 / 32
 default_grid_params.number_of_nodes_per_direction = [8, 8, 32]
 default_grid_params.blocks_per_direction = [1, 1, 4]
 default_grid_params.periodic_boundary_in_x1 = True
 default_grid_params.periodic_boundary_in_x2 = True
 
-default_physical_params = PhysicalParameters()
+default_physical_params = cpu.parameters.PhysicalParameters()
 default_physical_params.lattice_viscosity = 0.005
 
-default_runtime_params = RuntimeParameters()
+default_runtime_params = cpu.parameters.RuntimeParameters()
 default_runtime_params.number_of_threads = 4
 default_runtime_params.number_of_timesteps = 10000
 default_runtime_params.timestep_log_interval = 1000
 
-default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
 default_kernel.use_forcing = True
 default_kernel.forcing_in_x1 = 1e-8
 
-default_writer = Writer()
+default_writer = cpu.writer.Writer()
 default_writer.output_path = "./output"
-default_writer.output_format = OutputFormat.BINARY
+default_writer.output_format = cpu.writer.OutputFormat.BINARY
 
 
-default_kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+default_kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
 default_kernel.use_forcing = True
 default_kernel.forcing_in_x1 = 1e-8
 
@@ -39,7 +35,7 @@ def run_simulation(physical_params=default_physical_params,
                    runtime_params=default_runtime_params,
                    kernel=default_kernel,
                    writer=default_writer):
-    simulation = Simulation()
+    simulation = cpu.Simulation()
 
     simulation.set_kernel_config(kernel)
     simulation.set_physical_parameters(physical_params)
@@ -47,11 +43,11 @@ def run_simulation(physical_params=default_physical_params,
     simulation.set_runtime_parameters(runtime_params)
     simulation.set_writer(writer)
 
-    no_slip_bc = NoSlipBoundaryCondition()
+    no_slip_bc = cpu.boundaryconditions.NoSlipBoundaryCondition()
 
     block_thickness = 3 * grid_params.node_distance
     simulation.add_object(
-        GbCuboid3D(
+        cpu.geometry.GbCuboid3D(
             grid_params.bounding_box.min_x1 - block_thickness,
             grid_params.bounding_box.min_x2 - block_thickness,
             grid_params.bounding_box.min_x3 - block_thickness,
@@ -59,10 +55,10 @@ def run_simulation(physical_params=default_physical_params,
             grid_params.bounding_box.max_x2 + block_thickness,
             grid_params.bounding_box.min_x3),
         no_slip_bc,
-        State.SOLID, "/geo/addWallZMin")
+        cpu.geometry.State.SOLID, "/geo/addWallZMin")
 
     simulation.add_object(
-        GbCuboid3D(
+        cpu.geometry.GbCuboid3D(
             grid_params.bounding_box.min_x1 - block_thickness,
             grid_params.bounding_box.min_x2 - block_thickness,
             grid_params.bounding_box.max_x3,
@@ -70,7 +66,7 @@ def run_simulation(physical_params=default_physical_params,
             grid_params.bounding_box.max_x2 + block_thickness,
             grid_params.bounding_box.max_x3 + block_thickness),
         no_slip_bc,
-        State.SOLID, "/geo/addWallZMax")
+        cpu.geometry.State.SOLID, "/geo/addWallZMax")
 
     simulation.run_simulation()
 
diff --git a/Python/poiseuille/test_poiseuille_l2.py b/Python/poiseuille/test_poiseuille_l2.py
index 93aa2600d5260dea7e72f3aa98db7334fe5285c6..1c6cdbd0116301db7552aa10092a2da97d1c81e8 100644
--- a/Python/poiseuille/test_poiseuille_l2.py
+++ b/Python/poiseuille/test_poiseuille_l2.py
@@ -5,8 +5,7 @@ import unittest
 import matplotlib.pyplot as plt
 import numpy as np
 import pyvista as pv
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import GridParameters, PhysicalParameters, RuntimeParameters
+from pyfluids import cpu
 from scipy import stats
 
 from errors import normalized_l2_error
@@ -33,13 +32,13 @@ class TestPoiseuilleFlow(unittest.TestCase):
         self.skipTest("This test is not implemented correctly yet")
         plt.ion()
 
-        physical_params = PhysicalParameters()
+        physical_params = cpu.parameters.PhysicalParameters()
 
-        runtime_params = RuntimeParameters()
+        runtime_params = cpu.parameters.RuntimeParameters()
         runtime_params.number_of_threads = os.cpu_count()
         runtime_params.timestep_log_interval = 10000
 
-        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
         kernel.use_forcing = True
 
         normalized_l2_errors = []
@@ -140,7 +139,7 @@ def get_heights_from_indices(mesh, indices):
 
 
 def create_grid_params_with_nodes_in_column(nodes_in_column, delta_x):
-    grid_params = GridParameters()
+    grid_params = cpu.parameters.GridParameters()
     grid_params.node_distance = delta_x
     grid_params.number_of_nodes_per_direction = [1, 1, nodes_in_column]
     grid_params.blocks_per_direction = [1, 1, 8]
diff --git a/Python/tests/test_acousticscaling.py b/Python/tests/test_acousticscaling.py
index 6413123a80db8c5882fcf1dbe6f72a1f5438736c..e37b57c98e625f642ec000ec7a8ea83fe0df204c 100644
--- a/Python/tests/test_acousticscaling.py
+++ b/Python/tests/test_acousticscaling.py
@@ -1,9 +1,7 @@
 import unittest
 from typing import List
 
-from pyfluids.cpu.kernel import LBMKernel, KernelType
-from pyfluids.cpu.parameters import GridParameters, PhysicalParameters, RuntimeParameters
-
+from pyfluids import cpu
 from acousticscaling import OneDirectionalAcousticScaling
 
 
@@ -58,18 +56,18 @@ class OneDirectionalAcousticScalingTest(unittest.TestCase):
         self.assertEqual(self.grid_params.periodic_boundary_in_x2, actual_grid_params.periodic_boundary_in_x2)
         self.assertEqual(self.grid_params.periodic_boundary_in_x3, actual_grid_params.periodic_boundary_in_x3)
 
-    def assert_physical_params_scaled_by_factor(self, actual_params: PhysicalParameters, factor: int):
+    def assert_physical_params_scaled_by_factor(self, actual_params: cpu.parameters.PhysicalParameters, factor: int):
         self.assertEqual(self.physical_params.lattice_viscosity * factor, actual_params.lattice_viscosity)
         self.assertEqual(self.physical_params.bulk_viscosity_factor, actual_params.bulk_viscosity_factor)
 
-    def assert_runtime_params_scaled_by_factor(self, actual_params: RuntimeParameters, factor: int):
+    def assert_runtime_params_scaled_by_factor(self, actual_params: cpu.parameters.RuntimeParameters, factor: int):
         self.assertEqual(self.runtime_params.number_of_timesteps * factor, actual_params.number_of_timesteps)
         self.assertEqual(self.runtime_params.number_of_threads, actual_params.number_of_threads)
         self.assertEqual(self.runtime_params.timestep_log_interval, actual_params.timestep_log_interval)
 
-    def assert_kernel_forcing_scaled_by_factor(self, actual_kernel: LBMKernel, factor: int):
+    def assert_kernel_forcing_scaled_by_factor(self, actual_kernel: cpu.kernel.LBMKernel, factor: int):
         self.assertEqual(self.kernel.type, actual_kernel.type)
-        self.assertEqual(self.kernel.use_forcing, actual_kernel.use_forcing)
+        self.assertEqual(self.kernel.use_forcing, actual_kernel.cpu.parameters.use_forcing)
         self.assertAlmostEqual(self.kernel.forcing_in_x1 / factor, actual_kernel.forcing_in_x1)
         self.assertAlmostEqual(self.kernel.forcing_in_x2, actual_kernel.forcing_in_x2)
         self.assertAlmostEqual(self.kernel.forcing_in_x3, actual_kernel.forcing_in_x3)
@@ -80,14 +78,14 @@ class OneDirectionalAcousticScalingTest(unittest.TestCase):
 
     @staticmethod
     def make_kernel():
-        kernel = LBMKernel(KernelType.CompressibleCumulantFourthOrderViscosity)
+        kernel = cpu.kernel.LBMKernel(cpu.kernel.KernelType.CompressibleCumulantFourthOrderViscosity)
         kernel.use_forcing = True
         kernel.forcing_in_x1 = 5e-10
         return kernel
 
     @staticmethod
     def make_runtime_params():
-        runtime_params = RuntimeParameters()
+        runtime_params = cpu.parameters.RuntimeParameters()
         runtime_params.number_of_threads = 4
         runtime_params.number_of_timesteps = 4_000_000
         runtime_params.timestep_log_interval = 1_000_000
@@ -95,13 +93,13 @@ class OneDirectionalAcousticScalingTest(unittest.TestCase):
 
     @staticmethod
     def make_physical_params():
-        physical_params = PhysicalParameters()
+        physical_params = cpu.parameters.PhysicalParameters()
         physical_params.lattice_viscosity = 1e-4
         return physical_params
 
     @staticmethod
     def make_grid_params():
-        grid_params = GridParameters()
+        grid_params = cpu.parameters.GridParameters()
         grid_params.node_distance = 1
         grid_params.number_of_nodes_per_direction = [1, 1, 16]
         grid_params.blocks_per_direction = [1, 1, 16]
diff --git a/Python/tests/test_boundaryconditions.py b/Python/tests/test_boundaryconditions.py
index e004ddfa21c78ea3d63a89f5dbc3bd7438a18ff1..c2f6c494febe0708754f686bb804fa6a3baec6f1 100644
--- a/Python/tests/test_boundaryconditions.py
+++ b/Python/tests/test_boundaryconditions.py
@@ -1,5 +1,5 @@
 import unittest
-from pyfluids.cpu.boundaryconditions import *
+from pyfluids import cpu
 
 
 class BoundaryConditionsTest(unittest.TestCase):
@@ -8,13 +8,13 @@ class BoundaryConditionsTest(unittest.TestCase):
         """
         Should be able to create NoSlipBoundaryCondition
         """
-        sut = NoSlipBoundaryCondition()
+        sut = cpu.boundaryconditions.NoSlipBoundaryCondition()
 
     def test__can_create_velocity_bc(self):
         """
         Should be able to create VelocityBoundaryCondition
         """
-        sut = VelocityBoundaryCondition()
+        sut = cpu.boundaryconditions.VelocityBoundaryCondition()
 
     def test__can_create_velocity_bc_with_directions_function_and_time(self):
         """
@@ -24,7 +24,7 @@ class BoundaryConditionsTest(unittest.TestCase):
 
         parser = Parser()
         parser.expression = "1"
-        sut = VelocityBoundaryCondition(True, True, True, parser, 0, 1)
+        sut = cpu.boundaryconditions.VelocityBoundaryCondition(True, True, True, parser, 0, 1)
 
     def test__can_create_velocity_bc_with_directions__function_per_direction__and__time(self):
         """
@@ -40,7 +40,7 @@ class BoundaryConditionsTest(unittest.TestCase):
 
         f3 = Parser()
         f3.expression = "1"
-        sut = VelocityBoundaryCondition(True, True, True, f1, f2, f3, 0, 1)
+        sut = cpu.boundaryconditions.VelocityBoundaryCondition(True, True, True, f1, f2, f3, 0, 1)
 
     def test__can_create_velocity_bc_with_speeds_and_times_per_direction(self):
         """
@@ -51,11 +51,11 @@ class BoundaryConditionsTest(unittest.TestCase):
         start2, end2 = 1, 2
         start3, end3 = 2, 3
 
-        sut = VelocityBoundaryCondition(vx1, start1, end1, vx2, start2, end2, vx3, start3, end3)
+        sut = cpu.boundaryconditions.VelocityBoundaryCondition(vx1, start1, end1, vx2, start2, end2, vx3, start3, end3)
 
     def test__can_create_non_reflecting_outflow(self):
         """
         Should be able to create NonReflectingOutflow
         """
 
-        sut = NonReflectingOutflow()
+        sut = cpu.boundaryconditions.NonReflectingOutflow()
diff --git a/Python/tests/test_geometry.py b/Python/tests/test_geometry.py
index 5bb89eb245b6055653b78fde381da050d402b0cc..b0ac2f64b19c4002a208a9e85e952727b4ff2ee0 100644
--- a/Python/tests/test_geometry.py
+++ b/Python/tests/test_geometry.py
@@ -1,6 +1,6 @@
 import unittest
 
-from pyfluids.cpu.geometry import *
+from pyfluids import cpu
 
 
 class TestGeometry(unittest.TestCase):
@@ -9,7 +9,7 @@ class TestGeometry(unittest.TestCase):
         """
         WHEN setting point coordinates in constructor THEN point should have coordinates
         """
-        sut = GbPoint3D(4, 8, 3)
+        sut = cpu.geometry.GbPoint3D(4, 8, 3)
 
         self.assertEqual(sut.x1, 4)
         self.assertEqual(sut.x2, 8)
@@ -19,7 +19,7 @@ class TestGeometry(unittest.TestCase):
         """
         WHEN setting point coordinates THEN point should have coordinates
         """
-        sut = GbPoint3D()
+        sut = cpu.geometry.GbPoint3D()
 
         sut.x1 = 4
         sut.x2 = 8
@@ -33,10 +33,10 @@ class TestGeometry(unittest.TestCase):
         """
         WHEN setting line points THEN line should have points
         """
-        sut = GbLine3D()
+        sut = cpu.geometry.GbLine3D()
 
-        point1 = GbPoint3D()
-        point2 = GbPoint3D()
+        point1 = cpu.geometry.GbPoint3D()
+        point2 = cpu.geometry.GbPoint3D()
         sut.point1 = point1
         sut.point2 = point2
 
diff --git a/Python/tests/test_kernel.py b/Python/tests/test_kernel.py
index 8f58a1c869f9e292856268d43245a75f1dcfe213..4736f4ebe41f3c93ca433c1577bf38d9d95ce6cb 100644
--- a/Python/tests/test_kernel.py
+++ b/Python/tests/test_kernel.py
@@ -1,12 +1,12 @@
 import unittest
 
-from pyfluids.cpu.kernel import LBMKernel, KernelType
+from pyfluids import cpu
 
 
 class TestLBMKernel(unittest.TestCase):
 
     def setUp(self) -> None:
-        self.sut = LBMKernel(KernelType.BGK)
+        self.sut = cpu.kernel.LBMKernel(cpu.kernel.KernelType.BGK)
 
     def test_lbm_kernel__when_use_forcing_set_to_true__use_forcing_should_be_true(self) -> None:
         """
@@ -57,4 +57,4 @@ class TestLBMKernel(unittest.TestCase):
         """
 
         actual = self.sut.type
-        self.assertEqual(KernelType.BGK, actual)
+        self.assertEqual(cpu.kernel.KernelType.BGK, actual)