Add 2-particle TEOS PE-CVD model (#87)

* Add 2-particle PE-CVD model

* Add option for material specific rates in psSingleParticleProcess

* Correct order in psSingleParticleProcess constructor

.gitignore

@@ -1,5 +1,4 @@
 build/
-venv/
 
 .vscode
 .vs

include/viennaps/models/psSingleParticleProcess.hpp

@@ -13,13 +13,14 @@ using namespace viennacore;
 namespace impl {
 template <typename NumericType, int D>
 class SingleParticleSurfaceModel : public viennaps::SurfaceModel<NumericType> {
-  const NumericType rateFactor_;
-  const std::vector<Material> maskMaterials_;
+  const NumericType rate_;
+  const std::unordered_map<Material, NumericType> materialRates_;
 
 public:
-  SingleParticleSurfaceModel(NumericType rate,
-                             const std::vector<Material> &mask)
-      : rateFactor_(rate), maskMaterials_(mask) {}
+  SingleParticleSurfaceModel(
+      const NumericType rate,
+      const std::unordered_map<Material, NumericType> &mask)
+      : rate_(rate), materialRates_(mask) {}
 
   SmartPointer<std::vector<NumericType>> calculateVelocities(
       SmartPointer<viennals::PointData<NumericType>> rates,
@@ -31,22 +32,17 @@ class SingleParticleSurfaceModel : public viennaps::SurfaceModel<NumericType> {
     auto flux = rates->getScalarData("particleFlux");
 
     for (std::size_t i = 0; i < velocity->size(); i++) {
-      if (!isMaskMaterial(materialIds[i])) {
-        velocity->at(i) = flux->at(i) * rateFactor_;
+      if (auto matRate =
+              materialRates_.find(MaterialMap::mapToMaterial(materialIds[i]));
+          matRate == materialRates_.end()) {
+        velocity->at(i) = flux->at(i) * rate_;
+      } else {
+        velocity->at(i) = flux->at(i) * matRate->second;
       }
     }
 
     return velocity;
   }
-
-private:
-  bool isMaskMaterial(const NumericType &material) const {
-    for (const auto &mat : maskMaterials_) {
-      if (MaterialMap::isMaterial(material, mat))
-        return true;
-    }
-    return false;
-  }
 };
 
 template <typename NumericType, int D>
@@ -97,22 +93,34 @@ class SingleParticleProcess : public ProcessModel<NumericType, D> {
                         NumericType stickingProbability = 1.,
                         NumericType sourceDistributionPower = 1.,
                         Material maskMaterial = Material::None) {
-    std::vector<Material> maskMaterialVec = {maskMaterial};
+    std::unordered_map<Material, NumericType> maskMaterialMap = {
+        {maskMaterial, 0.}};
     initialize(rate, stickingProbability, sourceDistributionPower,
-               std::move(maskMaterialVec));
+               std::move(maskMaterialMap));
   }
 
   SingleParticleProcess(NumericType rate, NumericType stickingProbability,
                         NumericType sourceDistributionPower,
                         std::vector<Material> maskMaterial) {
+    std::unordered_map<Material, NumericType> maskMaterialMap;
+    for (auto &mat : maskMaterial) {
+      maskMaterialMap[mat] = 0.;
+    }
     initialize(rate, stickingProbability, sourceDistributionPower,
-               std::move(maskMaterial));
+               std::move(maskMaterialMap));
+  }
+
+  SingleParticleProcess(std::unordered_map<Material, NumericType> materialRates,
+                        NumericType stickingProbability,
+                        NumericType sourceDistributionPower) {
+    initialize(0., stickingProbability, sourceDistributionPower,
+               std::move(materialRates));
   }
 
 private:
   void initialize(NumericType rate, NumericType stickingProbability,
                   NumericType sourceDistributionPower,
-                  std::vector<Material> &&maskMaterial) {
+                  std::unordered_map<Material, NumericType> &&maskMaterial) {
     // particles
     auto particle = std::make_unique<impl::SingleParticle<NumericType, D>>(
         stickingProbability, sourceDistributionPower);

include/viennaps/models/psTEOSPECVD.hpp

@@ -0,0 +1,170 @@
+#pragma once
+
+#include "../psProcessModel.hpp"
+
+#include <rayParticle.hpp>
+#include <rayReflection.hpp>
+#include <rayUtil.hpp>
+
+namespace viennaps {
+
+namespace impl {
+template <class NumericType>
+class PECVDSurfaceModel : public SurfaceModel<NumericType> {
+  const NumericType radicalRate;
+  const NumericType radicalReactionOrder;
+  const NumericType ionRate;
+  const NumericType ionReactionOrder;
+
+public:
+  PECVDSurfaceModel(const NumericType passedRadicalRate,
+                    const NumericType passedRadicalReactionOrder,
+                    const NumericType passedIonRate,
+                    const NumericType passedIonReactionOrder)
+      : radicalRate(passedRadicalRate),
+        radicalReactionOrder(passedRadicalReactionOrder),
+        ionRate(passedIonRate), ionReactionOrder(passedIonReactionOrder) {}
+
+  SmartPointer<std::vector<NumericType>> calculateVelocities(
+      SmartPointer<viennals::PointData<NumericType>> rates,
+      const std::vector<std::array<NumericType, 3>> &coordinates,
+      const std::vector<NumericType> &materialIDs) override {
+    // define the surface reaction here
+    auto particleFluxRadical = rates->getScalarData("radicalFlux");
+    auto particleFluxIon = rates->getScalarData("ionFlux");
+
+    std::vector<NumericType> velocity(particleFluxRadical->size(), 0.);
+
+    for (std::size_t i = 0; i < velocity.size(); i++) {
+      // calculate surface velocity based on particle fluxes
+      velocity[i] =
+          radicalRate *
+              std::pow(particleFluxRadical->at(i), radicalReactionOrder) +
+          ionRate * std::pow(particleFluxIon->at(i), ionReactionOrder);
+    }
+
+    return SmartPointer<std::vector<NumericType>>::New(velocity);
+  }
+};
+
+// Particle type (modify at you own risk)
+template <class NumericType, int D>
+class Radical
+    : public viennaray::Particle<Radical<NumericType, D>, NumericType> {
+public:
+  Radical(const NumericType pStickingProbability,
+          const std::string pDataLabel = "radicalFlux")
+      : stickingProbability(pStickingProbability), dataLabel(pDataLabel) {}
+  std::pair<NumericType, Vec3D<NumericType>>
+  surfaceReflection(NumericType rayWeight, const Vec3D<NumericType> &rayDir,
+                    const Vec3D<NumericType> &geomNormal,
+                    const unsigned int primID, const int materialId,
+                    const viennaray::TracingData<NumericType> *globalData,
+                    RNG &Rng) override final {
+    auto direction =
+        viennaray::ReflectionDiffuse<NumericType, D>(geomNormal, Rng);
+    return std::pair<NumericType, Vec3D<NumericType>>{stickingProbability,
+                                                      direction};
+  }
+  void surfaceCollision(NumericType rayWeight, const Vec3D<NumericType> &rayDir,
+                        const Vec3D<NumericType> &geomNormal,
+                        const unsigned int primID, const int materialId,
+                        viennaray::TracingData<NumericType> &localData,
+                        const viennaray::TracingData<NumericType> *globalData,
+                        RNG &Rng) override final {
+    localData.getVectorData(0)[primID] += rayWeight;
+  }
+  NumericType getSourceDistributionPower() const override final { return 1; }
+  std::vector<std::string> getLocalDataLabels() const override final {
+    return {dataLabel};
+  }
+
+private:
+  const NumericType stickingProbability;
+  const std::string dataLabel = "radicalFlux";
+};
+
+template <typename NumericType, int D>
+class Ion : public viennaray::Particle<Ion<NumericType, D>, NumericType> {
+public:
+  Ion(const NumericType pStickingProbability, const NumericType pExponent,
+      const NumericType pMinAngle, const std::string pDataLabel = "ionFlux")
+      : stickingProbability(pStickingProbability), exponent(pExponent),
+        minAngle(pMinAngle), dataLabel(pDataLabel) {}
+
+  std::pair<NumericType, Vec3D<NumericType>>
+  surfaceReflection(NumericType rayWeight, const Vec3D<NumericType> &rayDir,
+                    const Vec3D<NumericType> &geomNormal,
+                    const unsigned int primId, const int materialId,
+                    const viennaray::TracingData<NumericType> *globalData,
+                    RNG &Rng) override final {
+    auto cosTheta = -rayInternal::DotProduct(rayDir, geomNormal);
+
+    assert(cosTheta >= 0 && "Hit backside of disc");
+    assert(cosTheta <= 1 + 1e-6 && "Error in calculating cos theta");
+
+    NumericType incAngle =
+        std::acos(std::max(std::min(cosTheta, static_cast<NumericType>(1.)),
+                           static_cast<NumericType>(0.)));
+
+    auto direction = viennaray::ReflectionConedCosine<NumericType, D>(
+        rayDir, geomNormal, Rng, std::max(incAngle, minAngle));
+    return std::pair<NumericType, Vec3D<NumericType>>{stickingProbability,
+                                                      direction};
+  }
+
+  void surfaceCollision(NumericType rayWeight, const Vec3D<NumericType> &rayDir,
+                        const Vec3D<NumericType> &geomNormal,
+                        const unsigned int primID, const int materialId,
+                        viennaray::TracingData<NumericType> &localData,
+                        const viennaray::TracingData<NumericType> *globalData,
+                        RNG &Rng) override final {
+    localData.getVectorData(0)[primID] += rayWeight;
+  }
+
+  NumericType getSourceDistributionPower() const override final {
+    return exponent;
+  }
+  std::vector<std::string> getLocalDataLabels() const override final {
+    return {dataLabel};
+  }
+
+private:
+  const NumericType stickingProbability;
+  const NumericType exponent;
+  const NumericType minAngle;
+  const std::string dataLabel = "ionFlux";
+};
+} // namespace impl
+
+template <class NumericType, int D>
+class TEOSPECVD : public ProcessModel<NumericType, D> {
+public:
+  TEOSPECVD(const NumericType pRadicalSticking, const NumericType pRadicalRate,
+            const NumericType pIonRate, const NumericType pIonExponent,
+            const NumericType pIonSticking = 1.,
+            const NumericType pRadicalOrder = 1.,
+            const NumericType pIonOrder = 1.,
+            const NumericType pIonMinAngle = 0.) {
+    // velocity field
+    auto velField = SmartPointer<DefaultVelocityField<NumericType>>::New(2);
+    this->setVelocityField(velField);
+
+    // particles
+    auto radical = std::make_unique<impl::Radical<NumericType, D>>(
+        pRadicalSticking, "radicalFlux");
+    auto ion = std::make_unique<impl::Ion<NumericType, D>>(
+        pIonSticking, pIonExponent, pIonMinAngle, "ionFlux");
+
+    // surface model
+    auto surfModel = SmartPointer<impl::PECVDSurfaceModel<NumericType>>::New(
+        pRadicalRate, pRadicalOrder, pIonRate, pIonOrder);
+
+    this->setSurfaceModel(surfModel);
+    this->insertNextParticleType(radical);
+    this->insertNextParticleType(ion);
+    this->setProcessName("TEOSPECVD");
+  }
+};
+
+} // namespace viennaps

python/pyWrap.cpp

@@ -57,6 +57,7 @@
 #include <models/psSingleParticleALD.hpp>
 #include <models/psSingleParticleProcess.hpp>
 #include <models/psTEOSDeposition.hpp>
+#include <models/psTEOSPECVD.hpp>
 
 // visualization
 #include <psToDiskMesh.hpp>
@@ -574,7 +575,10 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
                  rate, sticking, power, mask);
            }),
            pybind11::arg("rate"), pybind11::arg("stickingProbability"),
