testparticle question?

[kaiwu-astro]

In the simple 2-body simulation, rebound gives very different result when I set the small mass companion as testparticle.

A central body with a small mass companion

import rebound
s1 = rebound.Simulation()
s1.add(m=1.)
s1.add(m=1.0e-10, a=10)

s1.integrate(100)
print(s1.particles[1].a)

it prints 10.0, which is good

then I set the companion as a testparticle, but give it a huge mass.

s2 = rebound.Simulation()
s2.add(m=1.)
s2.add(m=0.1, a=10)

s2.testparticle_type = 0
s2.N_active = 1

s2.integrate(100)
print(s2.particles[1].a)

it gives 10.3302539, which is 3.3% more

I was expecting the same or very similar result, because the central body does not feel the testparticle, and the companion feels the same acceleration. Actually in the 2nd simulation, data of the companion (a, x, y, z) changes a lot when I change its mass (1e-1, 1e-5, 1e-10).

In practice, I can set the mass of the testparticle to a very small value (1e-10) to get the same result. I just wonder why would this happen. Maybe it would cause potential problems when users forget to set the testparticle mass to very small.

[Findus23]

This might be a bit related to #505 where I accidentally set the wrong testparticle type and had massive testparticles

[hannorein]

If you want to have a massless particle, you need to set the mass to zero. N_active exists to speed up the calculation with many massless particles.

For everything else, e.g. "semi-active particles", it's really important to understand the details because most likely the simulation is always "unphysical" in one way or the other.

One could add a warning message for incorrect use (testparticles with mass), but this would require REBOUND to check all the particles every timestep. I'm not sure that's worth it.

Also note that calculating orbital elements might depend on the masses of particles but the actual integration might not. If you're trying to understand what's going on, it might be helpful to look at the cartesian coordinates, not only the orbital parameters.

[Findus23]

Maybe one could just add a single check at the start of the simulation. It wouldn't add much overhead, but might catch a large fraction of people accidently misusing rebound?

[hannorein]

I'm a bit hesitant because there are cases where one might want this behaviour. For example, I once ran a simulation of Saturn rings where each particle had a mass which was used to resolve collisions. However, N_active was set to 1 so that only the gravity from a moonlet was calculated.

[hannorein]

Maybe improving the documentation should be the first goal...

[kaiwu-astro]

Thanks for replying. Besides updating the documentation, I would like to suggest adding a check at the beginning of the simulation as well, especially when a user use testparticle type 0 and assign non-zero mass, and issue a warning, (maybe with an option to suppress the warning), so misusing users can be warned, and those who do this by purpose will not be influenced.

[hannorein]

Fair enough. Looks like I got overruled ;-)

[mbeutel]

Having run into the "At least one test particle (type 0) has finite mass" warning, I'll follow up on this issue with a question:

From @hannorein's explanations, my understanding is that giving type-0 test particles a mass should not affect the results of the integration, but orbital elements may turn out differently because the coordinate conversion depends on the particle mass. In other words, for particles that have a mass but whose gravitational impact is to be kept out of the N-body calculation (e.g. because it is accounted for by different means), using type-0 test particles with finite masses is the right thing to do. Is that correct?

I'm glad that the warning exists; but perhaps the documentation could explain why the behavior might be unexpected, and which behavior is expected for massive type-0 test particles. I feel bad about suppressing the warning if I don't know what exactly it warns about.

[hannorein]

Hi Moritz!

The code tries to optimize various operations if you turn the test-particle setting on. Ignoring optimizations, you should get the same results if you just have normal particles and set the mass to zero.

If you give test particles a mass, the coordinate conversion are definitely a potentially problematic area as you point out. To be honest, I don't know out of my head which other parts of the code would be affected by this (maybe someone else remembers scenarios where this would be an issue). The best I can do is call it "undefined behaviour". I guess the question is why do you give your particles a mass? If you just use it in your own parts of the code (e.g. to identify various species of particles), then it might be best to store that data somewhere else, for example by using the particles' hash or ap properties. REBOUNDx can also help with that sort of thing.

It's probably not the answer you were hoping for, but I hope it helps. Feel free to reach out again if anything is unclear.

Hanno

[mbeutel]

Hi Hanno,

Thanks for your quick reply and for your guidance. I admit that it was mainly convenience which made me give test particles a mass. As I prefer to avoid undefined behavior, it'll be best if I keep my test particles massless and store the physical mass separately.

Cheers,
Moritz

[hannorein]

Great. Then I would recommend to keep the masses at zero and attach any additional properties using either REBOUNDx or the ap pointer of each particle. Feel free to ask if it's unclear how to do that.