Fix directional etch model

examples/boschProcess/boschProcessEmulate.cpp

@@ -13,7 +13,7 @@ void etch(ps::SmartPointer<ps::Domain<NumericType, D>> domain,
   direction[D - 1] = -1.;
   auto etchModel =
       ps::SmartPointer<ps::DirectionalEtching<NumericType, D>>::New(
-          direction, params.get("ionRate"), params.get("neutralRate"), false,
+          direction, params.get("ionRate"), params.get("neutralRate"),
           std::vector<ps::Material>{ps::Material::Mask, ps::Material::Polymer});
   ps::Process<NumericType, D>(domain, etchModel, params.get("etchTime"))
       .apply();
@@ -30,7 +30,7 @@ void punchThrough(ps::SmartPointer<ps::Domain<NumericType, D>> domain,
       ps::SmartPointer<ps::DirectionalEtching<NumericType, D>>::New(
           direction,
           params.get("depositionThickness") + params.get("gridDelta") / 2., 0.0,
-          true, ps::Material::Mask);
+          ps::Material::Mask);
   ps::Process<NumericType, D>(domain, depoRemoval, 1.).apply();
 }
 

examples/boschProcess/boschProcessEmulate.py

@@ -2,7 +2,7 @@
 
 # parse config file name and simulation dimension
 parser = ArgumentParser(
-    prog="boschProcess",
+    prog="boschProcessEmulate",
     description="Run a Bosch process on a trench geometry.",
 )
 parser.add_argument("-D", "-DIM", dest="dim", type=int, default=2)
@@ -44,14 +44,12 @@
     direction,
     params["depositionThickness"] + params["gridDelta"] / 2.0,
     0.0,
-    True,
     vps.Material.Mask,
 )
 etchModel = vps.DirectionalEtching(
     direction,
     params["ionRate"],
     params["neutralRate"],
-    False,
     [vps.Material.Mask, vps.Material.Polymer],
 )
 

examples/boschProcess/config.txt

@@ -7,11 +7,11 @@ maskHeight = 20.0
 
 
 neutralStickingProbability = 0.1
-neutralRate = 0.8
+neutralRate = 1.0
 ionSourceExponent = 200
 ionRate = 1.0
 
-etchTime = 6.0
+etchTime = 5.0
 
 depositionThickness = 1.0
 

include/viennaps/models/psDirectionalEtching.hpp

@@ -23,7 +23,7 @@ class DirectionalEtchVelocityField : public VelocityField<NumericType> {
                                const NumericType directionalVelocity,
                                const NumericType isotropicVelocity,
                                const std::vector<int> &mask,
-                               const bool useVisibilities = false)
+                               const bool useVisibilities)
       : direction_(direction), directionalVelocity_(directionalVelocity),
         isotropicVelocity_(isotropicVelocity), maskMaterials_(mask),
         useVisibilities_(useVisibilities) {}
