Thermal neutron physics drastic change from 10.5 to 10.6

sergio · December 4, 2020, 11:38am

Quick summary: the behaviour of thermal neutrons in my simulation is very different when using Geant4 10.5 compared to Geant4 10.6.

I am simulating a monoenergetic 10 keV neutron source on a water phantom, and scoring the fluence of photons and thermal, epithermal and fast neutrons on a voxel grid. The following figures show the thermal and epithermal neutron fluence profiles respectively along the beam axis, normalized so that the maximum is 1. Results are provided for two fully independent code bases. A result calculated with MCNP used as reference is also shown.

l thermal_neutron_fluence_compare_geant4_versions

epithermal_neutron_fluence_compare_geant4_versions

Although it might not be obvious from the images, for epithermal neutrons all calculations agree, while for thermal neutrons all calculations using Geant4 10.6.0 or later have a different profile (also different to the MCNP reference).

All results have been calculated using G4ENDL-4.5. G4ENDL-4.6 has shown similar results when tested.

Any ideas? Is this a bug or a feature?

sergio · December 10, 2020, 2:39pm

I looked at the neutron cross-sections and they are the same, in case this helps troubleshooting.

civanch · December 11, 2020, 10:42pm

Hello,

I would advice to compared initialization printout for hadronic physics in Geant4 10.5 and 10.6 and be sure, that in 10.6 thermal neutrons are enabled. I mean printout at the beginning of simulation.

VI

sergio · December 14, 2020, 8:41am

Thank you Vladimir, good point!

This is from 10.5:

====================================================================
                  HADRONIC PROCESSES SUMMARY (verbose level 1)

---------------------------------------------------
                           Hadronic Processes for gamma

  Process: photonNuclear
        Model:            BertiniCascade: 0 eV  ---> 3.5 GeV
     Cr_sctns:            PhotoNuclearXS: 0 eV  ---> 100 TeV

---------------------------------------------------
                           Hadronic Processes for neutron

  Process: hadElastic
        Model:          NeutronHPElastic: 4 eV  ---> 20 MeV
        Model: NeutronHPThermalScattering: 0 eV  ---> 4 eV 
     Cr_sctns: NeutronHPThermalScatteringData: 0 eV  ---> 4 eV 
     Cr_sctns:        NeutronHPElasticXS: 0 eV  ---> 20 MeV
     Cr_sctns:            GheishaElastic: 0 eV  ---> 100 TeV

  Process: neutronInelastic
        Model:        NeutronHPInelastic: 0 eV  ---> 20 MeV
     Cr_sctns:      NeutronHPInelasticXS: 0 eV  ---> 20 MeV
     Cr_sctns:          GheishaInelastic: 0 eV  ---> 100 TeV

  Process: nCapture
        Model:          NeutronHPCapture: 0 eV  ---> 20 MeV
     Cr_sctns:        NeutronHPCaptureXS: 0 eV  ---> 20 MeV
     Cr_sctns:        G4NeutronCaptureXS: 0 eV  ---> 100 TeV

================================================================

And this is from 10.6:

====================================================================
                  HADRONIC PROCESSES SUMMARY (verbose level 1)

---------------------------------------------------
                           Hadronic Processes for gamma

  Process: photonNuclear
        Model:            BertiniCascade: 0 eV  ---> 3.5 GeV
     Cr_sctns:            PhotoNuclearXS: 0 eV  ---> 100 TeV

---------------------------------------------------
                           Hadronic Processes for neutron

  Process: hadElastic
        Model:          NeutronHPElastic: 4 eV  ---> 20 MeV
        Model: NeutronHPThermalScattering: 0 eV  ---> 4 eV 
     Cr_sctns: NeutronHPThermalScatteringData: 0 eV  ---> 4 eV 
     Cr_sctns:        NeutronHPElasticXS: 0 eV  ---> 20 MeV

  Process: neutronInelastic
        Model:        NeutronHPInelastic: 0 eV  ---> 20 MeV
     Cr_sctns:      NeutronHPInelasticXS: 0 eV  ---> 20 MeV

  Process: nCapture
        Model:          NeutronHPCapture: 0 eV  ---> 20 MeV
     Cr_sctns:        NeutronHPCaptureXS: 0 eV  ---> 20 MeV
     Cr_sctns:        G4NeutronCaptureXS: 0 eV  ---> 100 TeV

