possible issue with gravity tree and non-leapfrog-integrators

[birnstiel]

As I ran into issues with a more complicated setup, I tried the very basic setup below, setting the sim.gravity to tree or basic and trying the integrators ias15, whfast, and leapfrog. While they all work with basic gravity, tree-gravity causes issues with ias15 and whfast (also: hermes crashes, but I know it's still in development). To summarize:

	basic	tree
ias15	ok	different
whfast	ok	very different
hermes	ok	Kernel restarts
leapfrog	ok	ok

Is that to be expected? Or does it have to do with the low number of particles? Still new to this -- sorry if it's an obvious thing.

I'm using rebound 3.5.2.

import rebound

sim = rebound.Simulation()

sim.integrator          = 'whfast'
sim.units               = ('Jyr','AU','Msun')
sim.dt                  = 1e-3
sim.collision_resolve   = 'merge'
sim.boundary            = 'open'
sim.gravity             = 'tree'
sim.collision_resolve_keep_sorted = 1
sim.configure_box(20.)

sim.add(m=1)

sim.add(m=1e-5,a=1)
sim.add(m=1e-5,a=4)

sim.move_to_com()

# integrate and plot, inner planet should be at 0, outer at 180 degree

sim.integrate(100.)

print('integrator = {}'.format(sim.integrator))
print('gravity    = {}'.format(sim.gravity))
rebound.OrbitPlot(sim,color=True)

[hannorein]

This is expected. The tree code is an approximate gravity solver. It is only useful when an approximate gravitational force is good enough for the problem at hand. If that is the case, then you also want to use a simple and fast integrator (leapfrog). I cannot think of a case where the combination of an approximate gravity solver and a specialized high accuracy integrator (IAS15, WHFast or Hermes) would be useful.

But maybe you have a use case that I was not thinking about. What do you want to use the tree code for?

[birnstiel]

It's a self-stirring swarm of planetesimals or protoplaneta where there should be close encounters. it's just for a lecture demonstration, but would leapfrog be a fair choice for such a scenario? I thought for collisions frequencies getting the encounters right might matter.

[hannorein]

No, that would not be a good choice. If the gravitational force for the planetesimals is only approximate, they will spiral into the star over secular timescales, hitting the planet in the process, and thus leading to un unphysical migration. Similarly, the planet will drift because of approximate gravitational forces from the planetesimals. What you can do is ignore the planetesimal-planetesimal interactions all together, then you have an O(N) scaling (set the Nactive variable for that). This would be a self-consistent and well posed problem. Having a self-gravitational disk of planetesimals and a planet all together is a very difficult task. REBOUND doesn't have an out-of-the box solution for that, unless you use the slow O(N^2) gravity solver.

[r-zolotarev]

Dear developers team, possibly I have the same issue when test WHFast integrator with tree gravity. I tried to test integrate simple two-body problem and with basic gravity everything is OK, leapfrog with tree gravity is OK, but WHFast with tree gravity shows different behavior, even if I set opening angle to zero in params. In my expectation Barnes-Hut algorithm should provide the same gravity as direct sum when set opening angle to zero or in case of two bodies. Is it possible to setup WHFast to work with tree? (I would like to try rebound with WHFast and BH-tree in future project)

I am new here, sorry if my question is incorrect

Thank you! Roman Zolotarev

here is the text of my setup in problem.c and some figures int main(int argc, char argv[]){ struct reb_simulation const r = reb_create_simulation(); r->integrator = REB_INTEGRATOR_WHFAST; r->gravity = REB_GRAVITY_TREE; r->boundary = REB_BOUNDARY_OPEN; r->opening_angle2 = 0.0; r->G = 1;
r->softening = 1.0e-2; r->dt = 1.0e-2; const double boxsize = 20.2; reb_configure_box(r,boxsize,1,1,1); struct reb_particle star = {0}; star.m = 1; star.x = 0.0; star.y = 0.0; star.z = 0.0; star.vx = 0.0; star.vy = 0.0; star.vz = 0.0; struct reb_particle pt = {0}; pt.m = 1.0e-3; pt.x = 5.0; pt.y = 0.0; pt.z = 0.0; pt.vx = 0.0; pt.vy = 1./sqrt(5.0); pt.vz = 0.0; reb_add(r, star); reb_add(r, pt); r->heartbeat = heartbeat; double t_fin = 1.0e2 2.0 M_PI; reb_integrate(r, t_fin); reb_free_simulation(r); }

[hannorein]

Hi Roman. As I mentioned earlier in the thread, I just don't see how using WHFast together with a tree code could possibly make sense. I'd encourage you to think about another solution (i.e. using active and non-active particles).

[r-zolotarev]

Thank you! I will try with test particles

Roman Zolotarev

[hannorein]

I'll close this for now, but feel free to reopen the issue!