@@ -76,45 +76,58 @@ class DirectionalEtching : public ProcessModel<NumericType, D> {
   DirectionalEtching(const Vec3D<NumericType> &direction,
                      const NumericType directionalVelocity = 1.,
                      const NumericType isotropicVelocity = 0.,
-                     const bool useVisibilities = false,
-                     const Material mask = Material::Mask) {
-    // default surface model
-    auto surfModel = SmartPointer<SurfaceModel<NumericType>>::New();
-
-    // velocity field
-    std::vector<int> maskMaterialsInt = {static_cast<int>(mask)};
-    auto velField =
-        SmartPointer<impl::DirectionalEtchVelocityField<NumericType, D>>::New(
-            direction, directionalVelocity, isotropicVelocity, maskMaterialsInt,
-            useVisibilities);
-
-    this->setSurfaceModel(surfModel);
-    this->setVelocityField(velField);
-    this->setProcessName("DirectionalEtching");
+                     const Material mask = Material::Mask)
+      : direction_(direction), directionalVelocity_(directionalVelocity),
+        isotropicVelocity_(isotropicVelocity) {
+    if (mask != Material::None)
+      maskMaterials_.push_back(static_cast<int>(mask));
+    initialize(direction_, directionalVelocity_, isotropicVelocity_, true,
+               maskMaterials_);
   }
 
   DirectionalEtching(const Vec3D<NumericType> &direction,
                      const NumericType directionalVelocity,
                      const NumericType isotropicVelocity,
-                     const bool useVisibilities,
-                     const std::vector<Material> maskMaterials) {
+                     const std::vector<Material> maskMaterials)
+      : direction_(direction), directionalVelocity_(directionalVelocity),
+        isotropicVelocity_(isotropicVelocity) {
+    for (const auto &mat : maskMaterials) {
+      maskMaterials_.push_back(static_cast<int>(mat));
+    }
+    initialize(direction_, directionalVelocity_, isotropicVelocity_, true,
+               maskMaterials_);
+  }
+
+  void disableVisibilityCheck() {
+    initialize(direction_, directionalVelocity_, isotropicVelocity_, false,
+               maskMaterials_);
+  }
+
+private:
+  void initialize(const Vec3D<NumericType> &direction,
+                  const NumericType directionalVelocity,
+                  const NumericType isotropicVelocity,
+                  const bool useVisibilities,
+                  const std::vector<int> &maskMaterials) {
     // default surface model
     auto surfModel = SmartPointer<SurfaceModel<NumericType>>::New();
 
-    std::vector<int> maskMaterialsInt;
-    for (const auto &mat : maskMaterials) {
-      maskMaterialsInt.push_back(static_cast<int>(mat));
-    }
     // velocity field
     auto velField =
         SmartPointer<impl::DirectionalEtchVelocityField<NumericType, D>>::New(
-            direction, directionalVelocity, isotropicVelocity, maskMaterialsInt,
+            direction, directionalVelocity, isotropicVelocity, maskMaterials,
             useVisibilities);
 
     this->setSurfaceModel(surfModel);
     this->setVelocityField(velField);
     this->setProcessName("DirectionalEtching");
   }
+
+private:
+  Vec3D<NumericType> direction_;
+  NumericType directionalVelocity_;
+  NumericType isotropicVelocity_;
+  std::vector<int> maskMaterials_;
 };
 
 } // namespace viennaps

include/viennaps/psProcess.hpp

@@ -99,6 +99,12 @@ template <typename NumericType, int D> class Process {
     integrationScheme = passedIntegrationScheme;
   }
 
+  // Enable the output of the advection velocities on the level-set mesh.
+  void enableAdvectionVelocityOutput() { lsVelocityOutput = true; }
+
+  // Disable the output of the advection velocities on the level-set mesh.
+  void disableAdvectionVelocityOutput() { lsVelocityOutput = false; }
+
   // Enable the use of random seeds for ray tracing. This is useful to
   // prevent the formation of artifacts in the flux calculation.
   void enableRandomSeeds() { useRandomSeeds_ = true; }
