Open DavidWalz opened 8 years ago
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 ?
I cannot reproduce this in a large simulation, so I assume that this is fixed by now.
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)
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,
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
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.
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'
The feature was added to CRPropa in October cc0baa12fafa07296cc4ef2a9e8e966b04e0eb93 .
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)
`
Sophia occasionally crashes with
sophiaevent: theta < -1.D0: -1.0000067840120266
reported by @DavidWit