Fatal Error trying Poincare Map
Bruntonian opened this issue · 11 comments
Hello,
I recently wanted to try creating a Poincare Map and ran into issues. I then attempted to follow the Poincare Map example verbatim, and got the following error(s) upon running this line:
Nsim = 10
with Pool() as pool:
res = pool.map(runone,range(Nsim))
Fatal error! Exiting now. Variational gravity calculation not yet implemented.
Fatal error! Exiting now. Variational gravity calculation not yet implemented.
Fatal error! Exiting now. Variational gravity calculation not yet implemented.
Fatal error! Exiting now. Variational gravity calculation not yet implemented.
---------------------------------------------------------------------------
KeyboardInterrupt Traceback (most recent call last)
Cell In[15], line 3
1 Nsim = 10
2 with Pool() as pool:
----> 3 res = pool.map(runone,range(Nsim))
File /Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/multiprocess/pool.py:367, in Pool.map(self, func, iterable, chunksize)
362 def map(self, func, iterable, chunksize=None):
363 '''
364 Apply `func` to each element in `iterable`, collecting the results
365 in a list that is returned.
366 '''
--> 367 return self._map_async(func, iterable, mapstar, chunksize).get()
File /Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/multiprocess/pool.py:768, in ApplyResult.get(self, timeout)
767 def get(self, timeout=None):
--> 768 self.wait(timeout)
769 if not self.ready():
770 raise TimeoutError
File /Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/site-packages/multiprocess/pool.py:765, in ApplyResult.wait(self, timeout)
764 def wait(self, timeout=None):
--> 765 self._event.wait(timeout)
File /Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/threading.py:655, in Event.wait(self, timeout)
653 signaled = self._flag
654 if not signaled:
--> 655 signaled = self._cond.wait(timeout)
656 return signaled
File /Library/Frameworks/Python.framework/Versions/3.12/lib/python3.12/threading.py:355, in Condition.wait(self, timeout)
353 try: # restore state no matter what (e.g., KeyboardInterrupt)
354 if timeout is None:
--> 355 waiter.acquire()
356 gotit = True
357 else:
KeyboardInterrupt:
Rebound Version: 4.4.1
Python 3.12.2
macOS sonoma 14.6.1
Any help would be greatly appreciated.
Thanks,
Ian
Hi Ian. That's strange. I've tried it just now (latest REBOUND version, macOS sonoma 14.6.1, but Python 3.9.6). It works for me. It looks like there was a Keyboard Interrupt, If you didn't do that then I'm not sure what could be the problem. Try running the example without the multiprocessing part, just by doing it in serial. It'll be a bit slower, but it should help to narrow down the problem.
Hi Hanno, thanks for the quick reply. I did a keyboard interrupt because I was getting the fatal error message one-by-one so stopped it. Like below:
Just also tried without the multiprocessing part (I think I did it correctly) and my kernel died and had to restart, so not really sure.
Thanks,
-Ian
Hm. And just to confirm, you didn’t make any changes to the notebook?
Yup, I just copy-pasted from the example page to start.
Did you install Rebound from source? If so it could be that you did an upgrade and still need to do another pip install.
I did, have since tried reinstalling but still no luck so I'm wondering it may just be something else on my end. I'm going to keep trying some different things and will update if I find a solution.
Strange. Maybe you can post the full code here?
Sorry do you mean the full code for the Poincare Map? Because for right now I have just copied and pasted verbatim from the Poincare Map Example on the rebound site:
import rebound
import numpy as np
import warnings
sim = rebound.Simulation()
sim.add(m=1.)
sim.add(m=1e-3,a=1,e=0.001)
sim.add(m=0.,a=1.65)
sim.move_to_com()
def migrationForce(reb_sim):
tau = 40000.
ps[2].ax -= ps[2].vx/tau
ps[2].ay -= ps[2].vy/tau
ps[2].az -= ps[2].vz/tau
sim.additional_forces = migrationForce
ps = sim.particles
sim.integrate(3000.)
sim.save_to_file("resonant_system.bin")
def hyper(sim):
dp = sim.particles[2]-sim.particles[0]
return dp.x*dp.vx + dp.y*dp.vy
def mod2pi(x):
if x>np.pi:
return mod2pi(x-2.*np.pi)
if x<-np.pi:
return mod2pi(x+2.*np.pi)
return x
def runone(args):
i = args # integer numbering the run
N_points_max = 2000 # maximum number of point in our Poincare Section
N_points = 0
poincare_map = np.zeros((N_points_max,2))
# setting up simulation from binary file
import warnings # ignore warning for function pointers
warnings.filterwarnings('ignore')
sim = rebound.Simulation("resonant_system.bin")
vx = 0.97+0.06*(float(i)/float(Nsim))
sim.particles[2].vx *= vx
sim.t = 0. # reset time to 0
# Integrate simulation in small intervals
# After each interval check if we crossed the
# hypersurface. If so, bisect until we hit the
# hypersurface exactly up to a precision
# of dt_epsilon
dt = 0.13
dt_epsilon = 0.001
sign = hyper(sim)
while sim.t<15000. and N_points < N_points_max:
oldt = sim.t
olddt = sim.dt
sim.integrate(oldt+dt)
nsign = hyper(sim)
if sign*nsign < 0.:
# Hyper surface crossed.
leftt = oldt
rightt = sim.t
sim.dt = -olddt
while (rightt-leftt > dt_epsilon):
# Bisection.
midt = (leftt+rightt)/2.
sim.integrate(midt)
msign = hyper(sim)
if msign*sign > 0.:
leftt = midt
sim.dt = 0.3*olddt
else:
rightt = midt
sim.dt = -0.3*olddt
# Hyper surface found up to precision of dt_epsilon.
# Calculate orbital elements
o = sim.orbits()
# Check if we cross hypersurface in one direction or the other.
if o[1].d<o[1].a:
# Calculate resonant angle phi and its time derivative
tp = np.pi*2.
phi = mod2pi(o[0].l-2.*o[1].l+o[1].omega+o[1].Omega)
phid = (tp/o[0].P-2.*tp/o[1].P)/(tp/o[0].P)
# Store value for map
poincare_map[N_points] = [phi,phid]
N_points += 1
sim.dt = olddt
sim.integrate(oldt+dt)
sign = nsign
# Rerun to calculate Megno
sim = rebound.Simulation("resonant_system.bin")
vx = 0.97+0.06*(float(i)/float(Nsim))
sim.particles[2].vx *= vx
sim.t = 0. # reset time to 0
sim.init_megno() # adds variational particles and initialized MEGNO
sim.integrate(15000.)
return (poincare_map, sim.megno(),vx)
from multiprocess import Pool
Nsim = 10
with Pool() as pool:
res = pool.map(runone,range(Nsim))
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(14,8))
ax = plt.subplot(111)
ax.set_xlabel("$\phi$"); ax.set_ylabel("$\dot{\phi}$")
ax.set_xlim([-np.pi,np.pi]); ax.set_ylim([-0.06,0.1])
for m, megno, vx in res:
c = np.empty(len(m[:,0])); c.fill(megno)
p = ax.scatter(m[:,0],m[:,1],marker=".",c=c, vmin=1.4, vmax=3, s=25,edgecolor='none', cmap="brg")
cb = plt.colorbar(p, ax=ax)
cb.set_label("MEGNO $<Y>$")
Was just trying to replicate that.
Thanks. I wast just trying to make sure nothing got accidentally changed. This works for me. I suspect there is something wrong with the installation on your end. All I can suggest is:
- create a new virtual environment
- install the latest version of rebound in the new virtual environment
Ok, thanks Hanno. So I just reinstalled everything from the ground up and restarted my computer this time and miraculously it now works. Not entirely sure what I did different than before, but alas, it is now working. Sorry for the trouble and thanks for the help.
No worries and glad it works now! It can be tricky sometimes to figure out which package python is actually using... (been there myself!)