Fix GPU buffer freeing (#180)

* Bump dependencies, update transmission benchmark

* Bump ViennaRay

* Free coverage buffer

* Remove ViennaRay dev branch

Author: tobre1

CMakeLists.txt

@@ -101,7 +101,7 @@ include("cmake/cpm.cmake")
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 1.6.3
+  VERSION 1.7.0
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}
   OPTIONS "VIENNACORE_USE_GPU ${VIENNAPS_USE_GPU}")
@@ -113,7 +113,7 @@ CPMAddPackage(
 
 CPMFindPackage(
   NAME ViennaRay
-  VERSION 3.8.0
+  VERSION 3.8.2
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaRay"
   EXCLUDE_FROM_ALL ${VIENNAPS_BUILD_PYTHON}
   OPTIONS "VIENNARAY_USE_GPU ${VIENNAPS_USE_GPU}")

gpu/benchmark/transmission.cpp

@@ -5,11 +5,16 @@
 #include <raygTraceDisk.hpp>
 #include <raygTraceTriangle.hpp>
 
-// #define useDiskMesh
+#define USE_DISK_MESH
+#ifdef USE_DISK_MESH
+#define FILE_POST "_disk"
+#else
+#define FILE_POST "_tri"
+#endif
 
 constexpr int D = 3;
 using NumericType = float;
-constexpr bool useGPU = true;
+constexpr bool useGPU = false;
 
 template <class T> auto createCylinderLS(T gridDelta, T radius, T length) {
   double bounds[2 * D] = {-1.0, 1.0, -1.0, 1.0, 0.0, length + gridDelta};
@@ -104,12 +109,12 @@ auto createCylinderDiskMesh(NumericType gridDelta, NumericType radius,
   lsMesh->insertNextNode(Vec3Df{0.0f, 0.0f, 0.0f});
   materialIds.push_back(1);
   normals->push_back(Vec3Df{0.0f, 0.0f, 1.0f});
-  radii.push_back(radius);
+  radii.push_back(radius * 1.5f); // make sure rays hit the plane
 
   lsMesh->insertNextNode(Vec3Df{0.0f, 0.0f, float(length)});
   materialIds.push_back(2);
   normals->push_back(Vec3Df{0.0f, 0.0f, -1.0f});
-  radii.push_back(radius);
+  radii.push_back(radius * 1.5f); // make sure rays hit the plane
 
   // viennals::VTKWriter<float>(lsMesh, "cylinder.vtp").apply();
 
@@ -209,16 +214,19 @@ std::array<NumericType, N> logSpace(NumericType start, NumericType end) {
 }
 
 int main() {
-  auto lengths = logSpace<20>(.01, 10000.0);
+  auto lengths = logSpace<25>(.01, 10000.0);
   // std::array<NumericType, 1> lengths = {100.0};
   constexpr NumericType radius = 10.0;
   constexpr NumericType gridDelta = .5;
   constexpr int nRays = int(1e7);
   std::vector<unsigned> counts(4, 0);
 
   if constexpr (useGPU) {
-    viennacore::CudaBuffer customDataBuffer;
-    customDataBuffer.allocUpload(counts);
+#ifdef USE_DISK_MESH
+    viennacore::Logger::getInstance()
+        .addError("No GPU disk mesh implementation")
+        .print();
+#endif
 
     auto context = viennacore::DeviceContext::createContext();
     viennaray::gpu::TraceTriangle<NumericType, D> tracer(context);
@@ -245,7 +253,10 @@ int main() {
     tracer.insertNextParticle(particle);
     tracer.prepareParticlePrograms();
 
-    std::ofstream logFile("transmission_log_tri_GPU.txt");
+    std::string fileName = "transmission_log_GPU";
+    fileName.append(FILE_POST);
+    fileName.append(".txt");
+    std::ofstream logFile(fileName);
     logFile << "LR\tTransmission\tReflection\tHitsToTransmission\tHitsToReflect"
                "ion\n";
 
@@ -281,10 +292,9 @@ int main() {
     }
 
     logFile.close();
-    customDataBuffer.free();
   } else {
 
-#ifdef useDiskMesh
+#ifdef USE_DISK_MESH
     viennaray::TraceDisk<NumericType, D> tracer;
 #else
     viennaray::TraceTriangle<NumericType, D> tracer;
@@ -302,14 +312,17 @@ int main() {
     auto particle = std::make_unique<TransmissionParticle>();
     tracer.setParticleType(particle);
 
-    std::ofstream logFile("transmission_log_disk_CPU.txt");
+    std::string fileName = "transmission_log_CPU";
+    fileName.append(FILE_POST);
+    fileName.append(".txt");
+    std::ofstream logFile(fileName);
     logFile << "LR\tTransmission\tReflection\tHitsToTransmission\tHitsToReflect"
                "ion\n";
 
     for (auto length : lengths) {
       std::cout << "Length: " << length << std::endl;
       std::cout << "L/R: " << length / radius << std::endl;
-#ifdef useDiskMesh
+#ifdef USE_DISK_MESH
       auto [mesh, materialIds] =
           createCylinderDiskMesh(gridDelta, radius, length);
 #else
@@ -323,6 +336,13 @@ int main() {
 
       tracer.apply();
 
+      // auto info = tracer.getRayTraceInfo();
+      // std::cout << "Non geo hits: " << info.nonGeometryHits << std::endl;
+      // std::cout << "Num rays: " << info.numRays << std::endl;
+      // std::cout << "Boundary hits: " << info.boundaryHits << std::endl;
+      // std::cout << "Reflections: " << info.reflections << std::endl;
+      // std::cout << "Rays terminated: " << info.raysTerminated << std::endl;
+
       auto &localData = tracer.getLocalData();
       counts[0] = static_cast<unsigned>(
           localData.getVectorData(0)[0]); // Transmission count

gpu/benchmark/transmission.ipynb

@@ -1,6 +1,6 @@
 #pragma once
 
-#include <raygRNG.hpp>
+#include <vcRNG.hpp>
 
 namespace viennaps::gpu::impl {
 struct IonParams {
@@ -64,7 +64,7 @@ __forceinline__ __device__ void updateEnergy(viennaray::gpu::PerRayData *prd,
 
   float newEnergy;
   do {
-    const float u = viennaray::gpu::getNormalDistRand(&prd->RNGstate);
+    const float u = viennacore::getNormalDistRand(&prd->RNGstate);
     newEnergy = prd->energy * (Eref_peak + sigma * u);
   } while (newEnergy < 0.f || newEnergy > prd->energy);
   prd->energy = newEnergy;
@@ -97,7 +97,7 @@ initNormalDistEnergy(viennaray::gpu::PerRayData *prd, const float mean,
     energy = mean;
   } else {
     do {
-      const float u = viennaray::gpu::getNormalDistRand(&prd->RNGstate);
+      const float u = viennacore::getNormalDistRand(&prd->RNGstate);
       energy = mean + sigma * u;
     } while (energy < 0.f);
   }

gpu/include/psgPipelineParameters.hpp

@@ -207,11 +207,9 @@ class GPUDiskEngine final : public FluxEngine<NumericType, D> {
           context.diskMesh, context.getProcessName() + "_flux_" +
                                 std::to_string(iterations++) + ".vtp")
           .apply();
-
-      if (context.flags.useCoverages) {
-        d_coverages.free();
-      }
     }
+
+    d_coverages.free();
     this->timer_.finish();
 
     return ProcessResult::SUCCESS;

include/viennaps/process/psGPUDiskEngine.hpp

@@ -114,8 +114,8 @@ class GPULineEngine final : public FluxEngine<NumericType, D> {
       elementKdTree_ = KDTreeType::New();
 
     CreateSurfaceMesh<NumericType, float, D>(
-        context.domain->getLevelSets().back(), surfaceMesh_, elementKdTree_,
-        1e-12, context.rayTracingParams.minNodeDistanceFactor)
+        context.domain->getSurface(), surfaceMesh_, elementKdTree_, 1e-12,
+        context.rayTracingParams.minNodeDistanceFactor)
         .apply();
 
     viennaray::LineMesh lineMesh(
@@ -258,11 +258,9 @@ class GPULineEngine final : public FluxEngine<NumericType, D> {
           context.diskMesh, context.getProcessName() + "_flux_" +
                                 std::to_string(iterations++) + ".vtp")
           .apply();
-
-      if (context.flags.useCoverages) {
-        d_coverages.free();
-      }
     }
+
+    d_coverages.free();
     this->timer_.finish();
 
     return ProcessResult::SUCCESS;

include/viennaps/process/psGPULineEngine.hpp

@@ -111,8 +111,8 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
       elementKdTree_ = KDTreeType::New();
 
     CreateSurfaceMesh<NumericType, float, D>(
-        context.domain->getLevelSets().back(), surfaceMesh_, elementKdTree_,
-        1e-12, context.rayTracingParams.minNodeDistanceFactor)
+        context.domain->getSurface(), surfaceMesh_, elementKdTree_, 1e-12,
+        context.rayTracingParams.minNodeDistanceFactor)
         .apply();
     assert(!surfaceMesh_->nodes.empty());
 
@@ -248,11 +248,9 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
                                  context.getProcessName() + "_flux_" +
                                      std::to_string(iterations++) + ".vtp")
           .apply();
-
-      if (context.flags.useCoverages) {
-        d_coverages.free();
-      }
     }
+
+    d_coverages.free();
     this->timer_.finish();
 
     return ProcessResult::SUCCESS;
@@ -285,7 +283,7 @@ class GPUTriangleEngine final : public FluxEngine<NumericType, D> {
     const auto numPoints = rayTracer_.getNumberOfElements();
     assert(numRates > 0);
     assert(numPoints == surfaceMesh_->getElements<3>().size());
-    auto valueBuffer = rayTracer_.getResults();
+    auto &valueBuffer = rayTracer_.getResults();
     std::vector<float> tmpBuffer(numRates * numPoints);
     valueBuffer.download(tmpBuffer.data(), numPoints * numRates);
     auto particles = rayTracer_.getParticles();

include/viennaps/process/psGPUTriangleEngine.hpp

@@ -193,6 +193,11 @@ template <typename NumericType, int D> class Process {
 
   // Factory method for creating flux engines
   std::unique_ptr<FluxEngine<NumericType, D>> createFluxEngine() {
+    assert(fluxEngineType_ != FluxEngineType::AUTO &&
+           "Flux engine type must be specified before creation.");
+    Logger::getInstance()
+        .addDebug("Creating flux engine of type: " + to_string(fluxEngineType_))
+        .print();
     // Create CPU engine
     if (fluxEngineType_ == FluxEngineType::CPU_DISK) {
       return std::make_unique<CPUDiskEngine<NumericType, D>>();

include/viennaps/process/psProcess.hpp

@@ -511,6 +511,8 @@ template <int D> void bindApi(py::module &module) {
       .def("generateCellSet", &Domain<T, D>::generateCellSet,
            "Generate the cell set.")
       .def("getLevelSets", &Domain<T, D>::getLevelSets)
+      .def("getSurface", &Domain<T, D>::getSurface,
+           "Get the surface level set.")
       .def("getCellSet", &Domain<T, D>::getCellSet, "Get the cell set.")
       .def("getGrid", &Domain<T, D>::getGrid, "Get the grid")
       .def("getGridDelta", &Domain<T, D>::getGridDelta, "Get the grid delta.")

python/pyWrapDimension.hpp

@@ -549,6 +549,11 @@ class Domain:
         Get the domain setup.
         """
 
+    def getSurface(self) -> viennals.d2.Domain:
+        """
+        Get the surface level set.
+        """
+
     def getSurfaceMesh(
         self,
         addInterfaces: bool = False,
@@ -1446,6 +1451,7 @@ class ToDiskMesh:
     def __init__(self, domain: Domain, mesh: viennals._core.Mesh) -> None: ...
     @typing.overload
     def __init__(self) -> None: ...
+    def apply(self) -> None: ...
     def setDomain(self, arg0: Domain) -> None:
         """
         Set the domain in the mesh converter.