taby said:
Is there an easy way to do double-sized floating point numbers in OpenGL 4.3 compute shaders?

In short: no.

All GPUs can do doubles natively, but it's 64 times slower than float.
Still, i don't think you could beat this.

An arch which is fast for doubles is AMDs CDNA for example. Maybe there are some ‘affordable’ GPUs of that for workstations.

Bummer. Well, I still have CPU cores to think of. Thanks for your input once again.

taby said:
Well, I still have CPU cores to think of.

Yeah, that's convenient to use and maintainable.

Ofc. you can also port the way high precision libraries work to GPU as well.

But what do you do? Motion of a single planet is not parallelizable at all?

Path integral… You calculate many paths and take the integrated version in the end. Integration of the integrated. This was your idea???

taby said:
Path integral… You calculate many paths and take the integrated version in the end. Integration of the integrated. This was your idea???

My idea? Can't be. Idk what a path integral is.

My idea would be to draw an ellipse locally, and rotate that globally while drawing. Problem solved. No Nobel price, though. ; )

Ok, so where the initial conditions are:

const MyBig initial_vel = 25000; // 38858.47;

custom_math::vector_3 mercury_pos(0, 69817079000.0, 0);
custom_math::vector_3 mercury_vel(-initial_vel, 0, 0);

I get the following bit to angle conversion table:

bits   angle
------------
20     17.04
21     17.04
22     17.04
23     47.8
24     38.95
25     34.8
26     34.96
27     34.9
28     35.07
29     35.08
30     35.06

For this setup, the analytical solution for the relativistic precession is 103.7 arcseconds per Earth century.

If you add bits 22, 23, and 24 together (e.g. 17.04 + 47.8 + 38.95) you get 103.79 arcseconds per Earth century.

It works for many different initial speeds.

Edit: There is a weighting function that weighs the bit contributions. I have a vague vision in my mind, but I'm not certain yet.

Must be more luck…….

So, Joe, there is double precision in GLSL, but not OpenGL? In that case, can I pass in two float textures and then decode them in the shader as one double texture?

taby said:
For this setup, the analytical solution for the relativistic precession is 103.7 arcseconds per Earth century.

Why don't you use the analytical method if it exists? It would allow huge timesteps and high accuracy. It would prevent piling up hacks out of thin air until the result becomes close to correct?

taby said:
If you add bits 22, 23, and 24 together (e.g. 17.04 + 47.8 + 38.95) you get 103.79 arcseconds per Earth century. It works for many different initial speeds.

Proposal:

Use the 30 bits value, which does your math most accurately, and multiply the result by 103.7 / 35.06.
That's just as good, just as wrong, much faster, and works for many different speeds as well, i guess.

Conclusion:

Either the math you use probably isn't totally wrong, but also not yet totally right. Something is missing or wrong. Otherwise there would be no need to introduce hacks.

Or the error you get is actually just integration error. If so it would get smaller with smaller timesteps.

Use the math of the analytical solution to figure out the mistake / make a specialized integrator to minimize error.

taby said:
So, Joe, there is double precision in GLSL, but not OpenGL? In that case, can I pass in two float textures and then decode them in the shader as one double texture?

If there are no double texture formats, use shader storage buffers instead.

Textures are specialized image formats which the texture filtering HW can process.

SSBO is for general data which can be anything. It's usually the first choice for general purpose programs not working with images.

The analytical method that I'm talking about is either (they're super close):

const MyBig delta = 6.0f * pi * grav_constant * sun_mass / (speed_of_light * speed_of_light * (1.0f - eccentricity * eccentricity) * semi_major_axis);

const MyBig delta2 = 24.0f * pi * pi * pi * semi_major_axis * semi_major_axis / (orbital_period*orbital_period * speed_of_light * speed_of_light * (1 - eccentricity*eccentricity));

I used 0.01 seconds for the timestep for Euler, and 0.1 for symplectic 4th-order. Both give practically the same result.

SSBO it is! Thanks again for your insight!

P.S. I am only quantizing the time dilation because I can. It’s interesting how the numbers pop out. Why do these numbers exist if not to be found out?