sophiaevent: theta < -1

[DavidWalz]

Sophia occasionally crashes with
sophiaevent: theta < -1.D0: -1.0000067840120266
reported by @davidwit

[TobiasWinchen]

Line 16822 in https://github.com/CRPropa/CRPropa3/blob/master/libs/sophia/sophia_interface.f ; Doesn't the fix by Gero already prevent the crash ?

[TobiasWinchen]

I cannot reproduce this in a large simulation, so I assume that this is fixed by now.

[avvliet]

I'm actually still running into this crash. Although it doesn't actually crash the simulation, it does say:

sophiaevent: theta < -1.D0: 
input energy below threshold for photopion production !
sqrt(s) =   0.93826997280120850

This is an example script that, at least sometimes, gives this crash:

from crpropa import *

Run = 'Test'

neutrinos = True
photons = False

# module setup
#EBL = IRB_Gilmore12
EBL = IRB_Dominguez11
m = ModuleList()
m.add(SimplePropagation(1*kpc, 10*Mpc))
m.add(Redshift())
m.add(PhotoPionProduction(CMB, photons, neutrinos))
m.add(PhotoPionProduction(EBL, photons, neutrinos))
m.add(PhotoDisintegration(CMB, photons))
m.add(PhotoDisintegration(EBL, photons))
m.add(NuclearDecay(photons, photons, neutrinos))
m.add(ElectronPairProduction(CMB, photons))
m.add(ElectronPairProduction(EBL, photons))
m.add(MinimumEnergy(10.**10 * eV))
#m.add(PhotonOutput1D('photon-input'+Run+'.txt'))

# observer for hadrons
obs1 = Observer()
obs1.add(ObserverPoint())
obs1.add(ObserverNeutrinoVeto())  # we don't want neutrinos here
output1 = TextOutput('nucleons-'+Run+'.gz', Output.Event1D)
obs1.onDetection( output1 )
m.add(obs1)
# observer for neutrinos
obs2 = Observer()
obs2.add(ObserverPoint())
obs2.add(ObserverNucleusVeto())  # we don't want hadrons here
output2 = TextOutput('neutrinos-'+Run+'.gz', Output.Event1D)
obs2.onDetection( output2 )
m.add(obs2)

# source
source = Source()
source.add(SourceUniform1D(0, redshift2ComovingDistance(6)))
#source.add(SourceUniform1D(0. * Mpc, 1000.* Mpc))
source.add(SourceRedshift1D())
source.add(SourcePowerLawSpectrum(10.**16 * eV, 10.**23 * eV, -1.0))
source.add(SourceParticleType(nucleusId(56, 26)))

# run simulation for 1000 primaries and propagate all secondaries
m.setShowProgress(True)
m.run(source, 100000, True)

[saikatxdas]

Thanks for reporting this, I had some problems reproducing the behavior. Could you maybe provide an example that reproduces the behavior and also post the random seed that is stored in the output file (in recent CRPropa versions)? Also, which CRPropa version are you using? Which operating system? Which compiler? And what is very frequently? To keep the discussion organized I propose to continue in he thread for issue #85.

Hi Tobias,

1. Here is a minimal script, that gives me the error - sim.txt. How can I store the random seed in the output file?
2. I am using CRPropa3-3.1. The latest release version available on github.
3. The operating system is Ubuntu 16.04 (the error appears also in Ubuntu 18.04 in another machine)
4. f95, gcc, g++ (version 7.3.0) compilers
5. For simulation with < 1e6 injected primaries, I do not get the message. But beyond that, it occurs at least once. Currently, I couldn't reproduce this message with primaries of Z<14. For 3e6 injected particles I get the message at least 10 times.

[TobiasWinchen]

Thank you for the infos. You can enable the storage of the random seeds in TextOutput using

t = crpropa.TextOutput("out.dat")
t.enableRandomSeeds()

The textoutput will then contain lines such as

#
# CRPropa version: 3.1-309-gf59cdc73
#
# Random seeds:
#   Thread 0: KRoPs5j6bliu6lyuZY1i4C+3kUr3+zS .... the rest of the base64 encoded random seed

