Adaptive time stepping feature in ViennaLS

CMakeLists.txt

@@ -120,7 +120,7 @@ CPMFindPackage(
 
 CPMFindPackage(
   NAME ViennaLS
-  VERSION 5.2.1
+  VERSION 5.3.0
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaLS"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON})
 

examples/boschProcess/boschProcessSimulate.cpp

@@ -8,9 +8,10 @@
 using namespace viennaps;
 constexpr int D = 2;
 using NumericType = double;
+const std::string name = "boschProcessSimulate_";
 
-void etch(SmartPointer<Domain<NumericType, D>> domain,
-          util::Parameters &params) {
+void etch(SmartPointer<Domain<NumericType, D>> domain, util::Parameters &params,
+          int &n) {
   std::cout << "  - Etching - " << std::endl;
   auto etchModel = SmartPointer<MultiParticleProcess<NumericType, D>>::New();
   etchModel->addNeutralParticle(params.get("neutralStickingProbability"));
@@ -29,10 +30,11 @@ void etch(SmartPointer<Domain<NumericType, D>> domain,
         return rate;
       });
   Process<NumericType, D>(domain, etchModel, params.get("etchTime")).apply();
+  domain->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
 }
 
 void punchThrough(SmartPointer<Domain<NumericType, D>> domain,
-                  util::Parameters &params) {
+                  util::Parameters &params, int &n) {
   std::cout << "  - Punching through - " << std::endl;
   NumericType depositionThickness = params.get("depositionThickness");
 
@@ -41,22 +43,25 @@ void punchThrough(SmartPointer<Domain<NumericType, D>> domain,
       -depositionThickness, 1. /* sticking */, params.get("ionSourceExponent"),
       Material::Mask);
   Process<NumericType, D>(domain, depoRemoval, 1.).apply();
+  domain->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
 }
 
 void deposit(SmartPointer<Domain<NumericType, D>> domain,
-             util::Parameters &params) {
+             util::Parameters &params, int &n) {
   std::cout << "  - Deposition - " << std::endl;
   NumericType depositionThickness = params.get("depositionThickness");
   NumericType depositionSticking = params.get("depositionStickingProbability");
   domain->duplicateTopLevelSet(Material::Polymer);
   auto model = SmartPointer<SingleParticleProcess<NumericType, D>>::New(
       depositionThickness, depositionSticking);
   Process<NumericType, D>(domain, model, 1.).apply();
+  domain->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
 }
 
-void ash(SmartPointer<Domain<NumericType, D>> domain) {
+void ash(SmartPointer<Domain<NumericType, D>> domain, int &n) {
   domain->removeTopLevelSet();
   domain->removeStrayPoints();
+  domain->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
 }
 
 int main(int argc, char **argv) {
@@ -87,25 +92,18 @@ int main(int argc, char **argv) {
                              params.get("maskHeight"))
       .apply();
 
-  const NumericType depositionThickness = params.get("depositionThickness");
   const int numCycles = params.get<int>("numCycles");
-  const std::string name = "boschProcessSimulate_";
 
   int n = 0;
   geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
-  etch(geometry, params);
-  geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
+  etch(geometry, params, n);
 
   for (int i = 0; i < numCycles; ++i) {
     std::cout << "Cycle " << i + 1 << std::endl;
-    deposit(geometry, params);
-    geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
-    punchThrough(geometry, params);
-    geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
-    etch(geometry, params);
-    geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
-    ash(geometry);
-    geometry->saveSurfaceMesh(name + std::to_string(n++) + ".vtp");
+    deposit(geometry, params, n);
+    punchThrough(geometry, params, n);
+    etch(geometry, params, n);
+    ash(geometry, n);
   }
 
   geometry->saveVolumeMesh(name + "final");

examples/boschProcess/boschProcessSimulate.py

@@ -65,23 +65,25 @@ def rateFunction(fluxes, material):
 
 etchModel.setRateFunction(rateFunction)
 etchTime = params["etchTime"]
-
 n = 0
 
 
-def runProcess(model, name, time=1.0):
+def saveGeometry(geometry):
     global n
+    geometry.saveSurfaceMesh("boschProcessSimulate_{}".format(n), addInterfaces=True)
+    n += 1
+
+
+def runProcess(model, name, time=1.0):
     print("  - {} - ".format(name))
     ps.Process(geometry, model, time).apply()
-    geometry.saveSurfaceMesh("boschProcessSimulate_{}".format(n))
-    n += 1
+    saveGeometry(geometry)
 
 
 numCycles = int(params["numCycles"])
 
 # Initial geometry
-geometry.saveSurfaceMesh("boschProcessSimulate_{}".format(n))
-n += 1
+saveGeometry(geometry)
 
 runProcess(etchModel, "Etching", etchTime)
 
@@ -101,8 +103,7 @@ def runProcess(model, name, time=1.0):
     # Ash (remove) the polymer
     geometry.removeTopLevelSet()
     geometry.removeStrayPoints()
-    geometry.saveSurfaceMesh("boschProcessSimulate_{}".format(n))
-    n += 1
+    saveGeometry(geometry)
 
 # save the final geometry
 geometry.saveVolumeMesh("boschProcessSimulate_final")

include/viennaps/process/psAdvectionHandler.hpp

@@ -29,7 +29,8 @@ template <typename NumericType, int D> class AdvectionHandler {
         (intSchem != IntegrationScheme::ENGQUIST_OSHER_1ST_ORDER &&
          intSchem != IntegrationScheme::ENGQUIST_OSHER_2ND_ORDER &&
          intSchem != IntegrationScheme::LOCAL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER &&
-         intSchem != IntegrationScheme::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER)) {
+         intSchem != IntegrationScheme::LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER &&
+         intSchem != IntegrationScheme::WENO_5TH_ORDER)) {
       VIENNACORE_LOG_WARNING(
           "Translation field method not supported in combination "
           "with integration scheme.");
@@ -50,6 +51,11 @@ template <typename NumericType, int D> class AdvectionHandler {
     advectionKernel_.setIgnoreVoids(context.advectionParams.ignoreVoids);
     advectionKernel_.setCheckDissipation(
         context.advectionParams.checkDissipation);
+    advectionKernel_.setAdaptiveTimeStepping(
+        context.advectionParams.adaptiveTimeStepping);
+    advectionKernel_.setAdaptiveTimeStepThreshold(
+        context.advectionParams.adaptiveTimeStepThreshold);
+
     // normals vectors are only necessary for analytical velocity fields
     if (translationMethod > 0)
       advectionKernel_.setCalculateNormalVectors(false);

include/viennaps/process/psProcess.hpp

@@ -41,10 +41,12 @@ template <typename NumericType, int D> class Process {
   Process(SmartPointer<Domain<NumericType, D>> domain) : context_{domain} {
     initializeStrategies();
   }
+  template <typename... ParamArgs>
   Process(SmartPointer<Domain<NumericType, D>> domain,
           SmartPointer<ProcessModelBase<NumericType, D>> model,
-          NumericType processDuration = 0.)
+          NumericType processDuration = 0., ParamArgs... params)
       : context_{domain, model, processDuration} {
+    (setParameters(params), ...);
     initializeStrategies();
   }
 

include/viennaps/process/psProcessParams.hpp

@@ -55,9 +55,11 @@ struct AdvectionParameters {
       IntegrationScheme::ENGQUIST_OSHER_1ST_ORDER;
   double timeStepRatio = 0.4999;
   double dissipationAlpha = 1.0;
+  double adaptiveTimeStepThreshold = 0.05;
   bool checkDissipation = true;
   bool velocityOutput = false;
   bool ignoreVoids = false;
+  bool adaptiveTimeStepping = false;
 
   auto toMetaData() const {
     std::unordered_map<std::string, std::vector<double>> metaData;

include/viennaps/psDomain.hpp

@@ -234,11 +234,7 @@ template <class NumericType, int D> class Domain {
 
     levelSets_.pop_back();
     if (materialMap_) {
-      auto newMatMap = MaterialMapType::New();
-      for (std::size_t i = 0; i < levelSets_.size(); i++) {
-        newMatMap->insertNextMaterial(materialMap_->getMaterialAtIdx(i));
-      }
-      materialMap_ = newMatMap;
+      materialMap_->removeMaterial();
     }
   }
 

include/viennaps/psMaterials.hpp

@@ -182,6 +182,12 @@ class MaterialMap {
     map_->insertNextMaterial(static_cast<int>(material));
   }
 
+  void removeMaterial() {
+    if (map_) {
+      map_->removeLastMaterial();
+    }
+  }
+
   // Returns the material at the given index. If the index is out of bounds, it
   // returns Material::GAS.
   [[nodiscard]] Material getMaterialAtIdx(std::size_t idx) const {

include/viennaps/psUtil.hpp

@@ -36,6 +36,8 @@ convertIntegrationScheme(const std::string &s) {
   if (s == "STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER" || s == "SLLF_1")
     return viennals::IntegrationSchemeEnum::
         STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER;
+  if (s == "WENO_5TH_ORDER" || s == "WENO_5")
+    return viennals::IntegrationSchemeEnum::WENO_5TH_ORDER;
   throw std::invalid_argument(
       "The value must be one of the following: "
       "ENGQUIST_OSHER_1ST_ORDER, ENGQUIST_OSHER_2ND_ORDER, "
@@ -45,7 +47,8 @@ convertIntegrationScheme(const std::string &s) {
       "LOCAL_LOCAL_LAX_FRIEDRICHS_2ND_ORDER, "
       "LOCAL_LAX_FRIEDRICHS_1ST_ORDER, "
       "LOCAL_LAX_FRIEDRICHS_2ND_ORDER, "
-      "STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER");
+      "STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER, "
+      "WENO_5TH_ORDER");
 }
 
 [[nodiscard]] inline std::string
@@ -72,6 +75,8 @@ convertIntegrationSchemeToString(viennals::IntegrationSchemeEnum scheme) {
     return "LOCAL_LAX_FRIEDRICHS_2ND_ORDER";
   case viennals::IntegrationSchemeEnum::STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER:
     return "STENCIL_LOCAL_LAX_FRIEDRICHS_1ST_ORDER";
+  case viennals::IntegrationSchemeEnum::WENO_5TH_ORDER:
+    return "WENO_5TH_ORDER";
   default:
     throw std::invalid_argument("Unknown integration scheme.");
   }

python/pyWrap.cpp

@@ -462,6 +462,10 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
       .def_readwrite("checkDissipation", &AdvectionParameters::checkDissipation)
       .def_readwrite("velocityOutput", &AdvectionParameters::velocityOutput)
       .def_readwrite("ignoreVoids", &AdvectionParameters::ignoreVoids)
+      .def_readwrite("adaptiveTimeStepping",
+                     &AdvectionParameters::adaptiveTimeStepping)
+      .def_readwrite("adaptiveTimeStepThreshold",
+                     &AdvectionParameters::adaptiveTimeStepThreshold)
       .def("toMetaData", &AdvectionParameters::toMetaData,
            "Convert the advection parameters to a metadata dict.")
       .def("toMetaDataString", &AdvectionParameters::toMetaDataString,

python/pyWrapDimension.hpp

@@ -1244,6 +1244,28 @@ template <int D> void bindApi(py::module &module) {
       .def(py::init<DomainType>(), py::arg("domain"))
       .def(py::init<DomainType, SmartPointer<ProcessModelBase<T, D>>, T>(),
            py::arg("domain"), py::arg("model"), py::arg("duration") = 0.)
+      .def(py::init([](DomainType domain,
+                       SmartPointer<ProcessModelBase<T, D>> model, T duration,
+                       py::args args) {
+             ProcessTD process(domain, model, duration);
+             for (auto arg : args) {
+               if (py::isinstance<AdvectionParameters>(arg)) {
+                 process.setParameters(arg.cast<AdvectionParameters>());
+               } else if (py::isinstance<RayTracingParameters>(arg)) {
+                 process.setParameters(arg.cast<RayTracingParameters>());
+               } else if (py::isinstance<CoverageParameters>(arg)) {
+                 process.setParameters(arg.cast<CoverageParameters>());
+               } else if (py::isinstance<AtomicLayerProcessParameters>(arg)) {
+                 process.setParameters(
+                     arg.cast<AtomicLayerProcessParameters>());
+               } else {
+                 throw py::type_error(
+                     "Unsupported parameter type for Process constructor");
+               }
+             }
+             return process;
+           }),
+           py::arg("domain"), py::arg("model"), py::arg("duration") = 0.)
       // methods
       .def("apply", &ProcessTD::apply,
            //  py::call_guard<py::gil_scoped_release>(),

python/scripts/install_ViennaPS_linux.py

@@ -9,7 +9,7 @@
 
 REQUIRED_GCC = "12"
 REQUIRED_NVCC_MAJOR = 12
-DEFAULT_VIENNALS_VERSION = "5.2.0"
+DEFAULT_VIENNALS_VERSION = "5.2.1"
 
 
 def run(cmd, **kwargs):

python/viennaps/_core/__init__.pyi

@@ -7,10 +7,10 @@ import collections.abc
 import enum
 import typing
 import viennals._core
-from viennaps import d2
 import viennaps.d2
-import viennaps.d3
+from viennaps import d2
 from viennaps import d3
+import viennaps.d3
 from . import constants
 from . import gpu
 from . import util
@@ -64,6 +64,7 @@ __all__: list[str] = [
 ]
 
 class AdvectionParameters:
+    adaptiveTimeStepping: bool
     checkDissipation: bool
     ignoreVoids: bool
     integrationScheme: viennals._core.IntegrationSchemeEnum
@@ -79,6 +80,10 @@ class AdvectionParameters:
         Convert the advection parameters to a metadata string.
         """
 
+    @property
+    def adaptiveTimeStepThreshold(self) -> float: ...
+    @adaptiveTimeStepThreshold.setter
+    def adaptiveTimeStepThreshold(self, arg0: typing.SupportsFloat) -> None: ...
     @property
     def dissipationAlpha(self) -> float: ...
     @dissipationAlpha.setter

python/viennaps/d2/__init__.pyi

@@ -1200,6 +1200,14 @@ class Process:
         model: ProcessModelBase,
         duration: typing.SupportsFloat = 0.0,
     ) -> None: ...
+    @typing.overload
+    def __init__(
+        self,
+        domain: Domain,
+        model: ProcessModelBase,
+        duration: typing.SupportsFloat = 0.0,
+        *args,
+    ) -> None: ...
     def apply(self) -> None:
         """
         Run the process.

python/viennaps/d3/__init__.pyi

@@ -1205,6 +1205,14 @@ class Process:
         model: ProcessModelBase,
         duration: typing.SupportsFloat = 0.0,
     ) -> None: ...
+    @typing.overload
+    def __init__(
+        self,
+        domain: Domain,
+        model: ProcessModelBase,
+        duration: typing.SupportsFloat = 0.0,
+        *args,
+    ) -> None: ...
     def apply(self) -> None:
         """
         Run the process.

tests/util/testUtil.cpp

@@ -11,6 +11,12 @@ void TestIntegrationSchemeConversion() {
   VC_TEST_ASSERT(convertIntegrationScheme("EO_1") ==
                  viennals::IntegrationSchemeEnum::ENGQUIST_OSHER_1ST_ORDER);
 
+  // Test string to enum
+  VC_TEST_ASSERT(convertIntegrationScheme("WENO_5TH_ORDER") ==
+                 viennals::IntegrationSchemeEnum::WENO_5TH_ORDER);
+  VC_TEST_ASSERT(convertIntegrationScheme("WENO_5") ==
+                 viennals::IntegrationSchemeEnum::WENO_5TH_ORDER);
+
   // Test enum to string
   VC_TEST_ASSERT(
       convertIntegrationSchemeToString(