feat: gpu surface source

.github/actions/setup/action.yml

@@ -51,9 +51,3 @@ runs:
       run: |
         echo "MACOSX_DEPLOYMENT_TARGET=15.0" >> $GITHUB_ENV
 
-    - name: Setup vcpkg (Windows)
-      shell: bash
-      if: ${{ inputs.os == 'windows-latest' }}
-      run: |
-        git clone https://github.com/microsoft/vcpkg.git
-        ./vcpkg/bootstrap-vcpkg.bat

.github/workflows/test.yml

@@ -33,42 +33,52 @@ jobs:
     container: ${{ matrix.container }}
 
     env:
+      VCPKG_COMMIT: c3867e714dd3a51c272826eea77267876517ed99
       VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/.vcpkg-bincache
 
     name: "🧪 Test on ${{ matrix.os }} (⚙️: ${{ matrix.config }}, 💿: ${{ matrix.container || matrix.os }})"
 
     steps:
       - name: 📥 Checkout
-        uses: actions/checkout@v4
+        uses: actions/checkout@v6
 
       - name: 🖥️ Setup Environment
         uses: ./.github/actions/setup
         with:
           os: ${{ matrix.os }}
 
+      - name: 🖥️ Setup vcpkg (Windows)
+        shell: pwsh
+        if: ${{ matrix.os == 'windows-latest' }}
+        run: |
+          Remove-Item Env:VCPKG_ROOT -ErrorAction SilentlyContinue
+
+          git clone https://github.com/microsoft/vcpkg.git
+          git -C vcpkg checkout $env:VCPKG_COMMIT
+          & .\vcpkg\bootstrap-vcpkg.bat
+
+          "VCPKG_ROOT=$($env:GITHUB_WORKSPACE)\vcpkg" | Out-File -FilePath $env:GITHUB_ENV -Encoding utf8 -Append
+
+          New-Item -ItemType Directory -Force -Path "${{ env.VCPKG_DEFAULT_BINARY_CACHE }}"
+
       - name: 📋 Install MacOS Dependencies
         shell: bash
         if: ${{ matrix.os == 'macos-latest' }}
         run: |
           brew install embree tbb
 
-      - name: Create vcpkg binary cache dir
+      - name: 🦥 Cache vcpkg binary (Windows)
         if: ${{ matrix.os == 'windows-latest' }}
-        shell: pwsh
-        run: |
-          New-Item -ItemType Directory -Force -Path "${{ env.VCPKG_DEFAULT_BINARY_CACHE }}"
-
-      - name: 🦥 Cache vcpkg build
-        if: ${{ matrix.os == 'windows-latest' }}
-        uses: actions/cache@v4
+        uses: actions/cache@v5
         with:
           path: ${{ env.VCPKG_DEFAULT_BINARY_CACHE }}
-          key: vcpkg-${{ matrix.os }}-${{ matrix.config }}
+          key: vcpkg-${{ matrix.os }}-${{ matrix.config }}-${{ env.VCPKG_COMMIT }}
 
-      - name: 🛠️ Build embree (Windows)
+      - name: 🛠️ Build Dependencies (Windows)
         if: ${{ matrix.os == 'windows-latest' }}
+        shell: pwsh
         run: |
-          ./vcpkg/vcpkg install embree
+          ./vcpkg/vcpkg install --triplet x64-windows
 
       - name: 🏗️ Compile (Windows)
         if: ${{ matrix.os == 'windows-latest' }}

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaRay
   LANGUAGES CXX
-  VERSION 4.1.2)
+  VERSION 4.2.0)
 
 # --------------------------------------------------------------------------------------------------------
 # Library switches
@@ -98,7 +98,7 @@ include("cmake/cpm.cmake")
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 2.1.1
+  VERSION 2.1.2
   GIT_REPOSITORY "https://github.com/ViennaTools/ViennaCore"
   OPTIONS "VIENNACORE_USE_GPU ${VIENNARAY_USE_GPU}")
 

README.md

