diff --git a/PythonAPI.md b/PythonAPI.md
index ac7c7ad6..002965d9 100644
--- a/PythonAPI.md
+++ b/PythonAPI.md
@@ -60,7 +60,7 @@ These are imported directly from the **root** `viennals` module.
 
 ### **Enums**
 - `LogLevel`
-- `IntegrationSchemeEnum`
+- `SpatialSchemeEnum`
 - `BooleanOperationEnum`
 - `CurvatureEnum`
 - `FeatureDetectionEnum`
diff --git a/examples/Deposition/Deposition.cpp b/examples/Deposition/Deposition.cpp
index 7b087f89..52cabd78 100644
--- a/examples/Deposition/Deposition.cpp
+++ b/examples/Deposition/Deposition.cpp
@@ -113,8 +113,7 @@ int main() {
   advectionKernel.setVelocityField(velocities);
   // advectionKernel.setAdvectionTime(4.);
   unsigned counter = 1;
-  advectionKernel.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::WENO_5TH_ORDER);
+  advectionKernel.setSpatialScheme(ls::SpatialSchemeEnum::WENO_5TH_ORDER);
   for (NumericType time = 0; time < 4.;
        time += advectionKernel.getAdvectedTime()) {
     advectionKernel.apply();
diff --git a/examples/Epitaxy/Epitaxy.cpp b/examples/Epitaxy/Epitaxy.cpp
index 5b8cb2a5..92ab5ede 100644
--- a/examples/Epitaxy/Epitaxy.cpp
+++ b/examples/Epitaxy/Epitaxy.cpp
@@ -21,7 +21,7 @@ class epitaxy final : public VelocityField<T> {
   static constexpr double high = 1.0;
 
 public:
-  epitaxy(std::vector<double> vel) : velocities(vel){};
+  epitaxy(std::vector<double> vel) : velocities(vel) {};
 
   double getScalarVelocity(const std::array<T, 3> & /*coordinate*/,
                            int material, const std::array<T, 3> &normal,
@@ -82,8 +82,8 @@ int main(int argc, char *argv[]) {
   auto velocityField = SmartPointer<epitaxy>::New(std::vector<double>{0., -.5});
 
   Advect<T, D> advectionKernel(levelSets, velocityField);
-  advectionKernel.setIntegrationScheme(
-      IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
+  advectionKernel.setSpatialScheme(
+      SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
   advectionKernel.setAdvectionTime(.5);
 
   Timer timer;
diff --git a/examples/Epitaxy/Epitaxy.py b/examples/Epitaxy/Epitaxy.py
index d5990c12..84f47922 100644
--- a/examples/Epitaxy/Epitaxy.py
+++ b/examples/Epitaxy/Epitaxy.py
@@ -5,6 +5,7 @@
 D = 3
 vls.setNumThreads(8)
 
+
 class EpitaxyVelocity(vls.VelocityField):
     def __init__(self, velocities):
         super().__init__()
@@ -30,13 +31,15 @@ def getScalarVelocity(self, coord, material, normal, point_id):
         return vel * mat_vel
 
     def getVectorVelocity(self, coord, material, normal, point_id):
-        return (0., 0., 0.)
+        return (0.0, 0.0, 0.0)
+
 
 def write_surface(domain, filename):
     mesh = vls.Mesh()
     vls.ToSurfaceMesh(domain, mesh).apply()
     vls.VTKWriter(mesh, filename).apply()
 
+
 def main():
     # Set dimension to 3D
     vls.setDimension(D)
@@ -50,7 +53,7 @@ def main():
     boundary_conditions = [
         vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY,
         vls.BoundaryConditionEnum.REFLECTIVE_BOUNDARY,
-        vls.BoundaryConditionEnum.INFINITE_BOUNDARY
+        vls.BoundaryConditionEnum.INFINITE_BOUNDARY,
     ]
 
     # Create Geometry
@@ -80,7 +83,9 @@ def main():
     for ls in level_sets:
         advection.insertNextLevelSet(ls)
     advection.setVelocityField(velocity_field)
-    advection.setIntegrationScheme(vls.IntegrationSchemeEnum.STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
+    advection.setSpatialScheme(
+        vls.SpatialSchemeEnum.STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER
+    )
     advection.setAdvectionTime(0.5)
 
     start_time = time.time()
@@ -92,5 +97,6 @@ def main():
 
     write_surface(substrate, "epitaxy.vtp")
 
+
 if __name__ == "__main__":
     main()
diff --git a/examples/PeriodicBoundary/PeriodicBoundary.cpp b/examples/PeriodicBoundary/PeriodicBoundary.cpp
index 6d47107d..c36d0078 100644
--- a/examples/PeriodicBoundary/PeriodicBoundary.cpp
+++ b/examples/PeriodicBoundary/PeriodicBoundary.cpp
@@ -90,8 +90,8 @@ int main() {
   ls::Advect<double, D> advectionKernel;
   advectionKernel.insertNextLevelSet(substrate);
   advectionKernel.setVelocityField(velocities);
-  advectionKernel.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::ENGQUIST_OSHER_2ND_ORDER);
+  advectionKernel.setSpatialScheme(
+      ls::SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER);
 
   // Now advect the level set 50 times, outputting every
   // advection step. Save the physical time that
diff --git a/examples/SquareEtch/SquareEtch.cpp b/examples/SquareEtch/SquareEtch.cpp
index a63d8780..1bfc6534 100644
--- a/examples/SquareEtch/SquareEtch.cpp
+++ b/examples/SquareEtch/SquareEtch.cpp
@@ -171,21 +171,21 @@ int main() {
   if (useAnalyticalVelocity) {
     advectionKernel.setVelocityField(analyticalVelocities);
     // Analytical velocity fields and dissipation coefficients
-    // can only be used with this integration scheme
-    advectionKernel.setIntegrationScheme(
-        // ls::IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
-        ls::IntegrationSchemeEnum::WENO_5TH_ORDER);
+    // can only be used with this spatial discretization scheme
+    advectionKernel.setSpatialScheme(
+        // ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
+        ls::SpatialSchemeEnum::WENO_5TH_ORDER);
   } else {
     // for numerical velocities, just use the default
-    // integration scheme, which is not accurate for certain
+    // spatial discretization scheme, which is not accurate for certain
     // velocity functions but very fast
     advectionKernel.setVelocityField(velocities);
 
     // For coordinate independent velocity functions
     // this numerical scheme is superior though.
     // However, it is slower.
-    // advectionKernel.setIntegrationScheme(
-    //     ls::IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
+    // advectionKernel.setSpatialScheme(
+    //     ls::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
   }
 
   // advect the level set until 50s have passed
diff --git a/examples/SquareEtch/SquareEtch.py b/examples/SquareEtch/SquareEtch.py
index c0f190af..347a50b5 100644
--- a/examples/SquareEtch/SquareEtch.py
+++ b/examples/SquareEtch/SquareEtch.py
@@ -3,6 +3,7 @@
 
 vls.Logger.setLogLevel(vls.LogLevel.INFO)
 
+
 # 1. Define the Velocity Field
 # We inherit from viennals.VelocityField to define custom logic in Python
 class EtchingField(vls.VelocityField):
@@ -10,12 +11,13 @@ def getScalarVelocity(self, coordinate, material, normalVector, pointId):
         if material == 2:
             return -0.1
         if material == 1:
-            return -1.
+            return -1.0
         return 0.0
 
     def getVectorVelocity(self, coordinate, material, normalVector, pointId):
         return [0.0] * len(coordinate)
 
+
 # 1. Define the Velocity Field
 # We inherit from viennals.VelocityField to define custom logic in Python
 class DepositionField(vls.VelocityField):
@@ -25,25 +27,29 @@ def getScalarVelocity(self, coordinate, material, normalVector, pointId):
     def getVectorVelocity(self, coordinate, material, normalVector, pointId):
         return [0.0] * len(coordinate)
 
+
 def main():
     # 1. Parse Arguments
-    parser = argparse.ArgumentParser(description="Run Square Etch simulation in 2D or 3D.")
+    parser = argparse.ArgumentParser(
+        description="Run Square Etch simulation in 2D or 3D."
+    )
     parser.add_argument(
-        "-D", "--dim", 
-        type=int, 
-        default=2, 
-        choices=[2, 3], 
-        help="Dimension of the simulation (2 or 3). Default is 2."
+        "-D",
+        "--dim",
+        type=int,
+        default=2,
+        choices=[2, 3],
+        help="Dimension of the simulation (2 or 3). Default is 2.",
     )
     args = parser.parse_args()
-    
+
     DIMENSION = args.dim
     vls.setDimension(DIMENSION)
     vls.setNumThreads(8)
 
     extent = 30.0
     gridDelta = 0.47
-    
+
     # Define bounds and boundary conditions
     trenchBottom = -2.0
     bounds = [-extent, extent, -extent, extent]
@@ -73,14 +79,14 @@ def main():
     # 3. Trench Geometry
     # --------------------------------------
     trench = vls.Domain(bounds, boundaryCons, gridDelta)
-    
+
     # Define Box Corners based on dimension
     minCorner = list(origin)
     maxCorner = list(origin)
-    originMask = list (origin)
+    originMask = list(origin)
     for i in range(DIMENSION - 1):
         minCorner[i] = -extent / 1.5
-        maxCorner[i] =  extent / 1.5
+        maxCorner[i] = extent / 1.5
     minCorner[DIMENSION - 1] = trenchBottom
     maxCorner[DIMENSION - 1] = 1.0
     originMask[DIMENSION - 1] = trenchBottom + 1e-9
@@ -89,7 +95,9 @@ def main():
     vls.MakeGeometry(trench, box).apply()
 
     # Subtract trench from substrate (Relative Complement)
-    vls.BooleanOperation(substrate, trench, vls.BooleanOperationEnum.RELATIVE_COMPLEMENT).apply()
+    vls.BooleanOperation(
+        substrate, trench, vls.BooleanOperationEnum.RELATIVE_COMPLEMENT
+    ).apply()
 
     # 4. Create the Mask Layer
     # The mask prevents etching at the bottom of the trench in specific configurations,
@@ -97,7 +105,7 @@ def main():
     mask = vls.Domain(bounds, boundaryCons, gridDelta)
 
     vls.MakeGeometry(mask, vls.Plane(originMask, downNormal)).apply()
-    
+
     # Intersect with substrate geometry
     vls.BooleanOperation(mask, substrate, vls.BooleanOperationEnum.INTERSECT).apply()
 
@@ -105,11 +113,11 @@ def main():
     print(f"Running in {DIMENSION}D.")
     print("Extracting initial meshes...")
     mesh = vls.Mesh()
-    
+
     # Save substrate
     vls.ToSurfaceMesh(substrate, mesh).apply()
     vls.VTKWriter(mesh, f"substrate-{DIMENSION}D.vtp").apply()
-    
+
     # Save mask
     vls.ToSurfaceMesh(mask, mesh).apply()
     vls.VTKWriter(mesh, f"mask-{DIMENSION}D.vtp").apply()
@@ -117,7 +125,6 @@ def main():
     print("Creating new layer...")
     polymer = vls.Domain(substrate)
 
-
     # 6. Setup Advection
     deposition = DepositionField()
     etching = EtchingField()
@@ -142,9 +149,9 @@ def main():
 
     advectionKernel.setSaveAdvectionVelocities(True)
     advectionKernel.setVelocityField(etching)
-    
-    # Use default integration scheme (Lax Friedrichs 1st order) as in the C++ else branch
-    # advectionKernel.setIntegrationScheme(viennals.IntegrationSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER)
+
+    # Use default spatial discretization scheme (Lax Friedrichs 1st order) as in the C++ else branch
+    # advectionKernel.setSpatialScheme(viennals.SpatialSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER)
 
     # 7. Time Loop
     finalTime = 10.0
@@ -153,7 +160,7 @@ def main():
 
     # The advect kernel calculates stable time steps automatically.
     # We call apply() repeatedly until we reach finalTime.
-    
+
     # Note: In the C++ example, the loop structure is:
     # for (double time = 0.; time < finalTime; time += advectionKernel.getAdvectedTime())
     # This implies one step is taken, then time is updated.
@@ -161,19 +168,21 @@ def main():
     # Save initial state
     vls.ToSurfaceMesh(polymer, mesh).apply()
     vls.VTKWriter(mesh, f"numerical-{DIMENSION}D-0.vtp").apply()
-    
+
     while currentTime < finalTime:
         advectionKernel.apply()
-        
+
         stepTime = advectionKernel.getAdvectedTime()
         currentTime += stepTime
-        
-        print(f"Advection step: {counter}, time: {stepTime:.4f} (Total: {currentTime:.4f})")
-        
+
+        print(
+            f"Advection step: {counter}, time: {stepTime:.4f} (Total: {currentTime:.4f})"
+        )
+
         # Export result
         vls.ToSurfaceMesh(polymer, mesh).apply()
         vls.VTKWriter(mesh, f"numerical-{DIMENSION}D-{counter}.vtp").apply()
-        
+
         counter += 1
 
     print(f"\nNumber of Advection steps taken: {counter}")
@@ -182,5 +191,6 @@ def main():
     vls.ToSurfaceMesh(polymer, mesh).apply()
     vls.VTKWriter(mesh, f"final-{DIMENSION}D.vtp").apply()
 
+
 if __name__ == "__main__":
-    main()
\ No newline at end of file
+    main()
diff --git a/include/viennals/lsAdvect.hpp b/include/viennals/lsAdvect.hpp
index 4aa2603b..3ab9d5b1 100644
--- a/include/viennals/lsAdvect.hpp
+++ b/include/viennals/lsAdvect.hpp
@@ -16,7 +16,7 @@
 #include <lsMarkVoidPoints.hpp>
 #include <lsReduce.hpp>
 
-// Integration schemes
+// Spatial discretization schemes
 #include <lsEngquistOsher.hpp>
 #include <lsLaxFriedrichs.hpp>
 #include <lsLocalLaxFriedrichs.hpp>
@@ -38,9 +38,9 @@ namespace viennals {
 
 using namespace viennacore;
 
-/// Enumeration for the different Integration schemes
+/// Enumeration for the different spatial discretization schemes
 /// used by the advection kernel
-enum struct IntegrationSchemeEnum : unsigned {
+enum struct SpatialSchemeEnum : unsigned {
   ENGQUIST_OSHER_1ST_ORDER = 0,
   ENGQUIST_OSHER_2ND_ORDER = 1,
   LAX_FRIEDRICHS_1ST_ORDER = 2,
@@ -54,6 +54,9 @@ enum struct IntegrationSchemeEnum : unsigned {
   WENO_5TH_ORDER = 10
 };
 
+// Legacy naming (deprecated, will be removed in future versions)
+using IntegrationSchemeEnum = SpatialSchemeEnum;
+
 /// This class is used to advance level sets over time.
 /// Level sets are passed to the constructor in a std::vector, with
 /// the last element being the level set to advect, or "top level set", while
@@ -71,8 +74,7 @@ template <class T, int D> class Advect {
 protected:
   std::vector<SmartPointer<Domain<T, D>>> levelSets;
   SmartPointer<VelocityField<T>> velocities = nullptr;
-  IntegrationSchemeEnum integrationScheme =
-      IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER;
+  SpatialSchemeEnum spatialScheme = SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER;
   double timeStepRatio = 0.4999;
   double dissipationAlpha = 1.0;
   bool calculateNormalVectors = true;
@@ -213,12 +215,12 @@ template <class T, int D> class Advect {
     auto &newDomain = newlsDomain->getDomain();
     auto &domain = levelSets.back()->getDomain();
 
-    // Determine cutoff and width based on integration scheme to avoid immediate
-    // re-expansion
+    // Determine cutoff and width based on discretization scheme to avoid
+    // immediate re-expansion
     T cutoff = 1.0;
     int finalWidth = 2;
-    if (integrationScheme ==
-        IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    if (spatialScheme ==
+        SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       cutoff = 1.5;
       finalWidth = 3;
     }
@@ -396,10 +398,11 @@ template <class T, int D> class Advect {
   }
 
   /// Internal function used to calculate the deltas to be applied to the LS
-  /// values from the given velocities and the integration scheme to be used.
-  /// This function fills up the storedRates to be used when moving the LS
-  template <class IntegrationSchemeType>
-  double integrateTime(IntegrationSchemeType IntegrationScheme,
+  /// values from the given velocities and the spatial discretization scheme to
+  /// be used. This function fills up the storedRates to be used when moving the
+  /// LS
+  template <class DiscretizationSchemeType>
+  double integrateTime(DiscretizationSchemeType spatialScheme,
                        double maxTimeStep) {
 
     auto &topDomain = levelSets.back()->getDomain();
@@ -457,7 +460,7 @@ template <class T, int D> class Advect {
         iterators.emplace_back(ls->getDomain());
       }
 
-      IntegrationSchemeType scheme(IntegrationScheme);
+      DiscretizationSchemeType scheme(spatialScheme);
 
       for (ConstSparseIterator it(topDomain, startVector);
            it.getStartIndices() < endVector; ++it) {
@@ -580,71 +583,67 @@ template <class T, int D> class Advect {
     return maxTimeStep;
   }
 
-  /// This function applies the integration scheme and calculates the rates and
-  /// the maximum time step, but it does **not** move the surface.
+  /// This function applies the discretization scheme and calculates the rates
+  /// and the maximum time step, but it does **not** move the surface.
   void computeRates(double maxTimeStep = std::numeric_limits<double>::max()) {
     prepareLS();
-    if (integrationScheme == IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
+    if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
       auto is = lsInternal::EngquistOsher<T, D, 1>(levelSets.back(), velocities,
                                                    calculateNormalVectors);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
       auto is = lsInternal::EngquistOsher<T, D, 2>(levelSets.back(), velocities,
                                                    calculateNormalVectors);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
       auto alphas = findGlobalAlphas();
       auto is = lsInternal::LaxFriedrichs<T, D, 1>(levelSets.back(), velocities,
                                                    dissipationAlpha, alphas,
                                                    calculateNormalVectors);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
       auto alphas = findGlobalAlphas();
       auto is = lsInternal::LaxFriedrichs<T, D, 2>(levelSets.back(), velocities,
                                                    dissipationAlpha, alphas,
                                                    calculateNormalVectors);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::
-                   LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
       auto is = lsInternal::LocalLaxFriedrichsAnalytical<T, D, 1>(
           levelSets.back(), velocities);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       auto is = lsInternal::LocalLocalLaxFriedrichs<T, D, 1>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
       auto is = lsInternal::LocalLocalLaxFriedrichs<T, D, 2>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       auto is = lsInternal::LocalLaxFriedrichs<T, D, 1>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
       auto is = lsInternal::LocalLaxFriedrichs<T, D, 2>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       auto is = lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
-    } else if (integrationScheme == IntegrationSchemeEnum::WENO_5TH_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
       // Instantiate WENO5 with order 3 (neighbors +/- 3)
       auto is = lsInternal::WENO5<T, D, 3>(levelSets.back(), velocities,
                                            dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
     } else {
-      VIENNACORE_LOG_ERROR("Advect: Integration scheme not found.");
+      VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
       currentTimeStep = -1.;
     }
   }
@@ -770,7 +769,7 @@ template <class T, int D> class Advect {
   }
 
   /// internal function used as a wrapper to call specialized integrateTime
-  /// with the chosen integrationScheme
+  /// with the chosen spatial discretization scheme
   virtual double advect(double maxTimeStep) {
     prepareLS();
 
@@ -782,8 +781,8 @@ template <class T, int D> class Advect {
     rebuildLS();
 
     // Adjust all level sets below the advected one
-    if (integrationScheme !=
-        IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    if (spatialScheme !=
+        SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       for (unsigned i = 0; i < levelSets.size() - 1; ++i) {
         BooleanOperation<T, D>(levelSets[i], levelSets.back(),
                                BooleanOperationEnum::INTERSECT)
@@ -896,10 +895,16 @@ template <class T, int D> class Advect {
   /// Get whether normal vectors were calculated.
   bool getCalculateNormalVectors() const { return calculateNormalVectors; }
 
-  /// Set which integration scheme should be used out of the ones specified
-  /// in IntegrationSchemeEnum.
+  /// Set which spatial discretization scheme should be used out of the ones
+  /// specified in SpatialSchemeEnum.
+  void setSpatialScheme(SpatialSchemeEnum scheme) { spatialScheme = scheme; }
+
+  // Deprecated, will be remove in future versions: use setSpatialScheme instead
   void setIntegrationScheme(IntegrationSchemeEnum scheme) {
-    integrationScheme = scheme;
+    VIENNACORE_LOG_WARNING(
+        "Advect::setIntegrationScheme is deprecated and will be removed in "
+        "future versions. Use setSpatialScheme instead.");
+    spatialScheme = scheme;
   }
 
   /// Set the alpha dissipation coefficient.
@@ -915,8 +920,8 @@ template <class T, int D> class Advect {
   /// be translated to the advected LS. Defaults to true.
   void setUpdatePointData(bool update) { updatePointData = update; }
 
-  // Prepare the levelset for advection, based on the provided integration
-  // scheme.
+  // Prepare the levelset for advection, based on the provided spatial
+  // discretization scheme.
   void prepareLS() {
     // check whether a level set and velocities have been given
     if (levelSets.empty()) {
@@ -925,46 +930,42 @@ template <class T, int D> class Advect {
       return;
     }
 
-    if (integrationScheme == IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
+    if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER) {
       lsInternal::EngquistOsher<T, D, 1>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER) {
       lsInternal::EngquistOsher<T, D, 2>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::LaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER) {
       lsInternal::LaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::
-                   LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER) {
       lsInternal::LocalLaxFriedrichsAnalytical<T, D, 1>::prepareLS(
           levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::LocalLocalLaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
       lsInternal::LocalLocalLaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::LocalLaxFriedrichs<T, D, 1>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER) {
       lsInternal::LocalLaxFriedrichs<T, D, 2>::prepareLS(levelSets.back());
-    } else if (integrationScheme ==
-               IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    } else if (spatialScheme ==
+               SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>::prepareLS(
           levelSets.back());
-    } else if (integrationScheme == IntegrationSchemeEnum::WENO_5TH_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
       // WENO5 requires a stencil radius of 3 (template parameter 3)
       lsInternal::WENO5<T, D, 3>::prepareLS(levelSets.back());
-    } else if (integrationScheme == IntegrationSchemeEnum::WENO_5TH_ORDER) {
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_5TH_ORDER) {
       // WENO5 requires a stencil radius of 3 (template parameter 3)
       lsInternal::WENO5<T, D, 3>::prepareLS(levelSets.back());
     } else {
-      VIENNACORE_LOG_ERROR("Advect: Integration scheme not found.");
+      VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
     }
   }
 
diff --git a/include/viennals/lsAdvectForwardEuler.hpp b/include/viennals/lsAdvectForwardEuler.hpp
index 4c608530..6530f64f 100644
--- a/include/viennals/lsAdvectForwardEuler.hpp
+++ b/include/viennals/lsAdvectForwardEuler.hpp
@@ -11,4 +11,4 @@ template <class T, int D> class AdvectForwardEuler : public Advect<T, D> {
 
 PRECOMPILE_PRECISION_DIMENSION(AdvectForwardEuler);
 
-} // namespace viennals
\ No newline at end of file
+} // namespace viennals
diff --git a/include/viennals/lsAdvectRungeKutta3.hpp b/include/viennals/lsAdvectRungeKutta3.hpp
index 34209380..964460f2 100644
--- a/include/viennals/lsAdvectRungeKutta3.hpp
+++ b/include/viennals/lsAdvectRungeKutta3.hpp
@@ -102,8 +102,8 @@ template <class T, int D> class AdvectRungeKutta3 : public Advect<T, D> {
     Base::rebuildLS();
 
     // Adjust lower layers
-    if (Base::integrationScheme !=
-        IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
+    if (Base::spatialScheme !=
+        SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       for (unsigned i = 0; i < levelSets.size() - 1; ++i) {
         BooleanOperation<T, D>(levelSets[i], levelSets.back(),
                                BooleanOperationEnum::INTERSECT)
diff --git a/include/viennals/lsEngquistOsher.hpp b/include/viennals/lsEngquistOsher.hpp
index f31bcfac..f61c5d0a 100644
--- a/include/viennals/lsEngquistOsher.hpp
+++ b/include/viennals/lsEngquistOsher.hpp
@@ -12,8 +12,8 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Engquist-Osher integration scheme based on the
-/// upwind integration scheme. Offers high performance
+/// Engquist-Osher spatial discretization scheme based on the
+/// upwind spatial discretization scheme. Offers high performance
 /// but lower accuracy for complex velocity fields.
 template <class T, int D, int order> class EngquistOsher {
   SmartPointer<viennals::Domain<T, D>> levelSet;
@@ -69,7 +69,7 @@ template <class T, int D, int order> class EngquistOsher {
       T diffPos = (phiPos - phi0) / deltaPos;
       T diffNeg = (phiNeg - phi0) / deltaNeg;
 
-      if (order == 2) { // if second order time integration scheme is used
+      if (order == 2) { // if second order spatial discretization scheme is used
         const T deltaPosPos = 2 * gridDelta;
         const T deltaNegNeg = -2 * gridDelta;
 
diff --git a/include/viennals/lsFiniteDifferences.hpp b/include/viennals/lsFiniteDifferences.hpp
index 79cc2ad6..e4258ef7 100644
--- a/include/viennals/lsFiniteDifferences.hpp
+++ b/include/viennals/lsFiniteDifferences.hpp
@@ -154,7 +154,7 @@ class FiniteDifferences {
     if constexpr (scheme == DifferentiationSchemeEnum::FIRST_ORDER) {
       return (values[1] - values[0]) / delta;
     } else if (scheme == DifferentiationSchemeEnum::SECOND_ORDER) {
-      // TODO: implement second order integration here
+      // TODO: implement second order differentiation here
       Logger::getInstance().addError("Second order scheme not implemented!");
       return 0;
     } else if (scheme == DifferentiationSchemeEnum::WENO3) {
@@ -177,8 +177,8 @@ class FiniteDifferences {
     if constexpr (scheme == DifferentiationSchemeEnum::FIRST_ORDER) {
       return (values[2] - values[1]) / delta;
     } else if (scheme == DifferentiationSchemeEnum::SECOND_ORDER) {
-      // TODO: implement second order integration here
-      Logger::getInstance().addError("Seconed order scheme not implemented!");
+      // TODO: implement second order differentiation here
+      Logger::getInstance().addError("Second order scheme not implemented!");
       return 0;
     } else if (scheme == DifferentiationSchemeEnum::WENO3) {
       return weno3(values, delta, true);
diff --git a/include/viennals/lsLaxFriedrichs.hpp b/include/viennals/lsLaxFriedrichs.hpp
index 45c9d929..7320b518 100644
--- a/include/viennals/lsLaxFriedrichs.hpp
+++ b/include/viennals/lsLaxFriedrichs.hpp
@@ -12,7 +12,7 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Lax Friedrichs integration scheme with constant alpha
+/// Lax Friedrichs spatial discretization scheme with constant alpha
 /// value for dissipation. This alpha value should be fitted
 /// based on the results of the advection and passed to the
 /// advection Kernel.
@@ -76,7 +76,7 @@ template <class T, int D, int order> class LaxFriedrichs {
       T diffPos = (phiPos - phi0) / deltaPos;
       T diffNeg = (phiNeg - phi0) / deltaNeg;
 
-      if (order == 2) { // if second order time integration scheme is used
+      if (order == 2) { // if second order spatial discretization scheme is used
         const T deltaPosPos = 2 * gridDelta;
         const T deltaNegNeg = -2 * gridDelta;
 
diff --git a/include/viennals/lsLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLaxFriedrichs.hpp
index 7cb64220..e73a9809 100644
--- a/include/viennals/lsLocalLaxFriedrichs.hpp
+++ b/include/viennals/lsLocalLaxFriedrichs.hpp
@@ -12,7 +12,7 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Lax Friedrichs integration scheme, which uses a first neighbour
+/// Lax Friedrichs spatial discretization scheme, which uses a first neighbour
 /// stencil to calculate the alpha values for all neighbours.
 /// The largest alpha value is then chosen for dissipation.
 /// Slower than lsLocalLocalLaxFriedrichs or lsEngquistOsher
@@ -105,7 +105,7 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
       T diffPos = (phiPos - phi0) / deltaPos;
       T diffNeg = (phiNeg - phi0) / deltaNeg;
 
-      if (order == 2) { // if second order time integration scheme is used
+      if (order == 2) { // if second order spatial discretization scheme is used
         posUnit[i] = 2;
         negUnit[i] = -2;
 
diff --git a/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp b/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
index 02815238..1e276fc4 100644
--- a/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
+++ b/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
@@ -12,11 +12,11 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Lax Friedrichs integration scheme, which uses alpha values
+/// Lax Friedrichs spatial discretization scheme, which uses alpha values
 /// provided by the user in getDissipationAlphas in lsVelocityField.
 /// If it is possible to derive analytical solutions for the velocityField
-/// and the alpha values, this integration scheme should be used and never
-/// otherwise.
+/// and the alpha values, this spatial discretization scheme should be used and
+/// never otherwise.
 template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
   SmartPointer<viennals::Domain<T, D>> levelSet;
   SmartPointer<viennals::VelocityField<T>> velocities;
@@ -104,7 +104,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
       T diffPos = (phiPos - phi0) / deltaPos;
       T diffNeg = (phiNeg - phi0) / deltaNeg;
 
-      if (order == 2) { // if second order time integration scheme is used
+      if (order == 2) { // if second order spatial discretization scheme is used
         posUnit[i] = 2;
         negUnit[i] = -2;
 
diff --git a/include/viennals/lsLocalLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLocalLaxFriedrichs.hpp
index 7d4493e5..4361172f 100644
--- a/include/viennals/lsLocalLocalLaxFriedrichs.hpp
+++ b/include/viennals/lsLocalLocalLaxFriedrichs.hpp
@@ -12,8 +12,8 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Lax Friedrichs integration scheme, which considers only the current
-/// point for alpha calculation. Faster than lsLocalLaxFriedrichs but
+/// Lax Friedrichs spatial discretization scheme, which considers only the
+/// current point for alpha calculation. Faster than lsLocalLaxFriedrichs but
 /// not as accurate.
 template <class T, int D, int order> class LocalLocalLaxFriedrichs {
   SmartPointer<viennals::Domain<T, D>> levelSet;
@@ -79,7 +79,7 @@ template <class T, int D, int order> class LocalLocalLaxFriedrichs {
       T diffPos = (phiPos - phi0) / deltaPos;
       T diffNeg = (phiNeg - phi0) / deltaNeg;
 
-      if (order == 2) { // if second order time integration scheme is used
+      if (order == 2) { // if second order spatial discretization scheme is used
         const T deltaPosPos = 2 * gridDelta;
         const T deltaNegNeg = -2 * gridDelta;
 
diff --git a/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp b/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
index 911701a8..1393cbb9 100644
--- a/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
+++ b/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
@@ -21,7 +21,7 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Stencil Local Lax Friedrichs Integration Scheme.
+/// Stencil Local Lax Friedrichs Discretization Scheme.
 /// It uses a stencil of order around active points, in order to
 /// evaluate dissipation values for each point, taking into account
 /// the mathematical nature of the speed function.
@@ -447,7 +447,7 @@ namespace viennals {
 using namespace viennacore;
 
 /// This function creates the specialized layer wrapping which
-/// produces better results for the SSLF integration scheme.
+/// produces better results for the SSLF spatial discretization scheme.
 /// isDepo must contain whether the corresponding level sets
 /// are used for deposition or not.
 /// This function assumes that the layers where deposition is
diff --git a/include/viennals/lsVelocityField.hpp b/include/viennals/lsVelocityField.hpp
index 50444a1c..3b853690 100644
--- a/include/viennals/lsVelocityField.hpp
+++ b/include/viennals/lsVelocityField.hpp
@@ -29,9 +29,9 @@ template <class T> class VelocityField {
     return Vec3D<T>{0, 0, 0};
   }
 
-  /// If lsLocalLaxFriedrichsAnalytical is used as the advection scheme,
-  /// this is called to provide the analytical solution for the alpha
-  /// values, needed for stable integration.
+  /// If lsLocalLaxFriedrichsAnalytical is used as the spatial discretization
+  /// scheme, this is called to provide the analytical solution for the alpha
+  /// values, needed for numerical stability.
   virtual T getDissipationAlpha(int /*direction*/, int /*material*/,
                                 const Vec3D<T> & /*centralDifferences*/) {
     return 0;
diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp
index 332cb05c..af6b0644 100644
--- a/python/pyWrap.cpp
+++ b/python/pyWrap.cpp
@@ -81,31 +81,34 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
       .def("print", [](Logger &instance) { instance.print(std::cout); });
 
   // ------ ENUMS ------
-  py::native_enum<IntegrationSchemeEnum>(module, "IntegrationSchemeEnum",
-                                         "enum.IntEnum")
+  py::native_enum<SpatialSchemeEnum>(module, "SpatialSchemeEnum",
+                                     "enum.IntEnum")
       .value("ENGQUIST_OSHER_1ST_ORDER",
-             IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER)
+             SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER)
       .value("ENGQUIST_OSHER_2ND_ORDER",
-             IntegrationSchemeEnum::ENGQUIST_OSHER_2ND_ORDER)
+             SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER)
       .value("LAX_FRIEDRICHS_1ST_ORDER",
-             IntegrationSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER)
+             SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER)
       .value("LAX_FRIEDRICHS_2ND_ORDER",
-             IntegrationSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER)
+             SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER)
       .value("LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER",
-             IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER)
+             SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER)
       .value("LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER",
-             IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
+             SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
       .value("LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER",
-             IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER)
+             SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER)
       .value("LOCAL_LAX_FRIEDRICHS_1ST_ORDER",
-             IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
+             SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
       .value("LOCAL_LAX_FRIEDRICHS_2ND_ORDER",
-             IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER)
+             SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER)
       .value("STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER",
-             IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
-      .value("WENO_5TH_ORDER", IntegrationSchemeEnum::WENO_5TH_ORDER)
+             SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
+      .value("WENO_5TH_ORDER", SpatialSchemeEnum::WENO_5TH_ORDER)
       .finalize();
 
+  module.attr("IntegrationSchemeEnum") =
+      module.attr("SpatialSchemeEnum"); // IntegrationSchemeEnum is deprecated
+
   py::native_enum<BooleanOperationEnum>(module, "BooleanOperationEnum",
                                         "enum.IntEnum")
       .value("INTERSECT", BooleanOperationEnum::INTERSECT)
@@ -163,7 +166,7 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
   py::class_<PointData<T>, SmartPointer<PointData<T>>>(module, "PointData")
       // constructors
       .def(py::init(&SmartPointer<PointData<T>>::New<>))
-      // methods
+      // methods    
       .def("insertNextScalarData",
            (void(PointData<T>::*)(const PointData<T>::ScalarDataType &,
                                   const std::string &)) &
diff --git a/python/pyWrap.hpp b/python/pyWrap.hpp
index 0f33df3d..9001ae7d 100644
--- a/python/pyWrap.hpp
+++ b/python/pyWrap.hpp
@@ -220,10 +220,14 @@ template <int D> void bindApi(py::module &module) {
       .def("getCalculateNormalVectors",
            &Advect<T, D>::getCalculateNormalVectors,
            "Get whether normal vectors are computed during advection.")
+      .def("setSpatialScheme", &Advect<T, D>::setSpatialScheme,
+           "Set the spatial discretization scheme to use during advection.")
       .def("setIntegrationScheme", &Advect<T, D>::setIntegrationScheme,
-           "Set the integration scheme to use during advection.")
+           "(DEPRECATED, use setSpatialScheme instead) Set the spatial "
+           "discretization scheme to use during advection.")
       .def("setDissipationAlpha", &Advect<T, D>::setDissipationAlpha,
-           "Set the dissipation value to use for Lax Friedrichs integration.")
+           "Set the dissipation value to use for Lax Friedrichs spatial "
+           "discretization.")
       .def("setUpdatePointData", &Advect<T, D>::setUpdatePointData,
            "Set whether the point data in the old LS should be translated to "
            "the advected LS. Defaults to true.")
@@ -620,26 +624,26 @@ template <int D> void bindApi(py::module &module) {
       .def("setLevelSet", &MakeGeometry<T, D>::setLevelSet,
            "Set the levelset in which to create the geometry.")
       .def("setGeometry",
-           (void(MakeGeometry<T, D>::*)(SmartPointer<Sphere<T, D>>)) &
-               MakeGeometry<T, D>::setGeometry)
+           (void (MakeGeometry<T, D>::*)(
+               SmartPointer<Sphere<T, D>>))&MakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void(MakeGeometry<T, D>::*)(SmartPointer<Plane<T, D>>)) &
-               MakeGeometry<T, D>::setGeometry)
+           (void (MakeGeometry<T, D>::*)(
+               SmartPointer<Plane<T, D>>))&MakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void(MakeGeometry<T, D>::*)(SmartPointer<Box<T, D>>)) &
-               MakeGeometry<T, D>::setGeometry)
+           (void (MakeGeometry<T, D>::*)(
+               SmartPointer<Box<T, D>>))&MakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void(MakeGeometry<T, D>::*)(SmartPointer<Cylinder<T, D>>)) &
-               MakeGeometry<T, D>::setGeometry)
+           (void (MakeGeometry<T, D>::*)(
+               SmartPointer<Cylinder<T, D>>))&MakeGeometry<T, D>::setGeometry)
       .def("setGeometry",
-           (void(MakeGeometry<T, D>::*)(SmartPointer<PointCloud<T, D>>)) &
-               MakeGeometry<T, D>::setGeometry)
+           (void (MakeGeometry<T, D>::*)(
+               SmartPointer<PointCloud<T, D>>))&MakeGeometry<T, D>::setGeometry)
       .def("setIgnoreBoundaryConditions",
-           (void(MakeGeometry<T, D>::*)(bool)) &
-               MakeGeometry<T, D>::setIgnoreBoundaryConditions)
+           (void (MakeGeometry<T, D>::*)(
+               bool))&MakeGeometry<T, D>::setIgnoreBoundaryConditions)
       .def("setIgnoreBoundaryConditions",
-           (void(MakeGeometry<T, D>::*)(std::array<bool, 3>)) &
-               MakeGeometry<T, D>::setIgnoreBoundaryConditions)
+           (void (MakeGeometry<T, D>::*)(std::array<bool, 3>))&MakeGeometry<
+               T, D>::setIgnoreBoundaryConditions)
       .def("apply", &MakeGeometry<T, D>::apply, "Generate the geometry.");
 
   // MarkVoidPoints
diff --git a/python/viennals/__init__.pyi b/python/viennals/__init__.pyi
index 43c7db13..9c32e40b 100644
--- a/python/viennals/__init__.pyi
+++ b/python/viennals/__init__.pyi
@@ -17,13 +17,14 @@ from viennals._core import CurvatureEnum
 from viennals._core import Extrude
 from viennals._core import FeatureDetectionEnum
 from viennals._core import FileFormatEnum
-from viennals._core import IntegrationSchemeEnum
 from viennals._core import LogLevel
 from viennals._core import Logger
 from viennals._core import MaterialMap
 from viennals._core import Mesh
 from viennals._core import PointData
 from viennals._core import Slice
+from viennals._core import SpatialSchemeEnum
+from viennals._core import SpatialSchemeEnum as IntegrationSchemeEnum
 from viennals._core import TransformEnum
 from viennals._core import TransformMesh
 from viennals._core import VTKReader
@@ -133,6 +134,7 @@ __all__: list[str] = [
     "Reduce",
     "RemoveStrayPoints",
     "Slice",
+    "SpatialSchemeEnum",
     "Sphere",
     "SphereDistribution",
     "StencilLocalLaxFriedrichsScalar",
diff --git a/python/viennals/_core.pyi b/python/viennals/_core.pyi
index f84b7449..81ff81af 100644
--- a/python/viennals/_core.pyi
+++ b/python/viennals/_core.pyi
@@ -8,8 +8,8 @@ import enum
 import typing
 from viennals import d2
 import viennals.d2
-from viennals import d3
 import viennals.d3
+from viennals import d3
 
 __all__: list[str] = [
     "BooleanOperationEnum",
@@ -25,6 +25,7 @@ __all__: list[str] = [
     "Mesh",
     "PointData",
     "Slice",
+    "SpatialSchemeEnum",
     "TransformEnum",
     "TransformMesh",
     "VTKReader",
@@ -187,47 +188,6 @@ class FileFormatEnum(enum.IntEnum):
         Convert to a string according to format_spec.
         """
 
-class IntegrationSchemeEnum(enum.IntEnum):
-    ENGQUIST_OSHER_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.ENGQUIST_OSHER_1ST_ORDER: 0>
-    ENGQUIST_OSHER_2ND_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.ENGQUIST_OSHER_2ND_ORDER: 1>
-    LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER: 2>
-    LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LAX_FRIEDRICHS_2ND_ORDER: 3>
-    LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 7>
-    LOCAL_LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LOCAL_LAX_FRIEDRICHS_2ND_ORDER: 8>
-    LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER: 4>
-    LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 5>
-    LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER: 6>
-    STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 9>
-    WENO_5TH_ORDER: typing.ClassVar[
-        IntegrationSchemeEnum
-    ]  # value = <IntegrationSchemeEnum.WENO_5TH_ORDER: 10>
-    @classmethod
-    def __new__(cls, value): ...
-    def __format__(self, format_spec):
-        """
-        Convert to a string according to format_spec.
-        """
-
 class LogLevel:
     """
     Members:
@@ -539,6 +499,47 @@ class Slice:
         Set the path where the slice should be written to.
         """
 
+class SpatialSchemeEnum(enum.IntEnum):
+    ENGQUIST_OSHER_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.ENGQUIST_OSHER_1ST_ORDER: 0>
+    ENGQUIST_OSHER_2ND_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.ENGQUIST_OSHER_2ND_ORDER: 1>
+    LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER: 2>
+    LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LAX_FRIEDRICHS_2ND_ORDER: 3>
+    LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 7>
+    LOCAL_LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LOCAL_LAX_FRIEDRICHS_2ND_ORDER: 8>
+    LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER: 4>
+    LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 5>
+    LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER: 6>
+    STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER: 9>
+    WENO_5TH_ORDER: typing.ClassVar[
+        SpatialSchemeEnum
+    ]  # value = <SpatialSchemeEnum.WENO_5TH_ORDER: 10>
+    @classmethod
+    def __new__(cls, value): ...
+    def __format__(self, format_spec):
+        """
+        Convert to a string according to format_spec.
+        """
+
 class TransformEnum(enum.IntEnum):
     ROTATION: typing.ClassVar[TransformEnum]  # value = <TransformEnum.ROTATION: 1>
     SCALE: typing.ClassVar[TransformEnum]  # value = <TransformEnum.SCALE: 2>
@@ -733,3 +734,4 @@ def setNumThreads(arg0: typing.SupportsInt) -> None: ...
 
 __version__: str = "5.3.0"
 version: str = "5.3.0"
+IntegrationSchemeEnum = SpatialSchemeEnum
diff --git a/python/viennals/d2.pyi b/python/viennals/d2.pyi
index 439ce917..d9d57727 100644
--- a/python/viennals/d2.pyi
+++ b/python/viennals/d2.pyi
@@ -133,7 +133,7 @@ class Advect:
 
     def setDissipationAlpha(self, arg0: typing.SupportsFloat) -> None:
         """
-        Set the dissipation value to use for Lax Friedrichs integration.
+        Set the dissipation value to use for Lax Friedrichs spatial discretization.
         """
 
     def setIgnoreVoids(self, arg0: bool) -> None:
@@ -141,9 +141,9 @@ class Advect:
         Set whether voids in the geometry should be ignored during advection or not.
         """
 
-    def setIntegrationScheme(self, arg0: viennals._core.IntegrationSchemeEnum) -> None:
+    def setIntegrationScheme(self, arg0: viennals._core.SpatialSchemeEnum) -> None:
         """
-        Set the integration scheme to use during advection.
+        (DEPRECATED, use setSpatialScheme instead) Set the spatial discretization scheme to use during advection.
         """
 
     def setSaveAdvectionVelocities(self, arg0: bool) -> None:
@@ -156,6 +156,11 @@ class Advect:
         Set whether only a single advection step should be performed.
         """
 
+    def setSpatialScheme(self, arg0: viennals._core.SpatialSchemeEnum) -> None:
+        """
+        Set the spatial discretization scheme to use during advection.
+        """
+
     def setTimeStepRatio(self, arg0: typing.SupportsFloat) -> None:
         """
         Set the maximum time step size relative to grid size. Advection is only stable for <0.5.
diff --git a/python/viennals/d3.pyi b/python/viennals/d3.pyi
index b6da65ef..de64dbbe 100644
--- a/python/viennals/d3.pyi
+++ b/python/viennals/d3.pyi
@@ -128,7 +128,7 @@ class Advect:
 
     def setDissipationAlpha(self, arg0: typing.SupportsFloat) -> None:
         """
-        Set the dissipation value to use for Lax Friedrichs integration.
+        Set the dissipation value to use for Lax Friedrichs spatial discretization.
         """
 
     def setIgnoreVoids(self, arg0: bool) -> None:
@@ -136,9 +136,9 @@ class Advect:
         Set whether voids in the geometry should be ignored during advection or not.
         """
 
-    def setIntegrationScheme(self, arg0: viennals._core.IntegrationSchemeEnum) -> None:
+    def setIntegrationScheme(self, arg0: viennals._core.SpatialSchemeEnum) -> None:
         """
-        Set the integration scheme to use during advection.
+        (DEPRECATED, use setSpatialScheme instead) Set the spatial discretization scheme to use during advection.
         """
 
     def setSaveAdvectionVelocities(self, arg0: bool) -> None:
@@ -151,6 +151,11 @@ class Advect:
         Set whether only a single advection step should be performed.
         """
 
+    def setSpatialScheme(self, arg0: viennals._core.SpatialSchemeEnum) -> None:
+        """
+        Set the spatial discretization scheme to use during advection.
+        """
+
     def setTimeStepRatio(self, arg0: typing.SupportsFloat) -> None:
         """
         Set the maximum time step size relative to grid size. Advection is only stable for <0.5.
diff --git a/tests/Advection/Advection.cpp b/tests/Advection/Advection.cpp
index 3806505c..4e0d4a2c 100644
--- a/tests/Advection/Advection.cpp
+++ b/tests/Advection/Advection.cpp
@@ -49,18 +49,18 @@ int main() {
 
   double gridDelta = 0.6999999;
 
-  std::vector<ls::IntegrationSchemeEnum> integrationSchemes = {
-      ls::IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER,
-      ls::IntegrationSchemeEnum::ENGQUIST_OSHER_2ND_ORDER,
-      ls::IntegrationSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER,
-      ls::IntegrationSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER,
-      ls::IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
-      ls::IntegrationSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
-      ls::IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
-      ls::IntegrationSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
-      ls::IntegrationSchemeEnum::WENO_5TH_ORDER};
-
-  for (auto integrationScheme : integrationSchemes) {
+  std::vector<ls::SpatialSchemeEnum> spatialSchemes = {
+      ls::SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER,
+      ls::SpatialSchemeEnum::ENGQUIST_OSHER_2ND_ORDER,
+      ls::SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER,
+      ls::SpatialSchemeEnum::LAX_FRIEDRICHS_2ND_ORDER,
+      ls::SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
+      ls::SpatialSchemeEnum::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
+      ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
+      ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
+      ls::SpatialSchemeEnum::WENO_5TH_ORDER};
+
+  for (auto scheme : spatialSchemes) {
     auto sphere1 = ls::Domain<double, D>::New(gridDelta);
 
     double origin[3] = {5., 0., 0.};
@@ -93,7 +93,7 @@ int main() {
     ls::AdvectForwardEuler<double, D> advectionKernel;
     advectionKernel.insertNextLevelSet(sphere1);
     advectionKernel.setVelocityField(velocities);
-    advectionKernel.setIntegrationScheme(integrationScheme);
+    advectionKernel.setSpatialScheme(scheme);
     advectionKernel.setSaveAdvectionVelocities(true);
 
     double time = 0.;
diff --git a/tests/Advection/Advection.py b/tests/Advection/Advection.py
index c7b66aa1..aebf3cb0 100644
--- a/tests/Advection/Advection.py
+++ b/tests/Advection/Advection.py
@@ -2,6 +2,7 @@
 
 vls.setDimension(3)
 
+
 # Implement own velocity field
 class VelocityField(vls.VelocityField):
     def __init__(self):
@@ -15,25 +16,26 @@ def getScalarVelocity(self, coord, material, normal, pointId):
     def getVectorVelocity(self, coord, material, normal, pointId):
         return [0.0, 0.0, 0.0]
 
+
 def main():
     # Set number of threads
     vls.setNumThreads(8)
 
     gridDelta = 0.6999999
 
-    integrationSchemes = [
-        vls.IntegrationSchemeEnum.ENGQUIST_OSHER_1ST_ORDER,
-        vls.IntegrationSchemeEnum.ENGQUIST_OSHER_2ND_ORDER,
-        vls.IntegrationSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER,
-        vls.IntegrationSchemeEnum.LAX_FRIEDRICHS_2ND_ORDER,
-        vls.IntegrationSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
-        vls.IntegrationSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
-        vls.IntegrationSchemeEnum.LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
-        vls.IntegrationSchemeEnum.LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
-        vls.IntegrationSchemeEnum.WENO_5TH_ORDER
+    discretizationSchemes = [
+        vls.SpatialSchemeEnum.ENGQUIST_OSHER_1ST_ORDER,
+        vls.SpatialSchemeEnum.ENGQUIST_OSHER_2ND_ORDER,
+        vls.SpatialSchemeEnum.LAX_FRIEDRICHS_1ST_ORDER,
+        vls.SpatialSchemeEnum.LAX_FRIEDRICHS_2ND_ORDER,
+        vls.SpatialSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
+        vls.SpatialSchemeEnum.LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
+        vls.SpatialSchemeEnum.LOCAL_LAX_FRIEDRICHS_1ST_ORDER,
+        vls.SpatialSchemeEnum.LOCAL_LAX_FRIEDRICHS_2ND_ORDER,
+        vls.SpatialSchemeEnum.WENO_5TH_ORDER,
     ]
 
-    for scheme in integrationSchemes:
+    for scheme in discretizationSchemes:
         sphere1 = vls.Domain(gridDelta)
 
         origin = [5.0, 0.0, 0.0]
@@ -47,7 +49,7 @@ def main():
         advectionKernel = vls.AdvectRungeKutta3()
         advectionKernel.insertNextLevelSet(sphere1)
         advectionKernel.setVelocityField(velocities)
-        advectionKernel.setIntegrationScheme(scheme)
+        advectionKernel.setSpatialScheme(scheme)
         advectionKernel.setSaveAdvectionVelocities(True)
 
         time = 0.0
@@ -65,5 +67,6 @@ def main():
             i += 1
         print(f"Done scheme {scheme}")
 
+
 if __name__ == "__main__":
     main()
diff --git a/tests/Advection/StencilLaxFriedrichsTest.cpp b/tests/Advection/StencilLaxFriedrichsTest.cpp
index 1f6e20b0..0cfef606 100644
--- a/tests/Advection/StencilLaxFriedrichsTest.cpp
+++ b/tests/Advection/StencilLaxFriedrichsTest.cpp
@@ -61,9 +61,9 @@ int main() {
   advectionKernel.setVelocityField(velocityField);
   advectionKernel.setAdvectionTime(0.5);
 
-  // Set the specific integration scheme
-  advectionKernel.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
+  // Set the specific spatial discretization scheme
+  advectionKernel.setSpatialScheme(
+      ls::SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER);
 
   // Run Advection
   std::cout << "Running Stencil Local Lax Friedrichs Advection..." << std::endl;
diff --git a/tests/Advection/TimeIntegrationComparison.cpp b/tests/Advection/TimeIntegrationComparison.cpp
index 5b1dbe21..f1f2d018 100644
--- a/tests/Advection/TimeIntegrationComparison.cpp
+++ b/tests/Advection/TimeIntegrationComparison.cpp
@@ -65,16 +65,14 @@ int main() {
   advectFE.insertNextLevelSet(sphereFE);
   advectFE.setVelocityField(velocityField);
   advectFE.setAdvectionTime(2.0);
-  advectFE.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
+  advectFE.setSpatialScheme(ls::SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
 
   // Setup Advection: Runge-Kutta 3
   ls::AdvectRungeKutta3<T, D> advectRK3;
   advectRK3.insertNextLevelSet(sphereRK3);
   advectRK3.setVelocityField(velocityField);
   advectRK3.setAdvectionTime(2.0);
-  advectRK3.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
+  advectRK3.setSpatialScheme(ls::SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
 
   // Run Advection
   std::cout << "Running Forward Euler Advection..." << std::endl;
diff --git a/tests/Advection/TimeIntegrationComparison.py b/tests/Advection/TimeIntegrationComparison.py
index 17b4c8e0..b179a245 100644
--- a/tests/Advection/TimeIntegrationComparison.py
+++ b/tests/Advection/TimeIntegrationComparison.py
@@ -1,6 +1,8 @@
 import viennals as vls
+
 vls.setDimension(3)
 
+
 # Define a constant velocity field
 class ConstantVelocity(vls.VelocityField):
     def __init__(self, vel):
@@ -13,6 +15,7 @@ def getScalarVelocity(self, coord, material, normal, pointId):
     def getVectorVelocity(self, coord, material, normal, pointId):
         return self.velocity
 
+
 def main():
     # Define grid and domain bounds
     gridDelta = 0.1
@@ -38,19 +41,19 @@ def main():
     advectFE.insertNextLevelSet(sphereFE)
     advectFE.setVelocityField(velocityField)
     advectFE.setAdvectionTime(2.0)
-    advectFE.setIntegrationScheme(vls.IntegrationSchemeEnum.ENGQUIST_OSHER_1ST_ORDER)
+    advectFE.setSpatialScheme(vls.SpatialSchemeEnum.ENGQUIST_OSHER_1ST_ORDER)
 
     # Setup Advection: Runge-Kutta 3
     advectRK3 = vls.AdvectRungeKutta3()
     advectRK3.insertNextLevelSet(sphereRK3)
     advectRK3.setVelocityField(velocityField)
     advectRK3.setAdvectionTime(2.0)
-    advectRK3.setIntegrationScheme(vls.IntegrationSchemeEnum.ENGQUIST_OSHER_1ST_ORDER)
+    advectRK3.setSpatialScheme(vls.SpatialSchemeEnum.ENGQUIST_OSHER_1ST_ORDER)
 
     # Run Advection
     print("Running Forward Euler Advection...")
     advectFE.apply()
-    
+
     checkFE = vls.Check(sphereFE)
     checkFE.apply()
 
@@ -60,7 +63,7 @@ def main():
 
     print("Running Runge-Kutta 3 Advection...")
     advectRK3.apply()
-    
+
     checkRK3 = vls.Check(sphereRK3)
     checkRK3.apply()
 
@@ -68,5 +71,6 @@ def main():
     vls.ToSurfaceMesh(sphereRK3, meshRK3).apply()
     vls.VTKWriter(meshRK3, "sphereRK3.vtp").apply()
 
+
 if __name__ == "__main__":
     main()
diff --git a/tests/Advection2D/Advection2D.cpp b/tests/Advection2D/Advection2D.cpp
index 51aa24dc..e49a9172 100644
--- a/tests/Advection2D/Advection2D.cpp
+++ b/tests/Advection2D/Advection2D.cpp
@@ -105,8 +105,8 @@ int main() {
   ls::Advect<double, D> advectionKernel;
   advectionKernel.insertNextLevelSet(sphere1);
   advectionKernel.setVelocityField(velocities);
-  advectionKernel.setIntegrationScheme(
-      ls::IntegrationSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER);
+  advectionKernel.setSpatialScheme(
+      ls::SpatialSchemeEnum::LAX_FRIEDRICHS_1ST_ORDER);
   advectionKernel.setAdvectionTime(20.);
   advectionKernel.apply();
 
diff --git a/tests/GeometricAdvectPerformance/GeometricAdvectPerformance.cpp b/tests/GeometricAdvectPerformance/GeometricAdvectPerformance.cpp
index 1229683d..64ebdb05 100644
--- a/tests/GeometricAdvectPerformance/GeometricAdvectPerformance.cpp
+++ b/tests/GeometricAdvectPerformance/GeometricAdvectPerformance.cpp
@@ -126,7 +126,7 @@ int main() {
   // ls::ToSurfaceMesh<double, D>(newLayer, mesh).apply();
   // ls::VTKWriter<double>(mesh, "finalSurface.vtk").apply();
 
-  // now rund lsAdvect for all other advection schemes
+  // now run lsAdvect for all other advection schemes
   // last scheme is SLLFS with i == 9
   for (unsigned i = 0; i < 10; ++i) {
     if (i == 4) {
@@ -139,8 +139,7 @@ int main() {
     auto velocities = ls::SmartPointer<velocityField>::New();
     advectionKernel.setVelocityField(velocities);
     advectionKernel.setAdvectionTime(depositionDistance);
-    advectionKernel.setIntegrationScheme(
-        static_cast<ls::IntegrationSchemeEnum>(i));
+    advectionKernel.setSpatialScheme(static_cast<ls::SpatialSchemeEnum>(i));
     {
       auto start = std::chrono::high_resolution_clock::now();
       advectionKernel.apply();