-           pybind11::arg("sourceExponent"), pybind11::arg("maskMaterials"));
+           pybind11::arg("sourceExponent"), pybind11::arg("maskMaterials"))
+      .def(pybind11::init<std::unordered_map<Material, T>, T, T>(),
+           pybind11::arg("materialRates"), pybind11::arg("stickingProbability"),
+           pybind11::arg("sourceExponent"));
 
   // TEOS Deposition
   pybind11::class_<TEOSDeposition<T, D>, SmartPointer<TEOSDeposition<T, D>>>(
@@ -588,6 +592,22 @@ PYBIND11_MODULE(VIENNAPS_MODULE_NAME, module) {
            pybind11::arg("stickingProbabilityP2") = 0.,
            pybind11::arg("rateP2") = 0., pybind11::arg("orderP2") = 0.);
 
+  // TEOS PE-CVD
+  pybind11::class_<TEOSPECVD<T, D>, SmartPointer<TEOSPECVD<T, D>>>(
+      module, "TEOSPECVD", processModel)
+      .def(
+          pybind11::init(&SmartPointer<TEOSPECVD<T, D>>::New<
+                         const T /*stR*/, const T /*rateR*/, const T /*orderR*/,
+                         const T /*stI*/, const T /*rateI*/, const T /*orderI*/,
+                         const T /*exponentI*/, const T /*minAngleIon*/>),
+          pybind11::arg("stickingProbabilityRadical"),
+          pybind11::arg("depositionRateRadical"),
+          pybind11::arg("depositionRateIon"), pybind11::arg("exponentIon"),
+          pybind11::arg("stickingProbabilityIon") = 1.,
+          pybind11::arg("reactionOrderRadical") = 1.,
+          pybind11::arg("reactionOrderIon") = 1.,
+          pybind11::arg("minAngleIon") = 0.);
+
   // SF6O2 Parameters
   pybind11::class_<SF6O2Parameters<T>::MaskType>(module, "SF6O2ParametersMask")
       .def(pybind11::init<>())