Alternative approach for spherical harmonic gravity / geopotential

[sblana]

spherical harmonic gravity can be implemented with an O(1) time complexity (with respect to the maximum degree/order) if a 3-dimensional map of gravity (multiplied by distance squared) is pre-calculated, stored in memory, and then during orbit propagation, the gravity*r² of the texel that corresponds closest to the spacecraft's position is fetched and divided by distance squared to get the gravity. the map could be computed via a compute shader, computed in several passes (in case of exceeding time limits), stored as a file, etc. i tested it with a simple compute shader, which takes like <5 seconds to compute a map with a maximum degree of 50

my testing shows improvement in performance compared to running the spherical harmonics equation every time you want gravity. it does result in way higher memory usage, but a resolution of e.g. 256x128y32z (longitudinal, latitudinal, radial), which still produces good results, is only 32 MiB.

MS is the approach i propose here. PSE is the same approach used in principia (though without the dampening i.e. with a fixed maximum degree, it being L). i did this in godot so the numbers might be inflated compared to this mod's implementation

idk how important the SH gravity is in terms of priority of development or if there even is a point in better performance, so this may be completely pointless even if it is several 100x faster. but i do wanna open discussion on this as a viable alternative approach.

[himisawww]

My test of geopotential summation in c++ (my code, but with similar implementation as pincipia (maybe, I can't really understand the source of Principia)) takes 0.084902/0.207465/0.623894us per gravity evaluation for degree && order 5/10/20 harmonics. (CPU: 7950x3d)

Note by default in my code, the maximum degree & order is 8, I changed static const int_t Max_N= to 20 to perform the test.

I think c++ is fast enough for reasonable number of degrees. Actually for integrating ephemeris of major bodies in Solar system, no more than 8th degree of harmonics is required. (However for space craft trajectory, as they can be arbitraily close to surface of a planet, higher degrees do make a more significant difference. Nevertheless, I don't think adding too much (more than 20) degrees in a game is necessary.)

[eggrobin]

Thanks for the suggestions.

As you wonder in your last paragraph, the performance of our spherical harmonics calculation is not a problem ; our benchmarks show that we can compute a field of degree and order 10 in less than 600 ns, so we are not interested in pursuing alternative approaches. Here are nevertheless some comments that you may or may not find useful for your own projects.

The approach you propose entails constantly fetching data from main memory (given the size of the precomputed data, especially with 30 bodies, this will not live in the cache), which you don’t get to do much of in that time. The numbers for the implementation labeled PSE in your graph are not just inflated: they are inflated by a factor of 1000, which is likely to reflect some structural issue with that reimplementation of the method used by Principia. There is no reason to think that the method labeled MS would conversely become a thousand times faster if reimplemented in Principia, so it is not clear at all that it would actually be faster.

The approach of getting the force from « the texel that corresponds closest to the spacecraft's position » is a degree-0 approximation, which entails discontinuities in the force. This is likely to lead to unphysical effects, and to break the conditions that are required for our numerical integration methods to have nice properties such as energy conservation. When you say that the 32 MiB « still produces good results », how have you tested that ? In Principia, have some fairly sensitive tests of the long-term behaviour of orbits about the Moon (in https://github.com/mockingbirdnest/Principia/blob/master/astronomy/lunar_orbit_test.cpp), and we chose the cutoff at degree and order 30 based on those, comparing with the results of papers that used a model of higher degree and order.

In general, open-ended discussion of numerical techniques is best had on the Discord, rather than in issues here; we try to use issues to track things we actually want to do.

[pleroy]

My test of geopotential summation in c++, but with similar implementation as pincipia (maybe, I can't really understand the source of Principia)) takes 0.084902/0.207465/0.623894us per gravity evaluation for degree && order 5/10/20 harmonics. (CPU: 7950x3d)

Interestingly, that's aligned with what we see in our benchmarks (185 ns/544 ns for degree and order 5/10, CPU 5965WX). Your implementation is faster but the conclusion remains that it's fast enough for practical purposes, and it beats fetching data from main memory.