================================================================

The difference between the two seems to be that GheishaElastic and GheishaInelastic are gone. I have no idea what these are… Do you think that might be the cause? I did not do anything to remove those though.

lauriw · December 16, 2020, 9:01am

Hi,

I think I found the source of the problem.

Going from 10.5 to 10.6, the random rotation around the z-axis was removed in /source/processes/hadronic/processes/src/G4HadronElasticProcess.cc (i.e. “outdir.rotate(phi, it);” and “pdir.rotate(phi, it);” were removed prior to applying “rotateUz(indir)”). This removal seems to cause the neutrons to change direction only in the transverse plane, example (10 keV neutrons in water, thermal neutron recipe in use, G4_10.06.p02, G4NDL4.6):
20201215_G4_10.6.2_G4NDL4.6_longitudinalview

Running the same user code with and older G4 version (G4_10.05.p01, G4NDL4.6) the longitudinal scattering looks normal, i.e. not confined to the transverse plane.

When the azimuthal rotation described above is reintroduced into 10.6, the longitudinal scattering looks again normal, similar to 10.5.

With this type of bugfix, the G4_10.6.p02 and MCNP cell fluxes again agree nicely for both epithermal and thermal neutrons as well as photons as a function of depth along the center of the beam axis.

/LW

lauriw · December 16, 2020, 9:03am

The visualization with the same code in 10.06.p02, but reintroducing the azimuthal rotation looks like the following:

sergio · January 10, 2021, 2:05pm

I reported a bug:

This figure summarizes the problem that @lauriw’s figures revealed:
angle_z_histogram

The x-axis is output angle of neutrons in hadronic ellastic collisions with respect to the Z-axis. The y axis is probability density.

jvavrek · April 20, 2021, 3:06am

Can anyone confirm whether this fix was added to 10.7 or its patches? The missing phi rotation does not appear to be in the 10.7 source code nor in the 10.7.p1 patch version I downloaded.

sergio · April 20, 2021, 6:23am

It should be fixed, the location in the code of the phi rotation has changed. What happened was that the phi rotation was copied from the general hadron elastic code to each and every specific hadron elastic process (some processes in the future will do something different), but unfortunately thermal scattering was forgotten.

vasilis · June 1, 2021, 2:47pm

I have a similar problem with the geant4-10-07-patch-01 version.
Setup: orb (R=150cm) of G4_C, starting with neutrons of 1eV in the center.
Compared with FLUKA (pointwise + groupwise and Mcnp)
FLUKA+Mcnp agree well, while G4 is lower, not only for the total fluence
but also on the peak at kT=10meV instead of kT=25meV.

sergio · June 1, 2021, 3:04pm

Have you tested the behavior with 10.5 and 10.6? I assume that you have followed the instructions for low-energy neutron simulations:
https://geant4-userdoc.web.cern.ch/UsersGuides/ForApplicationDeveloper/html/TrackingAndPhysics/physicsProcess.html#high-precision-neutron-interactions-neutronhp

vasilis · June 2, 2021, 10:21am

10.6 gives exactly the same result. 10.5 I didn’t try since on my setup I am using the external navigator and was not present in the 10.5.
Just as a note: I am not using the S(a,b) on C, I am using it as amorphous for all setups

sergio · June 23, 2021, 12:03pm

If 10.6 and 10.7 give the same result, I think that you are having a different problem. But did you have all the low-energy neutron physics set up correctly?