@@ -277,6 +283,7 @@ template <typename NumericType, int D> class Process {
     advectionKernel.setVelocityField(transField);
     advectionKernel.setIntegrationScheme(integrationScheme);
     advectionKernel.setTimeStepRatio(timeStepRatio);
+    advectionKernel.setSaveAdvectionVelocities(lsVelocityOutput);
 
     for (auto dom : domain->getLevelSets()) {
       meshConverter.insertNextLevelSet(dom);
@@ -459,6 +466,7 @@ template <typename NumericType, int D> class Process {
 
     double previousTimeStep = 0.;
     size_t counter = 0;
+    unsigned lsVelCounter = 0;
     Timer rtTimer;
     Timer callbackTimer;
     Timer advTimer;
@@ -623,6 +631,15 @@ template <typename NumericType, int D> class Process {
       advTimer.finish();
       Logger::getInstance().addTiming("Surface advection", advTimer).print();
 
+      if (lsVelocityOutput) {
+        auto lsMesh = SmartPointer<viennals::Mesh<NumericType>>::New();
+        viennals::ToMesh<NumericType, D>(domain->getLevelSets().back(), lsMesh)
+            .apply();
+        viennals::VTKWriter<NumericType>(
+            lsMesh, "ls_velocities_" + std::to_string(lsVelCounter++) + ".vtp")
+            .apply();
+      }
+
       // update the translator to retrieve the correct coverages from the LS
       meshConverter.apply();
       if (useCoverages)
@@ -788,6 +805,7 @@ template <typename NumericType, int D> class Process {
   bool smoothFlux = true;
   NumericType diskRadius = 0.;
   bool ignoreFluxBoundaries = false;
+  bool lsVelocityOutput = false;
   unsigned maxIterations = 20;
   bool coveragesInitialized_ = false;
   NumericType printTime = 0.;

python/pyWrap.cpp

@@ -844,17 +844,17 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
                    SmartPointer<DirectionalEtching<T, D>>>(
       module, "DirectionalEtching", processModel)
       .def(pybind11::init<const std::array<T, 3> &, const T, const T,
-                          const bool, const Material>(),
+                          const Material>(),
            pybind11::arg("direction"),
            pybind11::arg("directionalVelocity") = 1.,
            pybind11::arg("isotropicVelocity") = 0.,
-           pybind11::arg("useVisibilities") = false,
            pybind11::arg("maskMaterial") = Material::None)
       .def(pybind11::init<const std::array<T, 3> &, const T, const T,
-                          const bool, const std::vector<Material>>(),
+                          const std::vector<Material>>(),
            pybind11::arg("direction"), pybind11::arg("directionalVelocity"),
-           pybind11::arg("isotropicVelocity"), pybind11::arg("useVisibilities"),
-           pybind11::arg("maskMaterial"));
+           pybind11::arg("isotropicVelocity"), pybind11::arg("maskMaterial"))
+      .def("disableVisibilityCheck",
+           &DirectionalEtching<T, D>::disableVisibilityCheck);
 
   // Sphere Distribution
   pybind11::class_<SphereDistribution<T, D>,
@@ -1107,6 +1107,12 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
            "Set the integration scheme for solving the level-set equation. "
            "Possible integration schemes are specified in "
            "viennals::IntegrationSchemeEnum.")
+      .def("enableAdvectionVelocityOutput",
+           &Process<T, D>::enableAdvectionVelocityOutput,
+           "Enable the output of the advection velocity field on the ls-mesh.")
+      .def("disableAdvectionVelocityOutput",
+           &Process<T, D>::disableAdvectionVelocityOutput,
+           "Disable the output of the advection velocity field on the ls-mesh.")
       .def("setTimeStepRatio", &Process<T, D>::setTimeStepRatio,
            "Set the CFL condition to use during advection. The CFL condition "
            "sets the maximum distance a surface can be moved during one "

tests/directionalEtch/directionalEtch.cpp

@@ -21,7 +21,8 @@ template <class NumericType, int D> void RunTest() {
     Vec3D<NumericType> direction{0., 0., 0.};
     direction[D - 1] = -1.;
     auto model = SmartPointer<DirectionalEtching<NumericType, D>>::New(
-        direction, 1., 0., false, Material::Mask);
+        direction, 1., 0., Material::Mask);
+    model->disableVisibilityCheck();
 
     VC_TEST_ASSERT(model->getSurfaceModel());
     VC_TEST_ASSERT(model->getVelocityField());
@@ -45,7 +46,7 @@ template <class NumericType, int D> void RunTest() {
     Vec3D<NumericType> direction{0., 0., 0.};
     direction[D - 1] = -1.;
     auto model = SmartPointer<DirectionalEtching<NumericType, D>>::New(
-        direction, 1., 0., true, maskMaterials);
+        direction, 1., 0., maskMaterials);
 
     VC_TEST_ASSERT(model->getSurfaceModel());
     VC_TEST_ASSERT(model->getVelocityField());