To ease reproducibility and debugging it is best to disable multi threading with
OMP_NUM_THREADS=1 python sim.py

[TobiasWinchen]

I can reproduce the issue with Saikat's code. The issue appears to be related to particles created in a large distance. Using source.add(SourceUniform1D(4500*Mpc, 5000*Mpc)) the issue occurs with increased frequency.

[saikatxdas]

Hi Tobias,

You are very much right. I checked and too found that with increased source distances, the message shows up more frequently. Thank you for pointing this out. Does this affect the simulation output in any way?

On using t.enableRandomSeeds() I get the following error
AttributeError: 'crpropa.TextOutput' object has no attribute 'enableRandomSeeds'

[TobiasWinchen]

The feature was added to CRPropa in October cc0baa1 .

[JonPaulLundquist]

This error still appears with CRPropa v3.1.6. My sources are up to 215 Mpc and I'm generating 4,000,000 nitrogen particles.

This is the end of the output of my log file. It appears that the simulation does not finish according to the progress bar though the number of nucleons and neutrinos in my output file might be about what I would expect? The code itself does not crash as the DintElecaPropagation code runs afterwards. I am not sure if all 4,000,000 are being simulated.

Started Wed Jul 14 01:17:32 2021 : [=====> ] 57% Finish in: 08:56:38 ^M2021-07-14 13:11:31 [ ERROR ] Something went wrong in the PhotoDisentigration Please report this error on https://github.com/CRPropa/CRPropa3/issues including your simulation setup and the following random seed: 0usob7MX+uZ0ZCmEC9+BRvMEQtR0T4TIqaoBycqpssq1Nh3I311hz17ki66AFKhsKgibHox/pfPoDHSCzYpnLWk9nqhKoSK/MLPCIcw+B5WcsR01SMgF12hQq8fh2hnbUNHZZj1VN/jgaPDTxJOfC8qconFT2+aOr4PZTTNpFDAiLaLDq3bbx2yz597PsB/tw4poVdSG7pWTOcZLW8x1c5ngz+IjfaL5Hy2wNWu+2P22eLtOHPvZUEnz/+J/leZpDHzg5LGnXKQdjp60aas4wX2fWRI63nNvl3tZljbo0i0gsPXKFn23kvVBvXz/NrYBZD0DlI/WgHUiTkMmm5xXCZJU/qcVYSAO13GN+DPaO7EL/UXS5OwjSKUGg7zbAqBgS8cUg5MmCarTYzWbEsuf37mmJDUbcCMkTWcGI5tBw7Yqv0WNjISVV946iSX999dLkYwTIgfwMw2bbz61fTGZpBP7cEczJluFH5Hyj6jhLZ35bncsdnGyRTMEO9BXMXYIwaTqmV75kN6ZWZSnqof8wG1tOUESj4+Dwmy9c+1JFvDcWqM604NPKAr8dlsdbnictEbzLvO/CJNMzXWhTjEJ0cshpW8ws6jjCIlfwSJgocrI7AjihgJGAbLrQEkQmOw7Ab/bH3vQY8UTY7fYxP00xQCd8awu5AblOPVPlaxYmiSe0vWGBo/qG+HezKNLq0TRc9Hj9/hN1y4nb+53wy/fNEpBy1FBq11m67GzdMlXeSfCW6aeGx1P9ddoMJUIZ9VwUG0hW+O+4lmT41zjtOieDyO0G78HZJDWuaIw3RxfVQRo522jMjwqM2HG6VVqPACjyOswA1My4R8uDXXPCze0l1EatR9Aepy0l1oWM7fI/M6vPF62xYg6tGdSKiUDj/H8Cm6l4oWTXTjbxbIuiq+A9ZldzOrpMEPy2WoF1banRbt45Z3ZfRzQ0oW87qFpx9TZUThwzGZHcQmEkQbD/c/w8+FeL1QGgSjfAasdMP9JNkVC1hTrFQ1izS92y4JNjjy8zXLieGa0cF0tsmJFL4cVMDrmIwyv8NzNGzyQDxxhvchn/b/86GH7l6+tiVuK41VOWlJc3119rwl+nNshdsPIEs9h/kbCYMZhtxzBjXHWZk0tPX30XFS1PK722nzGT++ARAPQiWBSV7rcD5Wgd4m4I3k6p9RYIVy8/pJcpGGYWJYSAR5zu5WUI7nCCmqq1NC4cXHJLb3lNKZr8knsdPXArggjIXbv86bo6KGvmi2ekm436l74TodB+H071yOmvTjRtMzajM3zNBjQxePNE7XmnM5A8cNsx8vLEKnaIVg3LDI5jAswxCHT6PXhuRV635PeKoDyXd+LWrNJu9JCsHzclPW4PimEQ+2bnKXEHsmS7lm3wAhNgAPKjk8orZf/bLydpOjNnSDGATsBYWsTk6CGOCmora4TA8SlaKZbyD/kF7cIUiLpjxcSmt2VeT5S6vyW9YIKGbNOkoAE7AymRdYFcetIjgjaRkHcpctbIdoc0gEy04U3nLA64YjP35WV3RZauhJkIeKpZObsS604y72agR8Pvq11/EZSMONSChB/ByiDDkTZ22dUxPfPMpgIsyJZ/hJqPG+fjTKM1sfNAbTGyXfcL+n+8b7C2BgOV8zvcFLli+TIiGIMGNdlSv4JOO3JHY/WayN3G8BBjXg14O0SvKWZl25fkTzB9nRge47ZwMgFgzzzkQ/vStitdOJBEnQHCHfng0C9xnfkO1qFVaKbsUIAPbupYWoIQSL3jHhefd7Qy3wWIK7jgUrSbHGFDEzo2j2+3k6cdweu/d8xfJULqeZvvmL6hNgKosqHMmKfTFvcT/4IigSUYpWNUYNoY37A4u7tH4XXirq1Xo37LA9fnp6cvAXuGt/Zp29d//2BV5km7cZkchcAuVXD9BG3ecorl0Vglq+mEIRjef69IJMbJ6/DXPcc5+5iW8opaHTDt8sZELzNCDfc5+va/ZGZNR880/BV4fMERFqHIjA3cwUJfaECR/mFhgMJ93Ny5d47pd4YmqEYcLxodku2t7yxxNQi6/OgYu/pujE+GB6wGYoMmS2xjIQwIS+ldBMBoiYUTwGIql1cjrVqpLPqAWQeQwF605inBXYnIHA5r6WUcLKfDTegxU/pi7QRXeYoLuhxLU0OsAmME4YovZp42YDuyOXOxjeKSwe3qeN/eoJrWi8cxsYn7FGk845fOYkOqT6UOXuGMvGWEEY/5/a9z3XZDDzJJZaBRaXObvP1ROo+KfhbxB4q5w0L15IIqLylbX1r6V4/+9Jjl1kjqBj6CaNqP0NPZuuzizndV6wTNnEMhE7RIt2bOWwL3PwwLGhTK4XTi3WOwSRKo2FbnwSKNxanDJjoAwfdA9uA9TMZuCJyxACRdt0OTGDE7LquDFBxh13DrlO6Hip1OztNLk3uTJYSvzukL3wgtKEIE7UUxuModmA+LpV5Lr0/zHhgbxVVjfBF18/yXiZBnIFS9mcovwEYW0Y/IWGx+rksQHd9sEakdS0NeZwOaB0HusWKsPi++ad3GxFSSN8aiirXymjwvu4bdcSgXiGQrA9d+D9gMhJGCYcUwqV3TpDiYvwTrKGITeJMOt5J/7veZWrpc7lIbHG0v0jpULz7QXXkeDJ3DWOCP5D7Cmo/jMoAPw+c+xr6Igjj4ztu07DRtjRQOnGh2hNMQLDggkX008+fAFeX1D/VziLVPCsCQLVVezH7xL0qB6tYclehYZEB4BjM0Q9HGaVxTjuOnlJXiqnEzxPbM97USBqylNniukmKsaFitblVoC9+e/xaCvxsEueTY3TXnV9Hy4WoHZQxXdRn9e9Biw7HtRtF4dYHvzG7243S0RYBfUtscEOb8H4fCK0JrortvZPMd9VLKf8Gch/Vb51ZzqtMTB/8O1w8+grOdO7PeW8I6RWxgxyVpbtDEPCl+MjKzeFmUcxR4zDBtsbWBW4C9DW//dgibjdvpgYr5b1C0dkeaG2EvKI+44pApa0fDFguWSFKEm3U78ZR5t5yKeA02gmaaRZAJUQ0TGlk25gaa5dycgrAzQGfgOoCzgexuJSsLbKlPEQn3jCq5CNuXMCj5SgJG2u8KA7hcqjFefU2QhLQi7Y4RLpfe02Yao/Rae+Bpgqa3ILCTZYTTt2sMbYmGBhNSmUSIMTz0rBGDwlS7SQwOgIS43e+O13j6jyImLfw1e6j9WOVAchsZNAFdoePr+ITgMHDliqQbOEgwCrGcF2pXjHjsp5WcaaZnLPECV5Y1eoIZTJ2V/bCrG5KdaHwM0zk6sUD6YjmIn6tcYdhImoZYz1fmwSZSDopy8NgHhRD4RZGS9Tg \n Exception in crpropa::ModuleList::run: crpropa::Nucleus: no nucleus with Z < 0, A=-54 Z=-23 Run EleCa propagation Started Wed Jul 14 13:11:42 2021 : [=> ] 99% Finish in: 00:00:00 ^M Started Wed Jul 14 13:11:42 2021 : [^[[1;32m Finished ^[[0m] 100% Needed: 00:00:00 - Finished at Wed Jul 14 13:11:42 2021 sophiaevent: theta < -1.D0: -1.0000035726845624 input energy below threshold for photopion production ! sqrt(s) = 1.0791663303040897

