From 9987fa9a0785c9d225b55116c6cfa3258e6508e0 Mon Sep 17 00:00:00 2001 From: filipovic Date: Sat, 2 May 2026 00:17:34 +0200 Subject: [PATCH 1/8] Add GPU surface-source ray emission support for ViennaRay. --- include/viennaray/gpu/raygLaunchParams.hpp | 8 +++ include/viennaray/gpu/raygSource.hpp | 34 +++++++++ include/viennaray/gpu/raygTrace.hpp | 76 +++++++++++++++++++-- include/viennaray/gpu/raygTraceDisk.hpp | 4 +- include/viennaray/gpu/raygTraceTriangle.hpp | 4 +- 5 files changed, 119 insertions(+), 7 deletions(-) diff --git a/include/viennaray/gpu/raygLaunchParams.hpp b/include/viennaray/gpu/raygLaunchParams.hpp index 666df4b..3a357d0 100644 --- a/include/viennaray/gpu/raygLaunchParams.hpp +++ b/include/viennaray/gpu/raygLaunchParams.hpp @@ -53,6 +53,14 @@ struct LaunchParams { bool customDirectionBasis = false; } source; + bool useSurfaceSource = false; + viennacore::Vec3Df *surfaceSourcePositions = nullptr; + viennacore::Vec3Df *surfaceSourceNormals = nullptr; + float *surfaceSourceWeights = nullptr; + unsigned int surfaceSourceCount = 0; + float surfaceSourceArea = 0.f; + float surfaceSourceOffset = 0.f; + OptixTraversableHandle traversable; }; diff --git a/include/viennaray/gpu/raygSource.hpp b/include/viennaray/gpu/raygSource.hpp index 87cfcf7..ca20872 100644 --- a/include/viennaray/gpu/raygSource.hpp +++ b/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 &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,22 @@ 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, launchParams.cosineExponent, + getOrthonormalBasis(normal)); + return; + } + initializeRayPosition(prd, launchParams.source, launchParams.D); if (launchParams.source.customDirectionBasis) { initializeRayDirection(prd, launchParams.cosineExponent, diff --git a/include/viennaray/gpu/raygTrace.hpp b/include/viennaray/gpu/raygTrace.hpp index 5de79b1..8fdbbf0 100644 --- a/include/viennaray/gpu/raygTrace.hpp +++ b/include/viennaray/gpu/raygTrace.hpp @@ -114,20 +114,47 @@ template class Trace { launchParams_.maxReflections = config_.maxReflections; launchParams_.maxBoundaryHits = config_.maxBoundaryHits; + launchParams_.useSurfaceSource = surfaceSourceEnabled_; + if (surfaceSourceEnabled_) { + launchParams_.surfaceSourcePositions = + (Vec3Df *)surfaceSourcePositionsBuffer_.dPointer(); + launchParams_.surfaceSourceNormals = + (Vec3Df *)surfaceSourceNormalsBuffer_.dPointer(); + launchParams_.surfaceSourceWeights = + (float *)surfaceSourceWeightsBuffer_.dPointer(); + launchParams_.surfaceSourceCount = surfaceSourceCount_; + launchParams_.surfaceSourceArea = surfaceSourceArea_; + launchParams_.surfaceSourceOffset = surfaceSourceOffset_; + } + int numPointsPerDim = static_cast( std::sqrt(static_cast(launchParams_.numElements))); + unsigned int launchDimX = config_.numRaysPerPoint; + unsigned int launchDimY = numPointsPerDim; + unsigned int launchDimZ = numPointsPerDim; if (config_.numRaysFixed > 0) { numPointsPerDim = 1; config_.numRaysPerPoint = config_.numRaysFixed; + launchDimX = config_.numRaysPerPoint; + launchDimY = numPointsPerDim; + launchDimZ = numPointsPerDim; + } + + if (surfaceSourceEnabled_) { + launchDimY = surfaceSourceCount_; + launchDimZ = 1; } - numRays_ = numPointsPerDim * numPointsPerDim * config_.numRaysPerPoint; + numRays_ = static_cast(launchDimX) * launchDimY * launchDimZ; if (numRays_ > (1 << 29)) { VIENNACORE_LOG_WARNING("Too many rays for single launch: " + util::prettyDouble(numRays_)); - config_.numRaysPerPoint = (1 << 29) / (numPointsPerDim * numPointsPerDim); - numRays_ = numPointsPerDim * numPointsPerDim * config_.numRaysPerPoint; + const auto sourcePoints = static_cast(launchDimY) * launchDimZ; + config_.numRaysPerPoint = + std::max(1, (1 << 29) / sourcePoints); + launchDimX = config_.numRaysPerPoint; + numRays_ = static_cast(launchDimX) * sourcePoints; } VIENNACORE_LOG_DEBUG("Number of rays: " + util::prettyDouble(numRays_)); @@ -204,7 +231,7 @@ template class Trace { launchParamsBuffers_[i].dPointer(), launchParamsBuffers_[i].sizeInBytes, &shaderBindingTable_, /*! dimensions of the launch: */ - config_.numRaysPerPoint, numPointsPerDim, numPointsPerDim)); + launchDimX, launchDimY, launchDimZ)); } #else // Launch on multiple streams in release mode for (size_t i = 0; i < particles_.size(); i++) { @@ -214,7 +241,7 @@ template class Trace { launchParamsBuffers_[i].dPointer(), launchParamsBuffers_[i].sizeInBytes, &shaderBindingTable_, /*! dimensions of the launch: */ - config_.numRaysPerPoint, numPointsPerDim, numPointsPerDim)); + launchDimX, launchDimY, launchDimZ)); } #endif @@ -235,6 +262,37 @@ template class Trace { numCellData_ = numData; } + void setSurfaceSource(const std::vector &positions, + const std::vector &normals, + const std::vector &weights, + const float sourceArea, const float sourceOffset) { + if (positions.size() != normals.size() || positions.size() != weights.size()) { + VIENNACORE_LOG_ERROR("Surface source arrays must have matching sizes."); + } + if (positions.empty()) { + VIENNACORE_LOG_ERROR("Surface source must contain at least one point."); + } + + surfaceSourcePositionsBuffer_.allocUpload(positions); + surfaceSourceNormalsBuffer_.allocUpload(normals); + surfaceSourceWeightsBuffer_.allocUpload(weights); + surfaceSourceCount_ = static_cast(positions.size()); + surfaceSourceArea_ = sourceArea; + surfaceSourceOffset_ = sourceOffset; + surfaceSourceEnabled_ = true; + } + + void clearSurfaceSource() { + surfaceSourceEnabled_ = false; + surfaceSourceCount_ = 0; + surfaceSourceArea_ = 0.f; + surfaceSourceOffset_ = 0.f; + launchParams_.useSurfaceSource = false; + surfaceSourcePositionsBuffer_.free(); + surfaceSourceNormalsBuffer_.free(); + surfaceSourceWeightsBuffer_.free(); + } + template void setMaterialIds(const std::vector &materialIds, const bool mapToConsecutive = true) { @@ -361,6 +419,7 @@ template class Trace { for (auto &buffer : materialStickingBuffer_) { buffer.free(); } + clearSurfaceSource(); } void destroyMembers() { @@ -797,6 +856,9 @@ template class Trace { // sbt data CudaBuffer cellDataBuffer_; + CudaBuffer surfaceSourcePositionsBuffer_; + CudaBuffer surfaceSourceNormalsBuffer_; + CudaBuffer surfaceSourceWeightsBuffer_; OptixPipeline pipeline_{}; OptixPipelineCompileOptions pipelineCompileOptions_ = {}; @@ -828,11 +890,15 @@ template class Trace { rayInternal::KernelConfig config_; bool ignoreBoundary_ = false; + bool surfaceSourceEnabled_ = false; bool resultsDownloaded_ = false; bool isSynced_ = false; size_t numRays_ = 0; unsigned numCellData_ = 0; + unsigned int surfaceSourceCount_ = 0; + float surfaceSourceArea_ = 0.f; + float surfaceSourceOffset_ = 0.f; const std::string globalParamsName_ = "launchParams"; const std::string normModuleName_ = "normKernels.ptx"; diff --git a/include/viennaray/gpu/raygTraceDisk.hpp b/include/viennaray/gpu/raygTraceDisk.hpp index 1e9655f..e6cd995 100644 --- a/include/viennaray/gpu/raygTraceDisk.hpp +++ b/include/viennaray/gpu/raygTraceDisk.hpp @@ -87,7 +87,9 @@ template class TraceDisk final : public Trace { assert(this->resultBuffer_.sizeInBytes != 0 && "Normalization: Result buffer not initialized."); double sourceArea = 0.0; - if constexpr (D == 2) { + if (launchParams_.useSurfaceSource) { + sourceArea = launchParams_.surfaceSourceArea; + } else if constexpr (D == 2) { sourceArea = (launchParams_.source.maxPoint[0] - launchParams_.source.minPoint[0]); } else { diff --git a/include/viennaray/gpu/raygTraceTriangle.hpp b/include/viennaray/gpu/raygTraceTriangle.hpp index bda8f95..2c3eb34 100644 --- a/include/viennaray/gpu/raygTraceTriangle.hpp +++ b/include/viennaray/gpu/raygTraceTriangle.hpp @@ -42,7 +42,9 @@ template class TraceTriangle final : public Trace { void normalizeResults() override { double sourceArea = 0.0; - if constexpr (D == 2) { + if (launchParams_.useSurfaceSource) { + sourceArea = launchParams_.surfaceSourceArea; + } else if constexpr (D == 2) { sourceArea = (launchParams_.source.maxPoint[0] - launchParams_.source.minPoint[0]); } else { From 26956123b51f74473b0c4c0e89caf789e1899741 Mon Sep 17 00:00:00 2001 From: filipovic Date: Mon, 4 May 2026 09:48:13 +0200 Subject: [PATCH 2/8] match OptiX tracing for CPU to ignore backface triangle hits and continue instead of discarding the ray. --- include/viennaray/rayTraceKernel.hpp | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index 92b0f28..2d4d40c 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -241,10 +241,11 @@ template class TraceKernel { } } else { if (backfaceHit) { - // For triangle geometries, we simply discard backface hits as - // they are not considered valid geometry hits. - ++raysTerminated; - break; + // Match OptiX triangle tracing with back-face culling: ignore + // this hit and continue along the same ray. + reflect = true; + fillRayPosition(rayHit.ray, hitPoint); + continue; } } From 7faaf3d5e59c79d059023fa27ab3d680010c8c26 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Mon, 4 May 2026 10:50:41 +0200 Subject: [PATCH 3/8] refactor: remove unecessary launch params --- include/viennaray/gpu/raygLaunchParams.hpp | 6 +-- include/viennaray/gpu/raygSource.hpp | 3 +- include/viennaray/gpu/raygTrace.hpp | 41 ++++++++++----------- include/viennaray/gpu/raygTraceDisk.hpp | 2 +- include/viennaray/gpu/raygTraceTriangle.hpp | 2 +- 5 files changed, 25 insertions(+), 29 deletions(-) diff --git a/include/viennaray/gpu/raygLaunchParams.hpp b/include/viennaray/gpu/raygLaunchParams.hpp index 3a357d0..60b0856 100644 --- a/include/viennaray/gpu/raygLaunchParams.hpp +++ b/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; @@ -57,11 +59,7 @@ struct LaunchParams { viennacore::Vec3Df *surfaceSourcePositions = nullptr; viennacore::Vec3Df *surfaceSourceNormals = nullptr; float *surfaceSourceWeights = nullptr; - unsigned int surfaceSourceCount = 0; - float surfaceSourceArea = 0.f; float surfaceSourceOffset = 0.f; - - OptixTraversableHandle traversable; }; #ifdef __CUDACC__ diff --git a/include/viennaray/gpu/raygSource.hpp b/include/viennaray/gpu/raygSource.hpp index ca20872..381952f 100644 --- a/include/viennaray/gpu/raygSource.hpp +++ b/include/viennaray/gpu/raygSource.hpp @@ -113,8 +113,7 @@ initializeRayPositionAndDirection(PerRayData &prd, normal * launchParams.surfaceSourceOffset; prd.rayWeight = launchParams.surfaceSourceWeights[sourceIdx]; - initializeRayDirectionFromBasis(prd, launchParams.cosineExponent, - getOrthonormalBasis(normal)); + initializeRayDirectionFromBasis(prd, 1.f, getOrthonormalBasis(normal)); return; } diff --git a/include/viennaray/gpu/raygTrace.hpp b/include/viennaray/gpu/raygTrace.hpp index 8fdbbf0..4b6227d 100644 --- a/include/viennaray/gpu/raygTrace.hpp +++ b/include/viennaray/gpu/raygTrace.hpp @@ -122,8 +122,6 @@ template class Trace { (Vec3Df *)surfaceSourceNormalsBuffer_.dPointer(); launchParams_.surfaceSourceWeights = (float *)surfaceSourceWeightsBuffer_.dPointer(); - launchParams_.surfaceSourceCount = surfaceSourceCount_; - launchParams_.surfaceSourceArea = surfaceSourceArea_; launchParams_.surfaceSourceOffset = surfaceSourceOffset_; } @@ -142,7 +140,7 @@ template class Trace { } if (surfaceSourceEnabled_) { - launchDimY = surfaceSourceCount_; + launchDimY = surfaceSourcePositionsCount_; launchDimZ = 1; } @@ -225,23 +223,23 @@ template class Trace { #ifndef NDEBUG // Launch on single stream in debug mode for (size_t i = 0; i < particles_.size(); i++) { - OPTIX_CHECK(optixLaunch( - pipeline_, streams_[0], - /*! parameters and SBT */ - launchParamsBuffers_[i].dPointer(), - launchParamsBuffers_[i].sizeInBytes, &shaderBindingTable_, - /*! dimensions of the launch: */ - launchDimX, launchDimY, launchDimZ)); + OPTIX_CHECK(optixLaunch(pipeline_, streams_[0], + /*! parameters and SBT */ + launchParamsBuffers_[i].dPointer(), + launchParamsBuffers_[i].sizeInBytes, + &shaderBindingTable_, + /*! dimensions of the launch: */ + launchDimX, launchDimY, launchDimZ)); } #else // Launch on multiple streams in release mode for (size_t i = 0; i < particles_.size(); i++) { - OPTIX_CHECK(optixLaunch( - pipeline_, streams_[i], - /*! parameters and SBT */ - launchParamsBuffers_[i].dPointer(), - launchParamsBuffers_[i].sizeInBytes, &shaderBindingTable_, - /*! dimensions of the launch: */ - launchDimX, launchDimY, launchDimZ)); + OPTIX_CHECK(optixLaunch(pipeline_, streams_[i], + /*! parameters and SBT */ + launchParamsBuffers_[i].dPointer(), + launchParamsBuffers_[i].sizeInBytes, + &shaderBindingTable_, + /*! dimensions of the launch: */ + launchDimX, launchDimY, launchDimZ)); } #endif @@ -266,7 +264,8 @@ template class Trace { const std::vector &normals, const std::vector &weights, const float sourceArea, const float sourceOffset) { - if (positions.size() != normals.size() || positions.size() != weights.size()) { + if (positions.size() != normals.size() || + positions.size() != weights.size()) { VIENNACORE_LOG_ERROR("Surface source arrays must have matching sizes."); } if (positions.empty()) { @@ -276,7 +275,7 @@ template class Trace { surfaceSourcePositionsBuffer_.allocUpload(positions); surfaceSourceNormalsBuffer_.allocUpload(normals); surfaceSourceWeightsBuffer_.allocUpload(weights); - surfaceSourceCount_ = static_cast(positions.size()); + surfaceSourcePositionsCount_ = static_cast(positions.size()); surfaceSourceArea_ = sourceArea; surfaceSourceOffset_ = sourceOffset; surfaceSourceEnabled_ = true; @@ -284,7 +283,7 @@ template class Trace { void clearSurfaceSource() { surfaceSourceEnabled_ = false; - surfaceSourceCount_ = 0; + surfaceSourcePositionsCount_ = 0; surfaceSourceArea_ = 0.f; surfaceSourceOffset_ = 0.f; launchParams_.useSurfaceSource = false; @@ -896,7 +895,7 @@ template class Trace { size_t numRays_ = 0; unsigned numCellData_ = 0; - unsigned int surfaceSourceCount_ = 0; + unsigned int surfaceSourcePositionsCount_ = 0; float surfaceSourceArea_ = 0.f; float surfaceSourceOffset_ = 0.f; const std::string globalParamsName_ = "launchParams"; diff --git a/include/viennaray/gpu/raygTraceDisk.hpp b/include/viennaray/gpu/raygTraceDisk.hpp index e6cd995..dfebdd8 100644 --- a/include/viennaray/gpu/raygTraceDisk.hpp +++ b/include/viennaray/gpu/raygTraceDisk.hpp @@ -88,7 +88,7 @@ template class TraceDisk final : public Trace { "Normalization: Result buffer not initialized."); double sourceArea = 0.0; if (launchParams_.useSurfaceSource) { - sourceArea = launchParams_.surfaceSourceArea; + sourceArea = this->surfaceSourceArea_; } else if constexpr (D == 2) { sourceArea = (launchParams_.source.maxPoint[0] - launchParams_.source.minPoint[0]); diff --git a/include/viennaray/gpu/raygTraceTriangle.hpp b/include/viennaray/gpu/raygTraceTriangle.hpp index 2c3eb34..b37e00c 100644 --- a/include/viennaray/gpu/raygTraceTriangle.hpp +++ b/include/viennaray/gpu/raygTraceTriangle.hpp @@ -43,7 +43,7 @@ template class TraceTriangle final : public Trace { void normalizeResults() override { double sourceArea = 0.0; if (launchParams_.useSurfaceSource) { - sourceArea = launchParams_.surfaceSourceArea; + sourceArea = this->surfaceSourceArea_; } else if constexpr (D == 2) { sourceArea = (launchParams_.source.maxPoint[0] - launchParams_.source.minPoint[0]); From 7df44007116bdb606fc2068059e3004cafb45b1c Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Mon, 4 May 2026 10:51:01 +0200 Subject: [PATCH 4/8] feat: add surface source test --- gpu/examples/trenchTriangles.cpp | 46 +++++++++++++++++++++++++++++++- 1 file changed, 45 insertions(+), 1 deletion(-) diff --git a/gpu/examples/trenchTriangles.cpp b/gpu/examples/trenchTriangles.cpp index b4c0d18..a5f0af8 100644 --- a/gpu/examples/trenchTriangles.cpp +++ b/gpu/examples/trenchTriangles.cpp @@ -31,7 +31,7 @@ int main(int argc, char **argv) { particle.sticking = 1.f; particle.dataLabels = {"particleFlux"}; particle.materialSticking[7] = 0.1f; - particle.materialSticking[1] = 1.0f; + particle.materialSticking[1] = 0.1f; std::unordered_map pMap = {{"Particle", 0}}; std::vector cMap = { @@ -74,4 +74,48 @@ int main(int argc, char **argv) { rayCountBuffer.download(&rayCount, 1); std::cout << "Trace count: " << rayCount << std::endl; #endif + + // surface source test + gpu::TraceTriangle tracerSurface(context); + tracerSurface.setGeometry(mesh); + tracerSurface.setMaterialIds(materialIds); + tracerSurface.setCallables("ViennaRayCallableWrapper", context->modulePath); + tracerSurface.setParticleCallableMap({pMap, cMap}); + tracerSurface.setNumberOfRaysPerPoint(5000); + tracerSurface.insertNextParticle(particle); + tracerSurface.prepareParticlePrograms(); + + std::vector surfaceSourceWeights(mesh.triangles.size()); + assert(flux.size() == surfaceSourceWeights.size() && + "Flux size does not match surface source weights size."); + for (size_t i = 0; i < surfaceSourceWeights.size(); ++i) { + surfaceSourceWeights[i] = static_cast(flux[i]); + } + float sourceArea = 0.f; + std::vector surfaceSourcePosition(mesh.triangles.size()); + for (size_t i = 0; i < mesh.triangles.size(); ++i) { + sourceArea += + 0.5f * Norm(CrossProduct(mesh.nodes[mesh.triangles[i][1]] - + mesh.nodes[mesh.triangles[i][0]], + mesh.nodes[mesh.triangles[i][2]] - + mesh.nodes[mesh.triangles[i][0]])); + surfaceSourcePosition[i] = + (mesh.nodes[mesh.triangles[i][0]] + mesh.nodes[mesh.triangles[i][1]] + + mesh.nodes[mesh.triangles[i][2]]) / + 3.f; + } + + tracerSurface.setSurfaceSource(surfaceSourcePosition, mesh.normals, + surfaceSourceWeights, sourceArea, 1e-4f); + + tracerSurface.apply(); + tracerSurface.normalizeResults(); + + auto fluxSurface = tracerSurface.getFlux(0, 0); + + rayInternal::writeVTP( + "trenchTriangles_surfaceSource.vtp", mesh.nodes, mesh.triangles, + fluxSurface); + + return 0; } From 806460ca40410259307380115c0fb2f627674bab Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Wed, 6 May 2026 10:38:09 +0200 Subject: [PATCH 5/8] fix: use initial ray weight in rejection control --- gpu/examples/trenchTriangles.cpp | 72 ++++++++++++---------- gpu/pipelines/GeneralPipelineDisk.cu | 3 +- gpu/pipelines/GeneralPipelineLine.cu | 3 +- gpu/pipelines/GeneralPipelineTriangle.cu | 3 +- include/viennaray/gpu/raygLaunchParams.hpp | 7 ++- include/viennaray/gpu/raygTrace.hpp | 4 ++ 6 files changed, 52 insertions(+), 40 deletions(-) diff --git a/gpu/examples/trenchTriangles.cpp b/gpu/examples/trenchTriangles.cpp index a5f0af8..93a8320 100644 --- a/gpu/examples/trenchTriangles.cpp +++ b/gpu/examples/trenchTriangles.cpp @@ -26,12 +26,13 @@ int main(int argc, char **argv) { materialIds[i] = 1; } + float sticking = .5f; gpu::Particle particle; particle.name = "Particle"; - particle.sticking = 1.f; + particle.sticking = sticking; particle.dataLabels = {"particleFlux"}; - particle.materialSticking[7] = 0.1f; - particle.materialSticking[1] = 0.1f; + particle.materialSticking[7] = sticking; + particle.materialSticking[1] = sticking; std::unordered_map pMap = {{"Particle", 0}}; std::vector 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,7 +69,7 @@ int main(int argc, char **argv) { << std::endl; rayInternal::writeVTP( - "trenchTriangles_triMesh.vtp", mesh.nodes, mesh.triangles, flux); + "trenchTriangles_flux.vtp", mesh.nodes, mesh.triangles, flux); #ifdef COUNT_RAYS rayCountBuffer.download(&rayCount, 1); @@ -76,42 +77,45 @@ int main(int argc, char **argv) { #endif // surface source test - gpu::TraceTriangle tracerSurface(context); - tracerSurface.setGeometry(mesh); - tracerSurface.setMaterialIds(materialIds); - tracerSurface.setCallables("ViennaRayCallableWrapper", context->modulePath); - tracerSurface.setParticleCallableMap({pMap, cMap}); - tracerSurface.setNumberOfRaysPerPoint(5000); - tracerSurface.insertNextParticle(particle); - tracerSurface.prepareParticlePrograms(); - - std::vector surfaceSourceWeights(mesh.triangles.size()); - assert(flux.size() == surfaceSourceWeights.size() && - "Flux size does not match surface source weights size."); - for (size_t i = 0; i < surfaceSourceWeights.size(); ++i) { - surfaceSourceWeights[i] = static_cast(flux[i]); - } + double totalFlux = 0.0; float sourceArea = 0.f; std::vector surfaceSourcePosition(mesh.triangles.size()); + std::vector surfaceSourceWeights(mesh.triangles.size()); + std::vector areas(mesh.triangles.size()); for (size_t i = 0; i < mesh.triangles.size(); ++i) { - sourceArea += - 0.5f * Norm(CrossProduct(mesh.nodes[mesh.triangles[i][1]] - - mesh.nodes[mesh.triangles[i][0]], - mesh.nodes[mesh.triangles[i][2]] - - mesh.nodes[mesh.triangles[i][0]])); - surfaceSourcePosition[i] = - (mesh.nodes[mesh.triangles[i][0]] + mesh.nodes[mesh.triangles[i][1]] + - mesh.nodes[mesh.triangles[i][2]]) / - 3.f; + 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(flux[i]) * area * sticking; + totalFlux += flux[i] * area; + areas[i] = area; } - tracerSurface.setSurfaceSource(surfaceSourcePosition, mesh.normals, - surfaceSourceWeights, sourceArea, 1e-4f); + float averageArea = sourceArea / static_cast(mesh.triangles.size()); + for (size_t i = 0; i < surfaceSourceWeights.size(); ++i) { + surfaceSourceWeights[i] = surfaceSourceWeights[i] / averageArea; + } - tracerSurface.apply(); - tracerSurface.normalizeResults(); + std::cout << "Total flux from source plane: " << totalFlux * sticking + << std::endl; - auto fluxSurface = tracerSurface.getFlux(0, 0); + 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( "trenchTriangles_surfaceSource.vtp", mesh.nodes, mesh.triangles, diff --git a/gpu/pipelines/GeneralPipelineDisk.cu b/gpu/pipelines/GeneralPipelineDisk.cu index 8583feb..2f7134d 100644 --- a/gpu/pipelines/GeneralPipelineDisk.cu +++ b/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); diff --git a/gpu/pipelines/GeneralPipelineLine.cu b/gpu/pipelines/GeneralPipelineLine.cu index 51c9772..470738e 100644 --- a/gpu/pipelines/GeneralPipelineLine.cu +++ b/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); diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 997b4ee..6cffde6 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/gpu/pipelines/GeneralPipelineTriangle.cu @@ -86,6 +86,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 +98,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); diff --git a/include/viennaray/gpu/raygLaunchParams.hpp b/include/viennaray/gpu/raygLaunchParams.hpp index 60b0856..e1a29ae 100644 --- a/include/viennaray/gpu/raygLaunchParams.hpp +++ b/include/viennaray/gpu/raygLaunchParams.hpp @@ -82,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; @@ -93,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) { diff --git a/include/viennaray/gpu/raygTrace.hpp b/include/viennaray/gpu/raygTrace.hpp index 4b6227d..d2fc182 100644 --- a/include/viennaray/gpu/raygTrace.hpp +++ b/include/viennaray/gpu/raygTrace.hpp @@ -78,6 +78,10 @@ template class Trace { particles_.push_back(particle); } + void setRayWeightThreshold(float threshold) { + launchParams_.rayWeightThreshold = threshold; + } + void apply() { if (particles_.empty()) { VIENNACORE_LOG_ERROR( From 6c2c1e008101e0c8a509b2557e2befef6a6c70c7 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 7 May 2026 13:32:54 +0200 Subject: [PATCH 6/8] fix: github test CI --- .github/actions/setup/action.yml | 6 ------ .github/workflows/test.yml | 34 +++++++++++++++++++++----------- vcpkg.json | 8 ++++++++ 3 files changed, 30 insertions(+), 18 deletions(-) create mode 100644 vcpkg.json diff --git a/.github/actions/setup/action.yml b/.github/actions/setup/action.yml index 798ff42..eabde0a 100644 --- a/.github/actions/setup/action.yml +++ b/.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 diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml index 949735d..dce2937 100644 --- a/.github/workflows/test.yml +++ b/.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' }} diff --git a/vcpkg.json b/vcpkg.json new file mode 100644 index 0000000..2b6cb09 --- /dev/null +++ b/vcpkg.json @@ -0,0 +1,8 @@ +{ + "name": "viennaray-ci", + "version-string": "0.1.0", + "dependencies": [ + "embree" + ], + "builtin-baseline": "c3867e714dd3a51c272826eea77267876517ed99" +} \ No newline at end of file From 11c6a4625511f7f35a860d5a98f2c21395f04e65 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 7 May 2026 13:37:04 +0200 Subject: [PATCH 7/8] chore: bump version --- CMakeLists.txt | 4 ++-- README.md | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 18fa384..a645402 100644 --- a/CMakeLists.txt +++ b/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}") diff --git a/README.md b/README.md index 83d3e08..4ee04ec 100644 --- a/README.md +++ b/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 From 0fe89f4c19f4647bdd70613cf59e61ca79f90947 Mon Sep 17 00:00:00 2001 From: Tobias Reiter Date: Thu, 7 May 2026 15:17:02 +0200 Subject: [PATCH 8/8] refactor: discard backface hits always on triangles --- gpu/pipelines/GeneralPipelineTriangle.cu | 18 ++++++++++++++---- include/viennaray/rayTraceKernel.hpp | 16 +++++++--------- 2 files changed, 21 insertions(+), 13 deletions(-) diff --git a/gpu/pipelines/GeneralPipelineTriangle.cu b/gpu/pipelines/GeneralPipelineTriangle.cu index 6cffde6..6afedc5 100644 --- a/gpu/pipelines/GeneralPipelineTriangle.cu +++ b/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; @@ -110,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 diff --git a/include/viennaray/rayTraceKernel.hpp b/include/viennaray/rayTraceKernel.hpp index 2d4d40c..6473563 100644 --- a/include/viennaray/rayTraceKernel.hpp +++ b/include/viennaray/rayTraceKernel.hpp @@ -241,11 +241,10 @@ template class TraceKernel { } } else { if (backfaceHit) { - // Match OptiX triangle tracing with back-face culling: ignore - // this hit and continue along the same ray. - reflect = true; - fillRayPosition(rayHit.ray, hitPoint); - continue; + // For triangle geometries, we simply discard backface hits as + // they are not considered valid geometry hits. + ++raysTerminated; + break; } } @@ -290,11 +289,10 @@ template class TraceKernel { for (size_t diskId = 0; diskId < numDisksHit; ++diskId) { const auto matID = geometry_.getMaterialId(hitDiskIds[diskId]); const auto normal = geometry_.getPrimNormal(hitDiskIds[diskId]); + NumericType distRayWeight = rayWeight; #ifdef VIENNARAY_USE_WDIST - auto distRayWeight = rayWeight / impactDistances[diskId] / - invDistanceWeightSum * numDisksHit; -#else - auto distRayWeight = rayWeight; + distRayWeight = rayWeight / impactDistances[diskId] / + invDistanceWeightSum * numDisksHit; #endif particle->surfaceCollision(distRayWeight, rayDirection, normal, hitDiskIds[diskId], matID, myLocalData,