Difference in evolve_model between gadget2 and bridge

[fredt00]

Hi,

I have a problem that I think is related to issue #679. I am trying to couple a live NFW halo in gadget2 with an orbiting test mass, however I have noticed that no matter how long I pre-evolve the halo to allow it to reach equilibrium, as soon as I add it to bridge, the halo rapidly expands. I have tried the suggestions made in #679, in particular ensuring that max_size_timestep in gadget2 is a few factors of two less than the bridge timestep (dt), and using a multiple of two times the bridge timestep in the unit converter, but the issue persists.

Attached is an example script that demonstrates the problem by evolving the same halo first by calling evolve_model every bridge timestep, and then evolving with a bridge integrator. The resulting evolution is dramatically different. I ran the same test with BHTree and in this case the evolution is the same whether or not I call evolve_model directly or via bridge. The evolution with and without bridge does match for gadget2 if I add it to the bridge with add_code instead of add_system.

gadget_testing.py.txt

I'm also confused by the output file from gadget2, with every multiple of dt/max_size_timestep step having a systemstep=0, yet the time still increases:

Step= 93  t= 91  dt= 1 
Nf= 0000100000  total-Nf= 0011700000  ex-frac= 8.56557  iter= 1
work-load balance: 1.30849  max=0.104356 avg=0.079753 PE0=0.0842994
particle-load balance: 1.3368
max. nodes: 5556, filled: 0.694569
part/sec=83591.5 | 63883.9  ia/part=527.471 (0)

Step= 94  t= 92  dt= 0 
Nf= 0000100000  total-Nf= 0011800000  ex-frac= 14  iter= 1
work-load balance: 1.01866  max=9.64458 avg=9.46793 PE0=9.31209
particle-load balance: 1.0107
max. nodes: 4564, filled: 0.570557
part/sec=704.132 | 691.234  ia/part=100000 (0)

Step= 94  t= 92  dt= 0 
Nf= 0000100000  total-Nf= 0011900000  ex-frac= 8.62274  iter= 1
work-load balance: 1.39914  max=0.109665 avg=0.0783807 PE0=0.0941227
particle-load balance: 1.0107
max. nodes: 4564, filled: 0.570557
part/sec=85055 | 60790.9  ia/part=527.238 (0)

Step= 95  t= 93  dt= 1 
Nf= 0000100000  total-Nf= 0012000000  ex-frac= 8.56177  iter= 1
work-load balance: 1.30747  max=0.104649 avg=0.0800395 PE0=0.0844096
particle-load balance: 1.33185
max. nodes: 5513, filled: 0.689194
part/sec=83292.2 | 63704.8  ia/part=527.38 (0)

Am I implementing gadget2 incorrectly?

Thanks!
Fred

[fredt00]

I've narrowed down the problem to line 464 in bridge:
channel.copy_attributes(["vx","vy","vz"])
and more generally, copying to gadget2 between evolve_model calls seems to cause rounding errors in the timestep. The script test.py.txt demonstrates this with the output file for the case where data is copied to gadget2 showing the lines:

Begin Step 1, Time: 1, Systemstep: 0
Begin Step 2, Time: 2, Systemstep: 0
Begin Step 3, Time: 2.99994, Systemstep: 0
Begin Step 4, Time: 3.99994, Systemstep: 1
Begin Step 5, Time: 4, Systemstep: 0
Begin Step 6, Time: 5, Systemstep: 0

The two simulations diverge after step 5. Adding print(gravity.get_time_step().in_(units.Myr)) in the script shows a timestep of 0.999938964844 Myr followed by 6.103515625e-05 Myr at 4Myr and 5Myr respectively, yet the model_time still increases by 1Myr. I've tried different multiples of dt in the unit converter and parameters.max_size_timestep but the evolution never matches the case where no data is copied from the interface to gadget2.