Improved power cosine sampling

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 project(
   ViennaRay
   LANGUAGES CXX
-  VERSION 3.8.3)
+  VERSION 3.8.5)
 
 # --------------------------------------------------------------------------------------------------------
 # Library switches
@@ -92,7 +92,7 @@ include("cmake/cpm.cmake")
 
 CPMAddPackage(
   NAME ViennaCore
-  VERSION 1.7.0
+  VERSION 1.7.3
   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@3.8.2") # Use the latest release version
+  CPMAddPackage("gh:viennatools/viennaray@3.8.5") # Use the latest release version
   ```
 
 * With a local installation

examples/disk3D/disk3D.cpp

@@ -61,6 +61,7 @@ int main() {
   // Extract the normalized hit counts for each geometry point
   auto &flux = rayTracer.getLocalData().getVectorData("flux");
   rayTracer.normalizeFlux(flux, NormalizationType::SOURCE);
+  rayTracer.smoothFlux(flux);
 
   rayInternal::writeVTK<NumericType, D>("trenchResult.vtk", points, flux);
 

gpu/include/raygReflection.hpp

@@ -66,10 +66,10 @@ specularReflection(PerRayData *prd, const Vec3Df &geoNormal) {
 
 static __device__ Vec3Df PickRandomPointOnUnitSphere(CudaRNG *state) {
   const float4 u = curand_uniform4(state); // (0,1]
-  const float z = 1.0f - 2.0f * u.x;       // uniform in [-1,1]
-  const float r2 = fmaxf(0.0f, 1.0f - z * z);
+  const float z = 1.f - 2.f * u.x;         // uniform in [-1,1]
+  const float r2 = max(0.f, 1.f - z * z);
   const float r = sqrtf(r2);
-  const float phi = 2.0f * M_PIf * u.y;
+  const float phi = 2.f * M_PIf * u.y;
   float s, c;
   __sincosf(phi, &s, &c); // branch-free sin+cos
   return Vec3Df{r * c, r * s, z};

gpu/include/raygSource.hpp

@@ -29,19 +29,21 @@ __device__ __forceinline__ void
 initializeRayDirection(PerRayData *prd, const float power, const uint16_t D) {
   // source direction
   const float4 u = curand_uniform4(&prd->RNGstate); // (0,1]
-  const float tt = powf(u.w, 2.f / (power + 1.f));
-  float s, c;
-  __sincosf(2.f * M_PIf * u.x, &s, &c);
-  float sqrt1mtt = sqrtf(1 - tt);
-  prd->dir[0] = c * sqrt1mtt;
+  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);
+
+  prd->dir[0] = cosPhi * sinTheta;
   if (D == 2) {
-    prd->dir[1] = -1.f * sqrtf(tt);
+    prd->dir[1] = -cosTheta;
     prd->dir[2] = 0.f;
+    Normalize(prd->dir);
   } else {
-    prd->dir[1] = s * sqrt1mtt;
-    prd->dir[2] = -1.f * sqrtf(tt);
+    prd->dir[1] = sinPhi * sinTheta;
+    prd->dir[2] = -cosTheta;
+    // already normalized
   }
-  Normalize(prd->dir);
 }
 
 __device__ __forceinline__ void
@@ -50,20 +52,16 @@ initializeRayDirection(PerRayData *prd, const float power,
   // source direction
   do {
     const float4 u = curand_uniform4(&prd->RNGstate); // (0,1]
-    const float tt = powf(u.w, 2.f / (power + 1.f));
-    Vec3Df rndDirection;
-    rndDirection[0] = sqrtf(tt);
-    float s, c, sqrt1mtt = sqrtf(1 - tt);
-    __sincosf(2.f * M_PIf * u.x, &s, &c);
-    rndDirection[1] = c * sqrt1mtt;
-    rndDirection[2] = s * sqrt1mtt;
-
-    prd->dir[0] = basis[0][0] * rndDirection[0] +
-                  basis[1][0] * rndDirection[1] + basis[2][0] * rndDirection[2];
-    prd->dir[1] = basis[0][1] * rndDirection[0] +
-                  basis[1][1] * rndDirection[1] + basis[2][1] * rndDirection[2];
-    prd->dir[2] = basis[0][2] * rndDirection[0] +
-                  basis[1][2] * rndDirection[1] + basis[2][2] * rndDirection[2];
+    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);
+
+    float rx = cosTheta;
+    float ry = cosPhi * sinTheta;
+    float rz = sinPhi * sinTheta;
+
+    prd->dir = basis[0] * rx + basis[1] * ry + basis[2] * rz;
   } while (prd->dir[2] >= 0.f);
 
   if (D == 2)
@@ -102,14 +100,13 @@ initializeRayPositionAndDirection(PerRayData *prd,
                        launchParams->source.minPoint[1]);
   prd->pos[2] = launchParams->source.planeHeight;
 
-  const float tt = powf(u.w, 2.f / (launchParams->cosineExponent + 1.f));
+  const float tt = powf(u.w, 1.f / (launchParams->cosineExponent + 1.f));
   float s, c;
   __sincosf(2.f * M_PIf * u.z, &s, &c);
-  const float sqrt1mtt = sqrtf(1 - tt);
+  const float sqrt1mtt = sqrtf(1.f - tt * tt);
   prd->dir[0] = c * sqrt1mtt;
   prd->dir[1] = s * sqrt1mtt;
-  prd->dir[2] = -1.f * sqrtf(tt);
-  Normalize(prd->dir);
+  prd->dir[2] = -tt;
 }
 } // namespace viennaray::gpu
 #endif

gpu/include/raygTrace.hpp

@@ -39,11 +39,9 @@ template <class T, int D> class Trace {
       : geometryType_(std::move(geometryType)) {
     context_ = DeviceContext::getContextFromRegistry(deviceID);
     if (!context_) {
-      Logger::getInstance()
-          .addError("No context found for device ID " +
-                    std::to_string(deviceID) +
-                    ". Create and register a context first.")
-          .print();
+      VIENNACORE_LOG_ERROR("No context found for device ID " +
+                           std::to_string(deviceID) +
+                           ". Create and register a context first.");
     }
     initRayTracer();
   }
@@ -72,9 +70,7 @@ template <class T, int D> class Trace {
     callableFile_ = path / finalName;
 
     if (!std::filesystem::exists(callableFile_)) {
-      Logger::getInstance()
-          .addError("Callable file " + finalName + " not found.")
-          .print();
+      VIENNACORE_LOG_ERROR("Callable file " + finalName + " not found.");
     }
   }
 
@@ -84,16 +80,14 @@ template <class T, int D> class Trace {
 
   void apply() {
     if (particles_.empty()) {
-      Logger::getInstance()
-          .addError("No particles inserted. Use insertNextParticle first.")
-          .print();
+      VIENNACORE_LOG_ERROR(
+          "No particles inserted. Use insertNextParticle first.");
     }
 
     if (cellDataBuffer_.sizeInBytes / sizeof(float) !=
         numCellData * launchParams.numElements) {
-      Logger::getInstance()
-          .addError("Cell data buffer size does not match the expected size.")
-          .print();
+      VIENNACORE_LOG_ERROR(
+          "Cell data buffer size does not match the expected size.");
     }
 
     // Resize our cuda result buffer
@@ -127,26 +121,21 @@ template <class T, int D> class Trace {
 
     numRays = numPointsPerDim * numPointsPerDim * config_.numRaysPerPoint;
     if (numRays > (1 << 29)) {
-      Logger::getInstance()
-          .addWarning("Too many rays for single launch: " +
-                      util::prettyDouble(numRays))
-          .print();
+      VIENNACORE_LOG_WARNING("Too many rays for single launch: " +
+                             util::prettyDouble(numRays));
       config_.numRaysPerPoint = (1 << 29) / (numPointsPerDim * numPointsPerDim);
       numRays = numPointsPerDim * numPointsPerDim * config_.numRaysPerPoint;
     }
-    Logger::getInstance()
-        .addDebug("Number of rays: " + util::prettyDouble(numRays))
-        .print();
+    VIENNACORE_LOG_DEBUG("Number of rays: " + util::prettyDouble(numRays));
 
     // set up material specific sticking probabilities
     materialStickingBuffer_.resize(particles_.size());
     for (size_t i = 0; i < particles_.size(); i++) {
       if (!particles_[i].materialSticking.empty()) {
         if (uniqueMaterialIds_.empty() || materialIdsBuffer_.sizeInBytes == 0) {
-          Logger::getInstance()
-              .addError("Material IDs not set, when using material dependent "
-                        "sticking.")
-              .print();
+          VIENNACORE_LOG_ERROR(
+              "Material IDs not set, when using material dependent "
+              "sticking.");
         }
         std::vector<float> materialSticking(uniqueMaterialIds_.size());
         unsigned currentId = 0;
@@ -169,17 +158,13 @@ template <class T, int D> class Trace {
     launchParamsBuffers.resize(particles_.size());
 
     if (particleMap_.empty()) {
-      Logger::getInstance()
-          .addError("No particle name->particleType mapping provided.")
-          .print();
+      VIENNACORE_LOG_ERROR("No particle name->particleType mapping provided.");
     }
 
     for (size_t i = 0; i < particles_.size(); i++) {
       auto it = particleMap_.find(particles_[i].name);
       if (it == particleMap_.end()) {
-        Logger::getInstance()
-            .addError("Unknown particle name: " + particles_[i].name)
-            .print();
+        VIENNACORE_LOG_ERROR("Unknown particle name: " + particles_[i].name);
       }
       launchParams.particleType = it->second;
       launchParams.particleIdx = static_cast<unsigned>(i);
@@ -244,9 +229,8 @@ template <class T, int D> class Trace {
                       const unsigned numData) {
     if (passedCellDataBuffer.sizeInBytes / sizeof(float) / numData !=
         launchParams.numElements) {
-      Logger::getInstance()
-          .addWarning("Passed cell data does not match number of elements.")
-          .print();
+      VIENNACORE_LOG_WARNING(
+          "Passed cell data does not match number of elements.");
     }
     cellDataBuffer_ = passedCellDataBuffer;
 #ifndef NDEBUG
@@ -402,7 +386,7 @@ template <class T, int D> class Trace {
 
   unsigned int prepareParticlePrograms() {
     if (particles_.empty()) {
-      Logger::getInstance().addWarning("No particles defined.").print();
+      VIENNACORE_LOG_WARNING("No particles defined.");
       return 0;
     }
 
@@ -422,9 +406,8 @@ template <class T, int D> class Trace {
     dataPerParticleBuffer_.allocUpload(dataPerParticle);
     launchParams.dataPerParticle =
         (unsigned int *)dataPerParticleBuffer_.dPointer();
-    Logger::getInstance()
-        .addDebug("Number of flux arrays: " + std::to_string(numFluxes_))
-        .print();
+    VIENNACORE_LOG_DEBUG("Number of flux arrays: " +
+                         std::to_string(numFluxes_));
 
     return numFluxes_;
   }
@@ -503,16 +486,14 @@ template <class T, int D> class Trace {
     std::string pipelineFile = pipelineFileName + geometryType_ + ".optixir";
     std::filesystem::path pipelinePath = context_->modulePath / pipelineFile;
     if (!std::filesystem::exists(pipelinePath)) {
-      Logger::getInstance()
-          .addError("Pipeline file " + pipelinePath.string() + " not found.")
-          .print();
+      VIENNACORE_LOG_ERROR("Pipeline file " + pipelinePath.string() +
+                           " not found.");
     }
 
     auto pipelineInput = getInputData(pipelinePath.c_str(), inputSize);
     if (!pipelineInput) {
-      Logger::getInstance()
-          .addError("Pipeline file " + pipelinePath.string() + " not found.")
-          .print();
+      VIENNACORE_LOG_ERROR("Pipeline file " + pipelinePath.string() +
+                           " not found.");
     }
 
     OPTIX_CHECK(optixModuleCreate(context_->optix, &moduleCompileOptions_,
@@ -525,14 +506,13 @@ template <class T, int D> class Trace {
     size_t sizeof_log_callable = sizeof(logCallable);
 
     if (callableFile_.empty()) {
-      Logger::getInstance().addWarning("No callable file set.").print();
+      VIENNACORE_LOG_WARNING("No callable file set.");
       return;
     }
     auto callableInput = getInputData(callableFile_.c_str(), inputSize);
     if (!callableInput) {
-      Logger::getInstance()
-          .addError("Callable file " + callableFile_.string() + " not found.")
-          .print();
+      VIENNACORE_LOG_ERROR("Callable file " + callableFile_.string() +
+                           " not found.");
     }
 
     OPTIX_CHECK(optixModuleCreate(context_->optix, &moduleCompileOptions_,
@@ -606,9 +586,7 @@ template <class T, int D> class Trace {
   /// does all setup for the direct callables
   void createDirectCallablePrograms() {
     if (callableMap_.empty()) {
-      Logger::getInstance()
-          .addWarning("No particleType->callable mapping provided.")
-          .print();
+      VIENNACORE_LOG_WARNING("No particleType->callable mapping provided.");
       return;
     }
     unsigned maxParticleTypeId = 0;

gpu/include/raygTraceDisk.hpp

@@ -22,9 +22,7 @@ template <class T, int D> class TraceDisk final : public Trace<T, D> {
     assert(context_ && "Context not initialized.");
     diskMesh = passedMesh;
     if (diskMesh.gridDelta <= 0.f) {
-      Logger::getInstance()
-          .addError("DiskMesh gridDelta must be positive and non-zero.")
-          .print();
+      VIENNACORE_LOG_ERROR("DiskMesh gridDelta must be positive and non-zero.");
     }
     if (diskMesh.radius <= 0.f) {
       diskMesh.radius = rayInternal::DiskFactor<3> * diskMesh.gridDelta;

gpu/pipelines/CallableWrapper.cu

@@ -33,6 +33,11 @@ __direct_callable__particleReflection(const void *sbtData,
   particleReflection(sbtData, prd);
 }
 
+extern "C" __device__ void __direct_callable__particleReflectionConstSticking(
+    const void *sbtData, viennaray::gpu::PerRayData *prd) {
+  particleReflectionConstSticking(sbtData, prd);
+}
+
 extern "C" __device__ void
 __direct_callable__particleInit(const void *sbtData,
                                 viennaray::gpu::PerRayData *prd) {

gpu/pipelines/Particle.cuh

@@ -30,6 +30,14 @@ particleReflection(const void *sbtData, viennaray::gpu::PerRayData *prd) {
   viennaray::gpu::diffuseReflection(prd, geoNormal, launchParams.D);
 }
 
+__forceinline__ __device__ void
+particleReflectionConstSticking(const void *sbtData,
+                                viennaray::gpu::PerRayData *prd) {
+  prd->rayWeight -= prd->rayWeight * launchParams.sticking;
+  auto geoNormal = viennaray::gpu::computeNormal(sbtData, prd->primID);
+  viennaray::gpu::diffuseReflection(prd, geoNormal, launchParams.D);
+}
+
 __forceinline__ __device__ void particleInit(viennaray::gpu::PerRayData *prd) {
   // Optional
 }

gpu/tests/rngSeed/CMakeLists.txt

@@ -0,0 +1,9 @@
+project(GPU_rngSeed LANGUAGES CXX)
+
+add_gpu_executable(${PROJECT_NAME} target_name rngSeed.cpp)
+
+add_dependencies(ViennaRay-GPU_Tests ${PROJECT_NAME})
+add_test(NAME ${PROJECT_NAME} COMMAND $<TARGET_FILE:${PROJECT_NAME}>)
+
+configure_file(${CMAKE_SOURCE_DIR}/examples/triangle3D/trenchMesh.dat
+               ${CMAKE_CURRENT_BINARY_DIR}/trenchMesh.dat COPYONLY)