fix: adaptive time stepping in advection (#145)

* Limit time step at material interfaces and test with python SquareEtch example

* Update logger

* Reverted material interface behavior and added removeLastMaterial functionality

* format

* Fix material interface stability with adaptive sub-stepping of thin layers (soft landing)

* Added WENO integration scheme

* Add option to enable adaptive time stepping during advection

* Add python bindings

* Bump version

* Add option to set the adaptive time stepping threshold

* Fix MultiSurfaceMesh ctors

* Added Runge-Kutta 3 time integration

* Added Runge Kutta time integration scheme

* Fixed errors from previous commit and format

* Clean up

* Update Python stubs

* Fixed time step differences between FE and RK3

* Format AirGapDeposition example

* rename integration -> discretization where appropriate

* fix formatting

* Updated AirGapDeposition, Deposition, and SquareEtch examples to work in 3D. Added VoidEtching python example. Removed redundant prepareLS(); call in computeRates during advection. Updated RK3 substepping.

* format

* fixed SquareEtch example

* prepareLS(); in computeRates was added again which made it redundant in advect.

* Rename to SpatialScheme, add legacy IntegrationSchemeEnum

* Fixed RK3 time integration scheme to only need velocities on the sparse field

* Fixed RK3 (Strong-Stability preserving Runge-Kutta 3rd order) temporal scheme. The scheme uses a constant velocity during the entire time step, while the gradient is averaged according to SSP-RK3.

* Added compile-time "IntegrationScheme" naming depreciation narning.

* format

* Refactor Advect to use TemporalSchemeEnum and consolidate time integration

* Allow for calculating intermediate velocities during higher order time integration

* Added python bindings to the velocity calculation callback function

* Small fixes in FromMesh and MarkVoidPoints, add velocityUpdateCallback to support velocity calculation during stages in RK integration schemes.

* format

* remove last velocity calculation when in single step mode

* fixed failing PatternedSubstrate and VolumeToLevelSets examples, added 3D functionality to CompareSparseField and fixed ConvexHull generation for 3D.

* format

* remove GIT_TAG for ViennaHRLE

* Added 3D options and tests for comparing domains

* format

* use CompareVolume for 3D and CompareArea for 2D, fix minor errors, improve python tests, add changes to python wrapper and API

* Removed CompareArea from top level Python API, now it is only in d2

* Update python stubs

* Small fixes

* Correct inconsistency in custom z-increments.

* Small fixes and improvements to tests.

* Minor fix

* Enhance test for lsCompareCriticalDimensions

* Fixed RK2/RK3 workflow

* fix faulty multiMaterial averaging in combineLevelSets

* format

* add multi material etch test

* Fix adaptive time stepping logic for FE

* format

* refactor: use emplace_back, reduce test gridDelta

---------

Co-authored-by: Tobias Reiter <t.reiter1@live.de>
Co-authored-by: Roman Kostal <roman.kstl@gmail.com>

Author: 3 people

examples/AirGapDeposition/AirGapDeposition.cpp

@@ -90,7 +90,7 @@ double runSimulation(AdvectKernelType &kernel,
 int main() {
 
   constexpr int D = 2;
-  omp_set_num_threads(8);
+  omp_set_num_threads(16);
 
   NumericType extent = 30;
   NumericType gridDelta = 0.5;

include/viennals/lsAdvect.hpp

@@ -78,7 +78,7 @@ template <class T, int D> class Advect {
   bool adaptiveTimeStepping = false;
   unsigned adaptiveTimeStepSubdivisions = 20;
   static constexpr double wrappingLayerEpsilon = 1e-4;
-  SmartPointer<Domain<T, D>> originalLevelSet = nullptr;
+  std::vector<SmartPointer<Domain<T, D>>> initialLevelSets;
   std::function<bool(SmartPointer<Domain<T, D>>)> velocityUpdateCallback =
       nullptr;
 
@@ -196,41 +196,43 @@ template <class T, int D> class Advect {
 
   // Helper function for linear combination:
   // target = wTarget * target + wSource * source
-  void combineLevelSets(double wTarget, double wSource) {
+  bool combineLevelSets(T wTarget, T wSource) {
+    // Calculate required expansion width based on CFL and RK steps
+    int steps = 1;
+    if (temporalScheme == TemporalSchemeEnum::RUNGE_KUTTA_2ND_ORDER) {
+      steps = 2;
+    } else if (temporalScheme == TemporalSchemeEnum::RUNGE_KUTTA_3RD_ORDER) {
+      steps = 3;
+    }
 
-    auto &domainDest = levelSets.back()->getDomain();
-    auto &grid = levelSets.back()->getGrid();
+    // Expand both level sets to ensure sufficient overlap
+    int expansionWidth = std::ceil(2.0 * steps * timeStepRatio + 1);
+    viennals::Expand<T, D>(levelSets.back(), expansionWidth).apply();
+    viennals::Expand<T, D>(initialLevelSets.back(), expansionWidth).apply();
+
+    bool movedDown = false;
+
+    viennals::BooleanOperation<T, D> op(levelSets.back(),
+                                        initialLevelSets.back(),
+                                        viennals::BooleanOperationEnum::CUSTOM);
+    op.setBooleanOperationComparator(
+        [wTarget, wSource, &movedDown](const T &a, const T &b) {
+          if (a != Domain<T, D>::POS_VALUE && a != Domain<T, D>::NEG_VALUE &&
+              b != Domain<T, D>::POS_VALUE && b != Domain<T, D>::NEG_VALUE) {
+            T res = wSource * a + wTarget * b;
+            if (res > b)
+              movedDown = true;
+            return std::make_pair(res, true);
+          }
+          if (a == Domain<T, D>::POS_VALUE || a == Domain<T, D>::NEG_VALUE)
+            return std::make_pair(a, false);
+          return std::make_pair(b, false);
+        });
+    op.apply();
 
-#pragma omp parallel num_threads(domainDest.getNumberOfSegments())
-    {
-      int p = 0;
-#ifdef _OPENMP
-      p = omp_get_thread_num();
-#endif
-      auto &segDest = domainDest.getDomainSegment(p);
-
-      viennahrle::Index<D> start = (p == 0)
-                                       ? grid.getMinGridPoint()
-                                       : domainDest.getSegmentation()[p - 1];
-      viennahrle::Index<D> end =
-          (p != static_cast<int>(domainDest.getNumberOfSegments()) - 1)
-              ? domainDest.getSegmentation()[p]
-              : grid.incrementIndices(grid.getMaxGridPoint());
+    rebuildLS();
 
-      ConstSparseIterator itDest(domainDest, start);
-      ConstSparseIterator itTarget(originalLevelSet->getDomain(), start);
-
-      unsigned definedValueIndex = 0;
-      for (; itDest.getStartIndices() < end; ++itDest) {
-        if (itDest.isDefined()) {
-          itTarget.goToIndicesSequential(itDest.getStartIndices());
-          T valSource = itDest.getValue();
-          T valTarget = itTarget.getValue();
-          segDest.definedValues[definedValueIndex++] =
-              wTarget * valTarget + wSource * valSource;
-        }
-      }
-    }
+    return movedDown;
   }
 
   void rebuildLS() {
@@ -450,6 +452,7 @@ template <class T, int D> class Advect {
     }
     const bool ignoreVoidPoints = ignoreVoids;
     const bool useAdaptiveTimeStepping = adaptiveTimeStepping;
+    const auto adaptiveFactor = 1.0 / adaptiveTimeStepSubdivisions;
 
     if (!storedRates.empty()) {
       VIENNACORE_LOG_WARNING("Advect: Overwriting previously stored rates.");
@@ -473,8 +476,9 @@ template <class T, int D> class Advect {
               : grid.incrementIndices(grid.getMaxGridPoint());
 
       double tempMaxTimeStep = maxTimeStep;
-      auto &tempRates = storedRates[p];
-      tempRates.reserve(
+      // store the rates and value of underneath LS for this segment
+      auto &rates = storedRates[p];
+      rates.reserve(
           topDomain.getNumberOfPoints() /
               static_cast<double>((levelSets.back())->getNumberOfSegments()) +
           10);
@@ -527,15 +531,14 @@ template <class T, int D> class Advect {
             // Case 1: Growth / Deposition (Velocity > 0)
             // Limit the time step based on the standard CFL condition.
             maxStepTime += cfl / velocity;
-            tempRates.push_back(std::make_pair(gradNDissipation,
-                                               -std::numeric_limits<T>::max()));
+            rates.emplace_back(gradNDissipation,
+                               -std::numeric_limits<T>::max());
             break;
           } else if (velocity == 0.) {
             // Case 2: Static (Velocity == 0)
             // No time step limit imposed by this point.
             maxStepTime = std::numeric_limits<T>::max();
-            tempRates.push_back(std::make_pair(gradNDissipation,
-                                               std::numeric_limits<T>::max()));
+            rates.emplace_back(gradNDissipation, std::numeric_limits<T>::max());
             break;
           } else {
             // Case 3: Etching (Velocity < 0)
@@ -555,30 +558,30 @@ template <class T, int D> class Advect {
               // Sub-case 3a: Standard Advection
               // Far from interface: Use full CFL time step.
               maxStepTime -= cfl / velocity;
-              tempRates.push_back(std::make_pair(
-                  gradNDissipation, std::numeric_limits<T>::max()));
+              rates.emplace_back(gradNDissipation,
+                                 std::numeric_limits<T>::max());
               break;
-
             } else {
               // Sub-case 3b: Interface Interaction
-              auto adaptiveFactor = 1.0 / adaptiveTimeStepSubdivisions;
-              if (useAdaptiveTimeStepping && difference > 0.2 * cfl) {
+              // Use adaptiveFactor as threshold.
+              if (useAdaptiveTimeStepping &&
+                  difference > adaptiveFactor * cfl) {
                 // Adaptive Sub-stepping:
                 // Approaching boundary: Force small steps to gather
                 // flux statistics and prevent numerical overshoot ("Soft
                 // Landing").
                 maxStepTime -= adaptiveFactor * cfl / velocity;
-                tempRates.push_back(std::make_pair(
-                    gradNDissipation, std::numeric_limits<T>::min()));
+                rates.emplace_back(gradNDissipation,
+                                   std::numeric_limits<T>::max());
+                break;
               } else {
                 // Terminal Step:
                 // Within tolerance: Snap to boundary, consume budget, and
                 // switch material.
-                tempRates.push_back(
-                    std::make_pair(gradNDissipation, valueBelow));
                 cfl -= difference;
                 value = valueBelow;
                 maxStepTime -= difference / velocity;
+                rates.emplace_back(gradNDissipation, valueBelow);
               }
             }
           }
@@ -732,7 +735,8 @@ template <class T, int D> class Advect {
         // set the LS value to the one below
         auto const [gradient, dissipation] = itRS->first;
         T velocity = gradient - dissipation;
-        // check if dissipation is too high
+        // check if dissipation is too high and would cause a change in
+        // direction of the velocity
         if (checkDiss && (gradient < 0 && velocity > 0) ||
             (gradient > 0 && velocity < 0)) {
           velocity = 0;

include/viennals/lsAdvectIntegrationSchemes.hpp

@@ -43,110 +43,128 @@ namespace lsInternal {
 template <class T, int D> struct AdvectTimeIntegration {
   using AdvectType = viennals::Advect<T, D>;
 
-  static double evolveForwardEuler(AdvectType &kernel, double maxTimeStep) {
+  static double evolveForwardEuler(AdvectType &kernel, double maxTimeStep,
+                                   bool updateLowerLayers = true) {
     if (kernel.currentTimeStep < 0. || kernel.storedRates.empty())
       kernel.computeRates(maxTimeStep);
 
     kernel.updateLevelSet(kernel.currentTimeStep);
 
     kernel.rebuildLS();
 
-    kernel.adjustLowerLayers();
+    if (updateLowerLayers)
+      kernel.adjustLowerLayers();
 
     return kernel.currentTimeStep;
   }
 
   static double evolveRungeKutta2(AdvectType &kernel, double maxTimeStep) {
     // TVD Runge-Kutta 2nd Order (Heun's Method)
-    // 1. Determine time step
-    kernel.computeRates(maxTimeStep);
-    const double dt = kernel.getCurrentTimeStep();
-
-    // 2. Save u^n
-    if (kernel.originalLevelSet == nullptr) {
-      kernel.originalLevelSet =
-          viennals::Domain<T, D>::New(kernel.levelSets.back()->getGrid());
-    }
-    kernel.originalLevelSet->deepCopy(kernel.levelSets.back());
 
-    if (dt <= 0)
-      return 0.;
+    // Save initial level sets
+    if (kernel.initialLevelSets.size() != kernel.levelSets.size()) {
+      kernel.initialLevelSets.resize(kernel.levelSets.size());
+      for (auto &ls : kernel.initialLevelSets)
+        ls = viennals::Domain<T, D>::New(kernel.levelSets[0]->getGrid());
+    }
+    kernel.initialLevelSets.back()->deepCopy(kernel.levelSets.back());
+
+    // Save initial lower level sets only if Stage 1 will modify them (via
+    // callback)
+    if (kernel.velocityUpdateCallback) {
+      for (size_t i = 0; i < kernel.levelSets.size() - 1; ++i) {
+        kernel.initialLevelSets[i]->deepCopy(kernel.levelSets[i]);
+      }
+    }
 
     // Stage 1: u^(1) = u^n + dt * L(u^n)
-    kernel.updateLevelSet(dt);
+    // Update lower layers only if we have a callback
+    double dt1 = evolveForwardEuler(kernel, maxTimeStep,
+                                    kernel.velocityUpdateCallback != nullptr);
+
+    if (dt1 <= 0.)
+      return 0.;
 
     if (kernel.velocityUpdateCallback)
       kernel.velocityUpdateCallback(kernel.levelSets.back());
 
     // Stage 2: u^(n+1) = 1/2 u^n + 1/2 (u^(1) + dt * L(u^(1)))
     // Current level set is u^(1). Compute L(u^(1)).
-    kernel.computeRates(dt);
     // Update to u* = u^(1) + dt * L(u^(1))
-    kernel.updateLevelSet(dt);
-    // Combine: u^(n+1) = 0.5 * u^n + 0.5 * u*
-    kernel.combineLevelSets(0.5, 0.5);
+    double dt2 = evolveForwardEuler(kernel, dt1, false);
 
-    // If we are in single step mode, the level set will be rebuilt immediately
-    // after this, invalidating the velocity field. Thus, we skip the update.
-    if (kernel.velocityUpdateCallback && !kernel.performOnlySingleStep)
-      kernel.velocityUpdateCallback(kernel.levelSets.back());
-
-    // Finalize
-    kernel.rebuildLS();
+    // Combine: u^(n+1) = 0.5 * u^n + 0.5 * u*
+    bool etched = kernel.combineLevelSets(0.5, 0.5);
 
+    // Restore lower level sets if etched.
+    if (etched && kernel.velocityUpdateCallback) {
+      for (size_t i = 0; i < kernel.levelSets.size() - 1; ++i) {
+        kernel.levelSets[i]->deepCopy(kernel.initialLevelSets[i]);
+      }
+    }
     kernel.adjustLowerLayers();
 
-    return dt;
+    return 0.5 * dt1 + 0.5 * dt2;
   }
 
   static double evolveRungeKutta3(AdvectType &kernel, double maxTimeStep) {
-    // 1. Determine the single time step 'dt' for all stages.
-    kernel.computeRates(maxTimeStep);
-    const double dt = kernel.getCurrentTimeStep();
-
-    // 2. Save u^n (Deep copy to preserve topology)
-    if (kernel.originalLevelSet == nullptr) {
-      kernel.originalLevelSet =
-          viennals::Domain<T, D>::New(kernel.levelSets.back()->getGrid());
+    // Save initial level sets
+    if (kernel.initialLevelSets.size() != kernel.levelSets.size()) {
+      kernel.initialLevelSets.resize(kernel.levelSets.size());
+      for (auto &ls : kernel.initialLevelSets)
+        ls = viennals::Domain<T, D>::New(kernel.levelSets[0]->getGrid());
+    }
+    for (size_t i = 0; i < kernel.levelSets.size(); ++i) {
+      kernel.initialLevelSets[i]->deepCopy(kernel.levelSets[i]);
     }
-    kernel.originalLevelSet->deepCopy(kernel.levelSets.back());
-
-    // If dt is 0 or negative, no advection is possible or needed.
-    if (dt <= 0)
-      return 0.;
 
     // Stage 1: u^(1) = u^n + dt * L(u^n)
-    kernel.updateLevelSet(dt);
+    // This calculates dt based on u^n and advances to u^1.
+    double dt1 = evolveForwardEuler(kernel, maxTimeStep,
+                                    kernel.velocityUpdateCallback != nullptr);
+
+    if (dt1 <= 0.)
+      return 0.;
 
     if (kernel.velocityUpdateCallback)
       kernel.velocityUpdateCallback(kernel.levelSets.back());
 
     // Stage 2: u^(2) = 3/4 u^n + 1/4 (u^(1) + dt * L(u^(1)))
-    kernel.computeRates(dt);
-    kernel.updateLevelSet(dt);
+    // u* = u^(1) + dt * L(u^(1))
+    double dt2 = evolveForwardEuler(kernel, dt1, false);
+
     // Combine to get u^(2) = 0.75 * u^n + 0.25 * u*.
-    kernel.combineLevelSets(0.75, 0.25);
+    bool etched1 = kernel.combineLevelSets(0.75, 0.25);
 
-    if (kernel.velocityUpdateCallback)
-      kernel.velocityUpdateCallback(kernel.levelSets.back());
+    // Restore lower level sets if etched
+    if (etched1 && kernel.velocityUpdateCallback) {
+      for (size_t i = 0; i < kernel.levelSets.size() - 1; ++i) {
+        kernel.levelSets[i]->deepCopy(kernel.initialLevelSets[i]);
+      }
+    }
 
-    // Stage 3: u^(n+1) = 1/3 u^n + 2/3 (u^(2) + dt * L(u^(2)))
-    kernel.computeRates(dt);
-    kernel.updateLevelSet(dt);
-    // Combine to get u^(n+1) = 1/3 * u^n + 2/3 * u**.
-    kernel.combineLevelSets(1.0 / 3.0, 2.0 / 3.0);
+    kernel.adjustLowerLayers();
 
-    // If we are in single step mode, the level set will be rebuilt immediately
-    // after this, invalidating the velocity field. Thus, we skip the update.
-    if (kernel.velocityUpdateCallback && !kernel.performOnlySingleStep)
+    if (kernel.velocityUpdateCallback) {
       kernel.velocityUpdateCallback(kernel.levelSets.back());
+    }
 
-    // Finalize: Re-segment and renormalize the final result.
-    kernel.rebuildLS();
+    // Stage 3: u^(n+1) = 1/3 u^n + 2/3 (u^(2) + dt * L(u^(2)))
+    // u** = u^(2) + dt * L(u^(2))
+    double dt3 = evolveForwardEuler(kernel, dt1, false);
 
+    // Combine to get u^(n+1) = 1/3 * u^n + 2/3 * u**.
+    bool etched2 = kernel.combineLevelSets(1.0 / 3.0, 2.0 / 3.0);
+
+    // Restore lower level sets if etched.
+    if (etched2) {
+      for (size_t i = 0; i < kernel.levelSets.size() - 1; ++i) {
+        kernel.levelSets[i]->deepCopy(kernel.initialLevelSets[i]);
+      }
+    }
     kernel.adjustLowerLayers();
 
-    return dt;
+    return (dt1 + dt2 + 4.0 * dt3) / 6.0;
   }
 };
 

include/viennals/lsBooleanOperation.hpp

@@ -11,6 +11,8 @@
 #include <vcSmartPointer.hpp>
 #include <vcVectorType.hpp>
 
+#include <functional>
+
 namespace viennals {
 
 using namespace viennacore;
@@ -42,7 +44,8 @@ enum struct BooleanOperationEnum : unsigned {
 ///  ComparatorType.
 template <class T, int D> class BooleanOperation {
 public:
-  using ComparatorType = std::pair<T, bool> (*)(const T &, const T &);
+  using ComparatorType =
+      std::function<std::pair<T, bool>(const T &, const T &)>;
 
 private:
   typedef typename Domain<T, D>::DomainType hrleDomainType;