search

gandalfcode / gandalf

GANDALF (Graphical Astrophysics code for N-body Dynamics And Lagrangian Fluids)
GNU General Public License v2.0
44 stars 12 forks source link

Sph LeapfrogDKD #86

Closed rbooth200 closed 7 years ago

rbooth200 commented 7 years ago

I think there may still be bugs in the SPH leapfrog DKD algorithm. The problem as I see it is related to the fact that the force calculations need to be done in an unsyncrhonized manner:

Steps  : x - - - x - - - x - - - x - - - x - - - x - - - x - - - x - - - x
Level 0: D - - - - - - - - - - - - - - - K - - - - - - - - - - - - - - - D
Level 1: D - - - - - - - K - - - - - - - D - - - - - - - K - - - - - - - D
Level 2: D - - - K - - - D - - - K - - - D - - - K - - - D - - - K - - - D
Level 3: D - K - D - K - D - K - D - K - D - K - D - K - D - K - D - K - D

The x's denote the timesteps on the lowest level (when n is incremented by 1), while D's and K's are drifts and kicks. This shows that two force calculations are needed per timestep since we have kicks (K) appearing, i.e. at both the start and middle of the step. This doesn't fit into our SphSimulation MainLoop. This could work, e.g. if we are integrating an extra dummy level per timestep when we use DKD integration, so that each K appears syncrhonised with a time-step. However, this doesn't seem to be what we are doing.

The main question I have is why do we have both KDK and DKD? They are equivalent when the time-steps are fixed, otherwise KDK is both more accurate (e.g. Springel 2005) and more efficient (since the force computations appear synchronised so are fewer tree builds, commuincations and loops over particles).

I'd propose that we do away with the DKD all together unless there is a good reason to keep it? The main advantage being that we are getting rid of a scheme that no one will ever use which means a smaller code base to maintain.

dhubber commented 7 years ago

Actually, when running with LFDKD, we are running two half-steps per full integration step. There's a variables called integration_step which is set to 1 for one-step integration schemes (e.g. KDK) and 2 for 2-step schemes (e.g. DKD). Then in the block timesteps code, the variable level_step has an extra level added for 2-step schemes allowing the DKD to be integrated correctly. (Actually, the way I've it is not strictly correct if we wanted to upgrade to a 4-step scheme for instance, but for 1 or 2 steps, it's fine). So, the extra 'dummy level' per step is in effect there!

You are right that they are both equivalent for global, constant timesteps and also that the KDK is the most efficient CPU-wise. Also the DKD doesn't work so well with the Saitoh & Makino timestep limiter. But the question of which is better was certainly not resolved by Springel (2005). That was just a simple binary orbit and might be the case for 2-body dynamics, but when we were developing SEREN we tried both and found that their performance and accuracy characteristics varied a lot. It seemed for example that DKD could perhaps be more accurate for angular momentum conservation with discs, particularly with block timesteps. However, I've also seen cases where the KDK seems better. If it were a clear-cut case of one being better, I'd agree with getting rid of one (in fact, there wouldn't have been two in the first place). But I don't think that's the case. It might require a few well-posed test cases if we really wish to delete the DKD.

rbooth200 commented 7 years ago

Ok, which case I agree it's not a bug, but it could certainly use some more documentation in the code as I wasn't able to work out what was going on! If its working I also see no reason to get rid of it. I'll close the issue...