diff --git a/app/celer-g4/EventAction.cc b/app/celer-g4/EventAction.cc index d87b739add..bd56cd7fbb 100644 --- a/app/celer-g4/EventAction.cc +++ b/app/celer-g4/EventAction.cc @@ -49,9 +49,9 @@ void EventAction::BeginOfEventAction(G4Event const* event) if (SharedParams::CeleritasDisabled()) return; - // Set event ID in local transporter + // Set event ID in local transporter and reseed Celerits RNG ExceptionConverter call_g4exception{"celer0002"}; - CELER_TRY_HANDLE(transport_->SetEventId(event->GetEventID()), + CELER_TRY_HANDLE(transport_->InitializeEvent(event->GetEventID()), call_g4exception); } diff --git a/app/celer-g4/GeantDiagnostics.cc b/app/celer-g4/GeantDiagnostics.cc index 42706d9058..79d08a561f 100644 --- a/app/celer-g4/GeantDiagnostics.cc +++ b/app/celer-g4/GeantDiagnostics.cc @@ -13,7 +13,9 @@ #include "corecel/sys/MemRegistry.hh" #include "corecel/sys/MultiExceptionHandler.hh" #include "celeritas/Types.hh" +#include "celeritas/global/ActionRegistry.hh" #include "celeritas/global/CoreParams.hh" +#include "celeritas/user/StepDiagnostic.hh" #include "GlobalSetup.hh" @@ -57,9 +59,22 @@ GeantDiagnostics::GeantDiagnostics(SharedParams const& params) if (global_setup.StepDiagnostic()) { // Create the track step diagnostic and add to output registry - step_diagnostic_ = std::make_shared( - global_setup.GetStepDiagnosticBins(), num_threads); + auto num_bins = GlobalSetup::Instance()->GetStepDiagnosticBins(); + step_diagnostic_ + = std::make_shared(num_bins, num_threads); output_reg->insert(step_diagnostic_); + + // Add the Celeritas step diagnostic if Celeritas offloading is enabled + if (params) + { + auto step_diagnostic = std::make_shared( + params.Params()->action_reg()->next_id(), + params.Params()->particle(), + num_bins, + num_threads); + params.Params()->action_reg()->insert(step_diagnostic); + output_reg->insert(step_diagnostic); + } } if (!params) diff --git a/src/accel/LocalTransporter.cc b/src/accel/LocalTransporter.cc index 9487ee951b..2c3539c525 100644 --- a/src/accel/LocalTransporter.cc +++ b/src/accel/LocalTransporter.cc @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -114,14 +115,24 @@ LocalTransporter::LocalTransporter(SetupOptions const& options, //---------------------------------------------------------------------------// /*! - * Set the event ID at the start of an event. + * Set the event ID and reseed the Celeritas RNG at the start of an event. */ -void LocalTransporter::SetEventId(int id) +void LocalTransporter::InitializeEvent(int id) { CELER_EXPECT(*this); CELER_EXPECT(id >= 0); + event_id_ = EventId(id); track_counter_ = 0; + + if (!(G4Threading::IsMultithreadedApplication() + && G4MTRunManager::SeedOncePerCommunication())) + { + // Since Geant4 schedules events dynamically, reseed the Celeritas RNGs + // using the Geant4 event ID for reproducibility. This guarantees that + // an event can be reproduced given the event ID. + step_->reseed(event_id_); + } } //---------------------------------------------------------------------------// diff --git a/src/accel/LocalTransporter.hh b/src/accel/LocalTransporter.hh index 21abbee2f9..393a4e2879 100644 --- a/src/accel/LocalTransporter.hh +++ b/src/accel/LocalTransporter.hh @@ -63,8 +63,11 @@ class LocalTransporter inline void Initialize(SetupOptions const& options, SharedParams const& params); - // Set the event ID - void SetEventId(int); + // Set the event ID and reseed the Celeritas RNG (remove in v1.0) + [[deprecated]] void SetEventId(int id) { this->InitializeEvent(id); } + + // Set the event ID and reseed the Celeritas RNG at the start of an event + void InitializeEvent(int); // Offload this track void Push(G4Track const&); diff --git a/src/accel/SimpleOffload.cc b/src/accel/SimpleOffload.cc index 255a8acc0f..be2936712a 100644 --- a/src/accel/SimpleOffload.cc +++ b/src/accel/SimpleOffload.cc @@ -87,16 +87,17 @@ void SimpleOffload::BeginOfRunAction(G4Run const*) //---------------------------------------------------------------------------// /*! - * Send Celeritas the event ID. + * Send Celeritas the event ID and reseed the Celeritas RNG. */ void SimpleOffload::BeginOfEventAction(G4Event const* event) { if (!*this) return; - // Set event ID in local transporter + // Set event ID in local transporter and reseed RNG for reproducibility ExceptionConverter call_g4exception{"celer0002"}; - CELER_TRY_HANDLE(local_->SetEventId(event->GetEventID()), call_g4exception); + CELER_TRY_HANDLE(local_->InitializeEvent(event->GetEventID()), + call_g4exception); } //---------------------------------------------------------------------------// diff --git a/src/celeritas/CMakeLists.txt b/src/celeritas/CMakeLists.txt index db11b869b3..4d867ab513 100644 --- a/src/celeritas/CMakeLists.txt +++ b/src/celeritas/CMakeLists.txt @@ -289,6 +289,7 @@ celeritas_polysource(global/alongstep/AlongStepUniformMscAction) celeritas_polysource(global/alongstep/AlongStepRZMapFieldMscAction) celeritas_polysource(phys/detail/DiscreteSelectAction) celeritas_polysource(phys/detail/PreStepAction) +celeritas_polysource(random/RngReseed) celeritas_polysource(random/detail/CuHipRngStateInit) celeritas_polysource(track/detail/TrackInitAlgorithms) celeritas_polysource(track/detail/TrackSortUtils) diff --git a/src/celeritas/global/Stepper.cc b/src/celeritas/global/Stepper.cc index 70a19b0765..951aeec9de 100644 --- a/src/celeritas/global/Stepper.cc +++ b/src/celeritas/global/Stepper.cc @@ -14,6 +14,8 @@ #include "corecel/sys/ScopedProfiling.hh" #include "orange/OrangeData.hh" #include "celeritas/Types.hh" +#include "celeritas/random/RngParams.hh" +#include "celeritas/random/RngReseed.hh" #include "celeritas/track/TrackInitParams.hh" #include "CoreParams.hh" @@ -103,6 +105,20 @@ auto Stepper::operator()(SpanConstPrimary primaries) -> result_type return (*this)(); } +//---------------------------------------------------------------------------// +/*! + * Reseed the RNGs at the start of an event for "strong" reproducibility. + * + * This reinitializes the RNG states using a single seed and unique subsequence + * for each thread. It ensures that given an event number, that event can be + * reproduced. + */ +template +void Stepper::reseed(EventId event_id) +{ + reseed_rng(get_ref(*params_->rng()), state_.ref().rng, event_id.get()); +} + //---------------------------------------------------------------------------// // EXPLICIT INSTANTIATION //---------------------------------------------------------------------------// diff --git a/src/celeritas/global/Stepper.hh b/src/celeritas/global/Stepper.hh index 0ba59f7fb2..a2cc0ff6ac 100644 --- a/src/celeritas/global/Stepper.hh +++ b/src/celeritas/global/Stepper.hh @@ -90,6 +90,9 @@ class StepperInterface // Transport existing states and these new primaries virtual StepperResult operator()(SpanConstPrimary primaries) = 0; + // Reseed the RNGs at the start of an event for reproducibility + virtual void reseed(EventId event_id) = 0; + //! Get action sequence for timing diagnostics virtual ActionSequence const& actions() const = 0; @@ -139,6 +142,9 @@ class Stepper final : public StepperInterface // Transport existing states and these new primaries StepperResult operator()(SpanConstPrimary primaries) final; + // Reseed the RNGs at the start of an event for reproducibility + void reseed(EventId event_id) final; + //! Get action sequence for timing diagnostics ActionSequence const& actions() const final { return *actions_; } diff --git a/src/celeritas/random/RngReseed.cc b/src/celeritas/random/RngReseed.cc new file mode 100644 index 0000000000..fe2ed3b7aa --- /dev/null +++ b/src/celeritas/random/RngReseed.cc @@ -0,0 +1,40 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/random/RngReseed.cc +//---------------------------------------------------------------------------// +#include "RngReseed.hh" + +#include "corecel/cont/Range.hh" +#include "corecel/sys/ThreadId.hh" + +#include "RngEngine.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +/*! + * Reinitialize the RNG states on host at the start of an event. + * + * Each thread's state is initialized using same seed and skipped ahead a + * different number of subsequences so the sequences on different threads will + * not have statistically correlated values. + */ +void reseed_rng(HostCRef const& params, + HostRef const& state, + size_type event_id) +{ + for (auto tid : range(TrackSlotId{state.size()})) + { + RngEngine::Initializer_t init; + init.seed = params.seed; + init.subsequence = event_id * state.size() + tid.get(); + RngEngine engine(params, state, tid); + engine = init; + } +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/random/RngReseed.cu b/src/celeritas/random/RngReseed.cu new file mode 100644 index 0000000000..2f5aee079f --- /dev/null +++ b/src/celeritas/random/RngReseed.cu @@ -0,0 +1,67 @@ +//---------------------------------*-CUDA-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/random/RngReseed.cu +//---------------------------------------------------------------------------// +#include "RngReseed.hh" + +#include "corecel/device_runtime_api.h" +#include "corecel/Assert.hh" +#include "corecel/sys/Device.hh" +#include "corecel/sys/KernelParamCalculator.device.hh" + +#include "RngEngine.hh" + +namespace celeritas +{ +namespace +{ +//---------------------------------------------------------------------------// +// KERNELS +//---------------------------------------------------------------------------// +/*! + * Reinitialize the RNG states on device at the start of an event. + */ +__global__ void reseed_rng_kernel(DeviceCRef const params, + DeviceRef const state, + size_type event_id) +{ + auto tid = TrackSlotId{ + celeritas::KernelParamCalculator::thread_id().unchecked_get()}; + if (tid.get() < state.size()) + { + TrackSlotId tsid{tid.unchecked_get()}; + RngEngine::Initializer_t init; + init.seed = params.seed; + init.subsequence = event_id * state.size() + tsid.get(); + RngEngine rng(params, state, tsid); + rng = init; + } +} + +//---------------------------------------------------------------------------// +} // namespace + +//---------------------------------------------------------------------------// +// KERNEL INTERFACE +//---------------------------------------------------------------------------// +/*! + * Reinitialize the RNG states on device at the start of an event. + * + * Each thread's state is initialized using same seed and skipped ahead a + * different number of subsequences so the sequences on different threads will + * not have statistically correlated values. + */ +void reseed_rng(DeviceCRef const& params, + DeviceRef const& state, + size_type event_id) +{ + CELER_EXPECT(state); + CELER_EXPECT(params); + CELER_LAUNCH_KERNEL(reseed_rng, state.size(), 0, params, state, event_id); +} + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/random/RngReseed.hh b/src/celeritas/random/RngReseed.hh new file mode 100644 index 0000000000..ccc94093e8 --- /dev/null +++ b/src/celeritas/random/RngReseed.hh @@ -0,0 +1,43 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/random/RngReseed.hh +//---------------------------------------------------------------------------// +#pragma once + +#include "corecel/Assert.hh" +#include "corecel/Macros.hh" +#include "corecel/Types.hh" +#include "corecel/data/Collection.hh" + +#include "RngData.hh" + +namespace celeritas +{ +//---------------------------------------------------------------------------// +// Reinitialize the RNG states on host/device at the start of an event +void reseed_rng(DeviceCRef const&, + DeviceRef const&, + size_type); + +void reseed_rng(HostCRef const&, + HostRef const&, + size_type); + +#if !CELER_USE_DEVICE +//---------------------------------------------------------------------------// +/*! + * Reinitialize the RNG states on device at the start of an event. + */ +inline void reseed_rng(DeviceCRef const&, + DeviceRef const&, + size_type) +{ + CELER_ASSERT_UNREACHABLE(); +} +#endif + +//---------------------------------------------------------------------------// +} // namespace celeritas diff --git a/src/celeritas/random/XorwowRngData.hh b/src/celeritas/random/XorwowRngData.hh index bedff0e2c5..8cdcbc6d03 100644 --- a/src/celeritas/random/XorwowRngData.hh +++ b/src/celeritas/random/XorwowRngData.hh @@ -70,7 +70,7 @@ struct XorwowRngParamsData */ struct XorwowRngInitializer { - ull_int seed{0}; + Array seed{0}; ull_int subsequence{0}; ull_int offset{0}; }; diff --git a/src/celeritas/random/XorwowRngEngine.hh b/src/celeritas/random/XorwowRngEngine.hh index e3bdbe16c7..bf051432dc 100644 --- a/src/celeritas/random/XorwowRngEngine.hh +++ b/src/celeritas/random/XorwowRngEngine.hh @@ -164,7 +164,7 @@ XorwowRngEngine::operator=(Initializer_t const& init) auto& s = state_->xorstate; // Initialize the state from the seed - SplitMix64 rng{init.seed}; + SplitMix64 rng{init.seed[0]}; uint64_t seed = rng(); s[0] = static_cast(seed); s[1] = static_cast(seed >> 32); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index b3503ff48b..4e0b11cebc 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -557,6 +557,7 @@ set(CELERITASTEST_PREFIX celeritas/random) celeritas_add_device_test(celeritas/random/RngEngine) celeritas_add_test(celeritas/random/Selector.test.cc) +celeritas_add_test(celeritas/random/RngReseed.test.cc) celeritas_add_test(celeritas/random/XorwowRngEngine.test.cc GPU) celeritas_add_test(celeritas/random/distribution/BernoulliDistribution.test.cc) diff --git a/test/celeritas/random/RngReseed.test.cc b/test/celeritas/random/RngReseed.test.cc new file mode 100644 index 0000000000..10efa1cc54 --- /dev/null +++ b/test/celeritas/random/RngReseed.test.cc @@ -0,0 +1,89 @@ +//----------------------------------*-C++-*----------------------------------// +// Copyright 2021-2023 UT-Battelle, LLC, and other Celeritas developers. +// See the top-level COPYRIGHT file for details. +// SPDX-License-Identifier: (Apache-2.0 OR MIT) +//---------------------------------------------------------------------------// +//! \file celeritas/random/RngReseed.test.cc +//---------------------------------------------------------------------------// +#include "celeritas/random/RngReseed.hh" + +#include "corecel/OpaqueId.hh" +#include "corecel/Types.hh" +#include "corecel/data/CollectionStateStore.hh" +#include "celeritas/random/RngEngine.hh" +#include "celeritas/random/RngParams.hh" + +#include "celeritas_test.hh" + +namespace celeritas +{ +namespace test +{ +//---------------------------------------------------------------------------// +class RngReseedTest : public Test +{ + public: + using RngHostStore = CollectionStateStore; + + void SetUp() override { params = std::make_shared(12345); } + + std::shared_ptr params; +}; + +TEST_F(RngReseedTest, reseed) +{ + // Create and initialize states + size_type size = 1024; + RngHostStore states(params->host_ref(), StreamId{0}, size); + + size_type id = 8; + reseed_rng(params->host_ref(), states.ref(), id); + + RngEngine::Initializer_t init; + init.seed = params->host_ref().seed; + init.offset = 0; + + for (size_type i = 1; i < states.size(); ++i) + { + // Check that the reseeded RNGs were skipped ahead the correct number + // of subsequences + RngEngine skip_rng(params->host_ref(), states.ref(), TrackSlotId{0}); + init.subsequence = id * states.size() + i; + skip_rng = init; + + RngEngine rng(params->host_ref(), states.ref(), TrackSlotId{i}); + ASSERT_EQ(skip_rng(), rng()); + } + + reseed_rng(params->host_ref(), states.ref(), id); + std::vector values; + for (auto i : range(states.size()).step(128u)) + { + RngEngine rng(params->host_ref(), states.ref(), TrackSlotId{i}); + values.push_back(rng()); + } +#if CELERITAS_CORE_RNG == CELERITAS_CORE_RNG_CURAND + static unsigned int const expected_values[] = {65145249u, + 4154590960u, + 2216085262u, + 241608182u, + 2278993841u, + 1643630301u, + 2759037535u, + 3550652068u}; +#elif CELERITAS_CORE_RNG == CELERITAS_CORE_RNG_XORWOW + static unsigned int const expected_values[] = {3522223652u, + 296995412u, + 1414776235u, + 1609101469u, + 363980503u, + 2861073075u, + 1771581540u, + 3600889717u}; +#endif + EXPECT_VEC_EQ(values, expected_values); +} + +//---------------------------------------------------------------------------// +} // namespace test +} // namespace celeritas diff --git a/test/celeritas/random/XorwowRngEngine.test.cc b/test/celeritas/random/XorwowRngEngine.test.cc index 788c616904..28bbe45628 100644 --- a/test/celeritas/random/XorwowRngEngine.test.cc +++ b/test/celeritas/random/XorwowRngEngine.test.cc @@ -214,7 +214,7 @@ TEST_F(XorwowRngEngineTest, jump) XorwowRngEngine skip_rng(params->host_ref(), states.ref(), TrackSlotId{1}); XorwowRngInitializer init; - init.seed = 12345; + init.seed = {12345}; init.subsequence = 0; init.offset = 0; rng = init;