search

CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
70 stars 68 forks source link

Nucleus with Z < 0 produced #188

Open saikatxdas opened 6 years ago

saikatxdas commented 6 years ago

I am running a 1D simulation with injected Helium particles from a source redshift evolution and with relevant interactions for energy loss. The simulation fails with following Exception error saying that particle with negative atomic and mass number is produced. 

crpropa::ModuleList: Number of Threads: 8
Run ModuleList
Exception in crpropa::ModuleList::run: =>       ]  25%    Finish in: 01:49:48 
crpropa::Nucleus: no nucleus with Z < 0, A=-30 Z=-14

Can someone please tell me, what are the possible causes for this?

rafaelab commented 6 years ago

I have never encountered this problem. Could you please provide a minimum working example that reproduces it?

saikatxdas commented 6 years ago

This error happened in a simple 1D simulation. It never occurred before on this machine. On running again, the simulation finished successfully. I don't know, whether this can be reproduced, but here is the minimal script that I have been using: script.txt

marcus-wirtz-snkeos commented 5 years ago

The script runs also completely fine for me. I guess this issue can be closed?

saikatxdas commented 5 years ago

This happened only once. I didn't get the error again.

TobiasWinchen commented 5 years ago

We have two other issues #85 and #99 where on rare occasions CRPropa crashes in interactions due to something strange happening with the secondaries. Since recently the random number seeds are stored which may help to get to the bottom of these issues. Personally I think these issue are not resolved as we have no idea what happened and thus should not be closed.

saikatxdas commented 5 years ago

I encounter the issue mentioned in #85 very frequently. For elements heavier than Nitrogen, there is this error:

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

Is there a way this can be stopped?

TobiasWinchen commented 5 years ago

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.

adundovi commented 5 years ago

@Froehliche-Kernschmelze this could be related to SOPHIA and #225.

lukasmerten commented 3 years ago

Since we did not have any activity since 11/2018 on the original issue (production of nuclei with Z<0) and nobody was able to reproduce this error, I guess we can close this issue for now.

JonPaulLundquist commented 3 years ago

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. The code itself does not crash as the DintElecaPropagation code runs afterwards.

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)

`