Added redeposition model and example (#35)

* First version of redeposition model

* Seperate psMakeStack file

* Added config file and animation

* Added description and check gif

* Added coordinates in calculateVelocities function

* Fixed typos and removed intermediate meshes in visualization mesh

Examples/ExampleProcess/SurfaceModel.hpp

@@ -20,9 +20,10 @@ class SurfaceModel : public psSurfaceModel<NumericType> {
     processParams->insertNextScalar(0., "processParameter");
   }
 
-  psSmartPointer<std::vector<NumericType>>
-  calculateVelocities(psSmartPointer<psPointData<NumericType>> Rates,
-                      const std::vector<NumericType> &materialIds) override {
+  psSmartPointer<std::vector<NumericType>> calculateVelocities(
+      psSmartPointer<psPointData<NumericType>> Rates,
+      const std::vector<std::array<NumericType, 3>> &coordinates,
+      const std::vector<NumericType> &materialIds) override {
     // use coverages and rates here to calculate the velocity here
     return psSmartPointer<std::vector<NumericType>>::New(
         *Rates->getScalarData("particleRate"));

Examples/Redeposition/CMakeLists.txt

@@ -0,0 +1,10 @@
+cmake_minimum_required(VERSION 3.4)
+
+project("Redeposition")
+
+add_executable(${PROJECT_NAME} ${PROJECT_NAME}.cpp)
+target_include_directories(${PROJECT_NAME} PUBLIC ${VIENNAPS_INCLUDE_DIRS})
+target_link_libraries(${PROJECT_NAME} PRIVATE ${VIENNAPS_LIBRARIES})
+configure_file(config.txt ${CMAKE_CURRENT_BINARY_DIR}/config.txt COPYONLY)
+
+add_dependencies(buildExamples ${PROJECT_NAME})

Examples/Redeposition/Parameters.hpp

@@ -0,0 +1,44 @@
+#pragma once
+
+#include <string>
+#include <unordered_map>
+
+#include <psUtils.hpp>
+
+template <typename T> struct Parameters {
+  // Domain
+  T gridDelta = 0.2;
+  T xExtent = 10.0;
+
+  // Geometry
+  int numLayers = 26;
+  T layerHeight = 2;
+  T substrateHeight = 4;
+  T holeRadius = 2;
+
+  // Process
+  T diffusionCoefficient = 10.; // diffusion cofficient
+  T sink = 0.001;               // sink strength
+  // convection velocity in the scallops towards the center
+  T scallopStreamVelocity = 10.;
+  // convection velocity in the center towards the sink on the top
+  T holeStreamVelocity = 10.;
+
+  Parameters() {}
+
+  void fromMap(std::unordered_map<std::string, std::string> &m) {
+    psUtils::AssignItems(                                              //
+        m,                                                             //
+        psUtils::Item{"gridDelta", gridDelta},                         //
+        psUtils::Item{"xExtent", xExtent},                             //
+        psUtils::Item{"numLayers", numLayers},                         //
+        psUtils::Item{"layerHeight", layerHeight},                     //
+        psUtils::Item{"substrateHeight", substrateHeight},             //
+        psUtils::Item{"holeRadius", holeRadius},                       //
+        psUtils::Item{"diffusionCoefficient", diffusionCoefficient},   //
+        psUtils::Item{"sink", sink},                                   //
+        psUtils::Item{"scallopStreamVelocity", scallopStreamVelocity}, //
+        psUtils::Item{"holeStreamVelocity", holeStreamVelocity}        //
+    );
+  }
+};

Examples/Redeposition/Redeposition.cpp

@@ -0,0 +1,97 @@
+#include <csDenseCellSet.hpp>
+#include <csTracing.hpp>
+
+#include <psMakeStack.hpp>
+#include <psProcess.hpp>
+#include <psProcessModel.hpp>
+#include <psWriteVisualizationMesh.hpp>
+
+#include <StackRedeposition.hpp>
+
+#include "Parameters.hpp"
+
+int main(int argc, char **argv) {
+  using NumericType = double;
+
+  omp_set_num_threads(12);
+  // Attention:
+  // This model/example currently only works in 2D mode
+  constexpr int D = 2;
+
+  Parameters<NumericType> params;
+  if (argc > 1) {
+    auto config = psUtils::readConfigFile(argv[1]);
+    if (config.empty()) {
+      std::cerr << "Empty config provided" << std::endl;
+      return -1;
+    }
+    params.fromMap(config);
+  }
+
+  auto domain = psSmartPointer<psDomain<NumericType, D>>::New();
+  psMakeStack<NumericType, D> builder(
+      domain, params.gridDelta, params.xExtent, 0., params.numLayers,
+      params.layerHeight, params.substrateHeight, params.holeRadius, false);
+  builder.apply();
+  // copy top layer for deposition
+  auto depoLayer = psSmartPointer<lsDomain<NumericType, D>>::New(
+      domain->getLevelSets()->back());
+  domain->insertNextLevelSet(depoLayer);
+  domain->generateCellSet(2., true /* true means cell set above surface */);
+  auto &cellSet = domain->getCellSet();
+  cellSet->writeVTU("initial.vtu");
+  // we need neighborhood information for solving the
+  // convection-diffusion equation on the cell set
+  cellSet->buildNeighborhood();
+
+  // The redeposition model captures byproducts from the selective etching
+  // process in the cell set. The byproducts are then distributed by solving a
+  // convection-diffusion equation on the cell set.
+  auto redepoModel = psSmartPointer<RedepositionDynamics<NumericType, D>>::New(
+      domain, params.diffusionCoefficient, params.sink,
+      params.scallopStreamVelocity, params.holeStreamVelocity,
+      builder.getTopLayer(), builder.getHeight(), builder.getHoleRadius());
+  auto velField = psSmartPointer<psDefaultVelocityField<NumericType>>::New();
+  auto etchSurfModel =
+      psSmartPointer<SelectiveEtchingSurfaceModel<NumericType>>::New();
+  auto depoSurfModel =
+      psSmartPointer<RedepositionSurfaceModel<NumericType, D>>::New(
+          cellSet, 1., builder.getHeight(), builder.getTopLayer());
+
+  // run the selective etching process
+  {
+    auto etchModel = psSmartPointer<psProcessModel<NumericType, D>>::New();
+    etchModel->setSurfaceModel(etchSurfModel);
+    etchModel->setVelocityField(velField);
+    etchModel->setAdvectionCallback(redepoModel);
+    etchModel->setProcessName("SelectiveEtching");
+
+    psProcess<NumericType, D> processEtch;
+    processEtch.setDomain(domain);
+    processEtch.setProcessModel(etchModel);
+    processEtch.setProcessDuration(50.);
+
+    processEtch.apply();
+  }
+
+  // run the redeposition based on the captured byproducts
+  {
+    auto depoModel = psSmartPointer<psProcessModel<NumericType, D>>::New();
+    depoModel->setSurfaceModel(depoSurfModel);
+    depoModel->setVelocityField(velField);
+    depoModel->setProcessName("Redeposition");
+
+    psProcess<NumericType, D> processRedepo;
+    processRedepo.setDomain(domain);
+    processRedepo.setProcessModel(depoModel);
+    processRedepo.setProcessDuration(0.2);
+
+    processRedepo.apply();
+    cellSet->updateMaterials();
+    cellSet->writeVTU("RedepositionCellSet.vtu");
+  }
+
+  psWriteVisualizationMesh<NumericType, D>(domain, "FinalStack").apply();
+
+  return 0;
+}

Examples/Redeposition/config.txt

@@ -0,0 +1,16 @@
+# Config file for a redeposition example process
+# Domain
+gridDelta=0.1
+xExtent=10.0
+
+# Geometry
+int numLayers=26
+layerHeight=2
+substrateHeight=4
+holeRadius=2
+
+# Process
+diffusionCoefficient=10. 
+sink=0.001               
+scallopStreamVelocity=2.
+holeStreamVelocity=2.

README.md

@@ -97,6 +97,13 @@ By changing the dimension of the hole etching example (_D = 2_), we can easily s
   <img src="https://raw.githubusercontent.com/ViennaTools/ViennaPS/master/data/images/sidewall_tapering.svg" width=700 style="background-color:white;">
 </div>
 
+### Redeposition During Selective Etching
+
+This example demonstrates capturing etching byproducts and the subsequent redeposition during a selective etching process in a Si<sub>3</sub>N<sub>4</sub>/SiO<sub>2</sub> stack. The etching byproducts are captured in a cell set description of the etching plasma. To model the dynamics of these etching byproducts, a convection-diffusion equation is solved on the cell set using finite differences. The redeposition is then captured by adding up the byproducts in every step and using this information to generate a velocity field on the etched surface. 
+<div align="center">
+  <img src="https://raw.githubusercontent.com/ViennaTools/ViennaPS/master/data/images/redeposition.gif" width=700 style="background-color:white;">
+</div>
+
 ## Application
 
 It is also possible to build an application which can parse a custom configuration file and execute pre-defined processes. The application can be built using CMake (make sure all dependencies are installed/ have been built previously):

include/CellSet/csDenseCellSet.hpp

@@ -14,20 +14,23 @@
 
 #include <rayUtil.hpp>
 
+#include <psSmartPointer.hpp>
+
 /**
   This class represents a cell-based voxel implementation of a volume. The
   depth of the cell set in z-direction can be specified.
 */
 template <class T, int D> class csDenseCellSet {
 private:
-  using gridType = lsSmartPointer<lsMesh<T>>;
+  using gridType = psSmartPointer<lsMesh<T>>;
   using levelSetsType =
-      lsSmartPointer<std::vector<lsSmartPointer<lsDomain<T, D>>>>;
+      psSmartPointer<std::vector<psSmartPointer<lsDomain<T, D>>>>;
 
   levelSetsType levelSets = nullptr;
   gridType cellGrid = nullptr;
-  lsSmartPointer<lsDomain<T, D>> surface = nullptr;
-  lsSmartPointer<csBVH<T, D>> BVH = nullptr;
+  psSmartPointer<lsDomain<T, D>> surface = nullptr;
+  psSmartPointer<csBVH<T, D>> BVH = nullptr;
+  std::vector<std::set<unsigned long>> neighborhood;
   T gridDelta;
   size_t numberOfCells;
   T depth = 0.;
@@ -49,10 +52,10 @@ template <class T, int D> class csDenseCellSet {
     levelSets = passedLevelSets;
 
     if (cellGrid == nullptr)
-      cellGrid = lsSmartPointer<lsMesh<T>>::New();
+      cellGrid = psSmartPointer<lsMesh<T>>::New();
 
     if (surface == nullptr)
-      surface = lsSmartPointer<lsDomain<T, D>>::New(levelSets->back());
+      surface = psSmartPointer<lsDomain<T, D>>::New(levelSets->back());
     else
       surface->deepCopy(levelSets->back());
 
@@ -83,14 +86,14 @@ template <class T, int D> class csDenseCellSet {
     }
 
     lsToVoxelMesh<T, D> voxelConverter(cellGrid);
-    auto plane = lsSmartPointer<lsDomain<T, D>>::New(surface->getGrid());
+    auto plane = psSmartPointer<lsDomain<T, D>>::New(surface->getGrid());
     if (depth > 0.) {
       T origin[D] = {0.};
       T normal[D] = {0.};
       origin[D - 1] = depthPlanePos;
       normal[D - 1] = 1.;
       lsMakeGeometry<T, D>(plane,
-                           lsSmartPointer<lsPlane<T, D>>::New(origin, normal))
+                           psSmartPointer<lsPlane<T, D>>::New(origin, normal))
           .apply();
     }
     if (!cellSetAboveSurface && depth > 0.)
@@ -113,7 +116,7 @@ template <class T, int D> class csDenseCellSet {
     fillingFractions = cellGrid->getCellData().getScalarData("fillingFraction");
 
     calculateBounds(minBounds, maxBounds);
-    BVH = lsSmartPointer<csBVH<T, D>>::New(getBoundingBox(), BVHlayers);
+    BVH = psSmartPointer<csBVH<T, D>>::New(getBoundingBox(), BVHlayers);
     buildBVH();
   }
 
@@ -134,7 +137,7 @@ template <class T, int D> class csDenseCellSet {
 
   gridType getCellGrid() { return cellGrid; }
 
-  lsSmartPointer<csBVH<T, D>> getBVH() const { return BVH; }
+  psSmartPointer<csBVH<T, D>> getBVH() const { return BVH; }
 
   T getDepth() const { return depth; }
 
@@ -144,11 +147,11 @@ template <class T, int D> class csDenseCellSet {
     return cellGrid->getNodes();
   }
 
-  std::vector<std::array<unsigned, (1 << D)>> getElements() {
+  std::vector<std::array<unsigned, (1 << D)>> &getElements() {
     return cellGrid->template getElements<(1 << D)>();
   }
 
-  lsSmartPointer<lsDomain<T, D>> getSurface() { return surface; }
+  psSmartPointer<lsDomain<T, D>> getSurface() { return surface; }
 
   levelSetsType getLevelSets() const { return levelSets; }
 
@@ -250,14 +253,14 @@ template <class T, int D> class csDenseCellSet {
 
     lsToVoxelMesh<T, D> voxelConverter(cellGrid);
     auto plane =
-        lsSmartPointer<lsDomain<T, D>>::New(levelSets->back()->getGrid());
+        psSmartPointer<lsDomain<T, D>>::New(levelSets->back()->getGrid());
     if (depth > 0.) {
       T origin[D] = {0.};
       T normal[D] = {0.};
       origin[D - 1] = depthPlanePos;
       normal[D - 1] = 1.;
       lsMakeGeometry<T, D>(plane,
-                           lsSmartPointer<lsPlane<T, D>>::New(origin, normal))
+                           psSmartPointer<lsPlane<T, D>>::New(origin, normal))
           .apply();
     }
     if (!cellSetAboveSurface && depth > 0.)
@@ -287,18 +290,18 @@ template <class T, int D> class csDenseCellSet {
   // Updates the surface of the cell set. The new surface should be below the
   // old surface as this function can only remove cells from the cell set.
   void updateSurface() {
-    auto updateCellGrid = lsSmartPointer<lsMesh<T>>::New();
+    auto updateCellGrid = psSmartPointer<lsMesh<T>>::New();
 
     lsToVoxelMesh<T, D> voxelConverter(updateCellGrid);
     if (depth != 0.) {
-      auto plane = lsSmartPointer<lsDomain<T, D>>::New(surface->getGrid());
+      auto plane = psSmartPointer<lsDomain<T, D>>::New(surface->getGrid());
       T origin[D] = {0.};
       T normal[D] = {0.};
       origin[D - 1] = depthPlanePos;
       normal[D - 1] = 1.;
 
       lsMakeGeometry<T, D>(plane,
-                           lsSmartPointer<lsPlane<T, D>>::New(origin, normal))
+                           psSmartPointer<lsPlane<T, D>>::New(origin, normal))
           .apply();
       voxelConverter.insertNextLevelSet(plane);
     }
@@ -346,6 +349,42 @@ template <class T, int D> class csDenseCellSet {
     }
   }
 
+  void buildNeighborhood() {
+    const auto &cells = cellGrid->template getElements<(1 << D)>();
+    unsigned const numNodes = cellGrid->getNodes().size();
+    unsigned const numCells = cells.size();
+
+    neighborhood.clear();
+    neighborhood.resize(numCells);
+
+    std::vector<std::vector<unsigned>> nodeCellConnections(numNodes);
+
+    // for each node, store which cells are connected with the node
+    for (unsigned cellIdx = 0; cellIdx < numCells; cellIdx++) {
+      for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) {
+        nodeCellConnections[cells[cellIdx][cellNodeIdx]].push_back(cellIdx);
+      }
+    }
+
+    for (unsigned cellIdx = 0; cellIdx < numCells; cellIdx++) {
+      for (unsigned cellNodeIdx = 0; cellNodeIdx < (1 << D); cellNodeIdx++) {
+        auto &cellsAtNode = nodeCellConnections[cells[cellIdx][cellNodeIdx]];
+        for (const auto &neighborCell : cellsAtNode) {
+          if (neighborCell != cellIdx) {
+            neighborhood[cellIdx].insert(neighborCell);
+          }
+        }
+      }
+    }
+  }
+
+  std::set<size_t> &getNeighbors(size_t cellIdx) {
+    assert(!neighborhood.empty() &&
+           "Querying neighbors without creating neighborhood structure");
+    assert(cellIdx < numberOfCells && "Cell idx out of bounds");
+    return neighborhood[cellIdx];
+  }
+
 private:
   int findIndex(const csTriple<T> &point) {
     const auto &elems = cellGrid->template getElements<(1 << D)>();
@@ -453,4 +492,4 @@ template <class T, int D> class csDenseCellSet {
   }
 };
 
-#endif
+#endif

include/Geometries/psMakeFin.hpp

@@ -5,6 +5,9 @@
 #include <lsMakeGeometry.hpp>
 #include <psDomain.hpp>
 
+/**
+ * Creates a fin gemeotry in in z(3D)/y(2D) direction.
+ */
 template <class NumericType, int D> class psMakeFin {
   using LSPtrType = psSmartPointer<lsDomain<NumericType, D>>;
   using PSPtrType = psSmartPointer<psDomain<NumericType, D>>;