My program crashes in a multi-threaded run (1 000 particles) when using "TS_H_of_Water" for hydrogen. - Single-threaded run did not crash for 1 000 particles. - Not using any thermal scattering material id for 10 000 particles did not crash. ---------------------------------------------------------------------- In one crashing run, the two threads that crashed terminated like this: Thread 2: DEBUG: % Geant4 (WT2): % ********************************************************************************************************* DEBUG: % Geant4 (WT2): % * G4Track Information: Particle = neutron, Track ID = 1, Parent ID = 0 DEBUG: % Geant4 (WT2): % ********************************************************************************************************* DEBUG: % Geant4 (WT2): % DEBUG: % Geant4 (WT2): % Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName ... DEBUG: % Geant4 (WT2): % 406 224 -38.2 10.8 2.56e-08 0 2.58 594 Target hadElastic DEBUG: % Geant4 (WT2): % 407 223 -37.5 10.8 2.45e-08 0 0.834 595 Target hadElastic DEBUG: % Geant4 (WT2): % 408 221 -38.1 10.9 2.45e-08 0 2.09 597 Target Transportation DEBUG: % Geant4 (WT2): % 409 218 -39 11 2.45e-08 0 3.14 600 Target Transportation DEBUG: % Geant4 (WT2): % 410 218 -39.1 11 2.45e-08 0 0.185 600 Target Transportation DEBUG: % Geant4 (WT2): % 411 215 -40 11.1 2.45e-08 0 2.95 603 Target Transportation DEBUG: % Geant4 (WT2): % 412 212 -40.9 11.3 2.45e-08 0 3.14 606 Target Transportation DEBUG: % Geant4 (WT2): % 413 212 -41 11.3 3.12e-08 0 0.622 607 Target hadElastic DEBUG: % Geant4 (WT2): % G4ParticleChange::CheckIt : the Momentum Change is not unit vector !! Difference: 0.01495329191411354 DEBUG: % Geant4 (WT2): % neutron E=2.547093463277114e-08 pos=0.2119182362779898, -0.04109790054147729, 0.01127957004405689 DEBUG: % Geant4 (WT2): % ----------------------------------------------- DEBUG: % Geant4 (WT2): % G4ParticleChange Information DEBUG: % Geant4 (WT2): % ----------------------------------------------- DEBUG: % Geant4 (WT2): % # of 2ndaries : 0 DEBUG: % Geant4 (WT2): % ----------------------------------------------- DEBUG: % Geant4 (WT2): % Energy Deposit (MeV): 0 DEBUG: % Geant4 (WT2): % Non-ionizing Energy Deposit (MeV): 0 DEBUG: % Geant4 (WT2): % Track Status : Alive Thread 10: DEBUG: % Geant4 (WT10): % ********************************************************************************************************* DEBUG: % Geant4 (WT10): % * G4Track Information: Particle = neutron, Track ID = 1, Parent ID = 0 DEBUG: % Geant4 (WT10): % ********************************************************************************************************* DEBUG: % Geant4 (WT10): % DEBUG: % Geant4 (WT10): % Step# X(mm) Y(mm) Z(mm) KinE(MeV) dE(MeV) StepLeng TrackLeng NextVolume ProcName ... DEBUG: % Geant4 (WT10): % 173 297 -22.9 -40.1 4.29e-08 0 0.212 268 Target hadElastic DEBUG: % Geant4 (WT10): % 174 296 -22.9 -40.1 4.29e-08 0 0.0578 268 Target Transportation DEBUG: % Geant4 (WT10): % 175 294 -21.9 -40.1 3.12e-08 0 2.46 271 Target hadElastic DEBUG: % Geant4 (WT10): % G4ParticleChange::CheckIt : the Momentum Change is not unit vector !! Difference: 0.01495329191411354 DEBUG: % Geant4 (WT10): % neutron E=2.547093463277114e-08 pos=0.2941121789144017, -0.02295185791223434, -0.04006280754507685 DEBUG: % Geant4 (WT10): % ----------------------------------------------- DEBUG: % Geant4 (WT10): % G4ParticleChange Information DEBUG: % Geant4 (WT10): % ----------------------------------------------- DEBUG: % Geant4 (WT10): % # of 2ndaries : 0 DEBUG: % Geant4 (WT10): % ----------------------------------------------- DEBUG: % Geant4 (WT10): % Energy Deposit (MeV): 0 DEBUG: % Geant4 (WT10): % Non-ionizing Energy Deposit (MeV): 0 DEBUG: % Geant4 (WT10): % Track Status : Alive DEBUG: % Geant4 (WT10): % True Path Length (mm) : 1.04 Notice how both neutrons had the same energy and the momentum change difference is the same as well... It looks to me like some kind of race condition?
Never mind! In an attempt to simplify the thermal library classes, I ended up having only one model for all threads. The race condition happened when one thread modified one of the momentum direction components while another thread was copying it. I was wondering, it seems that using G4ThermalNeutrons is a much simpler and safer way to enable the HP stuff. I mean, it seems that all of this: G4HadronElasticProcess* theNeutronElasticProcess = new G4HadronElasticProcess; // Cross Section Data set G4NeutronHPElasticData* theHPElasticData = new G4NeutronHPElasticData; theNeutronElasticProcess->AddDataSet(theHPElasticData); G4NeutronHPThermalScatteringData* theHPThermalScatteringData = new G4NeutronHPThermalScatteringData; theNeutronElasticProcess->AddDataSet(theHPThermalScatteringData); // Models G4NeutronHPElastic* theNeutronElasticModel = new G4NeutronHPElastic; theNeutronElasticModel->SetMinEnergy(4.0*eV); theNeutronElasticProcess->RegisterMe(theNeutronElasticModel); G4NeutronHPThermalScattering* theNeutronThermalElasticModel = new G4NeutronHPThermalScattering; theNeutronThermalElasticModel->SetMaxEnergy(4.0*eV); theNeutronElasticProcess->RegisterMe(theNeutronThermalElasticModel); // Apply Processes to Process Manager of Neutron G4ProcessManager* pmanager = G4Neutron::Neutron()->GetProcessManager(); pmanager->AddDiscreteProcess(theNeutronElasticProcess); can be replaced with RegisterPhysics(new G4HadronElasticPhysicsHP); RegisterPhysics(new G4ThermalNeutrons); Is this correct? In such case, it would be worth updating the documentation in http://geant4-userdoc.web.cern.ch/geant4-userdoc/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html#hadronic-interactions ?
Yes, using the physics builders is the better way to do it, as long as they contain all the physics you need.
