dtamayo / reboundx

A library for adding additional forces to the REBOUND N-body integration package
GNU General Public License v3.0
80 stars 60 forks source link

gr_source flag not working #78

Closed tamidodo closed 2 years ago

tamidodo commented 2 years ago

Hi,

When adding General Relativity effects using "gr", it says in the documentation that the "gr_source" flag can be applied to a specific particle in the simulation to indicate that the effects should be applied using that particle as the source. It also seems that this effect should be able to be applied using multiple particles as sources. My use case is that I have multiple particles where GR is significant but too many particles to run "gr_full" without significant slowing. Looking in the gr.c file, it seems that the source is hardcoded to be the particle with index 0. Was this a functionality that got removed for some reason?

hannorein commented 2 years ago

Good question. I can see that this was removed in e3795b930fb33c0741c1f5dd5d68fa849e3c8f89. That was over 5 years ago and I don't remember why. The only thing that comes to mind is speed. Maybe Dan's memory is better. We need to either correct the documentation or reimplement gr_source.

dtamayo commented 2 years ago

Thank you! I believe I need to update the documentation. I think the reason we took it out for 'gr' is that Saha and Tremaine, which it is based on, assume that the planet-star mass ratio is << 1 and the central body is at index 0. I need to double check whether my memory is correct tomorrow.

What is your particular setup?

On Mon, Jun 20, 2022, 11:24 PM Hanno Rein @.***> wrote:

Good question. I can see that this was removed in e3795b9 https://github.com/dtamayo/reboundx/commit/e3795b930fb33c0741c1f5dd5d68fa849e3c8f89. That was over 5 years ago and I don't remember why. The only thing that comes to mind is speed. Maybe Dan's memory is better. We need to either correct the documentation or reimplement gr_source.

— Reply to this email directly, view it on GitHub https://github.com/dtamayo/reboundx/issues/78#issuecomment-1160897841, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA2ABF6O35PGCQZ5XRQTHLTVQDVRLANCNFSM5ZKD3FBQ . You are receiving this because you are subscribed to this thread.Message ID: @.***>

tamidodo commented 2 years ago

Looking at Apophis close approach structure over 100 year integration - there's a significant shift in target plane location for the various close approaches if I use gr_full vs gr but gr_full takes an unreasonable amount of time to run for the sims I want to do. I'm fairly certain what makes the big difference is the GR from Earth, Moon and Jupiter so I was hoping to just use GR for the Sun and those extra particles instead of for all the planets and 25 other asteroids I've included in the sims.

dtamayo commented 2 years ago

Hi Tammy,

I'm afraid that looking at the source code, I don't see an easy way to fix this. gr.c relies on Jacobi coordinates, which simplify the equations of motion (if you're interested see Appendix B of the REBOUNDx paper). These are a combination of all particles' positions and velocities, and use routines to go back and forth with inertial coordinates, which would (at least currently) require us to make reordered copies of the particles array for each gr source we wanted to consider.

I'm somewhat surprised that this matters though. Isn't the size of this acceleration GM_earth/rc^2 smaller than the acceleration Apophis feels from the Earth? For a close encounter distance of the Earth-Moon distance, I think that would be ~10^11 times smaller. If perturbations at that level matter, I would naively think that opens a whole can of worms. Doesn't that mean that one should also include not only the J2 of the Earth, but also J4, J6 (maybe J3)? Perhaps also these terms for the Moon? Could it be that Apophis has several close encounters with Earth over 100 years and its orbit is chaotic, so that any perturbation at all causes a significant shift in the target plane location? One could test that by changing the initial x position by Apophis at the 15th decimal and seeing what happens.

Apologies if the answer is already obvious to you, but if not, it might be worth thinking through whether this is definitely the limiting ingredient, because I think making this work in the code would take a fair bit of effort, unless @hannorein sees a simpler way. I've updated the documentation to not mention gr_source for now.

hannorein commented 2 years ago

This is a chaotic system, so small difference can grow exponentially fast. But I tend to agree that GR from other bodies should really be a very small effect. Even if you see something change when you turn on/off gr_full, that doesn't mean that GR is important. I think, like Dan suggested, that you need to integrate an ensemble of simulations initialized with small differences that are consistent with the observational constraints. Then you can look for differences in statistical quantities (probability of Apophis colliding with earth, etc). Sorry if this is something you already thought about and we're missing something!

tamidodo commented 2 years ago

I think we're on the same page in terms of other factors making a difference in the target plane location - I'd considered the simulation setup that included J2/J4 considerations for the Earth, radiation forces, GR (with the Sun as the only source), and 25 other asteroids included in the sim as the "nominal simulation". I then reran it removing one thing at a time to see how much of an effect on the target plane location for the 2036 encounter there was. I was surprised to find that the difference in target plane location between using gr and gr_full was only an order of magnitude less than the difference between using gr and not including it at all. These differences were also the biggest of all the scenarios I tried by a few orders of magnitude (removing some of the asteroids, removing radiation forces, etc.).

That might not mean that GR for other bodies is important, I just took it to mean that it was more important than some of the other things I'd made sure to include in the sims like J2/J4 for the Earth. If it had made a difference but not as large as other effects had on the sims, I would have been confident that my "nominal" setup was fairly robust for what knowledge we have, just using the Sun as a GR source. My original thinking was that GR for other bodies wouldn't make a huge difference so even though it would be more precise to include it for all bodies, it wasn't that much of a compromise to just consider it for the Sun.

To explore the target plane for the close approach structure, I ran 20 000 sims, each with varied initial position and velocity to make the Apophis particle go through different target plane locations for the 2036 encounter. The statistical quantity I'm looking at is the minimum distance between Apophis and Earth over the 100 year integration, mapped to a given target plane location. I have not run a full simulation exploring the target plane using gr_full so I'm not confident that the structure remains the same and it's just shifted by the difference.

dtamayo commented 2 years ago

Hope things were sufficiently resolved. I know that Matt Holman has been working a lot to build an ephemeris quality integrator from REBOUND, and has also been running tests on Apophis to compare with JPL. He's talking about it at the DPS meeting this year--I'm sure he'd be interested in your work! I'll close this for now--feel free to reopen at any point.