Implement 2D triangle engines and transmission benchmark (#178)

* Implement 2D triangle tracing on CPU

* Working engine test for CPU engines

* Engine test working for all engines

* Line engine smoothing

* 2D benchmark

* Assertions

* Implement ALD GPU model

* Add transmission benchmark

* Transmission test on different surfaces

* Fix CPU disk transmission test

* Bump ViennaRay

* Format and add missing ignoreFluxBoundary setters

Author: tobre1

CMakeLists.txt

@@ -113,7 +113,7 @@ CPMAddPackage(
 
 CPMFindPackage(
   NAME ViennaRay
-  VERSION 3.7.2
+  VERSION 3.8.0
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaRay"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}
   OPTIONS "VIENNARAY_USE_GPU ${VIENNAPS_USE_GPU}")

cmake/generate_ptx.cmake

@@ -84,7 +84,6 @@ function(add_GPU_executable target_name_base target_name_var)
   cuda_get_sources_and_options(source_files cmake_options options ${ARGN})
 
   # Isolate OBJ target files. NVCC should only process these files and leave PTX targets for NVRTC
-  set(cu_obj_source_files)
   set(cu_optix_source_files)
   foreach(file ${source_files})
     get_filename_component(_file_extension ${file} EXT)

examples/atomicLayerDeposition/atomicLayerDeposition.py

@@ -1,6 +1,8 @@
 import viennaps as ps
 import viennals as ls
 
+ps.setDimension(2)
+ls.setDimension(2)
 
 params = ps.readConfigFile("config.txt")
 geometry = ps.Domain()
@@ -50,7 +52,7 @@
 )
 print("Mean free path: ", gasMFP, " um")
 
-model = ps.d2.SingleParticleALD(
+model = ps.SingleParticleALD(
     stickingProbability=params["stickingProbability"],
     numCycles=int(params["numCycles"]),
     growthPerCycle=params["growthPerCycle"],
@@ -70,7 +72,7 @@
 rayParams = ps.RayTracingParameters()
 rayParams.raysPerPoint = int(params["numRaysPerPoint"])
 
-ALP = ps.d2.Process(geometry, model)
+ALP = ps.Process(geometry, model)
 ALP.setParameters(rayParams)
 ALP.setParameters(alpParams)
 ALP.apply()

examples/ionBeamEtching/config.txt

@@ -20,4 +20,4 @@ redepositionRate=0.05
 planeWaferRate=1.0
 angle=35.0
 
-processTime=5
+processTime=2

gpu/benchmark/Benchmark.hpp

@@ -6,25 +6,26 @@
 #include <psDomain.hpp>
 #include <rayParticle.hpp>
 
-#define MAKE_GEO Hole
+#define MAKE_GEO Trench
 #define DEFAULT_GRID_DELTA 0.1
 #define DEFAULT_STICKING 0.1
-#define DIM 3
+#define DIM 2
+#define FIXED_RAYS false
 
-constexpr int particleType = 0;
+constexpr int particleType = 1;
 using TranslatorType = std::unordered_map<unsigned long, unsigned long>;
+using namespace viennaps;
 
 template <class NumericType>
 auto Trench(NumericType gridDelta = DEFAULT_GRID_DELTA) {
-  NumericType xExtent = 10.;
-  NumericType yExtent = 10.;
-  NumericType width = 5.;
-  NumericType depth = 25.;
+  NumericType xExtent = 20.;
+  NumericType yExtent = 20.;
+  NumericType width = 10.;
+  NumericType depth = 50.;
 
-  using namespace viennaps;
   auto domain = Domain<NumericType, DIM>::New(
       gridDelta, xExtent, yExtent, BoundaryType::REFLECTIVE_BOUNDARY);
-  MakeTrench<NumericType, DIM>(domain, width, depth).apply().apply();
+  MakeTrench<NumericType, DIM>(domain, width, depth).apply();
   return domain;
 }
 
@@ -35,7 +36,6 @@ auto Hole(NumericType gridDelta = DEFAULT_GRID_DELTA) {
   NumericType radius = 3.0;
   NumericType depth = 30.;
 
-  using namespace viennaps;
   auto domain = Domain<NumericType, DIM>::New(
       gridDelta, xExtent, yExtent, BoundaryType::REFLECTIVE_BOUNDARY);
   MakeHole<NumericType, DIM>(domain, radius, depth).apply();
@@ -159,4 +159,189 @@ auto getDeviceParams() {
   deviceParams.redepositionThreshold = params_.redepositionThreshold;
 
   return deviceParams;
+}
+
+template <class NumericType, int D, class TracerType>
+void setupTriangleGeometry(
+    SmartPointer<Domain<NumericType, D>> &domain,
+    SmartPointer<viennals::Mesh<float>> &surfaceMesh_,
+    SmartPointer<KDTree<NumericType, Vec3D<NumericType>>> &elementKdTree_,
+    TracerType &rayTracer_) {
+  CreateSurfaceMesh<NumericType, float, D>(
+      domain->getLevelSets().back(), surfaceMesh_, elementKdTree_, 1e-12, 0.05)
+      .apply();
+
+  viennaray::TriangleMesh triangleMesh;
+
+  if constexpr (D == 2) {
+    viennaray::LineMesh lineMesh;
+    lineMesh.gridDelta = static_cast<float>(domain->getGridDelta());
+    lineMesh.lines = surfaceMesh_->lines;
+    lineMesh.nodes = surfaceMesh_->nodes;
+    lineMesh.minimumExtent = surfaceMesh_->minimumExtent;
+    lineMesh.maximumExtent = surfaceMesh_->maximumExtent;
+
+    triangleMesh = convertLinesToTriangles(lineMesh);
+    assert(triangleMesh.triangles.size() > 0);
+
+    std::vector<Vec3D<NumericType>> triangleCenters;
+    triangleCenters.reserve(triangleMesh.triangles.size());
+    for (const auto &tri : triangleMesh.triangles) {
+      Vec3D<NumericType> center = {0, 0, 0};
+      for (int i = 0; i < 3; ++i) {
+        center[0] += triangleMesh.nodes[tri[i]][0];
+        center[1] += triangleMesh.nodes[tri[i]][1];
+        center[2] += triangleMesh.nodes[tri[i]][2];
+      }
+      triangleCenters.push_back(center / static_cast<NumericType>(3.0));
+    }
+    assert(triangleCenters.size() > 0);
+    elementKdTree_->setPoints(triangleCenters);
+    elementKdTree_->build();
+  } else {
+    triangleMesh = CreateTriangleMesh(
+        static_cast<float>(domain->getGridDelta()), surfaceMesh_);
+  }
+
+  rayTracer_.setGeometry(triangleMesh);
+
+  if constexpr (D == 2) {
+    surfaceMesh_->nodes = std::move(triangleMesh.nodes);
+    surfaceMesh_->triangles = std::move(triangleMesh.triangles);
+    surfaceMesh_->getCellData().insertReplaceVectorData(
+        std::move(triangleMesh.normals), "Normals");
+    surfaceMesh_->minimumExtent = triangleMesh.minimumExtent;
+    surfaceMesh_->maximumExtent = triangleMesh.maximumExtent;
+  }
+}
+
+template <class NumericType, int D, class TracerType>
+void setupLineGeometry(
+    SmartPointer<Domain<NumericType, D>> &domain,
+    SmartPointer<viennals::Mesh<float>> &surfaceMesh_,
+    SmartPointer<KDTree<NumericType, Vec3D<NumericType>>> &elementKdTree_,
+    TracerType &rayTracer_) {
+  CreateSurfaceMesh<NumericType, float, D>(
+      domain->getLevelSets().back(), surfaceMesh_, elementKdTree_, 1e-12, 0.05)
+      .apply();
+
+  viennaray::LineMesh lineMesh(surfaceMesh_->nodes, surfaceMesh_->lines,
+                               static_cast<float>(domain->getGridDelta()));
+  // lines might have changed, so we need to update the surfaceMesh_ later
+
+  std::vector<Vec3D<NumericType>> elementCenters(lineMesh.lines.size());
+  for (int i = 0; i < lineMesh.lines.size(); ++i) {
+    auto const &p0 = lineMesh.nodes[lineMesh.lines[i][0]];
+    auto const &p1 = lineMesh.nodes[lineMesh.lines[i][1]];
+    auto center = (p0 + p1) / 2.f;
+    elementCenters[i] = Vec3D<NumericType>{static_cast<NumericType>(center[0]),
+                                           static_cast<NumericType>(center[1]),
+                                           static_cast<NumericType>(center[2])};
+  }
+  elementKdTree_->setPoints(elementCenters);
+  elementKdTree_->build();
+
+  rayTracer_.setGeometry(lineMesh);
+
+  surfaceMesh_->nodes = std::move(lineMesh.nodes);
+  surfaceMesh_->lines = std::move(lineMesh.lines);
+  surfaceMesh_->getCellData().insertReplaceVectorData(
+      std::move(lineMesh.normals), "Normals");
+  surfaceMesh_->minimumExtent = lineMesh.minimumExtent;
+  surfaceMesh_->maximumExtent = lineMesh.maximumExtent;
+}
+
+template <typename NumericType, class TracerType>
+void postProcessLineData(
+    viennals::PointData<NumericType> &pointData,
+    SmartPointer<viennals::Mesh<NumericType>> diskMesh, int smoothingNeighbors,
+    NumericType gridDelta, TracerType &rayTracer_,
+    SmartPointer<KDTree<NumericType, Vec3D<NumericType>>> elementKdTree_,
+    SmartPointer<viennals::Mesh<float>> surfaceMesh_) {
+
+  const auto numRates = rayTracer_.getNumberOfRates();
+  const auto numPoints = rayTracer_.getNumberOfElements();
+  const auto numDisks = diskMesh->nodes.size();
+  assert(numRates > 0);
+  auto particles = rayTracer_.getParticles();
+  auto const &elementNormals =
+      *surfaceMesh_->getCellData().getVectorData("Normals");
+  auto const &normals = *diskMesh->getCellData().getVectorData("Normals");
+  const auto numElements = surfaceMesh_->lines.size();
+
+  NumericType conversionRadius = gridDelta * (smoothingNeighbors + 1);
+  conversionRadius *= conversionRadius; // use squared radius
+
+  std::vector<std::vector<std::pair<unsigned, NumericType>>> elementsToPoint;
+  elementsToPoint.reserve(numDisks);
+
+  for (int i = 0; i < numDisks; i++) {
+    auto closePoints =
+        elementKdTree_
+            ->findNearestWithinRadius(diskMesh->nodes[i], conversionRadius)
+            .value();
+
+    std::vector<std::pair<unsigned, NumericType>> closePointsArray;
+    std::vector<NumericType> weights(closePoints.size(), NumericType(0));
+
+    unsigned numClosePoints = 0;
+    for (std::size_t n = 0; n < closePoints.size(); ++n) {
+      const auto &p = closePoints[n];
+      assert(p.first < numElements);
+
+      NumericType weight = 0;
+      for (int k = 0; k < 3; ++k) {
+        weight += normals[i][k] * elementNormals[p.first][k];
+      }
+      weight = std::max(weight, NumericType(0));
+
+      if (weight > NumericType(1e-6) && !std::isnan(weight)) {
+        weights[n] = weight;
+        ++numClosePoints;
+      }
+    }
+
+    if (numClosePoints == 0) { // fallback to nearest point
+      auto nearestPoint = elementKdTree_->findNearest(diskMesh->nodes[i]);
+      closePointsArray.emplace_back(nearestPoint->first, NumericType(1));
+    }
+
+    // Compute weighted average
+    const NumericType sum =
+        std::accumulate(weights.begin(), weights.end(), NumericType(0));
+
+    if (sum > NumericType(0)) {
+      for (std::size_t k = 0; k < closePoints.size(); ++k) {
+        if (weights[k] > NumericType(0)) {
+          closePointsArray.emplace_back(closePoints[k].first, weights[k] / sum);
+        }
+      }
+    } else {
+      // Fallback if all weights were discarded
+      auto nearestPoint = elementKdTree_->findNearest(diskMesh->nodes[i]);
+      closePointsArray.emplace_back(nearestPoint->first, NumericType(1));
+    }
+
+    elementsToPoint.push_back(closePointsArray);
+  }
+  assert(elementsToPoint.size() == numDisks);
+
+  for (int pIdx = 0; pIdx < particles.size(); pIdx++) {
+    for (int dIdx = 0; dIdx < particles[pIdx].dataLabels.size(); dIdx++) {
+      auto elementFlux = rayTracer_.getFlux(pIdx, dIdx, smoothingNeighbors);
+      auto name = particles[pIdx].dataLabels[dIdx];
+
+      // convert line fluxes to disk fluxes
+      std::vector<NumericType> diskFlux(numDisks, 0.);
+
+      for (int i = 0; i < numDisks; i++) {
+        for (const auto &elemPair : elementsToPoint[i]) {
+          diskFlux[i] += static_cast<NumericType>(elementFlux[elemPair.first]) *
+                         elemPair.second;
+        }
+      }
+
+      pointData.insertReplaceScalarData(std::move(diskFlux), name);
+    }
+  }
 }

gpu/benchmark/CMakeLists.txt

@@ -18,3 +18,7 @@ add_dependencies(ViennaPS-GPU_Benchmark holeComparison)
 
 add_gpu_executable(holeComparisonGrid target_name ${CMAKE_CURRENT_SOURCE_DIR}/holeCompareGrid.cpp)
 add_dependencies(ViennaPS-GPU_Benchmark holeComparisonGrid)
+
+add_gpu_executable(transmissionTest target_name ${CMAKE_CURRENT_SOURCE_DIR}/transmission.cpp
+                   TransmissionTest.cu)
+add_dependencies(ViennaPS-GPU_Benchmark transmissionTest)