@@ -63,7 +63,7 @@ We recommend using [CPM.cmake](https://github.com/cpm-cmake/CPM.cmake) to consum
 * Installation with CPM
 
   ```cmake
-  CPMAddPackage("gh:viennatools/viennaray@4.1.0") # Use the latest release version
+  CPMAddPackage("gh:viennatools/viennaray@4.2.0") # Use the latest release version
   ```
 
 * With a local installation

gpu/examples/trenchTriangles.cpp

@@ -26,12 +26,13 @@ int main(int argc, char **argv) {
     materialIds[i] = 1;
   }
 
+  float sticking = .5f;
   gpu::Particle<NumericType> particle;
   particle.name = "Particle";
-  particle.sticking = 1.f;
+  particle.sticking = sticking;
   particle.dataLabels = {"particleFlux"};
-  particle.materialSticking[7] = 0.1f;
-  particle.materialSticking[1] = 1.0f;
+  particle.materialSticking[7] = sticking;
+  particle.materialSticking[1] = sticking;
 
   std::unordered_map<std::string, unsigned int> pMap = {{"Particle", 0}};
   std::vector<gpu::CallableConfig> cMap = {
@@ -44,7 +45,7 @@ int main(int argc, char **argv) {
   tracer.setMaterialIds(materialIds);
   tracer.setCallables("ViennaRayCallableWrapper", context->modulePath);
   tracer.setParticleCallableMap({pMap, cMap});
-  tracer.setNumberOfRaysPerPoint(5000);
+  tracer.setNumberOfRaysPerPoint(3000);
   tracer.insertNextParticle(particle);
   tracer.prepareParticlePrograms();
 
@@ -68,10 +69,57 @@ int main(int argc, char **argv) {
             << std::endl;
 
   rayInternal::writeVTP<float, D, gpu::ResultType>(
-      "trenchTriangles_triMesh.vtp", mesh.nodes, mesh.triangles, flux);
+      "trenchTriangles_flux.vtp", mesh.nodes, mesh.triangles, flux);
 
 #ifdef COUNT_RAYS
   rayCountBuffer.download(&rayCount, 1);
   std::cout << "Trace count: " << rayCount << std::endl;
 #endif
+
+  // surface source test
+  double totalFlux = 0.0;
+  float sourceArea = 0.f;
+  std::vector<Vec3Df> surfaceSourcePosition(mesh.triangles.size());
+  std::vector<float> surfaceSourceWeights(mesh.triangles.size());
+  std::vector<float> areas(mesh.triangles.size());
+  for (size_t i = 0; i < mesh.triangles.size(); ++i) {
+    const auto &A = mesh.nodes[mesh.triangles[i][0]];
+    const auto &B = mesh.nodes[mesh.triangles[i][1]];
+    const auto &C = mesh.nodes[mesh.triangles[i][2]];
+    float area = 0.5f * Norm(CrossProduct(B - A, C - A));
+    sourceArea += area;
+    surfaceSourcePosition[i] = (A + B + C) / 3.f;
+    surfaceSourceWeights[i] = static_cast<float>(flux[i]) * area * sticking;
+    totalFlux += flux[i] * area;
+    areas[i] = area;
+  }
+
+  float averageArea = sourceArea / static_cast<float>(mesh.triangles.size());
+  for (size_t i = 0; i < surfaceSourceWeights.size(); ++i) {
+    surfaceSourceWeights[i] = surfaceSourceWeights[i] / averageArea;
+  }
+
+  std::cout << "Total flux from source plane: " << totalFlux * sticking
+            << std::endl;
+
+  tracer.setSurfaceSource(surfaceSourcePosition, mesh.normals,
+                          surfaceSourceWeights, sourceArea, 1e-4f);
+
+  tracer.apply();
+  tracer.normalizeResults();
+
+  auto fluxSurface = tracer.getFlux(0, 0);
+
+  double totalSurfaceFlux = 0.0;
+  for (size_t i = 0; i < mesh.triangles.size(); ++i) {
+    totalSurfaceFlux += fluxSurface[i] * areas[i];
+  }
+  std::cout << "Total flux from surface desorption: "
+            << totalSurfaceFlux * sticking << std::endl;
+
+  rayInternal::writeVTP<float, D, gpu::ResultType>(
+      "trenchTriangles_surfaceSource.vtp", mesh.nodes, mesh.triangles,
+      fluxSurface);
+
+  return 0;
 }

gpu/pipelines/GeneralPipelineDisk.cu

@@ -164,6 +164,7 @@ extern "C" __global__ void __raygen__() {
 
   // initialize ray position and direction
   initializeRayPositionAndDirection(prd, launchParams);
+  const float initialRayWeight = prd.rayWeight;
 
   unsigned callIdx =
       callableIndex(launchParams.particleType, CallableSlot::INIT);
@@ -175,7 +176,7 @@ extern "C" __global__ void __raygen__() {
   packPointer((void *)&prd, u0, u1);
   unsigned int hintBitLength = 2;
 
-  while (continueRay(launchParams, prd)) {
+  while (continueRay(launchParams, prd, initialRayWeight)) {
     if (launchParams.D == 2) {
       prd.traceDir[2] = 0.f;
       viennacore::Normalize(prd.traceDir);

gpu/pipelines/GeneralPipelineLine.cu

@@ -107,6 +107,7 @@ extern "C" __global__ void __raygen__() {
 
   // initialize ray position and direction
   initializeRayPositionAndDirection(prd, launchParams);
+  const float initialRayWeight = prd.rayWeight;
 
   unsigned callIdx =
       callableIndex(launchParams.particleType, CallableSlot::INIT);
@@ -118,7 +119,7 @@ extern "C" __global__ void __raygen__() {
   packPointer((void *)&prd, u0, u1);
   unsigned int hintBitLength = 2;
 
-  while (continueRay(launchParams, prd)) {
+  while (continueRay(launchParams, prd, initialRayWeight)) {
     if (launchParams.D == 2) {
       prd.traceDir[2] = 0.f;
       viennacore::Normalize(prd.traceDir);

gpu/pipelines/GeneralPipelineTriangle.cu

@@ -17,9 +17,16 @@ using namespace viennaray::gpu;
 extern "C" __constant__ LaunchParams launchParams;
 
 extern "C" __global__ void __closesthit__() {
+  PerRayData *prd = getPRD();
+
+  if (optixIsTriangleBackFaceHit()) {
+    // Discard geometry back face hits for triangles
+    prd->rayWeight = 0.f;
+    return;
+  }
+
   const HitSBTDataTriangle *sbtData =
       (const HitSBTDataTriangle *)optixGetSbtDataPointer();
-  PerRayData *prd = getPRD();
 
   const unsigned int primID = optixGetPrimitiveIndex();
   prd->tMin = optixGetRayTmax();
@@ -48,6 +55,11 @@ extern "C" __global__ void __closesthit__boundary__() {
   // update ray position to hit point
   prd->pos = prd->pos + prd->traceDir * optixGetRayTmax();
 
+  if (optixIsTriangleBackFaceHit()) {
+    // Continue ray without any changes
+    return;
+  }
+
   const unsigned int primID = optixGetPrimitiveIndex();
   // 0-3: X axis (dim 0), 4-7: Y axis (dim 1)
   const unsigned int dim = primID / 4;
@@ -86,6 +98,7 @@ extern "C" __global__ void __raygen__() {
 
   // initialize ray position and direction
   initializeRayPositionAndDirection(prd, launchParams);
+  const float initialRayWeight = prd.rayWeight;
 
   unsigned callIdx =
       callableIndex(launchParams.particleType, CallableSlot::INIT);
@@ -97,7 +110,7 @@ extern "C" __global__ void __raygen__() {
   packPointer((void *)&prd, u0, u1);
   unsigned int hintBitLength = 2;
 
-  while (continueRay(launchParams, prd)) {
+  while (continueRay(launchParams, prd, initialRayWeight)) {
     if (launchParams.D == 2) {
       prd.traceDir[2] = 0.f;
       viennacore::Normalize(prd.traceDir);
@@ -109,9 +122,7 @@ extern "C" __global__ void __raygen__() {
                   1e-4f,                        // tmin
                   1e20f,                        // tmax
                   0.0f,                         // rayTime
-                  OptixVisibilityMask(255),
-                  OPTIX_RAY_FLAG_DISABLE_ANYHIT |
-                      OPTIX_RAY_FLAG_CULL_BACK_FACING_TRIANGLES,
+                  OptixVisibilityMask(255), OPTIX_RAY_FLAG_DISABLE_ANYHIT,
                   0,       // SBT offset
                   1,       // SBT stride
                   0,       // missSBTIndex

include/viennaray/gpu/raygLaunchParams.hpp

@@ -20,6 +20,8 @@ __both__ __forceinline__ unsigned callableIndex(unsigned p, CallableSlot s) {
 }
 
 struct LaunchParams {
+  OptixTraversableHandle traversable;
+
   ResultType *resultBuffer;
 
   float rayWeightThreshold = 0.1f;
@@ -53,7 +55,11 @@ struct LaunchParams {
     bool customDirectionBasis = false;
   } source;
 
-  OptixTraversableHandle traversable;
+  bool useSurfaceSource = false;
+  viennacore::Vec3Df *surfaceSourcePositions = nullptr;
+  viennacore::Vec3Df *surfaceSourceNormals = nullptr;
+  float *surfaceSourceWeights = nullptr;
+  float surfaceSourceOffset = 0.f;
 };
 
 #ifdef __CUDACC__
@@ -76,7 +82,8 @@ getIdxOffset(int dataIdx, const LaunchParams &launchParams) {
 }
 
 __device__ __forceinline__ bool continueRay(const LaunchParams &launchParams,
-                                            PerRayData &prd) {
+                                            PerRayData &prd,
+                                            const float &initialRayWeight) {
   if (prd.rayWeight <= 0.f || prd.energy < 0.f)
     return false;
 
@@ -87,13 +94,13 @@ __device__ __forceinline__ bool continueRay(const LaunchParams &launchParams,
   // If the weight of the ray is above a certain threshold, we always reflect.
   // If the weight of the ray is below the threshold, we randomly decide to
   // either kill the ray or increase its weight (in an unbiased way).
-  if (prd.rayWeight >= launchParams.rayWeightThreshold && prd.energy >= 0.f)
+  if (prd.rayWeight >= launchParams.rayWeightThreshold * initialRayWeight)
     return true;
 
   // We want to set the weight of (the reflection of) the ray to the value of
   // renewWeight. In order to stay unbiased we kill the reflection with a
   // probability of (1 - rayWeight / renewWeight).
-  float renewWeight = 0.3f;
+  float renewWeight = 0.3f * initialRayWeight;
   float rnd = getNextRand(&prd.RNGstate);
   float killProbability = 1.f - prd.rayWeight / renewWeight;
   if (rnd < killProbability) {

include/viennaray/gpu/raygSource.hpp

@@ -62,6 +62,24 @@ initializeRayDirection(PerRayData &prd, const float power,
   prd.traceDir = prd.dir;
 }
 
+__device__ __forceinline__ void
+initializeRayDirectionFromBasis(PerRayData &prd, const float power,
+                                const std::array<Vec3Df, 3> &basis) {
+  const float4 u = curand_uniform4(&prd.RNGstate); // (0,1]
+  const float cosTheta = powf(u.w, 1.f / (power + 1.f));
+  const float sinTheta = sqrtf(max(0.f, 1.f - cosTheta * cosTheta));
+  float sinPhi, cosPhi;
+  __sincosf(2.f * M_PIf * u.x, &sinPhi, &cosPhi);
+
+  const float rx = cosTheta;
+  const float ry = cosPhi * sinTheta;
+  const float rz = sinPhi * sinTheta;
+
+  prd.dir = basis[0] * rx + basis[1] * ry + basis[2] * rz;
+  viennacore::Normalize(prd.dir);
+  prd.traceDir = prd.dir;
+}
+
 __device__ __forceinline__ void
 initializeRayPosition(PerRayData &prd, const LaunchParams::SourcePlane &source,
                       const uint16_t D) {
@@ -84,6 +102,21 @@ initializeRayPosition(PerRayData &prd, const LaunchParams::SourcePlane &source,
 __device__ __forceinline__ void
 initializeRayPositionAndDirection(PerRayData &prd,
                                   const LaunchParams &launchParams) {
+  if (launchParams.useSurfaceSource) {
+    const uint3 launchIdx = optixGetLaunchIndex();
+    const unsigned sourceIdx = launchIdx.y;
+
+    Vec3Df normal = launchParams.surfaceSourceNormals[sourceIdx];
+    viennacore::Normalize(normal);
+
+    prd.pos = launchParams.surfaceSourcePositions[sourceIdx] +
+              normal * launchParams.surfaceSourceOffset;
+    prd.rayWeight = launchParams.surfaceSourceWeights[sourceIdx];
+
+    initializeRayDirectionFromBasis(prd, 1.f, getOrthonormalBasis(normal));
+    return;
+  }
+
   initializeRayPosition(prd, launchParams.source, launchParams.D);
   if (launchParams.source.customDirectionBasis) {
     initializeRayDirection(prd, launchParams.cosineExponent,