Updated visibility to perform directional check

examples/boschProcess/boschProcessCombinedEmulate.py

@@ -1,4 +1,5 @@
 from argparse import ArgumentParser
+import time
 
 # parse config file name and simulation dimension
 parser = ArgumentParser(
@@ -42,25 +43,29 @@
 
 # Define directional rate
 etchDir = vps.RateSet(
-    direction=direction,
+    direction=[-1., -1., 0.],
+    directionalVelocity=params["ionRate"],
+)
+
+# Define directional rate
+etchDir2 = vps.RateSet(
+    direction=[0.25, -1., 0.],
     directionalVelocity=params["ionRate"],
-    isotropicVelocity=0.0,
-    maskMaterials=[vps.Material.Mask]
 )
 
 # Define isotropic rate
 etchIso = vps.RateSet(
     direction=direction,
-    directionalVelocity=0.0,
     isotropicVelocity=params["neutralRate"],
-    maskMaterials=[vps.Material.Mask, vps.Material.Polymer]
+    maskMaterials=[vps.Material.Mask, vps.Material.Polymer],
 )
 
 # List of rate sets
-etchRatesSet = [etchDir, etchIso]
+etchRatesSet = [etchDir, etchDir2]
 
 # Create the DirectionalEtching process with multiple rate sets
 etchModel = vps.DirectionalEtching(rateSets=etchRatesSet)
+# etchModel = vps.DirectionalEtching(rateSet=etchDir2)
 
 numCycles = 4 #int(params["numCycles"])
 n = 0
@@ -69,29 +74,36 @@
 geometry.saveVolumeMesh(f"boschProcessC_{n}")
 n += 1
 
-vps.Process(geometry, etchModel, params["etchTime"]).apply()
-geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
-geometry.saveVolumeMesh(f"boschProcessC_{n}")
-n += 1
-
-for i in range(numCycles):
-    # Deposit a layer of polymer
-    geometry.duplicateTopLevelSet(vps.Material.Polymer)
-    vps.Process(geometry, depoModel, 1.).apply()
-    geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
-    geometry.saveVolumeMesh(f"boschProcessC_{n}")
-    n += 1
-
-    # Etch the trench
-    vps.Process(geometry, etchModel, params["etchTime"]).apply()
-    geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
-    geometry.saveVolumeMesh(f"boschProcessC_{n}")
-    n += 1
-
-    # Ash the polymer
-    geometry.removeTopLevelSet()
+start_time = time.perf_counter()
+tot = 0
+while tot < params["etchTime"]:
+    vps.Process(geometry, etchModel, 1).apply()
     geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
     geometry.saveVolumeMesh(f"boschProcessC_{n}")
     n += 1
+    tot += 1
+end_time = time.perf_counter()
+elapsed_time = end_time - start_time
+print(f"Checking visited cells - execution time: {elapsed_time} seconds")
+
+# for i in range(numCycles):
+#     # Deposit a layer of polymer
+#     geometry.duplicateTopLevelSet(vps.Material.Polymer)
+#     vps.Process(geometry, depoModel, 1.).apply()
+#     geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
+#     geometry.saveVolumeMesh(f"boschProcessC_{n}")
+#     n += 1
+
+#     # Etch the trench
+#     vps.Process(geometry, etchModel, params["etchTime"]).apply()
+#     geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
+#     geometry.saveVolumeMesh(f"boschProcessC_{n}")
+#     n += 1
+
+#     # Ash the polymer
+#     geometry.removeTopLevelSet()
+#     geometry.saveSurfaceMesh("boschProcessC_{}".format(n))
+#     geometry.saveVolumeMesh(f"boschProcessC_{n}")
+#     n += 1
 
 geometry.saveVolumeMesh("finalC")

examples/boschProcess/config.txt

@@ -1,17 +1,17 @@
 # Geometry
-gridDelta = 0.9
+gridDelta = 0.5
 xExtent = 150.0
 yExtent = 50.0
 trenchWidth = 40.0
 maskHeight = 20.0
 
 
 neutralStickingProbability = 0.1
-neutralRate = 1.0
+neutralRate = -1.0
 ionSourceExponent = 200
-ionRate = 1.0
+ionRate = -1.0
 
-etchTime = 5.0
+etchTime = 20.0
 
 depositionThickness = 1.0
 

include/viennaps/models/psDirectionalEtching.hpp

@@ -19,13 +19,17 @@ struct RateSet {
     NumericType directionalVelocity;
     NumericType isotropicVelocity;
     std::vector<int> maskMaterials;
-
-    RateSet(const Vec3D<NumericType>& dir,
-            NumericType dirVel,
-            NumericType isoVel,
-            const std::vector<int>& masks)
+    bool calculateVisibility;
+
+    RateSet(const Vec3D<NumericType>& dir = Vec3D<NumericType>{0., 0., 0.},
+            NumericType dirVel = 0.,
+            NumericType isoVel = 0.,
+            const std::vector<int>& masks = std::vector<int>{Material::Mask}, // Default to Material::Mask
+            // const std::vector<int>& masks = {},
+            bool calcVis = true)
         : direction(dir), directionalVelocity(dirVel),
-          isotropicVelocity(isoVel), maskMaterials(masks) {}
+          isotropicVelocity(isoVel), maskMaterials(masks),
+          calculateVisibility(calcVis && (dir[0] != 0. || dir[1] != 0. || dir[2] != 0.) && dirVel != 0) { }
 };
 
 template <class NumericType, int D>
@@ -49,7 +53,7 @@ class DirectionalEtchVelocityField : public VelocityField<NumericType> {
                 continue; // Skip this rate set if material is masked
             }
             // Accumulate isotropic velocities
-            scalarVelocity -= rateSet.isotropicVelocity;
+            scalarVelocity += std::min(0., rateSet.isotropicVelocity);
         }
 
         return scalarVelocity;
