diff --git a/.gitignore b/.gitignore
index 26732282..12a3b820 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,4 +23,5 @@ _generate/
 *env*
 .mypy_cache/
 .eggs/
-test.ipynb
\ No newline at end of file
+test.ipynb
+testall.sh
diff --git a/CMakeLists.txt b/CMakeLists.txt
index b61b8e7f..ceced024 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -112,7 +112,7 @@ include(cmake/vtk.cmake)
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 1.9.0
+  VERSION 1.9.3
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   OPTIONS "VIENNACORE_FORMAT_EXCLUDE build/"
   EXCLUDE_FROM_ALL ${VIENNALS_BUILD_PYTHON})
@@ -125,7 +125,7 @@ CPMAddPackage(
 
 CPMFindPackage(
   NAME ViennaHRLE
-  VERSION 0.7.0
+  VERSION 0.8.0
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaHRLE"
   EXCLUDE_FROM_ALL ${VIENNALS_BUILD_PYTHON})
 
diff --git a/examples/SquareEtch/SquareEtch.cpp b/examples/SquareEtch/SquareEtch.cpp
index b6bf661e..f7673913 100644
--- a/examples/SquareEtch/SquareEtch.cpp
+++ b/examples/SquareEtch/SquareEtch.cpp
@@ -200,8 +200,7 @@ int main() {
     // Analytical velocity fields and dissipation coefficients
     // can only be used with this spatial discretization scheme
     advectionKernel.setSpatialScheme(
-        // ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
-        ls::SpatialSchemeEnum::WENO_5TH_ORDER);
+        ls::SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_ANALYTICAL_1ST_ORDER);
   } else {
     // for numerical velocities, just use the default
     // spatial discretization scheme, which is not accurate for certain
@@ -211,8 +210,8 @@ int main() {
     // For coordinate independent velocity functions
     // this numerical scheme is superior though.
     // However, it is slower.
-    // advectionKernel.setSpatialScheme(
-    //     ls::SpatialSchemeEnum::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/include/viennals/lsAdvect.hpp b/include/viennals/lsAdvect.hpp
index 33309638..64ad12b6 100644
--- a/include/viennals/lsAdvect.hpp
+++ b/include/viennals/lsAdvect.hpp
@@ -24,7 +24,7 @@
 #include <lsLocalLaxFriedrichsAnalytical.hpp>
 #include <lsLocalLocalLaxFriedrichs.hpp>
 #include <lsStencilLocalLaxFriedrichsScalar.hpp>
-#include <lsWENO5.hpp>
+#include <lsWENO.hpp>
 
 // Include implementation of time integration schemes
 #include <lsAdvectIntegrationSchemes.hpp>
@@ -147,12 +147,12 @@ template <class T, int D> class Advect {
             // move neighborIterator to current position
             neighborIterator.goToIndicesSequential(indices);
 
-            Vec3D<T> coords;
+            Vec3D<T> coords{};
             for (unsigned i = 0; i < D; ++i) {
               coords[i] = indices[i] * gridDelta;
             }
 
-            Vec3D<T> normal = {};
+            Vec3D<T> normal{};
             T normalModulus = 0.;
             for (unsigned i = 0; i < D; ++i) {
               const T phiPos = neighborIterator.getNeighbor(i).getValue();
@@ -664,10 +664,15 @@ template <class T, int D> class Advect {
       auto is = lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>(
           levelSets.back(), velocities, dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
+      // Instantiate WENO with order 3
+      auto is = lsInternal::WENO<T, D, 3>(levelSets.back(), velocities,
+                                          dissipationAlpha);
+      currentTimeStep = integrateTime(is, maxTimeStep);
     } 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);
+      // Instantiate WENO with order 5
+      auto is = lsInternal::WENO<T, D, 5>(levelSets.back(), velocities,
+                                          dissipationAlpha);
       currentTimeStep = integrateTime(is, maxTimeStep);
     } else {
       VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
@@ -998,9 +1003,10 @@ template <class T, int D> class Advect {
                SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER) {
       lsInternal::StencilLocalLaxFriedrichsScalar<T, D, 1>::prepareLS(
           levelSets.back());
+    } else if (spatialScheme == SpatialSchemeEnum::WENO_3RD_ORDER) {
+      lsInternal::WENO<T, D, 3>::prepareLS(levelSets.back());
     } 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());
+      lsInternal::WENO<T, D, 5>::prepareLS(levelSets.back());
     } else {
       VIENNACORE_LOG_ERROR("Advect: Discretization scheme not found.");
     }
diff --git a/include/viennals/lsAdvectIntegrationSchemes.hpp b/include/viennals/lsAdvectIntegrationSchemes.hpp
index 691c3ba5..ca005ae5 100644
--- a/include/viennals/lsAdvectIntegrationSchemes.hpp
+++ b/include/viennals/lsAdvectIntegrationSchemes.hpp
@@ -18,7 +18,8 @@ enum class SpatialSchemeEnum : unsigned {
   LOCAL_LAX_FRIEDRICHS_1ST_ORDER = 7,
   LOCAL_LAX_FRIEDRICHS_2ND_ORDER = 8,
   STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER = 9,
-  WENO_5TH_ORDER = 10
+  WENO_3RD_ORDER = 10,
+  WENO_5TH_ORDER = 11
 };
 
 // Legacy naming (deprecated, will be removed in future versions)
diff --git a/include/viennals/lsCalculateNormalVectors.hpp b/include/viennals/lsCalculateNormalVectors.hpp
index d24c77ce..3afd14a0 100644
--- a/include/viennals/lsCalculateNormalVectors.hpp
+++ b/include/viennals/lsCalculateNormalVectors.hpp
@@ -134,12 +134,12 @@ template <class T, int D> class CalculateNormalVectors {
           continue;
         } else if (std::abs(center.getValue()) > maxValue) {
           // push an empty vector to keep ordering correct
-          Vec3D<T> tmp = {};
+          Vec3D<T> tmp{};
           normalVectors.push_back(tmp);
           continue;
         }
 
-        Vec3D<T> n;
+        Vec3D<T> n{};
 
         T denominator = 0;
         for (int i = 0; i < D; i++) {
diff --git a/include/viennals/lsCalculateVisibilities.hpp b/include/viennals/lsCalculateVisibilities.hpp
index 33f10fb9..45bee138 100644
--- a/include/viennals/lsCalculateVisibilities.hpp
+++ b/include/viennals/lsCalculateVisibilities.hpp
@@ -41,8 +41,8 @@ template <class T, int D> class CalculateVisibilities {
       std::fill(visibilities.begin(), visibilities.end(), static_cast<T>(1.0));
     } else {
       // *** Determine extents of domain ***
-      Vec3D<T> minDefinedPoint;
-      Vec3D<T> maxDefinedPoint;
+      Vec3D<T> minDefinedPoint{};
+      Vec3D<T> maxDefinedPoint{};
       // Initialize with extreme values
       for (int i = 0; i < D; ++i) {
         minDefinedPoint[i] = std::numeric_limits<T>::max();
@@ -87,7 +87,7 @@ template <class T, int D> class CalculateVisibilities {
             continue;
 
           // Starting position of the point
-          Vec3D<T> currentPos;
+          Vec3D<T> currentPos{};
           for (int i = 0; i < D; ++i) {
             currentPos[i] = it.getStartIndices(i);
           }
diff --git a/include/viennals/lsCompareCriticalDimensions.hpp b/include/viennals/lsCompareCriticalDimensions.hpp
index 98e8270d..d2f2e6c2 100644
--- a/include/viennals/lsCompareCriticalDimensions.hpp
+++ b/include/viennals/lsCompareCriticalDimensions.hpp
@@ -339,7 +339,8 @@ template <class T, int D = 2> class CompareCriticalDimensions {
         continue;
 
       // Create points for target and sample positions
-      Vec3D<T> coordTarget = {0.0, 0.0, 0.0}, coordSample = {0.0, 0.0, 0.0};
+      Vec3D<T> coordTarget{};
+      Vec3D<T> coordSample{};
 
       for (int i = 0; i < D; ++i) {
         if (i == result.measureDimension) {
diff --git a/include/viennals/lsDomain.hpp b/include/viennals/lsDomain.hpp
index 5277c282..28633455 100644
--- a/include/viennals/lsDomain.hpp
+++ b/include/viennals/lsDomain.hpp
@@ -2,6 +2,7 @@
 
 #include <lsPreCompileMacros.hpp>
 
+#include <cstring>
 #include <limits>
 
 #include <hrleDomain.hpp>
@@ -221,7 +222,7 @@ template <class T, int D> class Domain {
     // Check identifier
     char identifier[8];
     stream.read(identifier, 8);
-    if (std::string(identifier).compare(0, 8, "lsDomain")) {
+    if (std::memcmp(identifier, "lsDomain", 8) != 0) {
       Logger::getInstance()
           .addError(
               "Reading Domain from stream failed. Header could not be found.")
diff --git a/include/viennals/lsEngquistOsher.hpp b/include/viennals/lsEngquistOsher.hpp
index f61c5d0a..c8f019df 100644
--- a/include/viennals/lsEngquistOsher.hpp
+++ b/include/viennals/lsEngquistOsher.hpp
@@ -121,7 +121,7 @@ template <class T, int D, int order> class EngquistOsher {
 
     // Calculate normal vector for velocity calculation
     // use std::array since it will be exposed to interface
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     if (calculateNormalVectors) {
       T denominator = 0;
       for (int i = 0; i < D; i++) {
diff --git a/include/viennals/lsExtrude.hpp b/include/viennals/lsExtrude.hpp
index 41f74423..c2d20580 100644
--- a/include/viennals/lsExtrude.hpp
+++ b/include/viennals/lsExtrude.hpp
@@ -19,7 +19,7 @@ template <class T> class Extrude {
   SmartPointer<Domain<T, 3>> outputLevelSet = nullptr;
   Vec2D<T> extent = {0., 0.};
   int extrusionAxis = 0;
-  std::array<BoundaryConditionEnum, 3> boundaryConds = {};
+  std::array<BoundaryConditionEnum, 3> boundaryConds{};
 
 public:
   Extrude() = default;
diff --git a/include/viennals/lsGeometricAdvect.hpp b/include/viennals/lsGeometricAdvect.hpp
index 60abb04e..5af9fffc 100644
--- a/include/viennals/lsGeometricAdvect.hpp
+++ b/include/viennals/lsGeometricAdvect.hpp
@@ -482,7 +482,7 @@ template <class T, int D> class GeometricAdvect {
     {
       std::vector<T> scalarData;
       for (auto it = newPoints[0].begin(); it != newPoints[0].end(); ++it) {
-        Vec3D<T> node{0., 0., 0.};
+        Vec3D<T> node{};
         for (unsigned i = 0; i < D; ++i) {
           node[i] = T((it->first)[i]) * gridDelta;
         }
diff --git a/include/viennals/lsGeometricAdvectDistributions.hpp b/include/viennals/lsGeometricAdvectDistributions.hpp
index c6bab5e4..55b8dd8f 100644
--- a/include/viennals/lsGeometricAdvectDistributions.hpp
+++ b/include/viennals/lsGeometricAdvectDistributions.hpp
@@ -108,7 +108,7 @@ class SphereDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -radius;
       bounds[2 * i + 1] = radius;
@@ -157,7 +157,7 @@ class BoxDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -posExtent[i];
       bounds[2 * i + 1] = posExtent[i];
@@ -240,7 +240,7 @@ class CustomSphereDistribution : public GeometricAdvectDistribution<T, D> {
   }
 
   std::array<viennahrle::CoordType, 6> getBounds() const override {
-    std::array<viennahrle::CoordType, 6> bounds = {};
+    std::array<viennahrle::CoordType, 6> bounds{};
     for (unsigned i = 0; i < D; ++i) {
       bounds[2 * i] = -maxRadius_;
       bounds[2 * i + 1] = maxRadius_;
diff --git a/include/viennals/lsLaxFriedrichs.hpp b/include/viennals/lsLaxFriedrichs.hpp
index 7320b518..04daf038 100644
--- a/include/viennals/lsLaxFriedrichs.hpp
+++ b/include/viennals/lsLaxFriedrichs.hpp
@@ -61,7 +61,7 @@ template <class T, int D, int order> class LaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
diff --git a/include/viennals/lsLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLaxFriedrichs.hpp
index e73a9809..371e82b6 100644
--- a/include/viennals/lsLocalLaxFriedrichs.hpp
+++ b/include/viennals/lsLocalLaxFriedrichs.hpp
@@ -86,7 +86,7 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
@@ -183,7 +183,7 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
     }
 
     // calculate alphas
-    T alpha[D] = {};
+    T alpha[D]{};
     {
       // alpha calculation is always on order 1 stencil
       const viennahrle::IndexType minIndex = -1;
@@ -192,11 +192,11 @@ template <class T, int D, int order> class LocalLaxFriedrichs {
 
       viennahrle::Index<D> neighborIndex(minIndex);
       for (unsigned i = 0; i < numNeighbors; ++i) {
-        Vec3D<T> coords;
+        Vec3D<T> coords{};
         for (unsigned dir = 0; dir < D; ++dir) {
           coords[dir] = coordinate[dir] + neighborIndex[dir] * gridDelta;
         }
-        Vec3D<T> normal = {};
+        Vec3D<T> normal{};
         double normalModulus = 0.;
         auto center = neighborIterator.getNeighbor(neighborIndex).getValue();
         for (unsigned dir = 0; dir < D; ++dir) {
diff --git a/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp b/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
index 1e276fc4..12ce9168 100644
--- a/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
+++ b/include/viennals/lsLocalLaxFriedrichsAnalytical.hpp
@@ -85,7 +85,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
@@ -182,7 +182,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
     }
 
     // calculate alphas
-    T alpha[D] = {};
+    T alpha[D]{};
     {
       // alpha calculation is always on order 1 stencil
       constexpr viennahrle::IndexType minIndex = -1;
@@ -192,7 +192,7 @@ template <class T, int D, int order> class LocalLaxFriedrichsAnalytical {
       viennahrle::Index<D> neighborIndex(minIndex);
       for (unsigned i = 0; i < numNeighbors; ++i) {
 
-        Vec3D<T> normal = {};
+        Vec3D<T> normal{};
         auto center = neighborIterator.getNeighbor(neighborIndex).getValue();
         for (unsigned dir = 0; dir < D; ++dir) {
           viennahrle::Index<D> unity(0);
diff --git a/include/viennals/lsLocalLocalLaxFriedrichs.hpp b/include/viennals/lsLocalLocalLaxFriedrichs.hpp
index 4361172f..c6da41dd 100644
--- a/include/viennals/lsLocalLocalLaxFriedrichs.hpp
+++ b/include/viennals/lsLocalLocalLaxFriedrichs.hpp
@@ -64,7 +64,7 @@ template <class T, int D, int order> class LocalLocalLaxFriedrichs {
     T grad = 0.;
     T dissipation = 0.;
 
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     T normalModulus = 0;
 
     for (int i = 0; i < D; i++) { // iterate over dimensions
diff --git a/include/viennals/lsMakeGeometry.hpp b/include/viennals/lsMakeGeometry.hpp
index 2a45bbd3..7b13dac5 100644
--- a/include/viennals/lsMakeGeometry.hpp
+++ b/include/viennals/lsMakeGeometry.hpp
@@ -441,10 +441,9 @@ template <class T, int D> class MakeGeometry {
 
         // create and insert points at base
         for (double angle = 0.; angle < limit; angle += smallAngle) {
-          Vec3D<T> point;
+          Vec3D<T> point{};
           point[0] = cylinder->radius * std::cos(angle);
           point[1] = cylinder->radius * std::sin(angle);
-          point[2] = 0.0;
           points.push_back(point);
           mesh->insertNextNode(point);
         }
diff --git a/include/viennals/lsMesh.hpp b/include/viennals/lsMesh.hpp
index d2ce80c1..007bc70b 100644
--- a/include/viennals/lsMesh.hpp
+++ b/include/viennals/lsMesh.hpp
@@ -29,8 +29,8 @@ template <class T = double> class Mesh {
   std::vector<std::array<unsigned, 8>> hexas;
   PointData<T> pointData;
   PointData<T> cellData;
-  Vec3D<T> minimumExtent;
-  Vec3D<T> maximumExtent;
+  Vec3D<T> minimumExtent{};
+  Vec3D<T> maximumExtent{};
 
 private:
   // iterator typedef
@@ -272,6 +272,8 @@ template <class T = double> class Mesh {
     hexas.clear();
     pointData.clear();
     cellData.clear();
+    minimumExtent = Vec3D<T>{};
+    maximumExtent = Vec3D<T>{};
   }
 
   void print() {
diff --git a/include/viennals/lsPointData.hpp b/include/viennals/lsPointData.hpp
index 1d603670..4bcadfa1 100644
--- a/include/viennals/lsPointData.hpp
+++ b/include/viennals/lsPointData.hpp
@@ -3,6 +3,7 @@
 #include <lsPreCompileMacros.hpp>
 
 #include <cstdint>
+#include <cstring>
 #include <utility>
 #include <vector>
 
@@ -385,7 +386,7 @@ class PointData {
   std::istream &deserialize(std::istream &stream) {
     char identifier[11];
     stream.read(identifier, 11);
-    if (std::string(identifier).compare(0, 11, "lsPointData")) {
+    if (std::memcmp(identifier, "lsPointData", 11) != 0) {
       Logger::getInstance()
           .addError("Reading PointData from stream failed. Header could "
                     "not be found.")
diff --git a/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp b/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
index 1393cbb9..503ecf27 100644
--- a/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
+++ b/include/viennals/lsStencilLocalLaxFriedrichsScalar.hpp
@@ -46,7 +46,7 @@ class StencilLocalLaxFriedrichsScalar {
 
   // Final dissipation coefficients that are used by the time integrator. If
   // D==2 last entries are 0.
-  Vec3D<T> finalAlphas;
+  Vec3D<T> finalAlphas{};
   static constexpr unsigned numStencilPoints = hrleUtil::pow(2 * order + 1, D);
   static double maxDissipation; // default: std::numeric_limits<double>::max();
 
@@ -124,7 +124,7 @@ class StencilLocalLaxFriedrichsScalar {
   }
 
   Vec3D<T> calculateNormal(const viennahrle::Index<D> &offset) {
-    Vec3D<T> normal = {0.0, 0.0, 0.0};
+    Vec3D<T> normal{};
     constexpr int startIndex = -1;
     T modulus = 0.;
 
@@ -204,11 +204,7 @@ class StencilLocalLaxFriedrichsScalar {
                                   SmartPointer<viennals::VelocityField<T>> vel,
                                   double a = 1.0)
       : levelSet(passedlsDomain), velocities(vel),
-        neighborIterator(levelSet->getDomain()), alphaFactor(a) {
-    for (int i = 0; i < 3; ++i) {
-      finalAlphas[i] = 0;
-    }
-  }
+        neighborIterator(levelSet->getDomain()), alphaFactor(a) {}
 
   static void setMaxDissipation(double maxDiss) { maxDissipation = maxDiss; }
 
@@ -232,7 +228,7 @@ class StencilLocalLaxFriedrichsScalar {
     auto &grid = levelSet->getGrid();
     double gridDelta = grid.getGridDelta();
 
-    Vec3D<T> coordinate{0., 0., 0.};
+    Vec3D<T> coordinate{};
     for (unsigned i = 0; i < D; ++i) {
       coordinate[i] = indices[i] * gridDelta;
     }
@@ -243,7 +239,7 @@ class StencilLocalLaxFriedrichsScalar {
     // if there is a vector velocity, we need to project it onto a scalar
     // velocity first using its normal vector
     // /*if (vectorVelocity != Vec3D<T>({}))*/ {
-    Vec3D<T> normalVector;
+    Vec3D<T> normalVector{};
     T denominator = 0; // normal modulus
     for (unsigned i = 0; i < D; i++) {
       viennahrle::Index<D> neighborIndex(0);
diff --git a/include/viennals/lsToDiskMesh.hpp b/include/viennals/lsToDiskMesh.hpp
index e2e67d49..baee36f7 100644
--- a/include/viennals/lsToDiskMesh.hpp
+++ b/include/viennals/lsToDiskMesh.hpp
@@ -116,8 +116,8 @@ template <class T, int D, class N = T> class ToDiskMesh {
     std::vector<N> materialIds;
 
     // save the extent of the resulting mesh
-    Vec3D<N> minimumExtent;
-    Vec3D<N> maximumExtent;
+    Vec3D<N> minimumExtent{};
+    Vec3D<N> maximumExtent{};
     for (unsigned i = 0; i < D; ++i) {
       minimumExtent[i] = std::numeric_limits<T>::max();
       maximumExtent[i] = std::numeric_limits<T>::lowest();
@@ -181,7 +181,7 @@ template <class T, int D, class N = T> class ToDiskMesh {
 
       // insert corresponding node shifted by ls value in direction of the
       // normal vector
-      Vec3D<N> node;
+      Vec3D<N> node{};
       node[2] = 0.;
       double max = 0.;
       for (unsigned i = 0; i < D; ++i) {
@@ -210,7 +210,7 @@ template <class T, int D, class N = T> class ToDiskMesh {
 
       // add data into mesh
       // copy normal
-      Vec3D<N> normal;
+      Vec3D<N> normal{};
       if (D == 2)
         normal[2] = 0.;
       for (unsigned i = 0; i < D; ++i) {
diff --git a/include/viennals/lsToMesh.hpp b/include/viennals/lsToMesh.hpp
index 2edc24a3..3669dc92 100644
--- a/include/viennals/lsToMesh.hpp
+++ b/include/viennals/lsToMesh.hpp
@@ -109,7 +109,7 @@ template <class T, int D> class ToMesh {
       mesh->insertNextVertex(vertex);
 
       // insert corresponding node
-      Vec3D<T> node;
+      Vec3D<T> node{};
       if (D == 2)
         node[2] = 0.;
       for (unsigned i = 0; i < D; ++i) {
diff --git a/include/viennals/lsToSurfaceMesh.hpp b/include/viennals/lsToSurfaceMesh.hpp
index 61fbb7ce..36e4f857 100644
--- a/include/viennals/lsToSurfaceMesh.hpp
+++ b/include/viennals/lsToSurfaceMesh.hpp
@@ -140,7 +140,7 @@ template <class T, int D> class ToSurfaceMesh {
           } else { // if node does not exist yet
 
             // calculate coordinate of new node
-            Vec3D<T> cc = {0., 0., 0.}; // initialise with zeros
+            Vec3D<T> cc{}; // initialise with zeros
             std::size_t currentPointId = 0;
             for (int z = 0; z < D; z++) {
               if (z != dir) {
diff --git a/include/viennals/lsToVoxelMesh.hpp b/include/viennals/lsToVoxelMesh.hpp
index f6befa7b..97b36978 100644
--- a/include/viennals/lsToVoxelMesh.hpp
+++ b/include/viennals/lsToVoxelMesh.hpp
@@ -198,7 +198,7 @@ template <class T, int D> class ToVoxelMesh {
     double gridDelta = grid.getGridDelta();
     mesh->nodes.resize(pointIdMapping.size());
     for (auto it = pointIdMapping.begin(); it != pointIdMapping.end(); ++it) {
-      Vec3D<T> coords;
+      Vec3D<T> coords{};
       for (unsigned i = 0; i < D; ++i) {
         coords[i] = gridDelta * it->first[i];
 
diff --git a/include/viennals/lsWENO5.hpp b/include/viennals/lsWENO.hpp
similarity index 68%
rename from include/viennals/lsWENO5.hpp
rename to include/viennals/lsWENO.hpp
index 5e6d5a6b..1d7ca139 100644
--- a/include/viennals/lsWENO5.hpp
+++ b/include/viennals/lsWENO.hpp
@@ -14,34 +14,38 @@ namespace lsInternal {
 
 using namespace viennacore;
 
-/// Fifth-order Weighted Essentially Non-Oscillatory (WENO5) scheme.
+/// Weighted Essentially Non-Oscillatory (WENO) scheme.
 /// This kernel acts as the grid-interface for the mathematical logic
 /// defined in lsFiniteDifferences.hpp.
-template <class T, int D, int order = 3> class WENO5 {
+template <class T, int D, int order> class WENO {
+  static_assert(order == 3 || order == 5, "WENO order must be 3 or 5.");
+
   SmartPointer<viennals::Domain<T, D>> levelSet;
   SmartPointer<viennals::VelocityField<T>> velocities;
 
-  // Iterator depth: WENO5 needs 3 neighbors on each side.
-  viennahrle::SparseStarIterator<viennahrle::Domain<T, D>, order>
+  static constexpr int stencilRadius = (order + 1) / 2;
+
+  // Iterator depth: WENO needs stencilRadius neighbors on each side.
+  viennahrle::SparseStarIterator<viennahrle::Domain<T, D>, stencilRadius>
       neighborIterator;
 
   const bool calculateNormalVectors = true;
 
-  // Use the existing math engine with WENO5 scheme
-  using MathScheme = FiniteDifferences<T, DifferentiationSchemeEnum::WENO5>;
+  // Use the existing math engine with WENO scheme
+  using MathScheme =
+      FiniteDifferences<T, (order == 3) ? DifferentiationSchemeEnum::WENO3
+                                        : DifferentiationSchemeEnum::WENO5>;
 
   static T pow2(const T &value) { return value * value; }
 
 public:
   static void prepareLS(SmartPointer<viennals::Domain<T, D>> passedlsDomain) {
-    // WENO5 uses a 7-point stencil (radius 3).
-    // Ensure we expand enough layers to access neighbors at +/- 3.
-    static_assert(order >= 3, "WENO5 requires an iterator order of at least 3");
-    viennals::Expand<T, D>(passedlsDomain, 2 * order + 1).apply();
+    // Ensure we expand enough layers to access neighbors.
+    viennals::Expand<T, D>(passedlsDomain, 2 * stencilRadius + 1).apply();
   }
 
-  WENO5(SmartPointer<viennals::Domain<T, D>> passedlsDomain,
-        SmartPointer<viennals::VelocityField<T>> vel, bool calcNormal = true)
+  WENO(SmartPointer<viennals::Domain<T, D>> passedlsDomain,
+       SmartPointer<viennals::VelocityField<T>> vel, bool calcNormal = true)
       : levelSet(passedlsDomain), velocities(vel),
         neighborIterator(levelSet->getDomain()),
         calculateNormalVectors(calcNormal) {}
@@ -66,29 +70,23 @@ template <class T, int D, int order = 3> class WENO5 {
     T wenoGradMinus[D]; // Approximates derivative from left (phi_x^-)
     T wenoGradPlus[D];  // Approximates derivative from right (phi_x^+)
 
-    // Array to hold the stencil values: [x-3, x-2, x-1, Center, x+1, x+2, x+3]
-    T stencil[7];
+    // Array to hold the stencil values
+    T stencil[2 * stencilRadius + 1];
 
     for (int i = 0; i < D; i++) {
       // 1. GATHER STENCIL
       // We map the SparseStarIterator (which uses encoded directions)
       // to the flat array expected by FiniteDifferences.
 
-      // Center (Index 3 in the size-7 array)
-      stencil[3] = neighborIterator.getCenter().getValue();
-
-      // Neighbors +/- 1
-      stencil[4] = neighborIterator.getNeighbor(i).getValue();     // +1
-      stencil[2] = neighborIterator.getNeighbor(i + D).getValue(); // -1
+      // Center
+      stencil[stencilRadius] = neighborIterator.getCenter().getValue();
 
-      // Neighbors +/- 2
-      // Note: SparseStarIterator encodes higher distances sequentially
-      stencil[5] = neighborIterator.getNeighbor(D * 2 + i).getValue();     // +2
-      stencil[1] = neighborIterator.getNeighbor(D * 2 + D + i).getValue(); // -2
-
-      // Neighbors +/- 3
-      stencil[6] = neighborIterator.getNeighbor(D * 4 + i).getValue();     // +3
-      stencil[0] = neighborIterator.getNeighbor(D * 4 + D + i).getValue(); // -3
+      for (int k = 1; k <= stencilRadius; ++k) {
+        stencil[stencilRadius + k] =
+            neighborIterator.getNeighbor((k - 1) * 2 * D + i).getValue();
+        stencil[stencilRadius - k] =
+            neighborIterator.getNeighbor((k - 1) * 2 * D + D + i).getValue();
+      }
 
       // 2. COMPUTE DERIVATIVES
       // Delegate the math to your existing library and store results
@@ -111,7 +109,7 @@ template <class T, int D, int order = 3> class WENO5 {
     T vel_grad = 0.;
 
     // --- Standard Normal Vector Calculation (for velocity lookup) ---
-    Vec3D<T> normalVector = {};
+    Vec3D<T> normalVector{};
     if (calculateNormalVectors) {
       T denominator = 0;
       for (int i = 0; i < D; i++) {
@@ -163,14 +161,16 @@ template <class T, int D, int order = 3> class WENO5 {
   }
 
   void reduceTimeStepHamiltonJacobi(double &MaxTimeStep,
-                                    double gridDelta) const {
-    // --- STABILITY IMPROVEMENT ---
-    // High-order schemes like WENO5 combined with simple time integration (like
-    // the likely Forward Euler used in Advect) can be less stable at CFL=0.5.
-    // We enforce a safety factor here to ensure robustness.
-    constexpr double wenoSafetyFactor = 0.5;
-    MaxTimeStep *= wenoSafetyFactor;
-  }
+                                    double gridDelta) const {}
+  //   // --- STABILITY IMPROVEMENT ---
+  //   // High-order schemes like WENO5 combined with simple time integration
+  //   (like
+  //   // the likely Forward Euler used in Advect) can be less stable at
+  //   CFL=0.5.
+  //   // We enforce a safety factor here to ensure robustness.
+  //   constexpr double wenoSafetyFactor = 0.5;
+  //   MaxTimeStep *= wenoSafetyFactor;
+  // }
 };
 
-} // namespace lsInternal
\ No newline at end of file
+} // namespace lsInternal
diff --git a/python/pyWrap.cpp b/python/pyWrap.cpp
index f231b5c7..6b7420d8 100644
--- a/python/pyWrap.cpp
+++ b/python/pyWrap.cpp
@@ -103,6 +103,7 @@ PYBIND11_MODULE(VIENNALS_MODULE_NAME, module) {
              SpatialSchemeEnum::LOCAL_LAX_FRIEDRICHS_2ND_ORDER)
       .value("STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER",
              SpatialSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER)
+      .value("WENO_3RD_ORDER", SpatialSchemeEnum::WENO_3RD_ORDER)
       .value("WENO_5TH_ORDER", SpatialSchemeEnum::WENO_5TH_ORDER)
       .finalize();
 
diff --git a/python/viennals/_core.pyi b/python/viennals/_core.pyi
index 3adc3a1c..b049d566 100644
--- a/python/viennals/_core.pyi
+++ b/python/viennals/_core.pyi
@@ -376,7 +376,8 @@ class SpatialSchemeEnum(enum.IntEnum):
     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>
+    WENO_3RD_ORDER: typing.ClassVar[SpatialSchemeEnum]  # value = <SpatialSchemeEnum.WENO_3RD_ORDER: 10>
+    WENO_5TH_ORDER: typing.ClassVar[SpatialSchemeEnum]  # value = <SpatialSchemeEnum.WENO_5TH_ORDER: 11>
     @classmethod
     def __new__(cls, value):
         ...
diff --git a/tests/Advection/Advection.cpp b/tests/Advection/Advection.cpp
index 876fdae3..7a3ea988 100644
--- a/tests/Advection/Advection.cpp
+++ b/tests/Advection/Advection.cpp
@@ -2,6 +2,7 @@
 #include <numeric>
 
 #include <lsAdvect.hpp>
+#include <lsCompareCriticalDimensions.hpp>
 #include <lsDomain.hpp>
 #include <lsExpand.hpp>
 #include <lsMakeGeometry.hpp>
@@ -45,7 +46,7 @@ class velocityField : public ls::VelocityField<double> {
 int main() {
 
   constexpr int D = 3;
-  omp_set_num_threads(1);
+  omp_set_num_threads(8);
 
   double gridDelta = 0.6999999;
 
@@ -58,6 +59,7 @@ int main() {
       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_3RD_ORDER,
       ls::SpatialSchemeEnum::WENO_5TH_ORDER};
 
   for (auto scheme : spatialSchemes) {
@@ -97,19 +99,47 @@ int main() {
     advectionKernel.setSaveAdvectionVelocities(true);
 
     double time = 0.;
-    for (unsigned i = 0; time < 1.0 && i < 1e2; ++i) {
+    for (unsigned i = 0; time < 1.0 && i < 50; ++i) {
       advectionKernel.apply();
       time += advectionKernel.getAdvectedTime();
 
-      std::string fileName = std::to_string(i) + ".vtp";
-      auto mesh = ls::Mesh<>::New();
-      ls::ToMesh<double, D>(sphere1, mesh).apply();
-      ls::VTKWriter<double>(mesh, "points_" + fileName).apply();
-      ls::ToSurfaceMesh<double, D>(sphere1, mesh).apply();
-      ls::VTKWriter(mesh, "surface_" + fileName).apply();
+      // std::string fileName = std::to_string(i) + ".vtp";
+      // auto mesh = ls::Mesh<>::New();
+      // ls::ToMesh<double, D>(sphere1, mesh).apply();
+      // ls::VTKWriter<double>(mesh, "points_" + fileName).apply();
+      // ls::ToSurfaceMesh<double, D>(sphere1, mesh).apply();
+      // ls::VTKWriter(mesh, "surface_" + fileName).apply();
     }
 
     LSTEST_ASSERT_VALID_LS(sphere1, double, D)
+
+    // Check critical dimensions against initial state
+    auto sphereRef = ls::Domain<double, D>::New(gridDelta);
+    ls::MakeGeometry<double, D>(sphereRef,
+                                ls::Sphere<double, D>::New(origin, radius))
+        .apply();
+
+    ls::CompareCriticalDimensions<double, D> compare(sphereRef, sphere1);
+    std::array<double, D> minB, maxB;
+    minB.fill(std::numeric_limits<double>::lowest());
+    maxB.fill(std::numeric_limits<double>::max());
+
+    // Measure Max X (Growth velocity ~3.3)
+    compare.addRange(0, minB, maxB, true);
+    // Measure Min X (Growth velocity ~1.5)
+    compare.addRange(0, minB, maxB, false);
+    // Measure Max Y (Growth velocity ~1.0)
+    compare.addRange(1, minB, maxB, true);
+
+    compare.apply();
+
+    double posRef, posSample, diffx, diffy, diffz;
+
+    compare.getCriticalDimensionResult(0, posRef, posSample, diffx);
+    compare.getCriticalDimensionResult(1, posRef, posSample, diffy);
+    compare.getCriticalDimensionResult(2, posRef, posSample, diffz);
+    VC_TEST_ASSERT(diffx > 3.45 && diffx < 3.55 && diffy > 1.45 &&
+                   diffy < 1.65 && diffz > 0.85 && diffz < 1.10);
   }
 
   // std::cout << sphere1->getNumberOfPoints() << std::endl;
diff --git a/tests/Advection/StencilLaxFriedrichsTest.cpp b/tests/Advection/StencilLaxFriedrichsTest.cpp
index 24429b1b..4efc2fef 100644
--- a/tests/Advection/StencilLaxFriedrichsTest.cpp
+++ b/tests/Advection/StencilLaxFriedrichsTest.cpp
@@ -2,6 +2,7 @@
 #include <iostream>
 
 #include <lsAdvect.hpp>
+#include <lsCompareChamfer.hpp>
 #include <lsDomain.hpp>
 #include <lsMakeGeometry.hpp>
 #include <lsTestAsserts.hpp>
@@ -43,14 +44,14 @@ int main() {
   // Create initial level set (Sphere)
   auto sphere =
       ls::SmartPointer<ls::Domain<T, D>>::New(bounds, boundaryCons, gridDelta);
-  T origin[3] = {0.0, 0.0, 0.0};
+  T origin[3]{};
   T radius = 1.0;
   ls::MakeGeometry<T, D>(sphere, ls::Sphere<T, D>::New(origin, radius)).apply();
 
   // Output initial geometry
-  auto mesh = ls::SmartPointer<ls::Mesh<T>>::New();
-  ls::ToSurfaceMesh<T, D>(sphere, mesh).apply();
-  ls::VTKWriter<T>(mesh, "sphere_initial.vtp").apply();
+  // auto mesh = ls::SmartPointer<ls::Mesh<T>>::New();
+  // ls::ToSurfaceMesh<T, D>(sphere, mesh).apply();
+  // ls::VTKWriter<T>(mesh, "sphere_initial.vtp").apply();
 
   // Define constant velocity field
   auto velocityField = ls::SmartPointer<ConstantScalarVelocity<T>>::New();
@@ -72,9 +73,19 @@ int main() {
   // Verify the result is a valid level set
   LSTEST_ASSERT_VALID_LS(sphere, T, D);
 
+  // Check result against analytical solution
+  auto sphereRef =
+      ls::SmartPointer<ls::Domain<T, D>>::New(bounds, boundaryCons, gridDelta);
+  ls::MakeGeometry<T, D>(sphereRef, ls::Sphere<T, D>::New(origin, radius + 0.5))
+      .apply();
+
+  auto chamfer = ls::CompareChamfer<T, D>(sphereRef, sphere);
+  chamfer.apply();
+  VC_TEST_ASSERT(chamfer.getChamferDistance() < 0.035);
+
   // Output final geometry
-  ls::ToSurfaceMesh<T, D>(sphere, mesh).apply();
-  ls::VTKWriter<T>(mesh, "sphere_final.vtp").apply();
+  // ls::ToSurfaceMesh<T, D>(sphere, mesh).apply();
+  // ls::VTKWriter<T>(mesh, "sphere_final.vtp").apply();
 
   std::cout << "Test passed!" << std::endl;
 
diff --git a/tests/Advection/TimeIntegrationComparison.cpp b/tests/Advection/TimeIntegrationComparison.cpp
index 60f713d9..c3d86d0d 100644
--- a/tests/Advection/TimeIntegrationComparison.cpp
+++ b/tests/Advection/TimeIntegrationComparison.cpp
@@ -2,12 +2,14 @@
 #include <iostream>
 
 #include <lsAdvect.hpp>
+#include <lsCompareChamfer.hpp>
 #include <lsDomain.hpp>
 #include <lsMakeGeometry.hpp>
 #include <lsTestAsserts.hpp>
 #include <lsToSurfaceMesh.hpp>
 #include <lsVTKWriter.hpp>
 #include <lsVelocityField.hpp>
+#include <vcTimer.hpp>
 
 namespace ls = viennals;
 
@@ -36,6 +38,7 @@ int main() {
   // Define dimension and type
   constexpr int D = 3;
   using T = double;
+  omp_set_num_threads(8);
 
   // Define grid and domain bounds
   double gridDelta = 0.1;
@@ -47,14 +50,32 @@ int main() {
   // Create initial level set (Sphere)
   auto sphere =
       ls::SmartPointer<ls::Domain<T, D>>::New(bounds, boundaryCons, gridDelta);
-  T origin[3] = {0.0, 0.0, 0.0};
+  T origin[3]{};
   T radius = 1.5;
   ls::MakeGeometry<T, D>(sphere, ls::Sphere<T, D>::New(origin, radius)).apply();
 
+  // Create analytical reference for error calculation at t=2.0.
+  // The sphere starts at x=0 and moves with v=1 for t=2, so it ends at x=2.
+  auto sphereRef =
+      ls::SmartPointer<ls::Domain<T, D>>::New(bounds, boundaryCons, gridDelta);
+  T originRef[3] = {2.0, 0.0, 0.0};
+  ls::MakeGeometry<T, D>(sphereRef, ls::Sphere<T, D>::New(originRef, radius))
+      .apply();
+
   // Create copies for Forward Euler, RK2 and RK3
   auto sphereFE = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
   auto sphereRK2 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
   auto sphereRK3 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO3_FE = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO3_RK2 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO3_RK3 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO5_FE = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO5_RK2 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+  auto sphereWENO5_RK3 = ls::SmartPointer<ls::Domain<T, D>>::New(sphere);
+
+  auto meshInit = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphere, meshInit).apply();
+  ls::VTKWriter<T>(meshInit, "sphereInit.vtp").apply();
 
   // Define constant velocity field (moving in x-direction)
   std::array<T, 3> vel = {1.0, 0.0, 0.0};
@@ -84,30 +105,212 @@ int main() {
   advectRK3.setSpatialScheme(ls::SpatialSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
   advectRK3.setTemporalScheme(ls::TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER);
 
+  // Setup Advection: WENO3 Forward Euler
+  ls::Advect<T, D> advectWENO3_FE;
+  advectWENO3_FE.insertNextLevelSet(sphereWENO3_FE);
+  advectWENO3_FE.setVelocityField(velocityField);
+  advectWENO3_FE.setAdvectionTime(2.0);
+  advectWENO3_FE.setSpatialScheme(ls::SpatialSchemeEnum::WENO_3RD_ORDER);
+  advectWENO3_FE.setTemporalScheme(ls::TemporalSchemeEnum::FORWARD_EULER);
+
+  // Setup Advection: WENO3 Runge-Kutta 2
+  ls::Advect<T, D> advectWENO3_RK2;
+  advectWENO3_RK2.insertNextLevelSet(sphereWENO3_RK2);
+  advectWENO3_RK2.setVelocityField(velocityField);
+  advectWENO3_RK2.setAdvectionTime(2.0);
+  advectWENO3_RK2.setSpatialScheme(ls::SpatialSchemeEnum::WENO_3RD_ORDER);
+  advectWENO3_RK2.setTemporalScheme(
+      ls::TemporalSchemeEnum::RUNGE_KUTTA_2ND_ORDER);
+
+  // Setup Advection: WENO3 Runge-Kutta 3
+  ls::Advect<T, D> advectWENO3_RK3;
+  advectWENO3_RK3.insertNextLevelSet(sphereWENO3_RK3);
+  advectWENO3_RK3.setVelocityField(velocityField);
+  advectWENO3_RK3.setAdvectionTime(2.0);
+  advectWENO3_RK3.setSpatialScheme(ls::SpatialSchemeEnum::WENO_3RD_ORDER);
+  advectWENO3_RK3.setTemporalScheme(
+      ls::TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER);
+
+  // Setup Advection: WENO5 Forward Euler
+  ls::Advect<T, D> advectWENO5_FE;
+  advectWENO5_FE.insertNextLevelSet(sphereWENO5_FE);
+  advectWENO5_FE.setVelocityField(velocityField);
+  advectWENO5_FE.setAdvectionTime(2.0);
+  advectWENO5_FE.setSpatialScheme(ls::SpatialSchemeEnum::WENO_5TH_ORDER);
+  advectWENO5_FE.setTemporalScheme(ls::TemporalSchemeEnum::FORWARD_EULER);
+
+  // Setup Advection: WENO5 Runge-Kutta 2
+  ls::Advect<T, D> advectWENO5_RK2;
+  advectWENO5_RK2.insertNextLevelSet(sphereWENO5_RK2);
+  advectWENO5_RK2.setVelocityField(velocityField);
+  advectWENO5_RK2.setAdvectionTime(2.0);
+  advectWENO5_RK2.setSpatialScheme(ls::SpatialSchemeEnum::WENO_5TH_ORDER);
+  advectWENO5_RK2.setTemporalScheme(
+      ls::TemporalSchemeEnum::RUNGE_KUTTA_2ND_ORDER);
+
+  // Setup Advection: WENO5 Runge-Kutta 3
+  ls::Advect<T, D> advectWENO5_RK3;
+  advectWENO5_RK3.insertNextLevelSet(sphereWENO5_RK3);
+  advectWENO5_RK3.setVelocityField(velocityField);
+  advectWENO5_RK3.setAdvectionTime(2.0);
+  advectWENO5_RK3.setSpatialScheme(ls::SpatialSchemeEnum::WENO_5TH_ORDER);
+  advectWENO5_RK3.setTemporalScheme(
+      ls::TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER);
+
   // Run Advection
   std::cout << "Running Forward Euler Advection..." << std::endl;
+  viennacore::Timer timer;
+  timer.start();
   advectFE.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
   LSTEST_ASSERT_VALID_LS(sphereFE, T, D);
 
   auto meshFE = ls::Mesh<T>::New();
   ls::ToSurfaceMesh<T, D>(sphereFE, meshFE).apply();
   ls::VTKWriter<T>(meshFE, "sphereFE.vtp").apply();
 
+  auto chamferFE = ls::CompareChamfer<T, D>(sphereRef, sphereFE);
+  chamferFE.apply();
+  std::cout << "Chamfer distance: " << chamferFE.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferFE.getChamferDistance() < 0.038);
+
   std::cout << "Running Runge-Kutta 2 Advection..." << std::endl;
+  timer.start();
   advectRK2.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
   LSTEST_ASSERT_VALID_LS(sphereRK2, T, D);
 
   auto meshRK2 = ls::Mesh<T>::New();
   ls::ToSurfaceMesh<T, D>(sphereRK2, meshRK2).apply();
   ls::VTKWriter<T>(meshRK2, "sphereRK2.vtp").apply();
 
+  auto chamferRK2 = ls::CompareChamfer<T, D>(sphereRef, sphereRK2);
+  chamferRK2.apply();
+  std::cout << "Chamfer distance: " << chamferRK2.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferRK2.getChamferDistance() < 0.067);
+
   std::cout << "Running Runge-Kutta 3 Advection..." << std::endl;
+  timer.start();
   advectRK3.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
   LSTEST_ASSERT_VALID_LS(sphereRK3, T, D);
 
   auto meshRK3 = ls::Mesh<T>::New();
   ls::ToSurfaceMesh<T, D>(sphereRK3, meshRK3).apply();
   ls::VTKWriter<T>(meshRK3, "sphereRK3.vtp").apply();
 
+  auto chamferRK3 = ls::CompareChamfer<T, D>(sphereRef, sphereRK3);
+  chamferRK3.apply();
+  std::cout << "Chamfer distance: " << chamferRK3.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferRK3.getChamferDistance() < 0.068);
+
+  std::cout << "Running WENO3 Forward Euler Advection..." << std::endl;
+  timer.start();
+  advectWENO3_FE.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO3_FE, T, D);
+
+  auto meshWENO3_FE = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO3_FE, meshWENO3_FE).apply();
+  ls::VTKWriter<T>(meshWENO3_FE, "sphereWENO3_FE.vtp").apply();
+
+  auto chamferWENO3_FE = ls::CompareChamfer<T, D>(sphereRef, sphereWENO3_FE);
+  chamferWENO3_FE.apply();
+  std::cout << "Chamfer distance: " << chamferWENO3_FE.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO3_FE.getChamferDistance() < 0.025);
+
+  std::cout << "Running WENO3 Runge-Kutta 2 Advection..." << std::endl;
+  timer.start();
+  advectWENO3_RK2.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO3_RK2, T, D);
+
+  auto meshWENO3_RK2 = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO3_RK2, meshWENO3_RK2).apply();
+  ls::VTKWriter<T>(meshWENO3_RK2, "sphereWENO3_RK2.vtp").apply();
+
+  auto chamferWENO3_RK2 = ls::CompareChamfer<T, D>(sphereRef, sphereWENO3_RK2);
+  chamferWENO3_RK2.apply();
+  std::cout << "Chamfer distance: " << chamferWENO3_RK2.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO3_RK2.getChamferDistance() < 0.01);
+
+  std::cout << "Running WENO3 Runge-Kutta 3 Advection..." << std::endl;
+  timer.start();
+  advectWENO3_RK3.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO3_RK3, T, D);
+
+  auto meshWENO3_RK3 = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO3_RK3, meshWENO3_RK3).apply();
+  ls::VTKWriter<T>(meshWENO3_RK3, "sphereWENO3_RK3.vtp").apply();
+
+  auto chamferWENO3_RK3 = ls::CompareChamfer<T, D>(sphereRef, sphereWENO3_RK3);
+  chamferWENO3_RK3.apply();
+  std::cout << "Chamfer distance: " << chamferWENO3_RK3.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO3_RK3.getChamferDistance() < 0.01);
+
+  std::cout << "Running WENO5 Forward Euler Advection..." << std::endl;
+  timer.start();
+  advectWENO5_FE.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO5_FE, T, D);
+
+  auto meshWENO5_FE = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO5_FE, meshWENO5_FE).apply();
+  ls::VTKWriter<T>(meshWENO5_FE, "sphereWENO5_FE.vtp").apply();
+
+  auto chamferWENO5_FE = ls::CompareChamfer<T, D>(sphereRef, sphereWENO5_FE);
+  chamferWENO5_FE.apply();
+  std::cout << "Chamfer distance: " << chamferWENO5_FE.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO5_FE.getChamferDistance() < 0.018);
+
+  std::cout << "Running WENO5 Runge-Kutta 2 Advection..." << std::endl;
+  timer.start();
+  advectWENO5_RK2.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO5_RK2, T, D);
+
+  auto meshWENO5_RK2 = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO5_RK2, meshWENO5_RK2).apply();
+  ls::VTKWriter<T>(meshWENO5_RK2, "sphereWENO5_RK2.vtp").apply();
+
+  auto chamferWENO5_RK2 = ls::CompareChamfer<T, D>(sphereRef, sphereWENO5_RK2);
+  chamferWENO5_RK2.apply();
+  std::cout << "Chamfer distance: " << chamferWENO5_RK2.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO5_RK2.getChamferDistance() < 0.004);
+
+  std::cout << "Running WENO5 Runge-Kutta 3 Advection..." << std::endl;
+  timer.start();
+  advectWENO5_RK3.apply();
+  timer.finish();
+  std::cout << "Time: " << timer.currentDuration / 1e9 << "s" << std::endl;
+  LSTEST_ASSERT_VALID_LS(sphereWENO5_RK3, T, D);
+
+  auto meshWENO5_RK3 = ls::Mesh<T>::New();
+  ls::ToSurfaceMesh<T, D>(sphereWENO5_RK3, meshWENO5_RK3).apply();
+  ls::VTKWriter<T>(meshWENO5_RK3, "sphereWENO5_RK3.vtp").apply();
+
+  auto chamferWENO5_RK3 = ls::CompareChamfer<T, D>(sphereRef, sphereWENO5_RK3);
+  chamferWENO5_RK3.apply();
+  std::cout << "Chamfer distance: " << chamferWENO5_RK3.getChamferDistance()
+            << std::endl;
+  VC_TEST_ASSERT(chamferWENO5_RK3.getChamferDistance() < 0.0034);
+
   return 0;
 }
diff --git a/tests/BooleanOperationExactZero/BooleanOperationExactZero.cpp b/tests/BooleanOperationExactZero/BooleanOperationExactZero.cpp
index dd77d52f..91b4aa0c 100644
--- a/tests/BooleanOperationExactZero/BooleanOperationExactZero.cpp
+++ b/tests/BooleanOperationExactZero/BooleanOperationExactZero.cpp
@@ -76,8 +76,8 @@ int main() {
   // now make substrate ( plane ) at the same height as the bottom of the mask
   auto substrate = LSType::New(mask->getGrid());
   {
-    double origin[D] = {};
-    double normal[D] = {};
+    double origin[D]{};
+    double normal[D]{};
     origin[D - 1] = 1.;
     normal[D - 1] = 1.;
     ls::MakeGeometry(
diff --git a/tests/GeometricAdvectMask/GeometricAdvectMask.cpp b/tests/GeometricAdvectMask/GeometricAdvectMask.cpp
index ba104199..6d3382cb 100644
--- a/tests/GeometricAdvectMask/GeometricAdvectMask.cpp
+++ b/tests/GeometricAdvectMask/GeometricAdvectMask.cpp
@@ -43,7 +43,7 @@ int main() {
   // create mask
   {
     NumericType normal[3] = {0., (D == 2) ? 1. : 0., (D == 3) ? 1. : 0.};
-    NumericType origin[3] = {};
+    NumericType origin[3]{};
     ls::MakeGeometry<NumericType, D>(
         mask, ls::SmartPointer<ls::Plane<NumericType, D>>::New(origin, normal))
         .apply();