Skip to content
This repository has been archived by the owner on Oct 14, 2023. It is now read-only.

Refactor perturbations relying on external bodies or satellites to allow acceleration #1485

Open
astrojuanlu opened this issue Feb 19, 2022 · 5 comments

Comments

@astrojuanlu
Copy link
Member

At the moment, we have 2 perturbations that make use of the position of external objects: third_body and radiation_pressure. Both perturbations take a callable f(t) that returns the position.

These are the only two perturbation functions that are not accelerated, and taking a callable is the reason for it:

def third_body(t0, state, k, k_third, perturbation_body):
r"""Calculate third body acceleration (km/s2).

def radiation_pressure(t0, state, k, R, C_R, A_over_m, Wdivc_s, star):
r"""Calculates radiation pressure acceleration (km/s2)

In #1434 it became apparent that being able to quickly interpolate the position of external objects is crucial to achieving good performance. At the moment, the "slow" tests are slow mainly because of these functions. While discussing adding an event that also would depend on the position of an external satellite (#1350) I started thinking again about how to accelerate them.

I see basically two ways:

  • Offer some tooling so that users can create fast interpolants using @jit, and pass those as callables. The disadvantage is that passing functions as arguments is still slow in numba Function accepting njitted functions as arguments is slow numba/numba#2952 - however, if the functions involved are slower than the overhead that numba introduces, this might still be a viable option.
  • Modify these perturbations so that they take an array of (times, positions). The disadvantage is that the user would need to precompute and sample those positions in a predefined time interval, and propagation outside of it would fail.

Both options would require numba-compatible interpolators. I found a couple of such projects: https://github.com/dbstein/fast_interp and https://pypi.org/project/interpolation/. The former is not on PyPI though dbstein/fast_interp#3 and seems more focused on grids, so perhaps our only acceptable option is the latter. The cool thing is that it supports extrapolation, although I'm not sure how much we should rely on it.

Adding numba-compatible interpolators would essentially mean replacing all of https://github.com/poliastro/poliastro/blob/38a28e4/src/poliastro/_math/interpolate.py by something else. Or perhaps we could add @jit to sinc_interp?

Related wiki pages: https://github.com/poliastro/poliastro/wiki/Orbit-interpolation, https://github.com/poliastro/poliastro/wiki/Ideas-for-fast-interpolation

@astrojuanlu
Copy link
Member Author

Well, pypi/interpolation is not very well maintained EconForge/interpolation.py#92. Perhaps our only realistic option is to try adding @jit to sinc_interp to start, and if that isn't enough (because we also need linear or even spline interpolation), we'll have to write our own...

@astrojuanlu
Copy link
Member Author

By the way, even if Numba @jitclass is not perfect numba/numba#7372, it is already powerful enough for some purposes, see for example https://github.com/mikegrudic/pytreegrav/blob/352ca273/src/pytreegrav/octree.py.

@astrojuanlu
Copy link
Member Author

Related: #1029

@astrojuanlu
Copy link
Member Author

Raising the priority of this.

@astrojuanlu astrojuanlu added this to the 0.17 milestone Jun 3, 2022
@astrojuanlu
Copy link
Member Author

Ran out of time to do this before 0.17, but I'll definitely tackle this in 0.18.

@astrojuanlu astrojuanlu modified the milestones: 0.17, 0.18 Jul 7, 2022
Sign up for free to subscribe to this conversation on GitHub. Already have an account? Sign in.
Projects
None yet
Development

No branches or pull requests

1 participant