diff --git a/README.md b/README.md index b665724f7..cec5db69f 100644 --- a/README.md +++ b/README.md @@ -69,6 +69,11 @@ docker run --rm --gpus all ghcr.io/bnlnpps/simphony simg4ox -g tests/geom/raindr apptainer exec --nv docker://ghcr.io/bnlnpps/simphony simg4ox -g /workspaces/simphony/tests/geom/raindrop.gdml -m /workspaces/simphony/tests/run.mac ``` +These commands use `simg4ox`'s default serial Geant4 run manager. The published +images include multithreaded Geant4; use `tests/run_mt.mac` and add +`--threads N` to run the same torch photons on `N` Geant4 CPU workers. Opticks +GPU launches remain serialized because the GPU event context is process-wide. + ### Install with Spack Simphony is also available through the BNLNPPS Spack repository: diff --git a/docs/assets/simg4ox-event-processing.svg b/docs/assets/simg4ox-event-processing.svg new file mode 100644 index 000000000..83772a1fb --- /dev/null +++ b/docs/assets/simg4ox-event-processing.svg @@ -0,0 +1 @@ + \ No newline at end of file diff --git a/docs/assets/simg4ox-event-processing.typ b/docs/assets/simg4ox-event-processing.typ new file mode 100644 index 000000000..55dda818a --- /dev/null +++ b/docs/assets/simg4ox-event-processing.typ @@ -0,0 +1,162 @@ +// Regenerate from the repository root with: +// typst compile --format svg docs/assets/simg4ox-event-processing.typ \ +// docs/assets/simg4ox-event-processing.svg + +#let ink = rgb("#17324d") +#let muted = rgb("#5d6975") +#let panel-fill = rgb("#f8fafc") +#let panel-stroke = rgb("#cbd5df") +#let cpu-fill = rgb("#dceeff") +#let cpu-stroke = rgb("#3978a8") +#let gpu-fill = rgb("#eee3ff") +#let gpu-stroke = rgb("#7651a8") +#let wait-fill = rgb("#f0f2f4") +#let wait-stroke = rgb("#8a929b") +#let output-fill = rgb("#e3f5e8") +#let output-stroke = rgb("#3d8b59") + +#set page(width: 1050pt, height: auto, margin: 22pt, fill: white) +#set text(font: "DejaVu Sans", size: 10pt, fill: ink) +#set par(leading: 0.7em) + +#let timeline-columns = ( + 90pt, + 1fr, 1fr, 1fr, 1fr, 1fr, 1fr, + 1fr, 1fr, 1fr, 1fr, 1fr, 1fr, +) + +#let lane-label(body) = align( + right + horizon, + text(size: 8.5pt, weight: "bold", fill: ink, body), +) + +#let segment(fill-color, border-color, body) = box( + width: 100%, + height: 29pt, + fill: fill-color, + stroke: 0.8pt + border-color, + radius: 4pt, + inset: (x: 3pt, y: 4pt), + align(center + horizon, text(size: 8pt, weight: "bold", fill: ink, body)), +) + +#let panel(title, subtitle, body) = block( + width: 100%, + fill: panel-fill, + stroke: 0.8pt + panel-stroke, + radius: 7pt, + inset: 12pt, +)[ + #text(size: 12pt, weight: "bold")[#title] + #h(8pt) + #text(size: 8.5pt, fill: muted)[#subtitle] + #v(8pt) + #body +] + +#align(center)[ + #text(size: 17pt, weight: "bold")[Current `simg4ox` event processing] + #v(2pt) + #text(size: 9pt, fill: muted)[E0–E2 denote complete Geant4 events, each potentially containing many G4 tracks. The GPU event context is process-wide in both modes.] +] + +#v(12pt) + +#panel( + [Serial (`--threads 1`)], + [Each Geant4 event and its GPU photon transport finish before the next event starts.], + grid( + columns: timeline-columns, + column-gutter: 2pt, + row-gutter: 5pt, + align: horizon, + + [], + grid.cell(colspan: 12)[#align(right)[#text(size: 7.5pt, fill: muted)[time →]]], + + [#lane-label[Geant4 CPU]], + grid.cell(colspan: 2)[#segment(cpu-fill, cpu-stroke)[process event E0]], + grid.cell(colspan: 2)[], + grid.cell(colspan: 2)[#segment(cpu-fill, cpu-stroke)[process event E1]], + grid.cell(colspan: 2)[], + grid.cell(colspan: 2)[#segment(cpu-fill, cpu-stroke)[process event E2]], + grid.cell(colspan: 2)[], + + [#lane-label[Shared GPU]], + grid.cell(colspan: 2)[], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E0 photons]], + grid.cell(colspan: 2)[], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E1 photons]], + grid.cell(colspan: 2)[], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E2 photons]], + ), +) + +#v(10pt) + +#panel( + [Multithreaded (`--threads 3`)], + [Geant4 events are processed concurrently; their photons enter GPU transport strictly in event-ID order.], + grid( + columns: timeline-columns, + column-gutter: 2pt, + row-gutter: 5pt, + align: horizon, + + [], + grid.cell(colspan: 12)[#align(right)[#text(size: 7.5pt, fill: muted)[time →]]], + + [#lane-label[Worker 0]], + grid.cell(colspan: 3)[#segment(cpu-fill, cpu-stroke)[process event E0]], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E0 photons]], + grid.cell(colspan: 7)[], + + [#lane-label[Worker 1]], + grid.cell(colspan: 2)[#segment(cpu-fill, cpu-stroke)[process event E1]], + grid.cell(colspan: 3)[#segment(wait-fill, wait-stroke)[wait for E0]], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E1 photons]], + grid.cell(colspan: 5)[], + + [#lane-label[Worker 2]], + grid.cell(colspan: 4)[#segment(cpu-fill, cpu-stroke)[process event E2]], + grid.cell(colspan: 3)[#segment(wait-fill, wait-stroke)[wait for E0–E1]], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E2 photons]], + grid.cell(colspan: 3)[], + + [#lane-label[Shared GPU]], + grid.cell(colspan: 3)[], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E0 photons]], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E1 photons]], + grid.cell(colspan: 2)[#segment(gpu-fill, gpu-stroke)[transport E2 photons]], + grid.cell(colspan: 3)[], + ), +) + +#v(10pt) + +#align(center)[ + #box(width: 11pt, height: 11pt, fill: cpu-fill, stroke: 0.8pt + cpu-stroke, radius: 2pt) + #h(4pt) Geant4 event processing + #h(18pt) + #box(width: 11pt, height: 11pt, fill: wait-fill, stroke: 0.8pt + wait-stroke, radius: 2pt) + #h(4pt) Waiting for event-ID turn + #h(18pt) + #box(width: 11pt, height: 11pt, fill: gpu-fill, stroke: 0.8pt + gpu-stroke, radius: 2pt) + #h(4pt) Serialized GPU photon transport + #h(18pt) + #box(width: 11pt, height: 11pt, fill: output-fill, stroke: 0.8pt + output-stroke, radius: 2pt) + #h(4pt) Run-end merge and save +] + +#v(8pt) + +#align(center)[ + #box( + fill: output-fill, + stroke: 0.8pt + output-stroke, + radius: 4pt, + inset: (x: 10pt, y: 5pt), + )[ + After all events: merge hits by event ID → `s_hits.npy` + `g_hits.npy` + ] +] \ No newline at end of file diff --git a/docs/getting-started.md b/docs/getting-started.md index 1a491b0d1..ed89b5042 100644 --- a/docs/getting-started.md +++ b/docs/getting-started.md @@ -126,8 +126,14 @@ only the target you are changing, and run a relevant test: ctest --test-dir build -N cmake --build build --target simg4ox ctest --test-dir build -R raindrop +ctest --test-dir build -R '^Integration\.simg4ox_multithread$' ``` +The `simg4ox_multithread` integration test runs five optical-photon events +with two Geant4 CPU workers and validates the merged CPU/GPU hit arrays. See +[the simg4ox example](../examples/README.md#example-4-simg4ox-g4--gpu-validation) +for serial and MT command lines and the supplied `tests/run_mt.mac` macro. + The full suite includes GPU-backed tests. You can build without a GPU, but running those tests requires a compatible NVIDIA driver and GPU access from the container. @@ -142,6 +148,10 @@ If you prefer to manage the toolchain yourself, install: - CMake 3.22+ - Python 3.10+ +Build Geant4 with multithreading enabled to use `simg4ox --threads N` with +`N>1`. The project container images already use a multithreaded Geant4 build; +`--threads 1` remains available for serial validation. + With those dependencies available, clone, build, and test the project: ```shell @@ -182,6 +192,9 @@ To try the latest published release and confirm that GPU access works: docker run --rm --gpus all ghcr.io/bnlnpps/simphony simg4ox -g tests/geom/raindrop.gdml -m tests/run.mac ``` +This uses the default serial run manager. To exercise Geant4 MT in the +container, use `tests/run_mt.mac` and append `--threads N`. + To test an image built from your current checkout instead: ```shell @@ -189,6 +202,9 @@ docker build -t simphony:develop . docker run --rm --gpus all simphony:develop simg4ox -g tests/geom/raindrop.gdml -m tests/run.mac ``` +As with the published image, select `tests/run_mt.mac` and pass `--threads N` +for an MT run. + For day-to-day development, the Dev Container is more convenient because source edits are available immediately without rebuilding the image. @@ -200,6 +216,9 @@ On systems that provide Apptainer, run the same published release with: apptainer exec --nv docker://ghcr.io/bnlnpps/simphony simg4ox -g /workspaces/simphony/tests/geom/raindrop.gdml -m /workspaces/simphony/tests/run.mac ``` +For MT, use the absolute macro path +`/workspaces/simphony/tests/run_mt.mac` and add `--threads N`. + Use `singularity` in place of `apptainer` on systems that provide the older command name. diff --git a/docs/inputs-outputs.md b/docs/inputs-outputs.md index 0967c4c2b..7fce5c398 100644 --- a/docs/inputs-outputs.md +++ b/docs/inputs-outputs.md @@ -57,6 +57,28 @@ tracking is also run for validation. | `numphoton` | Number of photons to generate | | `wavelength` | Photon wavelength (nm) | +### `simg4ox` execution model + +`simg4ox` gives the same generated torch photons to Geant4 CPU tracking and +Simphony GPU tracking. + +Select the Geant4 run manager with `--threads N`. `N=1`, the default, uses the +serial run manager. `N>1` requires a multithreaded Geant4 build and creates `N` +CPU workers. This choice must be made before the macro is read, so +`/run/numberOfThreads` is not the thread-selection interface for `simg4ox`. +Use a macro such as `tests/run_mt.mac` for run initialization and event count: + +```bash +simg4ox -g tests/geom/opticks_raindrop.gdml -c dev \ + -m tests/run_mt.mac -s 42 --threads 4 +``` + +In MT mode, actions and sensitive detectors are worker-local. Completed events +share one process-wide Opticks GPU context, so GPU launches are serialized in +event-ID order while Geant4 CPU tracking remains multithreaded. Run-wide CPU +and GPU hits are merged in event-ID order. Full CPU-side Opticks photon-history +recording is available only with the serial run manager. + ## Defining primary particles For charged-particle examples, the user or developer defines the primary @@ -150,6 +172,13 @@ without the final `A000` or `B000` event-index subdirectory. The app-level `Config` interface does not expose a separate save-mode field; prefer `Config::output_dir` for controlling where event folders are written. +`simg4ox` also writes application-level run aggregates directly under +`Config::output_dir`: `s_hits.npy` contains Simphony GPU hits and `g_hits.npy` +contains Geant4 sensitive-detector hits. Both have shape `(H, 4, 4)`, dtype +`float32`, and the `sphoton` layout described below. They are written in both +serial and MT modes. In MT mode, CPU-side history folders are not produced, but +`g_hits.npy` is still written from the merged Geant4 hit collections. + ### File schemas The table assumes the default full photon representation and unmerged hits: diff --git a/docs/performance-and-debugging.md b/docs/performance-and-debugging.md index 3f1ab2127..9c4d92aa6 100644 --- a/docs/performance-and-debugging.md +++ b/docs/performance-and-debugging.md @@ -26,6 +26,20 @@ docker run --rm -t -v /tmp/out:/tmp/out simphony:release \ run-performance -g tests/geom/opticks_raindrop.gdml -o /tmp/out/release ``` +### Interpreting `simg4ox` MT timings + +`simg4ox --threads N` parallelizes Geant4 CPU tracking of its configured torch +photons. It does not run multiple Opticks launches concurrently: the current +process-wide GPU event context requires launches to be serialized in event-ID +order. Consequently, an end-to-end `simg4ox` wall-clock measurement includes +parallel CPU work plus serialized GPU work and should not be interpreted as a +pure GPU speed-up measurement. + +![Serial and multithreaded simg4ox event-processing timelines](assets/simg4ox-event-processing.svg) + +The supplied `tests/run_mt.mac` contains only five low-statistics events and is +an integration check, not a benchmark. + ## Debug analysis with `optiphy/ana/photon_history_summary.py` The script analyzes GPU optical photon simulation output to debug where diff --git a/examples/README.md b/examples/README.md index 9610e0f35..47a42a5b3 100644 --- a/examples/README.md +++ b/examples/README.md @@ -7,7 +7,7 @@ Simphony provides several examples demonstrating GPU-accelerated optical photon | `simphox` | Optical photons (torch) | None | External project build and CPU/GPU photon generation smoke test | | `GPUCerenkov` | Cerenkov only | Simple nested boxes (raindrop) | Basic Cerenkov testing | | `GPURaytrace` | Cerenkov + Scintillation | 8x8 CsI crystal + SiPM array | Realistic detector simulation | -| `simg4ox` | Optical photons (torch) | Any GDML | G4 + GPU side-by-side validation | +| `simg4ox` | Optical photons (torch) | Any GDML | Serial/MT G4 + GPU side-by-side validation | | `GPUPhotonSourceMinimal` | Optical photons (torch) | Any GDML | GPU-only test | | `GPUPhotonFileSource` | Optical photons (text file) | Any GDML | GPU-only, user-defined photons from file | | WLS test | Wavelength shifting | WLS sphere + detector shell | Validate GPU WLS physics | @@ -26,7 +26,7 @@ Geant4 optical-photon tracking is also run for validation. | Photon input from text file | No | No | No | No | Yes | | G4 optical photon tracking | Yes | Yes | Yes | No | No | | GPU simulation (Simphony) | Yes | Yes | Yes | Yes | Yes | -| Multi-threaded | Yes | Yes | No | No | No | +| Geant4 CPU multithreading | Yes | Yes | Yes | No | No | `GPUCerenkov` and `GPURaytrace` collect gensteps from charged-particle interactions and pass them to Simphony for GPU photon generation and tracing. @@ -35,6 +35,10 @@ a torch configuration. `simg4ox` runs both G4 and GPU tracking for validation, while `GPUPhotonSourceMinimal` keeps only the GPU path. `GPUPhotonFileSource` reads user-defined photons from a text file. +For `simg4ox`, multithreading applies to Geant4 CPU tracking. Its process-wide +Opticks event context is protected by serializing GPU launches in event-ID +order. + ### Example 1: simphox (External build smoke test) Assuming Simphony is properly installed on the system, compile and run this @@ -109,13 +113,14 @@ grep -c "CreationProcessID=1" opticks_hits_output.txt # Scintillation ### Example 4: simg4ox (G4 + GPU Validation) -`simg4ox` generates optical photons from a configurable torch source and runs -both Geant4 and Simphony GPU simulation in parallel on the same input photons. This +`simg4ox` generates optical photons from a configurable torch source and +tracks the same input with Geant4 on the CPU and Simphony on the GPU. This enables direct comparison of hit counts and positions between the two engines. -Both engines detect photons using the same mechanism: border surface physics. On the G4 -side the `SteppingAction` records a hit when `G4OpBoundaryProcess` reports Detection at -the optical surface, matching how Simphony detects photons on the GPU. +Both engines detect photons using border-surface physics. On the Geant4 side, +the optical boundary process invokes the worker-local `PhotonSD` at configured +sensitive surfaces; the detector records and terminates the photon. Simphony +selects GPU hits from the corresponding optical-boundary result. | Argument | Description | Default | |----------|-------------|---------| @@ -124,12 +129,32 @@ the optical surface, matching how Simphony detects photons on the GPU. | `-m, --macro` | Path to G4 macro | `run.mac` | | `-i, --interactive` | Open interactive viewer | off | | `-s, --seed` | Fixed random seed | Geant4 default | +| `-t, --threads` | Geant4 CPU threads; `1` selects the serial run manager | `1` | + +The same executable supports serial and multi-threaded Geant4 runs: ```bash -simg4ox -g tests/geom/opticks_raindrop.gdml -c dev -m run.mac -s 42 +# Serial Geant4 (default) +simg4ox -g tests/geom/opticks_raindrop.gdml -c dev -m tests/run_5evt.mac -s 42 + +# Geant4 MT with four CPU workers +simg4ox -g tests/geom/opticks_raindrop.gdml -c dev -m tests/run_mt.mac -s 42 --threads 4 ``` +The thread count is a command-line option because the executable must select +the serial or MT run-manager type before Geant4 reads a macro. The example +`tests/run_mt.mac` contains the remaining run initialization and event-count +commands. + +In MT mode, Geant4 actions and sensitive detectors are worker-local and the +run-wide results are merged in event-ID order. Opticks currently exposes one +process-wide GPU event context, so GPU launches are serialized in that same +order while Geant4 CPU tracking remains multi-threaded. Full CPU-side Opticks +photon-history recording remains available in serial mode; both modes write +the run-wide hit arrays below. + **Output:** + - `s_hits.npy` — Simphony GPU hits - `g_hits.npy` — Geant4 hits diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 451988831..e2ad170b0 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -14,7 +14,7 @@ target_include_directories(consgeo PRIVATE ) target_link_libraries(consgeo simphony_argparse) -# simg4ox runs Geant4 and OptiX simulations +# simg4ox validates torch photons with serial/MT Geant4 and serialized OptiX launches add_executable(simg4ox simg4ox.cpp g4app.h) target_link_libraries(simg4ox SysRap simphony_g4_deps simphony_argparse) diff --git a/src/g4app.h b/src/g4app.h index 69dc3ee7a..fa6a09dad 100644 --- a/src/g4app.h +++ b/src/g4app.h @@ -1,12 +1,22 @@ +#pragma once + +#include +#include #include #include +#include +#include +#include #include #include #include #include +#include + #include "G4BooleanSolid.hh" #include "G4Event.hh" +#include "G4Exception.hh" #include "G4GDMLParser.hh" #include "G4LogicalVolumeStore.hh" #include "G4OpBoundaryProcess.hh" @@ -14,9 +24,12 @@ #include "G4PhysicalConstants.hh" #include "G4PrimaryParticle.hh" #include "G4PrimaryVertex.hh" +#include "G4Run.hh" +#include "G4RunManager.hh" #include "G4SDManager.hh" #include "G4SubtractionSolid.hh" #include "G4SystemOfUnits.hh" +#include "G4Threading.hh" #include "G4ThreeVector.hh" #include "G4Track.hh" #include "G4TrackStatus.hh" @@ -28,7 +41,9 @@ #include "G4UserTrackingAction.hh" #include "G4VPhysicalVolume.hh" #include "G4VProcess.hh" +#include "G4VUserActionInitialization.hh" #include "G4VUserDetectorConstruction.hh" +#include "G4VUserEventInformation.hh" #include "G4VUserPrimaryGeneratorAction.hh" #include "g4cx/G4CXOpticks.hh" @@ -77,6 +92,15 @@ using PhotonHitsCollection = G4THitsCollection; static_assert(sizeof(sphoton) == 16 * sizeof(float)); static_assert(std::is_trivially_copyable_v); +inline std::unique_ptr MakePhotonArray(const std::vector& photons) +{ + const size_t num_floats = photons.size() * 4 * 4; + const float* data = reinterpret_cast(photons.data()); + std::unique_ptr array(NP::MakeFromValues(data, num_floats)); + array->reshape({static_cast(photons.size()), 4, 4}); + return array; +} + struct PhotonSD : public G4VSensitiveDetector { PhotonHitsCollection* photon_hit_collection{nullptr}; @@ -171,10 +195,25 @@ struct DetectorConstruction : G4VUserDetectorConstruction } }; +struct PrimaryPhotonInfo : G4VUserEventInformation +{ + explicit PrimaryPhotonInfo(std::vector photons) : + photons(std::move(photons)) + { + } + + void Print() const override + { + } + + std::vector photons; +}; + struct PrimaryGenerator : G4VUserPrimaryGeneratorAction { - simphony::Config cfg; - SEvt* sev; + simphony::Config cfg; + SEvt* sev; + std::unique_ptr input_photon; PrimaryGenerator(const simphony::Config& cfg, SEvt* sev) : cfg(cfg), @@ -186,12 +225,6 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction { std::vector sphotons = generate_photons(cfg.torch); - size_t num_floats = sphotons.size() * 4 * 4; - float* data = reinterpret_cast(sphotons.data()); - NP* photons = NP::MakeFromValues(data, num_floats); - - photons->reshape({static_cast(sphotons.size()), 4, 4}); - for (const sphoton& p : sphotons) { G4ThreeVector position_mm(p.pos.x, p.pos.y, p.pos.z); @@ -212,47 +245,121 @@ struct PrimaryGenerator : G4VUserPrimaryGeneratorAction event->AddPrimaryVertex(vertex); } - sev->SetInputPhoton(photons); + // The Opticks CPU recorder is intentionally used only by the serial + // run manager. Its event instance is process-global and cannot be + // shared safely by Geant4 worker threads. Preserve MT event input on + // the G4Event until the event reaches the serialized GPU section. + if (sev) + { + input_photon = MakePhotonArray(sphotons); + SEvt::SetInputPhoton(input_photon.get()); + } + else + event->SetUserInformation(new PrimaryPhotonInfo(std::move(sphotons))); } }; -struct EventAction : G4UserEventAction +struct Simg4oxRun : G4Run { - simphony::Config cfg; - SEvt* sev; - std::vector g4_hits; - std::vector gpu_hits; + struct EventHits + { + std::vector gpu; + std::vector g4; + }; - EventAction(const simphony::Config& cfg, SEvt* sev) : - cfg(cfg), - sev(sev) + std::map hits_by_event; + + void AddEvent(G4int event_id, std::vector gpu_hits, std::vector g4_hits) { + hits_by_event.insert_or_assign(event_id, EventHits{std::move(gpu_hits), std::move(g4_hits)}); } - void BeginOfEventAction(const G4Event* event) override + void Merge(const G4Run* run) override { - sev->beginOfEvent(event->GetEventID()); + const auto* local_run = static_cast(run); + for (const auto& [event_id, hits] : local_run->hits_by_event) + hits_by_event.emplace(event_id, hits); + + G4Run::Merge(run); } - void ClearHits() + static std::vector Flatten( + const std::map& events, + const std::vector EventHits::* member) { - g4_hits.clear(); - gpu_hits.clear(); + size_t total = 0; + for (const auto& [event_id, hits] : events) + { + (void)event_id; + total += (hits.*member).size(); + } + + std::vector flattened; + flattened.reserve(total); + for (const auto& [event_id, hits] : events) + { + (void)event_id; + const auto& event_hits = hits.*member; + flattened.insert(flattened.end(), event_hits.begin(), event_hits.end()); + } + return flattened; + } + + std::vector GPUHits() const + { + return Flatten(hits_by_event, &EventHits::gpu); } + std::vector G4Hits() const + { + return Flatten(hits_by_event, &EventHits::g4); + } +}; - size_t CollectGPUHits(SEvt* sev_gpu) +struct Simg4oxSharedState +{ + std::mutex gpu_mutex; + std::condition_variable gpu_turn; + std::unique_ptr input_photon; + G4int next_gpu_event{0}; + + void BeginRun() { - const size_t num_gpu_hits = sev_gpu->getNumHit(); - const size_t offset = gpu_hits.size(); - gpu_hits.resize(offset + num_gpu_hits); + std::lock_guard lock(gpu_mutex); + next_gpu_event = 0; + } +}; + +struct EventAction : G4UserEventAction +{ + SEvt* sev; + std::shared_ptr shared_state; + bool order_gpu_events; + + EventAction(SEvt* sev, std::shared_ptr shared_state, bool order_gpu_events) : + sev(sev), + shared_state(std::move(shared_state)), + order_gpu_events(order_gpu_events) + { + } + + void BeginOfEventAction(const G4Event* event) override + { + if (sev) + sev->beginOfEvent(event->GetEventID()); + } + + static std::vector CollectGPUHits(SEvt* sev_gpu) + { + const size_t num_gpu_hits = sev_gpu->getNumHit(); + std::vector gpu_hits(num_gpu_hits); for (size_t idx = 0; idx < num_gpu_hits; idx++) - sev_gpu->getHit(gpu_hits[offset + idx], idx); + sev_gpu->getHit(gpu_hits[idx], idx); - return num_gpu_hits; + return gpu_hits; } - size_t CollectG4Hits(const G4Event* event) + static std::vector CollectG4Hits(const G4Event* event) { G4HCofThisEvent* hce = event->GetHCofThisEvent(); size_t num_g4_hits = 0; @@ -267,7 +374,8 @@ struct EventAction : G4UserEventAction } } - g4_hits.reserve(g4_hits.size() + num_g4_hits); + std::vector g4_hits; + g4_hits.reserve(num_g4_hits); if (hce) { @@ -282,53 +390,73 @@ struct EventAction : G4UserEventAction } } - return num_g4_hits; + return g4_hits; } - void SaveHits(const std::vector& source, const char* name) const + std::vector SimulateOnGPU(const G4Event* event) { - NP* hits = NP::Make(source.size(), 4, 4); - if (!source.empty()) - std::memcpy(hits->bytes(), source.data(), source.size() * sizeof(sphoton)); - - hits->save(cfg.output_dir.string().c_str(), name); - delete hits; - } + const PrimaryPhotonInfo* primary_info = nullptr; + if (order_gpu_events) + { + primary_info = dynamic_cast(event->GetUserInformation()); + if (!primary_info) + { + G4Exception("EventAction::SimulateOnGPU", "MissingPrimaryPhotonInfo", FatalException, + "MT event is missing its generated photons for GPU processing"); + return {}; + } + } - void SaveRunHits() const - { - SaveHits(gpu_hits, "s_hits.npy"); - SaveHits(g4_hits, "g_hits.npy"); - } + const G4int event_id = event->GetEventID(); + std::unique_lock lock(shared_state->gpu_mutex); + if (order_gpu_events) + shared_state->gpu_turn.wait(lock, [&] { return event_id == shared_state->next_gpu_event; }); - void EndOfEventAction(const G4Event* event) override - { - int eventID = event->GetEventID(); - sev->addEventConfigArray(); - sev->gather(); - sev->endOfEvent(eventID); + if (order_gpu_events) + { + shared_state->input_photon = MakePhotonArray(primary_info->photons); + SEvt::SetInputPhoton(shared_state->input_photon.get()); + } - // GPU-based simulation G4CXOpticks* gx = G4CXOpticks::Get(); - - gx->simulate(eventID, false); + gx->simulate(event_id, false); cudaDeviceSynchronize(); - SEvt* sev_gpu = SEvt::Get_EGPU(); - size_t num_hits_gpu = sev_gpu->getNumHit(); - size_t num_hits_cpu = sev->getNumHit(); + SEvt* sev_gpu = SEvt::Get_EGPU(); + auto gpu_hits = CollectGPUHits(sev_gpu); + gx->reset(event_id); + + if (order_gpu_events) + { + ++shared_state->next_gpu_event; + lock.unlock(); + shared_state->gpu_turn.notify_all(); + } + + return gpu_hits; + } + + void EndOfEventAction(const G4Event* event) override + { + const G4int event_id = event->GetEventID(); + if (sev) + { + sev->addEventConfigArray(); + sev->gather(); + sev->endOfEvent(event_id); + G4cout << "EventAction::EndOfEventAction: CPU hits: " << sev->getNumHit() << G4endl; + } - G4cout << "EventAction::EndOfEventAction: GPU hits: " << num_hits_gpu << G4endl; - G4cout << "EventAction::EndOfEventAction: CPU hits: " << num_hits_cpu << G4endl; + auto g4_hits = CollectG4Hits(event); + auto gpu_hits = SimulateOnGPU(event); - // Append the event-wide GPU buffer and all Geant4 hit collections to - // run-scoped buffers before the event data is reset by either backend. - size_t collected_gpu_hits = CollectGPUHits(sev_gpu); - size_t collected_g4_hits = CollectG4Hits(event); - G4cout << "EventAction::EndOfEventAction: Collected GPU hits: " << collected_gpu_hits << G4endl; - G4cout << "EventAction::EndOfEventAction: Collected G4 hits: " << collected_g4_hits << G4endl; + G4cout << "EventAction::EndOfEventAction: Event " << event_id + << ": Collected GPU hits: " << gpu_hits.size() << G4endl; + G4cout << "EventAction::EndOfEventAction: Event " << event_id + << ": Collected G4 hits: " << g4_hits.size() << G4endl; - gx->reset(eventID); + auto* run = static_cast(G4RunManager::GetRunManager()->GetNonConstCurrentRun()); + run->AddEvent(event_id, std::move(gpu_hits), std::move(g4_hits)); } }; @@ -453,6 +581,9 @@ struct TrackingAction : G4UserTrackingAction mutable_track->UseGivenVelocity(true); } + if (!sev) + return; + if (!STrackInfo::Exists(track)) PreUserTrackingAction_Optical_FabricateLabel(track); @@ -482,6 +613,9 @@ struct TrackingAction : G4UserTrackingAction void PostUserTrackingAction(const G4Track* track) override { + if (!sev) + return; + G4TrackStatus tstat = track->GetTrackStatus(); bool is_stop_and_kill = tstat == fStopAndKill; @@ -511,45 +645,88 @@ struct TrackingAction : G4UserTrackingAction struct RunAction : G4UserRunAction { - EventAction* event_action; + simphony::Config cfg; + std::shared_ptr shared_state; - RunAction(EventAction* eventAction) : - event_action(eventAction) + RunAction(const simphony::Config& cfg, std::shared_ptr shared_state) : + cfg(cfg), + shared_state(std::move(shared_state)) { } + G4Run* GenerateRun() override + { + return new Simg4oxRun; + } + void BeginOfRunAction(const G4Run*) override { - event_action->ClearHits(); + if (!G4Threading::IsWorkerThread()) + shared_state->BeginRun(); } - void EndOfRunAction(const G4Run*) override + void SaveHits(const std::vector& source, const char* name) const { - event_action->SaveRunHits(); - G4cout << "RunAction::EndOfRunAction: Total GPU hits: " << event_action->gpu_hits.size() << G4endl; - G4cout << "RunAction::EndOfRunAction: Total G4 hits: " << event_action->g4_hits.size() << G4endl; + NP* hits = NP::Make(source.size(), 4, 4); + if (!source.empty()) + std::memcpy(hits->bytes(), source.data(), source.size() * sizeof(sphoton)); + + hits->save(cfg.output_dir.string().c_str(), name); + delete hits; + } + + void EndOfRunAction(const G4Run* run) override + { + // Worker runs are merged by G4MTRunManager after their EndOfRunAction. + // Only the serial/master run owns the complete result. + if (G4Threading::IsWorkerThread()) + return; + + const auto* simg4ox_run = static_cast(run); + auto gpu_hits = simg4ox_run->GPUHits(); + auto g4_hits = simg4ox_run->G4Hits(); + + SaveHits(gpu_hits, "s_hits.npy"); + SaveHits(g4_hits, "g_hits.npy"); + G4cout << "RunAction::EndOfRunAction: Total GPU hits: " << gpu_hits.size() << G4endl; + G4cout << "RunAction::EndOfRunAction: Total G4 hits: " << g4_hits.size() << G4endl; } }; -struct G4App +struct ActionInitialization : G4VUserActionInitialization { - G4App(const simphony::Config& cfg, std::filesystem::path gdml_file) : - sev(SEvt::CreateOrReuse_ECPU()), - det_cons_(new DetectorConstruction(gdml_file)), - prim_gen_(new PrimaryGenerator(cfg, sev)), - event_act_(new EventAction(cfg, sev)), - run_act_(new RunAction(event_act_)), - stepping_(new SteppingAction(sev)), - tracking_(new TrackingAction(sev)) + simphony::Config cfg; + std::shared_ptr shared_state; + bool multithreaded; + + ActionInitialization( + const simphony::Config& cfg, + std::shared_ptr shared_state, + bool multithreaded) : + cfg(cfg), + shared_state(std::move(shared_state)), + multithreaded(multithreaded) { } - SEvt* sev; + void BuildForMaster() const override + { + SetUserAction(new RunAction(cfg, shared_state)); + } - G4VUserDetectorConstruction* det_cons_; - G4VUserPrimaryGeneratorAction* prim_gen_; - EventAction* event_act_; - RunAction* run_act_; - SteppingAction* stepping_; - TrackingAction* tracking_; + void Build() const override + { + // SEvt's CPU instance and its profiling/persistence helpers are + // process-global. Keep the full CPU history recorder in serial mode; + // MT workers still perform normal Geant4 tracking and collect SD hits. + SEvt* sev = multithreaded ? nullptr : SEvt::CreateOrReuse_ECPU(); + + SetUserAction(new PrimaryGenerator(cfg, sev)); + SetUserAction(new RunAction(cfg, shared_state)); + SetUserAction(new EventAction(sev, shared_state, multithreaded)); + SetUserAction(new TrackingAction(sev)); + + if (sev) + SetUserAction(new SteppingAction(sev)); + } }; diff --git a/src/simg4ox.cpp b/src/simg4ox.cpp index cb3068e54..becfce780 100644 --- a/src/simg4ox.cpp +++ b/src/simg4ox.cpp @@ -4,15 +4,17 @@ #include "FTFP_BERT.hh" #include "G4OpticalPhysics.hh" -#include "Randomize.hh" #include "G4RunManager.hh" +#include "G4RunManagerFactory.hh" #include "G4VModularPhysicsList.hh" +#include "Randomize.hh" #include "G4UIExecutive.hh" #include "G4UImanager.hh" #include "G4VisExecutive.hh" #include "sysrap/OPTICKS_LOG.hh" +#include "sysrap/SEventConfig.hh" #include "config.h" #include "g4app.h" @@ -53,6 +55,11 @@ int main(int argc, char** argv) program.add_argument("-s", "--seed").help("fixed random seed").scan<'i', long>(); + program.add_argument("-t", "--threads") + .help("number of Geant4 CPU worker threads (1 selects the serial run manager)") + .default_value(1) + .scan<'i', int>(); + try { program.parse_args(argc, argv); @@ -64,8 +71,28 @@ int main(int argc, char** argv) exit(EXIT_FAILURE); } + const int num_threads = program.get("--threads"); + if (num_threads < 1) + { + cerr << "--threads must be a positive integer" << endl; + return EXIT_FAILURE; + } + +#ifndef G4MULTITHREADED + if (num_threads > 1) + { + cerr << "This Geant4 installation was built without multithreading support" << endl; + return EXIT_FAILURE; + } +#endif + simphony::Config cfg(config_name); + // Device discovery used to happen as a side effect of constructing the + // serial CPU SEvt. MT workers deliberately do not share that instance, so + // initialize global event/device metadata explicitly on the master. + SEventConfig::Initialize(); + if (program.is_used("--seed")) { const long seed = program.get("--seed"); @@ -78,17 +105,29 @@ int main(int argc, char** argv) G4VModularPhysicsList* physics = new FTFP_BERT; physics->RegisterPhysics(new G4OpticalPhysics); - G4RunManager run_mgr; - run_mgr.SetUserInitialization(physics); + const bool multithreaded = num_threads > 1; + G4RunManager* run_mgr = multithreaded + ? G4RunManagerFactory::CreateRunManager(G4RunManagerType::MTOnly, true, num_threads) + : G4RunManagerFactory::CreateRunManager(G4RunManagerType::SerialOnly); + + const G4int configured_threads = multithreaded ? run_mgr->GetNumberOfThreads() : 1; + if (multithreaded && configured_threads != num_threads) + { + cerr << "Requested " << num_threads << " Geant4 CPU threads, but the run manager configured " + << configured_threads << endl; + delete run_mgr; + delete physics; + return EXIT_FAILURE; + } + + G4cout << "simg4ox: Geant4 run manager: " << (multithreaded ? "MT" : "serial") + << ", CPU threads: " << configured_threads << G4endl; + + run_mgr->SetUserInitialization(physics); + run_mgr->SetUserInitialization(new DetectorConstruction(gdml_file)); - G4App* g4app = new G4App(cfg, gdml_file); - run_mgr.SetUserInitialization(g4app->det_cons_); - run_mgr.SetUserAction(g4app->prim_gen_); - run_mgr.SetUserAction(g4app->run_act_); - run_mgr.SetUserAction(g4app->event_act_); - run_mgr.SetUserAction(g4app->tracking_); - run_mgr.SetUserAction(g4app->stepping_); - run_mgr.Initialize(); + auto shared_state = std::make_shared(); + run_mgr->SetUserInitialization(new ActionInitialization(cfg, shared_state, multithreaded)); G4UIExecutive* uix = nullptr; G4VisManager* vis = nullptr; @@ -109,6 +148,8 @@ int main(int argc, char** argv) } delete uix; + delete vis; + delete run_mgr; return EXIT_SUCCESS; } diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 725d4167f..06eb28498 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -1,4 +1,5 @@ find_package(Python3 REQUIRED COMPONENTS Interpreter) +find_package(Geant4 REQUIRED OPTIONAL_COMPONENTS multithreaded) find_program(BASH_EXECUTABLE bash REQUIRED) set(SIMPHONY_SIMG4OX_TEST_WORKDIR "${CMAKE_CURRENT_BINARY_DIR}/simg4ox") @@ -48,3 +49,25 @@ set_tests_properties(Integration.simg4ox_multievent PROPERTIES LABELS "integration" TIMEOUT 120 ) + +if(Geant4_multithreaded_FOUND) + set(SIMPHONY_SIMG4OX_MULTITHREAD_TEST_WORKDIR "${CMAKE_CURRENT_BINARY_DIR}/simg4ox_multithread") + file(MAKE_DIRECTORY "${SIMPHONY_SIMG4OX_MULTITHREAD_TEST_WORKDIR}") + + add_test( + NAME Integration.simg4ox_multithread + COMMAND ${CMAKE_COMMAND} -E env + REPO_DIR=${PROJECT_SOURCE_DIR} + PYTHON=${Python3_EXECUTABLE} + SIMG4OX_BIN=$ + THREADS=2 + ${BASH_EXECUTABLE} ${CMAKE_CURRENT_SOURCE_DIR}/test_simg4ox_multithread.sh + ) + + set_tests_properties(Integration.simg4ox_multithread PROPERTIES + WORKING_DIRECTORY "${SIMPHONY_SIMG4OX_MULTITHREAD_TEST_WORKDIR}" + LABELS "integration" + PROCESSORS 2 + TIMEOUT 120 + ) +endif() diff --git a/tests/run_mt.mac b/tests/run_mt.mac new file mode 100644 index 000000000..05ba1ee8a --- /dev/null +++ b/tests/run_mt.mac @@ -0,0 +1,7 @@ +# Multi-event batch run for the simg4ox Geant4 MT example. +# Select the worker count before this macro is read with: +# simg4ox --threads N ... -m tests/run_mt.mac +/run/verbose 1 +/process/optical/boundary/setInvokeSD true +/run/initialize +/run/beamOn 5 diff --git a/tests/test_simg4ox_multithread.sh b/tests/test_simg4ox_multithread.sh new file mode 100755 index 000000000..33700d391 --- /dev/null +++ b/tests/test_simg4ox_multithread.sh @@ -0,0 +1,29 @@ +#!/usr/bin/env bash + +set -euo pipefail + +SCRIPT_DIR=$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd) +REPO_DIR=${REPO_DIR:-$(cd "${SCRIPT_DIR}/.." && pwd)} + +SIMG4OX_BIN=${SIMG4OX_BIN:-simg4ox} +PYTHON=${PYTHON:-python3} +THREADS=${THREADS:-2} +RUN_LOG="${PWD}/simg4ox-multithread.log" + +export OPTICKS_HOME="${REPO_DIR}" +export SIMPHONY_CONFIG_DIR="${SIMPHONY_CONFIG_DIR:-${REPO_DIR}/config}" +export PYTHONPATH="${REPO_DIR}${PYTHONPATH:+:${PYTHONPATH}}" + +rm -f "${PWD}/g_hits.npy" "${PWD}/s_hits.npy" "${RUN_LOG}" + +"${SIMG4OX_BIN}" \ + -g "${REPO_DIR}/tests/geom/opticks_raindrop.gdml" \ + -m "${REPO_DIR}/tests/run_mt.mac" \ + -c dev \ + -s 42 \ + --threads "${THREADS}" 2>&1 | tee "${RUN_LOG}" + +"${PYTHON}" "${REPO_DIR}/tests/check_simg4ox_multievent.py" \ + --log "${RUN_LOG}" \ + --output-dir "${PWD}" \ + --events 5