diff --git a/docs/src/examples/G4example.cpp b/docs/src/examples/G4example.cpp index 1ec2656..1911c5a 100644 --- a/docs/src/examples/G4example.cpp +++ b/docs/src/examples/G4example.cpp @@ -14,6 +14,7 @@ //#include "G4SystemOfUnits.hh" #include "globals.hh" #include "G4VUserActionInitialization.hh" +#include "G4UserWorkerInitialization.hh" #include #include @@ -37,6 +38,7 @@ class DetectorConstruction : public G4VUserDetectorConstruction new G4PVPlacement(0, G4ThreeVector(), logicCore, "Core", logicWorld, false, 0); return physWorld; } + void ConstructSDandField() override {} }; //---Primary generator action class---------------------------------------------------------------- @@ -107,13 +109,32 @@ class RunAction : public G4UserRunAction runaction_f endrun_jl; }; + +class WorkerInitialization : public G4UserWorkerInitialization { + public: + WorkerInitialization() = default; + virtual ~WorkerInitialization() = default; + virtual void WorkerInitialize() const override { + if (jl_get_pgcstack() == NULL) jl_adopt_thread(); + } + virtual void WorkerStart() const override {} + virtual void WorkerRunStart() const override {} + virtual void WorkerRunEnd() const override { + jl_ptls_t ptls = jl_current_task->ptls; + jl_gc_safe_enter(ptls); + } + virtual void WorkerStop() const override {} +}; + //---Action initialization class------------------------------------------------------------------- class ActionInitialization : public G4VUserActionInitialization { public: ActionInitialization() = default; ~ActionInitialization() override = default; - void BuildForMaster() const override {} + void BuildForMaster() const override { + SetUserAction(new RunAction); + } void Build() const override { SetUserAction(new PrimaryGeneratorAction); SetUserAction(new RunAction); @@ -122,12 +143,17 @@ class ActionInitialization : public G4VUserActionInitialization } }; + + JULIA_DEFINE_FAST_TLS // only define this once, in an executable (not in a shared library) if you want fast code. //----Main program--------------------------------------------------------------------------------- int main(int, char**) { - //--- Required to setup the Julia context + int nthreads = 0; // Default number of threads + if (getenv("G4NUMTHREADS")) nthreads = atoi(getenv("G4NUMTHREADS")); + + //--- Required to setup the Julia context------------------------------------------------------ jl_init(); /* run Julia commands */ jl_eval_string("include(\"MyCode.jl\")"); @@ -135,10 +161,17 @@ int main(int, char**) std::cout << "=====> " << jl_typeof_str(jl_exception_occurred()) << std::endl; } - //---Construct the default run manager - auto runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial); + //---Construct the default run manager (taking into account the number of threads)------------- + G4RunManager* runManager; + if (nthreads > 0) { + runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::MT); + runManager->SetNumberOfThreads(nthreads); + runManager->SetUserInitialization(new WorkerInitialization()); + } else { + runManager = G4RunManagerFactory::CreateRunManager(G4RunManagerType::Serial); + } - //---Set mandatory initialization classes + //---Set mandatory initialization classes------------------------------------------------------ // Detector construction runManager->SetUserInitialization(new DetectorConstruction()); @@ -159,9 +192,14 @@ int main(int, char**) //UImanager->ApplyCommand("/tracking/verbose 1"); UImanager->ApplyCommand("/gun/particle e+"); UImanager->ApplyCommand("/gun/energy 100 MeV"); - UImanager->ApplyCommand("/run/beamOn 100000"); - // Job termination + // Start a run (we need to enter GC safe region here because the worker threads + // will enter a wait state while waiting for workers to finish and this can block GC + auto state = jl_gc_safe_enter(jl_current_task->ptls); + runManager->BeamOn(100000); + jl_gc_safe_leave(jl_current_task->ptls, state); + + // Job termination--------------------------------------------------------------------------- delete runManager; // strongly recommended: notify Julia that the program is about to terminate. diff --git a/docs/src/examples/JuliaAction_lit.jl b/docs/src/examples/JuliaAction_lit.jl index b20b2ed..0759c59 100644 --- a/docs/src/examples/JuliaAction_lit.jl +++ b/docs/src/examples/JuliaAction_lit.jl @@ -85,10 +85,11 @@ Base.run(`c++ -O2 -fPIC $cflags -I$jlprefix/include/julia $ldflags $g4libs $jlli # ## Run the application # We need to set the variable `JULIA_PROJECT` pointing to correctly setup Julia environment. -withenv("JULIA_PROJECT" => "@.") do - Base.run(`./G4example.exe`).exitcode == 0 || error("Execution failed"); +withenv("JULIA_PROJECT" => "@.", "G4NUMTHREADS" => "4") do + Base.run(`./G4example.exe`).exitcode == 0 || error("Execution failed") end # ## Display the results println("=====> The file edepHist.png should have been saved") #md # ![](edepHist.png) +#nb # ![](edepHist.png) diff --git a/docs/src/examples/MyCode.jl b/docs/src/examples/MyCode.jl index 4bafe50..36ee8cb 100644 --- a/docs/src/examples/MyCode.jl +++ b/docs/src/examples/MyCode.jl @@ -1,27 +1,45 @@ #---Load the needed Julia modules------------------------------------------------------------------ - using Geant4 -using GeometryBasics using FHist using Plots +using Parameters println("=====> Loading MyCode.jl") -edepHist = H1D("Event total Edep distribution", 100, 0., 110.) -edep = 0.0 +#---Simulation data struct------------------------------------------------------------------------- +# define a mutable struct to store the simulation data (counters, hitograms, etc.) +# make N+1 instances of this struct, where N is the number of threads +@with_kw mutable struct MyData + edep = 0.0 + edepHist = H1D("Event total Edep distribution", 100, 0., 110.) +end +add!(d::MyData, d2::MyData) = (d.edep += d2.edep; merge!(d.edepHist, d2.edepHist)) + +const nthreads = ENV["G4NUMTHREADS"] == nothing ? 0 : parse(Int, ENV["G4NUMTHREADS"]) +const simdata = [MyData() for i in 1:nthreads+1] # #workers + 1(master) + +function getMyData() + tid = G4Threading!G4GetThreadId() + tid < 0 && (tid = -1) # master thread (-2 for without multi-threading support) + simdata[tid+2] +end +#---Actions---------------------------------------------------------------------------------------- function end_of_event_action(event) - push!(edepHist, edep) + data = getMyData() + push!(data.edepHist, data.edep) return # This is mandatory to force to return nothing end function begin_of_event_action(event) - global edep = 0.0 + data = getMyData() + data.edep = 0.0 return # This is mandatory to force to return nothing end function stepping_action(step) - global edep += GetTotalEnergyDeposit(step) + data = getMyData() + data.edep += GetTotalEnergyDeposit(step) return # This is mandatory to force to return nothing end @@ -29,11 +47,17 @@ function begin_of_run_action(run) return # This is mandatory to force to return nothing end -function end_of_run_action(run) - println("=====> End of run") - h = edepHist - img = plot(h.hist, title=h.title) - savefig(img, "edepHist.png") - println("=====> edepHist.png saved") +function end_of_run_action(run) + if G4Threading!G4GetThreadId() < 0 # only for the master thread + println("=====> End of run") + data = simdata[1] + for d in simdata[2:end] + add!(data, d) # merge all thread data to the master data + end + h = data.edepHist + img = plot(h.hist, title=h.title) + savefig(img, "edepHist.png") + println("=====> edepHist.png saved") + end return # This is mandatory to force to return nothing end