Added MT to JuliaAction example

docs/src/examples/G4example.cpp

@@ -14,6 +14,7 @@
 //#include "G4SystemOfUnits.hh"
 #include "globals.hh"
 #include "G4VUserActionInitialization.hh"
+#include "G4UserWorkerInitialization.hh"
 #include <julia.h>
 #include <iostream>
 
@@ -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,23 +143,35 @@ 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\")");
     if (jl_exception_occurred()) {
         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. 

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)

docs/src/examples/MyCode.jl

@@ -1,39 +1,63 @@
 #---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
 
 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