refactor: use concepts and minor performance improvements

examples/holeEtching/config.txt

@@ -30,7 +30,7 @@ A_Si=7         # Si yield constant
 etchStopDepth=-10   # maximum etching depth
 spatialScheme=EO_1
 temporalScheme=RK3
-fluxEngine=CD
+fluxEngine=GT
 
 raysPerPoint=1000
 outputFile=final.vtp

examples/holeEtching/holeEtching.cpp

@@ -10,7 +10,7 @@ int main(int argc, char *argv[]) {
   using NumericType = double;
   constexpr int D = 2;
 
-  Logger::setLogLevel(LogLevel::ERROR);
+  Logger::setLogLevel(LogLevel::WARNING);
   omp_set_num_threads(8);
 
   // Parse the parameters

examples/holeEtching/holeEtching.py

@@ -20,23 +20,22 @@
 ps.Length.setUnit(params["lengthUnit"])
 ps.Time.setUnit(params["timeUnit"])
 
+
 def run_simulation(intermediate_velocities, suffix):
     # geometry setup, all units in um
     geometry = ps.Domain(
         gridDelta=params["gridDelta"],
         xExtent=params["xExtent"],
         yExtent=params["yExtent"],
     )
-
-    hole_shape = ps.HoleShape.QUARTER if args.dim == 3 else ps.HoleShape.HALF
 
     ps.MakeHole(
         domain=geometry,
         holeRadius=params["holeRadius"],
         holeDepth=0.0,
         maskHeight=params["maskHeight"],
         maskTaperAngle=params["taperAngle"],
-        holeShape=hole_shape,
+        holeShape=ps.HoleShape.QUARTER,
     ).apply()
 
     # use pre-defined model SF6O2 etching model
@@ -59,12 +58,8 @@ def run_simulation(intermediate_velocities, suffix):
     rayParams.raysPerPoint = int(params["raysPerPoint"])
 
     advParams = ps.AdvectionParameters()
-    advParams.spatialScheme = ps.util.convertSpatialScheme(
-        params["spatialScheme"]
-    )
-    advParams.temporalScheme = ps.util.convertTemporalScheme(
-        params["temporalScheme"]
-    )
+    advParams.spatialScheme = ps.util.convertSpatialScheme(params["spatialScheme"])
+    advParams.temporalScheme = ps.util.convertTemporalScheme(params["temporalScheme"])
     advParams.calculateIntermediateVelocities = intermediate_velocities
 
     # process setup
@@ -93,7 +88,10 @@ def run_simulation(intermediate_velocities, suffix):
     else:
         output_file += suffix
 
-    geometry.saveSurfaceMesh(filename=output_file, addInterfaces=True, boolMaterials=True)
+    geometry.saveSurfaceMesh(
+        filename=output_file, addInterfaces=True, boolMaterials=True
+    )
+
 
 print("Running simulation without intermediate velocity calculation...")
 run_simulation(False, "_noIntermediate")

examples/simpleEtching/simpleEtching.py

@@ -1,18 +1,17 @@
-import viennaps
-import viennaps.d2 as vps2
-import viennaps.d3 as vps3
+import viennaps as ps
+
+
+ps.setDimension(3)  # Set to 3 for hole etching in 3D
+ps.Logger.setLogLevel(ps.LogLevel.INFO)
+
 
 def main():
-    # Dimension of the domain: 2 for Trench, 3 for Hole
-    D = 2
-
-    viennaps.Logger.setLogLevel(viennaps.LogLevel.INFO)
 
     # Geometry parameters
     grid_delta = 0.05
     x_extent = 3.0
     y_extent = 3.0
-    feature_width = 1.0 # Diameter for hole, Width for trench
+    feature_width = 1.0  # Diameter for hole, Width for trench
     mask_height = 1.0
     taper_angle = 0.0
 
@@ -24,47 +23,43 @@ def main():
     process_time = 1.0
 
     def run_simulation(temporal_scheme, calc_intermediate):
-        suffix = str(temporal_scheme).split(".")[-1]
+        suffix = temporal_scheme.name
         if calc_intermediate:
             suffix += "_recalc"
-        if D == 3:
-            # Create a Hole in 3D
-            domain = vps3.Domain(grid_delta, x_extent, y_extent)
-            vps3.MakeHole(domain, feature_width / 2.0, 0.0, 0.0,
-                          mask_height, taper_angle,
-                          viennaps.HoleShape.QUARTER).apply()
-
-            model = vps3.SingleParticleProcess(rate, sticking_probability, source_power, viennaps.Material.Mask)
-            process = vps3.Process(domain, model, process_time)
-        else:
-            # Create a Trench in 2D
-            domain = vps2.Domain(grid_delta, x_extent, y_extent)
-            vps2.MakeTrench(domain, feature_width, 0.0, 0.0, mask_height,
-                            taper_angle, False).apply()
-
-            model = vps2.SingleParticleProcess(rate, sticking_probability, source_power, viennaps.Material.Mask)
-            process = vps2.Process(domain, model, process_time)
 
-        advection_params = viennaps.AdvectionParameters()
-        advection_params.spatialScheme = viennaps.util.convertSpatialScheme("WENO_5TH_ORDER")
+        # Create a Hole in 3D (defaults to trench in 2D)
+        domain = ps.Domain(grid_delta, x_extent, y_extent)
+        ps.MakeHole(
+            domain,
+            feature_width / 2.0,
+            0.0,
+            0.0,
+            mask_height,
+            taper_angle,
+        ).apply()
+
+        model = ps.SingleParticleProcess(
+            rate, sticking_probability, source_power, ps.Material.Mask
+        )
+        process = ps.Process(domain, model, process_time)
+
+        advection_params = ps.AdvectionParameters()
+        advection_params.spatialScheme = ps.SpatialScheme.WENO_5TH_ORDER
         advection_params.temporalScheme = temporal_scheme
         advection_params.calculateIntermediateVelocities = calc_intermediate
         process.setParameters(advection_params)
 
-        fluxEngine = viennaps.util.convertFluxEngineType("GT")
-        process.setFluxEngineType(fluxEngine)
-
-        viennaps.Logger.getInstance().addInfo(f"Running simulation: {suffix}").print()
+        ps.Logger.getInstance().addInfo(f"Running simulation: {suffix}").print()
         process.apply()
 
         domain.saveSurfaceMesh(f"simpleEtching_{suffix}.vtp")
 
-    run_simulation(viennaps.util.convertTemporalScheme("FORWARD_EULER"), False)
-    run_simulation(viennaps.util.convertTemporalScheme("RUNGE_KUTTA_2ND_ORDER"), False)
-    run_simulation(viennaps.util.convertTemporalScheme("RUNGE_KUTTA_2ND_ORDER"), True)
-    run_simulation(viennaps.util.convertTemporalScheme("RUNGE_KUTTA_3RD_ORDER"), False)
-    run_simulation(viennaps.util.convertTemporalScheme("RUNGE_KUTTA_3RD_ORDER"), True)
+    run_simulation(ps.util.convertTemporalScheme("FORWARD_EULER"), False)
+    run_simulation(ps.util.convertTemporalScheme("RUNGE_KUTTA_2ND_ORDER"), False)
+    run_simulation(ps.util.convertTemporalScheme("RUNGE_KUTTA_2ND_ORDER"), True)
+    run_simulation(ps.util.convertTemporalScheme("RUNGE_KUTTA_3RD_ORDER"), False)
+    run_simulation(ps.util.convertTemporalScheme("RUNGE_KUTTA_3RD_ORDER"), True)
+
 
 if __name__ == "__main__":
     main()
-

include/viennaps/gds/psGDSReader.hpp

@@ -1,8 +1,10 @@
 #pragma once
 
+#include <cstdint>
 #include <cstdio>
 #include <cstdlib>
 #include <string>
+#include <unordered_set>
 
 #include "psGDSGeometry.hpp"
 #include "psGDSUtils.hpp"
@@ -63,14 +65,17 @@ template <typename NumericType, int D> class GDSReader {
   float currentWidth = 0;
   char *tempStr = nullptr;
 
-  static bool contains(int32_t X, int32_t Y,
-                       const std::vector<std::array<int32_t, 2>> &uniPoints) {
-    for (const auto &p : uniPoints) {
-      if (p[0] == X && p[1] == Y)
-        return true;
-    }
+  static constexpr std::uint64_t packPointKey(int32_t x, int32_t y) noexcept {
+    return (static_cast<std::uint64_t>(static_cast<std::uint32_t>(x)) << 32) |
+           static_cast<std::uint32_t>(y);
+  }
 
-    return false;
+  // Returns true if the point already existed; inserts otherwise.
+  static bool contains(int32_t x, int32_t y,
+                       std::unordered_set<std::uint64_t> &uniquePointKeys) {
+    const auto key = packPointKey(x, y);
+    const auto [_, inserted] = uniquePointKeys.insert(key);
+    return !inserted;
   }
 
   void resetCurrentStructure() {
@@ -204,7 +209,8 @@ template <typename NumericType, int D> class GDSReader {
   void parseXYBoundary() {
     const unsigned int numPoints = currentRecordLen / 8;
     auto &currentElPointCloud = currentStructure.elements.back().pointCloud;
-    std::vector<std::array<int32_t, 2>> uniquePoints;
+    std::unordered_set<std::uint64_t> uniquePoints;
+    uniquePoints.reserve(numPoints);
 
     // do not include the last point since it
     // is just a copy of the first
@@ -213,7 +219,6 @@ template <typename NumericType, int D> class GDSReader {
       const auto pY = readFourByteSignedInt();
 
       if (!contains(pX, pY, uniquePoints)) {
-        uniquePoints.push_back({pX, pY});
 
         float X = units * static_cast<float>(pX);
         float Y = units * static_cast<float>(pY);

include/viennaps/models/psIonBeamEtching.hpp

@@ -40,6 +40,7 @@ class IBESurfaceModel : public SurfaceModel<NumericType> {
     auto velocity =
         SmartPointer<std::vector<NumericType>>::New(materialIds.size(), 0.);
     auto flux = rates->getScalarData(fluxLabel);
+    assert(flux && flux->size() == materialIds.size());
     std::vector<NumericType> redeposition(materialIds.size(), 0.);
     if (params_.redepositionRate > 0.) {
       redeposition = *rates->getScalarData(redepositionLabel);

include/viennaps/models/psSelectiveEpitaxy.hpp

@@ -30,7 +30,7 @@ class EpitaxyVelocityField : public VelocityField<NumericType, D> {
   NumericType getScalarVelocity(const Vec3D<NumericType> &coordinate,
                                 int material, const Vec3D<NumericType> &nv,
                                 unsigned long pointID) override {
-    for (auto epitaxyMaterial : materials) {
+    for (auto const &epitaxyMaterial : materials) {
       if (MaterialMap::isMaterial(material, epitaxyMaterial.first)) {
         double vel = std::max(std::abs(nv[0]), std::abs(nv[D - 1]));
         vel = (vel - low) * factor + R111;

include/viennaps/models/psSingleParticleALD.hpp

@@ -46,6 +46,7 @@ class SingleParticleALDSurfaceModel : public SurfaceModel<NumericType> {
     std::vector<NumericType> depoRate(numPoints, 0.);
 
     auto Coverage = coverages->getScalarData("Coverage");
+    assert(Coverage && Coverage->size() == numPoints);
 
     for (size_t i = 0; i < numPoints; ++i) {
       depoRate[i] = gpc_ * Coverage->at(i);

include/viennaps/models/psSingleParticleProcess.hpp

@@ -29,6 +29,7 @@ class SingleParticleSurfaceModel : public viennaps::SurfaceModel<NumericType> {
     auto velocity =
         SmartPointer<std::vector<NumericType>>::New(materialIds.size(), 0.);
     auto flux = rates->getScalarData("particleFlux");
+    assert(flux && velocity->size() == flux->size());
 
 #pragma omp parallel for
     for (size_t i = 0; i < velocity->size(); i++) {

include/viennaps/models/psTEOSDeposition.hpp

@@ -33,7 +33,7 @@ class SingleTEOSSurfaceModel : public SurfaceModel<NumericType> {
           depositionRate * std::pow(particleFlux->at(i), reactionOrder);
     }
 
-    return SmartPointer<std::vector<NumericType>>::New(velocity);
+    return SmartPointer<std::vector<NumericType>>::New(std::move(velocity));
   }
 
   void updateCoverages(SmartPointer<viennals::PointData<NumericType>> rates,

include/viennaps/models/psTEOSPECVD.hpp

@@ -41,7 +41,7 @@ class PECVDSurfaceModel : public SurfaceModel<NumericType> {
           ionRate_ * std::pow(particleFluxIon->at(i), ionReactionOrder_);
     }
 
-    return SmartPointer<std::vector<NumericType>>::New(velocity);
+    return SmartPointer<std::vector<NumericType>>::New(std::move(velocity));
   }
 };
 

include/viennaps/models/psWetEtching.hpp

@@ -47,7 +47,7 @@ class WetEtchingVelocityField : public VelocityField<NumericType, D> {
   NumericType getScalarVelocity(const Vec3D<NumericType> &coordinate,
                                 int material, const Vec3D<NumericType> &nv,
                                 unsigned long pointID) override {
-    for (auto etchingMaterial : materials) {
+    for (auto const &etchingMaterial : materials) {
       if (MaterialMap::isMaterial(material, etchingMaterial.first)) {
         if (std::abs(Norm(nv) - 1.) > 1e-4)
           return 0.;

include/viennaps/process/psALPStrategy.hpp

@@ -8,7 +8,7 @@
 
 namespace viennaps {
 
-template <typename NumericType, int D>
+VIENNAPS_TEMPLATE_ND(NumericType, D)
 class ALPStrategy final : public ProcessStrategy<NumericType, D> {
   using TranslatorType = std::unordered_map<unsigned long, unsigned long>;
 

include/viennaps/process/psAdvectionCallback.hpp

@@ -6,7 +6,7 @@ namespace viennaps {
 
 using namespace viennacore;
 
-template <typename NumericType, int D> class AdvectionCallback {
+VIENNAPS_TEMPLATE_ND(NumericType, D) class AdvectionCallback {
 protected:
   SmartPointer<Domain<NumericType, D>> domain = nullptr;
 

include/viennaps/process/psAdvectionHandler.hpp

@@ -8,7 +8,7 @@
 
 namespace viennaps {
 
-template <typename NumericType, int D> class AdvectionHandler {
+VIENNAPS_TEMPLATE_ND(NumericType, D) class AdvectionHandler {
   viennals::Advect<NumericType, D> advectionKernel_;
   viennacore::Timer<> timer_;
   unsigned lsVelOutputCounter = 0;
@@ -101,8 +101,7 @@ template <typename NumericType, int D> class AdvectionHandler {
 
     if (context.advectionParams.velocityOutput) {
       auto mesh = viennals::Mesh<NumericType>::New();
-      viennals::ToMesh<NumericType, D>(context.domain->getLevelSets().back(),
-                                       mesh)
+      viennals::ToMesh<NumericType, D>(context.domain->getSurface(), mesh)
           .apply();
       viennals::VTKWriter<NumericType>(
           mesh,
@@ -127,25 +126,31 @@ template <typename NumericType, int D> class AdvectionHandler {
       SmartPointer<std::unordered_map<unsigned long, unsigned long>> const
           &translator) {
     // Move coverages to the top level set
-    auto topLS = context.domain->getLevelSets().back();
+    auto topLS = context.domain->getSurface();
     auto coverages = context.model->getSurfaceModel()->getCoverages();
     assert(coverages != nullptr);
     assert(translator != nullptr);
 
-    for (size_t i = 0; i < coverages->getScalarDataSize(); i++) {
+    std::vector<std::vector<NumericType>> levelSetCoverages(
+        coverages->getScalarDataSize());
+
+#pragma omp parallel for
+    for (unsigned i = 0; i < levelSetCoverages.size(); i++) {
       auto covName = coverages->getScalarDataLabel(i);
       std::vector<NumericType> levelSetData(topLS->getNumberOfPoints(), 0);
       auto cov = coverages->getScalarData(covName);
-      for (const auto iter : *translator.get()) {
-        levelSetData[iter.first] = cov->at(iter.second);
-      }
-      if (auto data = topLS->getPointData().getScalarData(covName, true);
-          data != nullptr) {
-        *data = std::move(levelSetData);
-      } else {
-        topLS->getPointData().insertNextScalarData(std::move(levelSetData),
-                                                   covName);
+
+      for (const auto &[lsId, surfaceId] : *translator) {
+        levelSetData[lsId] = cov->at(surfaceId);
       }
+
+      levelSetCoverages[i] = std::move(levelSetData);
+    }
+
+    for (unsigned i = 0; i < levelSetCoverages.size(); i++) {
+      auto covName = coverages->getScalarDataLabel(i);
+      topLS->getPointData().insertReplaceScalarData(
+          std::move(levelSetCoverages[i]), covName);
     }
 
     return ProcessResult::SUCCESS;
@@ -156,17 +161,22 @@ template <typename NumericType, int D> class AdvectionHandler {
       SmartPointer<std::unordered_map<unsigned long, unsigned long>> const
           &translator) {
     // Update coverages from the advected surface
-    auto topLS = context.domain->getLevelSets().back();
+    auto topLS = context.domain->getSurface();
     auto coverages = context.model->getSurfaceModel()->getCoverages();
+    assert(coverages != nullptr);
+    assert(translator != nullptr);
+
     for (size_t i = 0; i < coverages->getScalarDataSize(); i++) {
       auto covName = coverages->getScalarDataLabel(i);
       auto levelSetData = topLS->getPointData().getScalarData(covName);
       auto covData = coverages->getScalarData(covName);
       covData->resize(translator->size());
-      for (const auto it : *translator.get()) {
-        covData->at(it.second) = levelSetData->at(it.first);
+
+      for (const auto &[lsId, surfaceId] : *translator) {
+        covData->at(surfaceId) = levelSetData->at(lsId);
       }
     }
+
     return ProcessResult::SUCCESS;
   }
   auto &getTimer() const { return timer_; }