This is the CRPropa component of my code:

`
obsPosition = Vector3d(0,0,0)
vgrid = Grid3f(obsPosition, 1024, 30 * kpc)
initTurbulence(vgrid, 1.09 * nG, 60 * kpc, 1000 * kpc, -11. / 3, 42)
Bfield = MagneticFieldGrid(vgrid)
aveField = (meanFieldStrength(vgrid) / nG) * nG

minStep = 15.*kpc                        ## minimum length of single step of calculation
maxStep = 4.*Mpc                         ## maximum length of single step of calculation
tolerance = 1e-2                         ## tolerance for error in iterative calculation of propagation step
    
sim = ModuleList()

neutrinos = True
photons = True
electrons = False
sim.add(PropagationCK( Bfield, tolerance, minStep, maxStep))
sim.add( Redshift() )
sim.add(PhotoPionProduction(CMB, photons, neutrinos, electrons))
sim.add(PhotoPionProduction(IRB_Gilmore12, photons, neutrinos, electrons))
sim.add(PhotoDisintegration(CMB)) 
sim.add(PhotoDisintegration(IRB_Gilmore12))
sim.add( ElectronPairProduction(CMB, electrons) )
sim.add( ElectronPairProduction(IRB_Gilmore12, electrons) )
sim.add(NuclearDecay(electrons, photons, neutrinos))
sim.add(EMPairProduction(CMB, electrons, thin))
sim.add(EMPairProduction(IRB_Gilmore12, electrons, thin))
sim.add(EMPairProduction(URB_Protheroe96, electrons, thin))
sim.add(EMDoublePairProduction(CMB, electrons, thin))
sim.add(EMDoublePairProduction(IRB_Gilmore12, electrons, thin))
sim.add(EMDoublePairProduction(URB_Protheroe96, electrons, thin))
sim.add(EMInverseComptonScattering(IRB_Gilmore12, photons, thin)) 
sim.add(EMInverseComptonScattering(CMB, photons, thin))
sim.add(EMInverseComptonScattering(URB_Protheroe96, photons, thin))
sim.add(EMTripletPairProduction(CMB, electrons, thin))
sim.add(EMTripletPairProduction(IRB_Gilmore12, electrons, thin))
sim.add(EMTripletPairProduction(URB_Protheroe96, electrons, thin))

minE = MinimumEnergyPerParticleId(3.98107*EeV)
minE.add(12, 100*TeV)
minE.add(-12, 100*TeV)
minE.add(14, 100*TeV)
minE.add(-14, 100*TeV)
minPhotonE = 100*MeV
minE.add(22, minPhotonE)
minE.add(11, minPhotonE*10)
minE.add(-11, minPhotonE*10)
sim.add(minE)

MaxTraj = MaximumTrajectoryLength( 2000*Mpc )
MaxTraj.addObserverPosition(obsPosition)
sim.add( MaxTraj )
sim.add( SphericalBoundary(obsPosition, 2000*Mpc) )

obs_nucleons = Observer()
obs_neutrinos = Observer()
obs_photons = Observer()

obsSize = 200*kpc  ## physical size of observer sphere

output_nucleons = TextOutput('FR0nucleons.txt', Output.Event3D)
output_neutrinos = TextOutput('FR0neutrinos.txt', Output.Event3D)
output_photons = TextOutput('FR0photons_init.txt', Output.Event1D)
output_photons.enable(Output.CreatedIdColumn) 
output_photons.enable(Output.CreatedEnergyColumn)
output_photons.enable(Output.CreatedPositionColumn)

obs_nucleons.add( ObserverSurface(Sphere( obsPosition, obsSize )) )
obs_nucleons.add(ObserverInactiveVeto())
obs_nucleons.add(ObserverElectronVeto())
obs_nucleons.add(ObserverPhotonVeto())
obs_nucleons.add(ObserverNeutrinoVeto())
obs_nucleons.onDetection(output_nucleons)
obs_nucleons.setDeactivateOnDetection(True)

obs_neutrinos.add(ObserverInactiveVeto())
obs_neutrinos.add(ObserverElectronVeto())
obs_neutrinos.add(ObserverPhotonVeto())
obs_neutrinos.add(ObserverNucleusVeto())
intObs = ObserverPathIntersect(obsPosition, obsSize)
intObs.addId(12)
intObs.addId(-12)
intObs.addId(14)
intObs.addId(-14)
intObs.addId(16)
intObs.addId(-16)
obs_neutrinos.add(intObs)
obs_neutrinos.onDetection(output_neutrinos)
obs_neutrinos.setDeactivateOnDetection(True)

obs_photons.add(ObserverInactiveVeto())
obs_photons.add(ObserverElectronVeto())
obs_photons.add(ObserverNeutrinoVeto())
obs_photons.add(ObserverNucleusVeto())
intObs = ObserverPathIntersect(obsPosition, obsSize)
intObs.addId(22)
obs_photons.add(intObs)
obs_photons.onDetection(output_photons)
obs_photons.setDeactivateOnDetection(True)

sim.add(obs_nucleons)
sim.add(obs_neutrinos)
sim.add(obs_photons)

#source information D2[i], theta[i], phi[i], sources, weight2 created here

sources = []

for i in range(0,D2.size):
    s = Source()
    direction = Vector3d()
    direction.setRThetaPhi(D2[i], theta[i], phi[i])
    s.add(SourcePosition(direction))
    s.add(SourceParticleType(nucleusId(14, 7)))
    s.add(SourcePowerLawSpectrum(3.98107*EeV, 1000 * EeV, -1))
    s.add(SourceIsotropicEmission())
    s.add(SourceRedshift1D())
    sources.append(s)

source_list = SourceList()
Ns = len(sources)
j = 0
for s in sources:
    source_list.add(s,weight2[j])
    j = j + 1

sim.setShowProgress(True)  # switch on the progress bar
sim.run(source_list, 4000000) #proton and nitrogen

output_nucleons.close()
output_neutrinos.close()
output_photons.close()

DintElecaPropagation('FR0photons_init.txt', 'FR0photons.txt', True, 0.0001*EeV, aveField)

`