Skip to content

Commit

Permalink
Change particle local data size evaluation and particle specific seed…
Browse files Browse the repository at this point in the history 
… in tracing kernel (#39)

* Change determination of particle local data size

* Fix typo

* Implement particle specific seed
  • Loading branch information
@tobre1
tobre1 authored Oct 25, 2023
1 parent 6bad1a0 commit b066b2d
Show file tree
Hide file tree
Showing 9 changed files with 65 additions and 79 deletions.
15 changes: 6 additions & 9 deletions Tests/createRay/createRay.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -33,9 +33,6 @@ int main() {
auto rng = rayRNG{};
unsigned seed = 31;
rayRNG rngstate1(seed + 0);
rayRNG rngstate2(seed + 1);
rayRNG rngstate3(seed + 2);
rayRNG rngstate4(seed + 3);

{
auto direction = rayTraceDirection::POS_Z;
Expand All @@ -49,7 +46,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_z < 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_z, (1. + 2 * gridDelta), eps)
}
Expand All @@ -67,7 +64,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_z > 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_z, (-1. - 2 * gridDelta), eps)
}
Expand All @@ -85,7 +82,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_x < 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_x, (1. + 2 * gridDelta), eps)
}
Expand All @@ -103,7 +100,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_x > 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_x, (-1. - 2 * gridDelta), eps)
}
Expand All @@ -121,7 +118,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_y < 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_y, (1. + 2 * gridDelta), eps)
}
Expand All @@ -139,7 +136,7 @@ int main() {
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < 10; ++i) {
source.fillRay(rayhit.ray, 0, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, 0, rngstate1);
RAYTEST_ASSERT(rayhit.ray.dir_y > 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_y, (-1. - 2 * gridDelta), eps)
}
Expand Down
7 changes: 2 additions & 5 deletions Tests/createSourceGrid/createSourceGrid.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,18 +26,15 @@ int main() {
auto grid = rayInternal::createSourceGrid<NumericType, D>(
boundingBox, points.size(), gridDelta, traceSettings);

rayRNG rngstate1(0);
rayRNG rngstate2(1);
rayRNG rngstate3(2);
rayRNG rngstate4(3);
rayRNG rngstate(0);
{
// build source in positive z direction;
auto source = raySourceGrid<NumericType, D>(grid, 1., traceSettings);
auto numGridPoints = source.getNumPoints();
alignas(128) auto rayhit =
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
for (size_t i = 0; i < numGridPoints; ++i) {
source.fillRay(rayhit.ray, i, rngstate1, rngstate2, rngstate3, rngstate4);
source.fillRay(rayhit.ray, i, rngstate);

RAYTEST_ASSERT(rayhit.ray.dir_z < 0.)
RAYTEST_ASSERT_ISCLOSE(rayhit.ray.org_z, (1. + 2 * gridDelta), eps)
Expand Down
14 changes: 5 additions & 9 deletions include/rayParticle.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,13 +53,12 @@ template <typename NumericType> class rayAbstractParticle {
const rayTracingData<NumericType> *globalData,
rayRNG &Rng) = 0;

/// Set the number of required data vectors for this particle to
/// collect data.
virtual int getRequiredLocalDataSize() const = 0;

/// Set the power of the cosine source distribution for this particle.
virtual NumericType getSourceDistributionPower() const = 0;

/// Set the number of required data vectors for this particle to
/// collect data. If an empty vector is returned, no local data will be
/// provided
virtual std::vector<std::string> getLocalDataLabels() const = 0;

virtual void logData(rayDataLog<NumericType> &log) = 0;
Expand Down Expand Up @@ -94,10 +93,9 @@ class rayParticle : public rayAbstractParticle<NumericType> {
const rayTracingData<NumericType> *globalData,
rayRNG &Rng) override { // collect data for this hit
}
virtual int getRequiredLocalDataSize() const override { return 0; }
virtual NumericType getSourceDistributionPower() const override { return 1.; }
virtual std::vector<std::string> getLocalDataLabels() const override {
return std::vector<std::string>(getRequiredLocalDataSize(), "localData");
return {};
}
virtual void logData(rayDataLog<NumericType> &log) override {}

Expand Down Expand Up @@ -133,12 +131,10 @@ class rayTestParticle
const rayTracingData<NumericType> *globalData,
rayRNG &Rng) override final {}

int getRequiredLocalDataSize() const override final { return 0; }

NumericType getSourceDistributionPower() const override final { return 1.; }

std::vector<std::string> getLocalDataLabels() const override final {
return std::vector<std::string>{};
return {};
}

void logData(rayDataLog<NumericType> &log) override final {}
Expand Down
21 changes: 20 additions & 1 deletion include/rayRNG.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,26 @@
#include <memory>
#include <random>

/// Use meresenne twister 19937 as random number generator.
/// Use mersenne twister 19937 as random number generator.
typedef std::mt19937_64 rayRNG;

namespace rayInternal {

// tiny encryption algortihm
template <unsigned int N>
static unsigned int tea(unsigned int val0, unsigned int val1) {
unsigned int v0 = val0;
unsigned int v1 = val1;
unsigned int s0 = 0;

for (unsigned int n = 0; n < N; n++) {
s0 += 0x9e3779b9;
v0 += ((v1 << 4) + 0xa341316c) ^ (v1 + s0) ^ ((v1 >> 5) + 0xc8013ea4);
v1 += ((v0 << 4) + 0xad90777d) ^ (v0 + s0) ^ ((v0 >> 5) + 0x7e95761e);
}

return v0;
}
} // namespace rayInternal

#endif // RAY_RNG_HPP
4 changes: 1 addition & 3 deletions include/raySource.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,9 +9,7 @@
template <typename NumericType, int D> class raySource {
public:
virtual ~raySource() {}
virtual void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState1,
rayRNG &RngState2, rayRNG &RngState3,
rayRNG &RngState4) {}
virtual void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) {}
virtual size_t getNumPoints() const { return 0; }
virtual void printIndexCounter(){};
};
Expand Down
12 changes: 5 additions & 7 deletions include/raySourceGrid.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,13 +17,11 @@ class raySourceGrid : public raySource<NumericType, D> {
ee(((NumericType)2) / (passedCosinePower + 1)),
indexCounter(sourceGrid.size(), 0) {}

void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState1,
rayRNG &RngState2, rayRNG &RngState3,
rayRNG &RngState4) override final {
void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) override final {
auto index = idx % mNumPoints;
indexCounter[index]++;
auto origin = mSourceGrid[idx % mNumPoints];
auto direction = getDirection(RngState3, RngState4);
auto direction = getDirection(RngState);

#ifdef ARCH_X86
reinterpret_cast<__m128 &>(ray) =
Expand Down Expand Up @@ -53,11 +51,11 @@ class raySourceGrid : public raySource<NumericType, D> {
}

private:
rayTriple<NumericType> getDirection(rayRNG &RngState1, rayRNG &RngState2) {
rayTriple<NumericType> getDirection(rayRNG &RngState) {
rayTriple<NumericType> direction{0., 0., 0.};
std::uniform_real_distribution<NumericType> uniDist;
auto r1 = uniDist(RngState1);
auto r2 = uniDist(RngState2);
auto r1 = uniDist(RngState);
auto r2 = uniDist(RngState);

NumericType tt = pow(r2, ee);
direction[rayDir] = posNeg * sqrtf(tt);
Expand Down
20 changes: 9 additions & 11 deletions include/raySourceRandom.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,11 +15,9 @@ class raySourceRandom : public raySource<NumericType, D> {
minMax(pTraceSettings[3]), posNeg(pTraceSettings[4]),
ee(((NumericType)2) / (pCosinePower + 1)), mNumPoints(pNumPoints) {}

void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState1,
rayRNG &RngState2, rayRNG &RngState3,
rayRNG &RngState4) override final {
auto origin = getOrigin(RngState1, RngState2);
auto direction = getDirection(RngState3, RngState4);
void fillRay(RTCRay &ray, const size_t idx, rayRNG &RngState) override final {
auto origin = getOrigin(RngState);
auto direction = getDirection(RngState);

#ifdef ARCH_X86
reinterpret_cast<__m128 &>(ray) =
Expand All @@ -43,10 +41,10 @@ class raySourceRandom : public raySource<NumericType, D> {
size_t getNumPoints() const override final { return mNumPoints; }

private:
rayTriple<NumericType> getOrigin(rayRNG &RngState1, rayRNG &RngState2) {
rayTriple<NumericType> getOrigin(rayRNG &RngState) {
rayTriple<NumericType> origin{0., 0., 0.};
std::uniform_real_distribution<NumericType> uniDist;
auto r1 = uniDist(RngState1);
auto r1 = uniDist(RngState);

origin[rayDir] = bdBox[minMax][rayDir];
origin[firstDir] =
Expand All @@ -55,19 +53,19 @@ class raySourceRandom : public raySource<NumericType, D> {
if constexpr (D == 2) {
origin[secondDir] = 0.;
} else {
auto r2 = uniDist(RngState2);
auto r2 = uniDist(RngState);
origin[secondDir] = bdBox[0][secondDir] +
(bdBox[1][secondDir] - bdBox[0][secondDir]) * r2;
}

return origin;
}

rayTriple<NumericType> getDirection(rayRNG &RngState1, rayRNG &RngState2) {
rayTriple<NumericType> getDirection(rayRNG &RngState) {
rayTriple<NumericType> direction{0., 0., 0.};
std::uniform_real_distribution<NumericType> uniDist;
auto r1 = uniDist(RngState1);
auto r2 = uniDist(RngState2);
auto r1 = uniDist(RngState);
auto r2 = uniDist(RngState);

const NumericType tt = pow(r2, ee);
direction[rayDir] = posNeg * sqrtf(tt);
Expand Down
12 changes: 6 additions & 6 deletions include/rayTrace.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -37,13 +37,13 @@ template <class NumericType, int D> class rayTrace {
boundingBox, mParticle->getSourceDistributionPower(), traceSettings,
mGeometry.getNumPoints());

auto numberOfLocalData = mParticle->getRequiredLocalDataSize();
if (numberOfLocalData) {
mLocalData.setNumberOfVectorData(numberOfLocalData);
auto localDataLabels = mParticle->getLocalDataLabels();
if (!localDataLabels.empty()) {
mLocalData.setNumberOfVectorData(localDataLabels.size());
auto numPoints = mGeometry.getNumPoints();
auto localDataLabes = mParticle->getLocalDataLabels();
for (int i = 0; i < numberOfLocalData; ++i) {
mLocalData.setVectorData(i, numPoints, 0., localDataLabes[i]);
auto localDataLabels = mParticle->getLocalDataLabels();
for (int i = 0; i < localDataLabels.size(); ++i) {
mLocalData.setVectorData(i, numPoints, 0., localDataLabels[i]);
}
}

Expand Down
39 changes: 11 additions & 28 deletions include/rayTraceKernel.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,30 +101,11 @@ template <typename NumericType, int D> class rayTraceKernel {
RTCRayHit{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};

const int threadID = omp_get_thread_num();
constexpr int numRngStates = 8;
unsigned int seeds[numRngStates];
unsigned int seed = mRunNumber;
if (mUseRandomSeeds) {
std::random_device rd;
for (size_t i = 0; i < numRngStates; ++i) {
seeds[i] = static_cast<unsigned int>(rd());
}
} else {
for (size_t i = 0; i < numRngStates; ++i) {
seeds[i] = static_cast<unsigned int>((omp_get_thread_num() + 1) * 31 +
i + mRunNumber);
}
seed = static_cast<unsigned int>(rd());
}
// It seems really important to use two separate seeds / states for
// sampling the source and sampling reflections. When we use only one
// state for both, then the variance is very high.
rayRNG RngState1(seeds[0]);
rayRNG RngState2(seeds[1]);
rayRNG RngState3(seeds[2]);
rayRNG RngState4(seeds[3]);
rayRNG RngState5(seeds[4]);
rayRNG RngState6(seeds[5]);
rayRNG RngState7(seeds[6]);
rayRNG RngState8(seeds[7]);

// thread-local particle object
auto particle = mParticle->clone();
Expand All @@ -143,12 +124,15 @@ template <typename NumericType, int D> class rayTraceKernel {

#pragma omp for schedule(dynamic)
for (long long idx = 0; idx < mNumRays; ++idx) {
particle->initNew(RngState8);
// particle specific RNG seed
auto particleSeed = rayInternal::tea<3>(idx, seed);
rayRNG RngState(particleSeed);

particle->initNew(RngState);
particle->logData(myDataLog);
NumericType rayWeight = initialRayWeight;

mSource.fillRay(rayHit.ray, idx, RngState1, RngState2, RngState3,
RngState4); // fills also tnear
mSource.fillRay(rayHit.ray, idx, RngState); // fills also tnear

#ifdef VIENNARAY_USE_RAY_MASKING
rayHit.ray.mask = -1;
Expand Down Expand Up @@ -266,14 +250,13 @@ template <typename NumericType, int D> class rayTraceKernel {
#endif
particle->surfaceCollision(distRayWeight, rayDir, normal,
hitDiskIds[diskId], matID, myLocalData,
globalData, RngState5);
globalData, RngState);
}

// get sticking probability and reflection direction
const auto stickingnDirection = particle->surfaceReflection(
rayWeight, rayDir, geomNormal, rayHit.hit.primID,
mGeometry.getMaterialId(rayHit.hit.primID), globalData,
RngState5);
mGeometry.getMaterialId(rayHit.hit.primID), globalData, RngState);
auto valueToDrop = rayWeight * stickingnDirection.first;

if (calcFlux) {
Expand All @@ -293,7 +276,7 @@ template <typename NumericType, int D> class rayTraceKernel {
if (rayWeight <= 0) {
break;
}
reflect = rejectionControl(rayWeight, initialRayWeight, RngState6);
reflect = rejectionControl(rayWeight, initialRayWeight, RngState);
if (!reflect) {
break;
}
Expand Down

0 comments on commit b066b2d

Please sign in to comment.