@@ -61,15 +65,32 @@ class DirectionalEtchVelocityField : public VelocityField<NumericType> {
                                          unsigned long pointId) override {
         Vec3D<NumericType> vectorVelocity = {0.0, 0.0, 0.0};
 
-        for (const auto& rateSet : rateSets_) {
+        // for (const auto& rateSet : rateSets_) {
+        for (int rateSetID = 0; rateSetID < rateSets_.size(); ++rateSetID) {
+            const auto& rateSet = rateSets_[rateSetID];
             if (isMaskMaterial(material, rateSet.maskMaterials)) {
                 continue; // Skip this rate set if material is masked
             }
-            if (useVisibilities_ && this->visibilities_->at(pointId) == 0.) {
+            // if (useVisibilities_ && this->visibilities_[rateSetID]->at(pointId) == 0.) {
+            if (rateSet.calculateVisibility && this->visibilities_[rateSetID]->at(pointId) == 0.) {
                 continue; // Skip if visibility check fails
             }
-            // Accumulate directional velocities
-            vectorVelocity = vectorVelocity + rateSet.direction * rateSet.directionalVelocity;
+
+            // Calculate the potential velocity vector for this rate set
+            Vec3D<NumericType> potentialVelocity = rateSet.direction * rateSet.directionalVelocity;
+
+            // Compute dot product manually for std::array
+            NumericType dotProduct = 0.0;
+            for (int i = 0; i < 3; ++i) {
+                dotProduct += potentialVelocity[i] * normalVector[i];
+            }
+
+            if (dotProduct > 0) { // Only include positive dot products (etching)
+                vectorVelocity = vectorVelocity - potentialVelocity;
+            }
+
+            // // Accumulate directional velocities
+            // vectorVelocity = vectorVelocity - rateSet.direction * rateSet.directionalVelocity;
         }
 
         return vectorVelocity;
@@ -79,8 +100,18 @@ class DirectionalEtchVelocityField : public VelocityField<NumericType> {
     // which only depends on an analytic velocity field
     int getTranslationFieldOptions() const override { return 0; }
 
+    // Return direction vector for a specific rateSet
+    Vec3D<NumericType> getDirection(const int rateSetId) const override {
+        const auto& rateSet = rateSets_[rateSetId];
+        return rateSet.direction;
+    }
+
+    bool useVisibilities(const int rateSetId) const override { return rateSets_[rateSetId].calculateVisibility; }
+
     bool useVisibilities() const override { return useVisibilities_; }
 
+    int numRates() const override { return rateSets_.size(); }
+
 private:
     bool isMaskMaterial(const int material, const std::vector<int>& maskMaterials) const {
         return std::find(maskMaterials.begin(), maskMaterials.end(), material) != maskMaterials.end();
@@ -98,19 +129,24 @@ class DirectionalEtching : public ProcessModel<NumericType, D> {
         NumericType directionalVelocity;
         NumericType isotropicVelocity;
         std::vector<Material> maskMaterials;
+        bool calculateVisibility;
 
         RateSet(const Vec3D<NumericType>& dir,
                 NumericType dirVel,
                 NumericType isoVel,
-                const std::vector<Material>& masks)
+                const std::vector<Material>& masks,
+                bool calcVis)
             : direction(dir), directionalVelocity(dirVel),
-              isotropicVelocity(isoVel), maskMaterials(masks) {}
+              isotropicVelocity(isoVel), maskMaterials(masks),
+              calculateVisibility(calcVis) {}
     };
 
-    /// Constructor accepting multiple rate sets
-    DirectionalEtching(const std::vector<RateSet>& rateSets,
-                      const bool useVisibilities = true)
-        : useVisibilities_(useVisibilities) {
+    // Constructor accepting single rate set
+    DirectionalEtching(const RateSet& rateSet) : DirectionalEtching(std::vector<RateSet>{rateSet}) { }
+
+    // Constructor accepting multiple rate sets
+    DirectionalEtching(const std::vector<RateSet>& rateSets) {
+        useVisibilities_ = false;
         // Convert RateSet materials from Material enum to int
         for (const auto& rs : rateSets) {
             std::vector<int> maskInts;
@@ -121,8 +157,10 @@ class DirectionalEtching : public ProcessModel<NumericType, D> {
                 rs.direction,
                 rs.directionalVelocity,
                 rs.isotropicVelocity,
-                maskInts
+                maskInts,
+                rs.calculateVisibility
             );
+            useVisibilities_ = useVisibilities_ || rs.calculateVisibility;
             rateSets_.push_back(internalRateSet);
         }
         initialize();

include/viennaps/psProcess.hpp

@@ -480,13 +480,17 @@ template <typename NumericType, int D> class Process {
 
       // check if visibilities are used
       if (model->getVelocityField()->useVisibilities()) {
-        viennals::CalculateVisibilities<NumericType, D>(
-            domain->getLevelSets().back())
-            .apply();
-        auto visibilities = SmartPointer<std::vector<NumericType>>::New(
-            *domain->getLevelSets().back()->getPointData().getScalarData(
-                "Visibilities"));
-        model->getVelocityField()->setVisibilities(visibilities);
+        for (int rateSetID = 0; rateSetID < model->getVelocityField()->numRates(); ++rateSetID) {
+          if (model->getVelocityField()->useVisibilities(rateSetID)) {
+            viennals::CalculateVisibilities<NumericType, D>(
+                domain->getLevelSets().back(), model->getVelocityField()->getDirection(rateSetID))
+                .apply();
+            auto visibilities = SmartPointer<std::vector<NumericType>>::New(
+                *domain->getLevelSets().back()->getPointData().getScalarData(
+                    "Visibilities"));
+            model->getVelocityField()->setVisibilities(visibilities, rateSetID);
+          }
+        }
       }
 
       auto rates = SmartPointer<viennals::PointData<NumericType>>::New();

include/viennaps/psVelocityField.hpp

@@ -36,8 +36,8 @@ template <typename NumericType> class VelocityField {
   virtual void
   setVelocities(SmartPointer<std::vector<NumericType>> velocities) {}
 
-  void setVisibilities(SmartPointer<std::vector<NumericType>> visibilities) {
-    visibilities_ = visibilities;
+  void setVisibilities(SmartPointer<std::vector<NumericType>> visibilities, int rateSetId = 0) {
+    visibilities_[rateSetId] = visibilities;
   }
 
   // translation field options
@@ -46,10 +46,19 @@ template <typename NumericType> class VelocityField {
   // 2: use kd-tree to translate level set ID to surface ID
   virtual int getTranslationFieldOptions() const { return 1; }
 
+  // Return direction vector for a specific rateSet
+  virtual Vec3D<NumericType> getDirection(const int rateSetId) const { return Vec3D<NumericType>{0., 0., 0.}; }
+
+  // Check if visibilities are defined for a specific rateSet
+  virtual bool useVisibilities(const int rateSetId) const { return false; }
+
   virtual bool useVisibilities() const { return false; }
 
+  virtual int numRates() const { return 0; }
+
 public:
-  SmartPointer<std::vector<NumericType>> visibilities_;
+  std::map<int, SmartPointer<std::vector<NumericType>>> visibilities_;
+  // SmartPointer<std::vector<NumericType>> visibilities_;
 };
 
 template <typename NumericType>

python/pyWrap.cpp

@@ -841,24 +841,27 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
 
 // Expose RateSet struct to Python
 pybind11::class_<DirectionalEtching<T, D>::RateSet>(module, "RateSet")
-    .def(pybind11::init<const std::array<T, 3>&, const T, const T, const std::vector<Material>&>(),
-         pybind11::arg("direction"),
-         pybind11::arg("directionalVelocity"),
-         pybind11::arg("isotropicVelocity"),
-         pybind11::arg("maskMaterials") = std::vector<Material>{Material::Mask})
+    .def(pybind11::init<const std::array<T, 3>&, const T, const T, const std::vector<Material>&, const bool>(),
+         pybind11::arg("direction") = std::array<T, 3>{0., 0., 0.},
+         pybind11::arg("directionalVelocity") = 0.,
+         pybind11::arg("isotropicVelocity") = 0.,
+         pybind11::arg("maskMaterials") = std::vector<Material>{Material::Mask},
+         pybind11::arg("calculateVisibility") = true)
     .def_readwrite("direction", &DirectionalEtching<T, D>::RateSet::direction)
     .def_readwrite("directionalVelocity", &DirectionalEtching<T, D>::RateSet::directionalVelocity)
     .def_readwrite("isotropicVelocity", &DirectionalEtching<T, D>::RateSet::isotropicVelocity)
-    .def_readwrite("maskMaterials", &DirectionalEtching<T, D>::RateSet::maskMaterials);
+    .def_readwrite("maskMaterials", &DirectionalEtching<T, D>::RateSet::maskMaterials)
+    .def_readwrite("calculateVisibility", &DirectionalEtching<T, D>::RateSet::calculateVisibility);
 
 // Expose DirectionalEtching class to Python
 pybind11::class_<DirectionalEtching<T, D>,
                  SmartPointer<DirectionalEtching<T, D>>>(
     module, "DirectionalEtching", processModel)
-    .def(pybind11::init<const std::vector<typename DirectionalEtching<T, D>::RateSet>&,
-                        const bool>(),
-         pybind11::arg("rateSets"),
-         pybind11::arg("useVisibilities") = true)
+    .def(pybind11::init<const std::vector<typename DirectionalEtching<T, D>::RateSet>&>(),
+         pybind11::arg("rateSets"))
+    // Constructor accepting a single rate set
+    .def(pybind11::init<const typename DirectionalEtching<T, D>::RateSet&>(),
+         pybind11::arg("rateSet"))
     .def("disableVisibilityCheck",
          &DirectionalEtching<T, D>::disableVisibilityCheck);