Problem 2214

Summary:	TRACK003 "the Momentum Change is not unit vector" error in multi-threaded run using thermal scattering libraries
Product:	Geant4	Reporter:	Sergio Losilla <sergio.losilla>
Component:	processes/hadronic/models/neutron_hp	Assignee:	dennis.herbert.wright
Status:	RESOLVED INVALID
Severity:	normal
Priority:	P4
Version:	10.6
Hardware:	All
OS:	All

Description Sergio Losilla 2020-01-17 15:28:44 CET

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?

Comment 1 Sergio Losilla 2020-01-24 14:47:38 CET

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 ?

Comment 2 dennis.herbert.wright 2020-01-24 17:25:22 CET

Yes, using the physics builders is the better way to do it, as long as they contain all the physics you need.