diff --git a/CSGOptiX/CMakeLists.txt b/CSGOptiX/CMakeLists.txt index d68ac98f7..4d0ea3540 100644 --- a/CSGOptiX/CMakeLists.txt +++ b/CSGOptiX/CMakeLists.txt @@ -22,6 +22,7 @@ list(APPEND SOURCES list(APPEND HEADERS CSGOptiX.h CSGOptiXService.h + CSGOptiXHelpers.h CSGOPTIX_API_EXPORT.hh CSGOPTIX_LOG.hh ) diff --git a/CSGOptiX/CSGOptiX.cc b/CSGOptiX/CSGOptiX.cc index 49c942e04..e95166a34 100644 --- a/CSGOptiX/CSGOptiX.cc +++ b/CSGOptiX/CSGOptiX.cc @@ -592,25 +592,14 @@ void CSGOptiX::initSimulate() params->PropagateEpsilon0Mask = SEventConfig::PropagateEpsilon0Mask(); // eg from CK|SI|TO|SC|RE params->PropagateRefine = SEventConfig::PropagateRefine(); - params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace - - params->tmin = SEventConfig::PropagateEpsilon() ; // eg 0.1 0.05 to avoid self-intersection off boundaries - params->tmax = 1000000.f ; - params->max_time = SEventConfig::MaxTime() ; + params->PropagateRefineDistance = SEventConfig::PropagateRefineDistance(); // approx distance beyond which to refine intersect with 2nd trace + params->tmin = SEventConfig::PropagateEpsilon(); // eg 0.1 0.05 to avoid self-intersection off boundaries + params->tmax = 1000000.f; + params->max_time = SEventConfig::MaxTime(); } - - - - - - - - - - /** CSGOptiX::initRender -------------------------- diff --git a/CSGOptiX/CSGOptiX7.cu b/CSGOptiX/CSGOptiX7.cu index 5821cbaf7..5c2d93b9e 100644 --- a/CSGOptiX/CSGOptiX7.cu +++ b/CSGOptiX/CSGOptiX7.cu @@ -77,6 +77,7 @@ __intersection__is #include "Binding.h" #include "Params.h" +#include "CSGOptiXHelpers.h" #ifdef WITH_PRD #include "scuda_pointer.h" @@ -439,6 +440,12 @@ static __forceinline__ __device__ void simulate( const uint3& launch_idx, const #endif sim->generate_photon(ctx.p, rng, gs, photon_idx, genstep_idx ); + // Only Cerenkov and scintillation gensteps store a source-material line in + // q0.u.z. TORCH and FRAME use that slot for a genstep id, so leave the + // carried line invalid for all other genstep types. + const unsigned gentype = gs.q0.u.x; + const bool carries_matline = CSGOptiX7_GenstepCarriesMaterialLine(gentype); + ctx.current_matline = carries_matline ? gs.q0.u.z : 0xFFFFFFFFu; FlowAction command = FlowAction::Start; int bounce = 0 ; @@ -902,8 +909,16 @@ extern "C" __global__ void __intersection__is() const unsigned boundary = node->boundary() ; // all CSGNode in the tree for one CSGPrim tree have same boundary const unsigned globalPrimIdx_boundary = (( globalPrimIdx & 0xffffu ) << 16 ) | ( boundary & 0xffffu ) ; + // Coincident sibling faces can produce equal entering and exiting hits. + // Move only the reported exiting distance to the next representable float + // during photon simulation so the entering sibling wins an exact tie. + // This one-ULP ordering change cannot overstep a representable physical + // gap. Keep the unbiased isect.w in the PRD for position advancement. + const float cosI = dot(ray_direction, make_float3(isect.x, isect.y, isect.z)); + const float t_report = CSGOptiX7_ReportedIntersectionDistance(isect.w, params.raygenmode, cosI); + #ifdef WITH_PRD - if(optixReportIntersection( isect.w, hitKind)) + if (optixReportIntersection(t_report, hitKind)) { quad2* prd = SOPTIX_getPRD(); // access prd addr from RG program prd->q0.f = isect ; // .w:distance and .xyz:normal which starts as the local frame one @@ -919,7 +934,7 @@ extern "C" __global__ void __intersection__is() a3 = __float_as_uint( isect.w ) ; a4 = globalPrimIdx_boundary ; a5 = __float_as_uint( lposcost ); - optixReportIntersection( isect.w, hitKind, a0, a1, a2, a3, a4, a5 ); + optixReportIntersection(t_report, hitKind, a0, a1, a2, a3, a4, a5); // IS:optixReportIntersection writes the attributes that can be read in CH and AH programs // max 8 attribute registers, see PIP::PIP, communicate to __closesthit__ch diff --git a/CSGOptiX/CSGOptiXHelpers.h b/CSGOptiX/CSGOptiXHelpers.h new file mode 100644 index 000000000..cc943db7b --- /dev/null +++ b/CSGOptiX/CSGOptiXHelpers.h @@ -0,0 +1,30 @@ +#pragma once + +#include + +#include "OpticksGenstep.h" +#include "SRG.h" + +#if defined(__CUDACC__) || defined(__CUDABE__) +#define CSGOPTIX7_HELPER_METHOD __device__ __forceinline__ +#else +#define CSGOPTIX7_HELPER_METHOD inline +#endif + +CSGOPTIX7_HELPER_METHOD bool CSGOptiX7_GenstepCarriesMaterialLine(unsigned gentype) +{ + return gentype == OpticksGenstep_CERENKOV || + gentype == OpticksGenstep_SCINTILLATION || + gentype == OpticksGenstep_G4Cerenkov_modified; +} + +CSGOPTIX7_HELPER_METHOD float CSGOptiX7_ReportedIntersectionDistance( + float distance, + unsigned raygenmode, + float direction_dot_normal) +{ + const bool defer_exit = raygenmode == SRG_SIMULATE && direction_dot_normal > 0.f; + return defer_exit ? nextafterf(distance, INFINITY) : distance; +} + +#undef CSGOPTIX7_HELPER_METHOD diff --git a/CSGOptiX/Params.cc b/CSGOptiX/Params.cc index cc3a8b7c5..3d388a2fe 100644 --- a/CSGOptiX/Params.cc +++ b/CSGOptiX/Params.cc @@ -123,12 +123,10 @@ std::string Params::detail() const << std::setw(20) << " evt " << std::setw(10) << evt << std::endl ; std::string s = ss.str(); - return s ; + return s; } - -Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) - : +Params::Params(int raygenmode_, unsigned width, unsigned height, unsigned depth) : raygenmode(SRG_RENDER), node(nullptr), plan(nullptr), diff --git a/CSGOptiX/Params.h b/CSGOptiX/Params.h index 54a93faea..7b927e1fa 100644 --- a/CSGOptiX/Params.h +++ b/CSGOptiX/Params.h @@ -61,12 +61,12 @@ struct Params float3 V ; float3 W ; float3 WNORM ; - float tmin ; + float tmin; - float tmin0 ; - unsigned PropagateEpsilon0Mask ; // default from SEventConfig TO,CK,SI,SC,RE - float PropagateRefineDistance ; - bool PropagateRefine ; + float tmin0; + unsigned PropagateEpsilon0Mask; // default from SEventConfig TO,CK,SI,SC,RE + float PropagateRefineDistance; + bool PropagateRefine; float tmax ; float4 ZPROJ ; diff --git a/CSGOptiX/tests/CMakeLists.txt b/CSGOptiX/tests/CMakeLists.txt index d5dbfe09b..3789957c1 100644 --- a/CSGOptiX/tests/CMakeLists.txt +++ b/CSGOptiX/tests/CMakeLists.txt @@ -10,6 +10,7 @@ or doing geometry specific things. set(SINGLE_SOURCES + CSGOptiXHelpersTest.cu CSGOptiXDescTest.cc CSGOptiXSimtraceTest.cc CSGOptiXSimulateTest.cc @@ -49,4 +50,6 @@ foreach(SRC ${TEST_SOURCES}) ) endforeach() +add_test(NAME ${name}.CSGOptiXHelpersTest COMMAND $) + set_tests_properties(${name}.CSGOptiXRenderTest PROPERTIES DISABLED TRUE) diff --git a/CSGOptiX/tests/CSGOptiXHelpersTest.cu b/CSGOptiX/tests/CSGOptiXHelpersTest.cu new file mode 100644 index 000000000..5e28e15ce --- /dev/null +++ b/CSGOptiX/tests/CSGOptiXHelpersTest.cu @@ -0,0 +1,92 @@ +#include +#include + +#include + +#include "CSGOptiXHelpers.h" + +namespace +{ + void require(bool condition, const char* expression, int line) + { + if (condition) return; + std::fprintf(stderr, "%s:%d requirement failed: %s\n", __FILE__, line, expression); + std::exit(EXIT_FAILURE); + } +} + +#define REQUIRE(expression) require((expression), #expression, __LINE__) + +struct CSGOptiXHelpersTestResult +{ + bool torch_carries_matline; + bool frame_carries_matline; + bool cerenkov_carries_matline; + bool scintillation_carries_matline; + bool modified_cerenkov_carries_matline; + float entering_distance; + float exiting_distance; + float non_simulation_exit_distance; +}; + +__global__ void CSGOptiXHelpersTest_kernel(CSGOptiXHelpersTestResult* result) +{ + result->torch_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_TORCH); + result->frame_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_FRAME); + result->cerenkov_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_CERENKOV); + result->scintillation_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_SCINTILLATION); + result->modified_cerenkov_carries_matline = + CSGOptiX7_GenstepCarriesMaterialLine(OpticksGenstep_G4Cerenkov_modified); + + const float coincident_distance = 10.f; + result->entering_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMULATE, + -1.f); + result->exiting_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMULATE, + 1.f); + result->non_simulation_exit_distance = CSGOptiX7_ReportedIntersectionDistance( + coincident_distance, + SRG_SIMTRACE, + 1.f); +} + +int main() +{ + CSGOptiXHelpersTestResult* device_result = nullptr; + REQUIRE(cudaMalloc(&device_result, sizeof(CSGOptiXHelpersTestResult)) == cudaSuccess); + + CSGOptiXHelpersTest_kernel<<<1, 1>>>(device_result); + REQUIRE(cudaGetLastError() == cudaSuccess); + REQUIRE(cudaDeviceSynchronize() == cudaSuccess); + + CSGOptiXHelpersTestResult result = {}; + REQUIRE(cudaMemcpy( + &result, + device_result, + sizeof(CSGOptiXHelpersTestResult), + cudaMemcpyDeviceToHost) == cudaSuccess); + REQUIRE(cudaFree(device_result) == cudaSuccess); + + REQUIRE(!result.torch_carries_matline); + REQUIRE(!result.frame_carries_matline); + REQUIRE(result.cerenkov_carries_matline); + REQUIRE(result.scintillation_carries_matline); + REQUIRE(result.modified_cerenkov_carries_matline); + + REQUIRE(result.entering_distance == 10.f); + REQUIRE(result.exiting_distance == nextafterf(10.f, INFINITY)); + REQUIRE(result.entering_distance < result.exiting_distance); + REQUIRE(result.non_simulation_exit_distance == 10.f); + + std::puts("CSGOptiXHelpersTest: PASS"); + return 0; +} + +#undef REQUIRE diff --git a/config/sibling_pair.json b/config/sibling_pair.json new file mode 100644 index 000000000..97ca98e8e --- /dev/null +++ b/config/sibling_pair.json @@ -0,0 +1,34 @@ +{ + "torch": { + "gentype": "TORCH", + "trackid": 0, + "matline": 0, + "numphoton": 2000, + + "pos": [-10.0, 0.0, -50.1], + "time": 0.0, + + "mom": [0.3, 0.0, 1.0], + "weight": 0.0, + + "pol": [0.0, 1.0, 0.0], + "wavelength": 420.0, + + "zenith": [-1.0, 1.0], + "azimuth": [0.25, 0.25], + + "radius": 0.05, + "distance": 0.0, + "mode": 255, + "type": "disc" + }, + + "event": { + "mode": "HitPhoton", + "maxslot": 10000, + "max_bounce": 8, + "propagate_epsilon": 0.05, + "propagate_epsilon0": 0.05, + "propagate_epsilon0_mask": "TO,CK,SI,SC,RE" + } +} diff --git a/qudarap/qbnd.h b/qudarap/qbnd.h index d6c5bdd9d..8dd2597fa 100644 --- a/qudarap/qbnd.h +++ b/qudarap/qbnd.h @@ -35,16 +35,13 @@ struct qbnd unsigned boundary_tex_MaterialLine_LS ; quad* optical ; -#if defined(__CUDACC__) || defined(__CUDABE__) || defined( MOCK_TEXTURE) || defined(MOCK_CUDA) - QBND_METHOD float4 boundary_lookup( unsigned ix, unsigned iy ); - QBND_METHOD float4 boundary_lookup( float nm, unsigned line, unsigned k ); - QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx ); +#if defined(__CUDACC__) || defined(__CUDABE__) || defined(MOCK_TEXTURE) || defined(MOCK_CUDA) + QBND_METHOD float4 boundary_lookup(unsigned ix, unsigned iy); + QBND_METHOD float4 boundary_lookup(float nm, unsigned line, unsigned k); + QBND_METHOD void fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline = 0xFFFFFFFFu); #endif - }; - - #if defined( MOCK_TEXTURE) || defined(MOCK_CUDA) #include "stexture.h" #endif @@ -194,30 +191,51 @@ s.optical.x **/ -inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx ) +inline QBND_METHOD void qbnd::fill_state(sstate& s, unsigned boundary, float wavelength, float cosTheta, unsigned long long idx, unsigned long long base_pidx, unsigned carried_matline) { - const int line = boundary*_BOUNDARY_NUM_MATSUR ; // now that are not signing boundary use 0-based - - const int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT ; - const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT ; - const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR ; - + const int line = boundary * _BOUNDARY_NUM_MATSUR; // now that are not signing boundary use 0-based + + int m1_line = cosTheta > 0.f ? line + IMAT : line + OMAT; + const int m2_line = cosTheta > 0.f ? line + OMAT : line + IMAT; + const int su_line = cosTheta > 0.f ? line + ISUR : line + OSUR; + + // A boundary normally describes a daughter volume relative to its parent. + // When two daughter volumes touch, however, a photon can move directly from + // one sibling to the other. The boundary selected by OptiX may then name the + // common parent instead of the material the photon is actually leaving. + // + // qsim carries the material-table row for the photon's current volume across + // interactions. Use that row as material1 when it is valid so absorption and + // scattering are sampled from the correct medium. Material2 and the surface + // still come from the selected boundary. + // + // UINT_MAX means that no current-material row is available. Line zero is a + // valid material row. Surface rows and rows beyond the uploaded boundary + // table must not be used for property or optical-buffer lookups. + const unsigned num_matline = boundary_meta->q0.u.y / _BOUNDARY_NUM_FLOAT4; + const unsigned matline_slot = carried_matline & 0x3u; + const bool carried_matline_valid = carried_matline < num_matline && (matline_slot == unsigned(OMAT) || matline_slot == unsigned(IMAT)); + + if (carried_matline_valid) + { + m1_line = int(carried_matline); + } - s.material1 = boundary_lookup( wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob - s.m1group2 = boundary_lookup( wavelength, m1_line, 1); // x: material1 group_velocity, y: material2 group_velocity, z/w unused - s.material2 = boundary_lookup( wavelength, m2_line, 0); // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used - s.surface = boundary_lookup( wavelength, su_line, 0); // detect, , absorb , (reflect_specular), reflect_diffuse [they add to 1. so one not used] + s.material1 = boundary_lookup(wavelength, m1_line, 0); // refractive_index, absorption_length, scattering_length, reemission_prob + s.m1group2 = boundary_lookup(wavelength, m1_line, 1); // x: material1 group_velocity, y: material2 group_velocity, z/w unused + s.material2 = boundary_lookup(wavelength, m2_line, 0); // refractive_index, (absorption_length, scattering_length, reemission_prob) only m2:refractive index actually used + s.surface = boundary_lookup(wavelength, su_line, 0); // detect, , absorb , (reflect_specular), reflect_diffuse [they add to 1. so one not used] s.set_material2_group_velocity(boundary_lookup(wavelength, m2_line, 1).x); - s.optical = optical[su_line].u ; // 1-based-surface-index-0-meaning-boundary/type/finish/value (type,finish,value not used currently) + s.optical = optical[su_line].u; // 1-based-surface-index-0-meaning-boundary/type/finish/value (type,finish,value not used currently) - s.index.x = optical[m1_line].u.x ; // m1 index (1-based, see sstandard::make_optical) - s.index.y = optical[m2_line].u.x ; // m2 index (1-based, see sstandard::make_optical) - s.index.z = optical[su_line].u.x ; // su index (1-based, see sstandard::make_optical) - s.index.w = 0u ; // avoid undefined memory comparison issues + s.index.x = optical[m1_line].u.x; // m1 index (1-based, see sstandard::make_optical) -- reflects override if any + s.index.y = optical[m2_line].u.x; // m2 index (1-based, see sstandard::make_optical) + s.index.z = optical[su_line].u.x; // su index (1-based, see sstandard::make_optical) + s.index.w = 0u; // avoid undefined memory comparison issues #if !defined(PRODUCTION) && defined(DEBUG_PIDX) - if( idx == base_pidx ) + if (idx == base_pidx) { printf("//qbnd.fill_state idx %7lld boundary %d line %d wavelength %10.4f m1_line %d m2_line %d su_line %d s.optical.x %d \n", idx, boundary, line, wavelength, m1_line, m2_line, su_line, s.optical.x ); printf("//qbnd.fill_state idx %7lld boundary %d [s.index.x-1](m1_index) %d [s.index.y-1](m2_index) %d [s.index.z-1](su_index) %d \n", idx, boundary, s.index.x-1u, s.index.y-1u, s.index.z-1u ); diff --git a/qudarap/qsim.h b/qudarap/qsim.h index 1e110505b..68de1fc34 100644 --- a/qudarap/qsim.h +++ b/qudarap/qsim.h @@ -1261,8 +1261,25 @@ inline QSIM_METHOD FlowAction qsim::propagate_at_boundary(unsigned& flag, RNG& r ctx.current_material_index = reflect ? s.index.x : s.index.y; ctx.current_group_velocity = reflect ? s.material1_group_velocity() : s.material2_group_velocity(); + // Remember the material entered by a transmitted photon so the next + // boundary lookup knows the photon's actual current medium. This matters + // when touching sibling volumes share a face because the selected boundary + // may describe one sibling relative to their common parent and omit the + // material the photon is leaving. + // + // Reflection does not change the current medium, so preserve the existing + // material row in that case. For transmission, _c1 identifies which side of + // the selected boundary is the material the photon enters. + if (!reflect) + { + const unsigned bnd_idx = ctx.prd->boundary(); + const unsigned imat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + IMAT; + const unsigned omat_line = bnd_idx * _BOUNDARY_NUM_MATSUR + OMAT; + ctx.current_matline = (_c1 < 0.f) ? omat_line : imat_line; // material entered by this transmission + } + #if !defined(PRODUCTION) && defined(DEBUG_TAG) - if( flag == BOUNDARY_REFLECT ) + if (flag == BOUNDARY_REFLECT) { const float u_br_align_0 = curand_uniform(&rng) ; const float u_br_align_1 = curand_uniform(&rng) ; @@ -2147,9 +2164,11 @@ QSIM_FORCEINLINE_METHOD FlowAction qsim::propagate(const int bounce, RNG& rng, s #endif // copy geometry info into the sphoton struct - ctx.p.set_prd(boundary, identity, cosTheta, iindex ); // HMM: lposcost not passed along + ctx.p.set_prd(boundary, identity, cosTheta, iindex); // HMM: lposcost not passed along - bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx ); + // Supply the current-material row so fill_state can correct material1 when + // the selected boundary does not describe the photon's actual current medium. + bnd->fill_state(ctx.s, boundary, ctx.p.wavelength, cosTheta, ctx.pidx, base->pidx, ctx.current_matline); if (ctx.current_material_index == 0u) { diff --git a/qudarap/tests/CMakeLists.txt b/qudarap/tests/CMakeLists.txt index 5deea0e09..2a5e173e2 100644 --- a/qudarap/tests/CMakeLists.txt +++ b/qudarap/tests/CMakeLists.txt @@ -18,6 +18,7 @@ set(TEST_SOURCES QTexMakerTest.cc QTexLookupTest.cc qbnd_boundary_lookup_test.cc + QSim_MaterialCarryTest.cc QMultiFilmTest.cc @@ -91,6 +92,10 @@ foreach(SRC ${TEST_SOURCES}) add_executable(${TGT} ${SRC}) target_link_libraries(${TGT} QUDARap CUDA::cudart) + if(TGT STREQUAL QSim_MaterialCarryTest) + target_compile_definitions(${TGT} PRIVATE MOCK_CUDA MOCK_CURAND) + endif() + if(TGT STREQUAL QEvt_Lifecycle_Test) foreach(INPUT_PHOTON_NAME ${QUDARAP_INPUT_PHOTON_NAMES}) set(TEST_NAME ${name}.${TGT}.${INPUT_PHOTON_NAME}) diff --git a/qudarap/tests/QSim_MaterialCarryTest.cc b/qudarap/tests/QSim_MaterialCarryTest.cc new file mode 100644 index 000000000..34b237b34 --- /dev/null +++ b/qudarap/tests/QSim_MaterialCarryTest.cc @@ -0,0 +1,121 @@ +#include +#include + +#include "scerenkov.h" +#include "scuda.h" +#include "smath.h" +#include "srec.h" +#include "srngcpu.h" +using RNG = srngcpu; + +#include "stexture.h" + +#include "OpticksPhoton.h" +#include "qsim.h" + +namespace +{ +void require(bool condition, const char* expression, int line) +{ + if (condition) + return; + std::fprintf(stderr, "%s:%d requirement failed: %s\n", __FILE__, line, expression); + std::exit(EXIT_FAILURE); +} + +#define REQUIRE(expression) require((expression), #expression, __LINE__) + +constexpr unsigned BOUNDARY_INDEX = 2u; +constexpr unsigned OUTER_MATERIAL_LINE = + BOUNDARY_INDEX * _BOUNDARY_NUM_MATSUR + OMAT; +constexpr unsigned INNER_MATERIAL_LINE = + BOUNDARY_INDEX * _BOUNDARY_NUM_MATSUR + IMAT; + +void initialize(sctx& ctx, quad2& prd) +{ + prd.zero(); + prd.q0.f.z = 1.f; + prd.q1.u.w = BOUNDARY_INDEX; + + ctx = {}; + ctx.prd = &prd; + ctx.pidx = 0u; + ctx.p.zero(); + ctx.p.mom = make_float3(0.f, 0.f, -1.f); + ctx.p.pol = make_float3(1.f, 0.f, 0.f); + ctx.p.wavelength = 420.f; + + ctx.s.material1 = make_float4(1.f, 0.f, 0.f, 0.f); + ctx.s.material2 = make_float4(1.f, 0.f, 0.f, 0.f); + ctx.s.m1group2 = make_float4(100.f, 200.f, 0.f, 0.f); + ctx.s.index = make_uint4(10u, 20u, 0u, 0u); +} + +void test_transmission_updates_carried_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.current_matline = 99u; + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 1.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_TRANSMIT); + REQUIRE(ctx.current_matline == INNER_MATERIAL_LINE); + REQUIRE(ctx.current_material_index == 20u); + REQUIRE(ctx.current_group_velocity == 200.f); +} + +void test_reflection_preserves_carried_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.current_matline = 99u; + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 0.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_REFLECT); + REQUIRE(ctx.current_matline == 99u); + REQUIRE(ctx.current_material_index == 10u); + REQUIRE(ctx.current_group_velocity == 100.f); +} + +void test_outward_transmission_selects_outer_material(const qsim& sim, RNG& rng) +{ + sctx ctx = {}; + quad2 prd = {}; + initialize(ctx, prd); + ctx.p.mom = make_float3(0.f, 0.f, 1.f); + + unsigned flag = 0u; + const FlowAction action = sim.propagate_at_boundary(flag, rng, ctx, 1.f); + + REQUIRE(action == FlowAction::Continue); + REQUIRE(flag == BOUNDARY_TRANSMIT); + REQUIRE(ctx.current_matline == OUTER_MATERIAL_LINE); +} +} // namespace + +int main() +{ + qbase base = {}; + base.pidx = ~0ull; + + qsim sim; + sim.base = &base; + + RNG rng; + rng.set_fake(0.5); + + test_transmission_updates_carried_material(sim, rng); + test_reflection_preserves_carried_material(sim, rng); + test_outward_transmission_selects_outer_material(sim, rng); + + std::puts("QSim_MaterialCarryTest: PASS"); + return 0; +} diff --git a/qudarap/tests/qbnd_boundary_lookup_test.cc b/qudarap/tests/qbnd_boundary_lookup_test.cc index 77cb99deb..e0e134854 100644 --- a/qudarap/tests/qbnd_boundary_lookup_test.cc +++ b/qudarap/tests/qbnd_boundary_lookup_test.cc @@ -102,6 +102,10 @@ int main() qbnd qb; qb.boundary_tex = texObj; qb.boundary_meta = &meta; + quad optical[ni * nj] = {}; + for (unsigned line = 0; line < ni * nj; line++) + optical[line].u.x = 1000u + line; + qb.optical = optical; printf("qbnd_boundary_lookup_test\n"); @@ -160,6 +164,34 @@ int main() check("row0 clamp above: .x == 30 (bin 3)", feq(p.x, 30.f)); } + // --- Test 4: Carried material-line validation --- + // Boundary 1 entered from outside normally uses line 4 (OMAT). + // A valid carry overrides only material1; invalid, surface, and out-of-range + // rows must fall back to line 4 without indexing outside the buffers. + { + const unsigned boundary = 1u; + const float cosTheta = -1.f; + sstate s = {}; + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 0u); + check("carried line 0 remains a valid OMAT override", feq(s.material1.x, 0.f)); + check("carried line 0 selects optical material index", s.index.x == 1000u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 7u); + check("carried line 7 remains a valid IMAT override", feq(s.material1.x, 1400.f)); + check("carried line 7 selects optical material index", s.index.x == 1007u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 0xFFFFFFFFu); + check("invalid sentinel falls back to boundary OMAT", feq(s.material1.x, 800.f)); + check("invalid sentinel keeps boundary material index", s.index.x == 1004u); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, 5u); + check("surface row carry falls back to boundary OMAT", feq(s.material1.x, 800.f)); + + qb.fill_state(s, boundary, 100.f, cosTheta, 0u, 0u, ni * nj); + check("out-of-range carry falls back to boundary OMAT", feq(s.material1.x, 800.f)); + } + printf("qbnd_boundary_lookup_test: %s (%d failure%s)\n", g_fail == 0 ? "PASS" : "FAIL", g_fail, g_fail == 1 ? "" : "s"); return g_fail == 0 ? 0 : 1; diff --git a/sysrap/sctx.h b/sysrap/sctx.h index 431975c54..88ef6650b 100644 --- a/sysrap/sctx.h +++ b/sysrap/sctx.h @@ -76,12 +76,13 @@ struct sctx #endif sphoton p ; - sstate s ; - float current_group_velocity = 0.f; + sstate s; + float current_group_velocity = 0.f; unsigned current_material_index = 0u; + unsigned current_matline = 0xFFFFFFFFu; // sibling-pair current-medium line; UINT_MAX means unavailable #ifndef PRODUCTION - srec rec ; + srec rec; sseq seq ; stagr tagr ; quad4 aux ; diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 6ad1f41ef..725d4167f 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -29,6 +29,7 @@ add_simg4ox_integration_test(dune_mock_wls_detector_box) add_simg4ox_integration_test(8x8SiPM_w_CSI_optial_grease) add_simg4ox_integration_test(opticks_two_spheres) add_simg4ox_integration_test(drich) +add_simg4ox_integration_test(sibling_pair) set(SIMPHONY_SIMG4OX_MULTIEVENT_TEST_WORKDIR "${CMAKE_CURRENT_BINARY_DIR}/simg4ox_multievent") file(MAKE_DIRECTORY "${SIMPHONY_SIMG4OX_MULTIEVENT_TEST_WORKDIR}") diff --git a/tests/geom/sibling_pair.gdml b/tests/geom/sibling_pair.gdml new file mode 100644 index 000000000..05c01350d --- /dev/null +++ b/tests/geom/sibling_pair.gdml @@ -0,0 +1,97 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/tests/test_simg4ox.sh b/tests/test_simg4ox.sh index 648a51fed..cf33613c5 100755 --- a/tests/test_simg4ox.sh +++ b/tests/test_simg4ox.sh @@ -22,6 +22,7 @@ usage() { echo " 8x8SiPM_w_CSI_optial_grease" >&2 echo " opticks_two_spheres" >&2 echo " drich" >&2 + echo " sibling_pair" >&2 } if [[ $# -gt 1 ]]; then @@ -91,6 +92,12 @@ case "${TEST_CASE}" in # Invoke Geant4 SDs for the EFFICIENCY=1 SiPM skin surfaces. run_hit_validation drich drich "${REPO_DIR}/tests/run_validate.mac" ;; + sibling_pair) + # The beam can reach the sensor only by crossing the exact shared face + # from the left sibling into the right sibling while retaining the + # correct current-medium properties for the segment ending at that face. + run_hit_validation sibling_pair sibling_pair "${REPO_DIR}/tests/run_validate.mac" + ;; *) echo "Unknown simg4ox test case: ${TEST_CASE}" >&2 usage