Gravity debug2

[dhubber]

This pull request contains various bug fixes, code improvements and other things for the Meshless FV gravity algorithm. These include :

• Various minor bug fixes in the general meshless gravity routine
• Correcting the variable smoothing length gravity (i.e. zeta) terms for the meshless
• Adding the analytical solution to the Jeans test
• Simplifying various OpenMP loops (mainly changing private variable to local thread variables)

The original gravity_debug branch got entangled with multiple commits due to rebase conflicts. So I created a new gravity_debug2 branch and cherry-picked the relevant commits. All should be there, but I'll double check soon before deleting it completely (from my local machine).

There are still a few issues with this branch which should be sorted before the pull request is accepted.

• Restarting simulations from snapshots does not appear to always work with the Meshless FV with gravity (in the Boss-Bodenheimer test) with various assertion errors being generated. This suggests either invalid values are being written to snapshots, or particle values are not being correctly converted to the various fluid vectors (e.g. Wprim, Qcons) when reading from snapshots. I tried also with hydro simulations and it did not appear to be an issue, so maybe this just involves gravitating simulations?
• The Meshless FV scheme will still crash with the masses of some particles dropping to zero. This may still be a simple fundamental flaw of the meshless scheme rather than an implementation bug. If that is the case, what should we do about it? We could throw and exception if self-gravity is enabled without the MFM scheme perhaps?

Once these two issues have been resolved/decided, then hopefully this will be ready for merging.

[rbooth200]

The compilation error for commit 5388781 is due to a MPI call in an OpenMP region (in block timesteps for the meshless, line 967). Do we do this anywhere else? Usually its a bad idea as it requires a higher level of OpenMP/MPI support. Since were are doing a synchronization here we should probably just exit/re-enter the parallel region.

[rbooth200]

I've added a fix in Ewald for the periodic gravity bug. However, I've noticed two things:

1) Ewald gravity requires compilation with GSL. I'm not convinced we check anywhere that the code has been compiled with GSL when using periodic boundaries.

Here's an off the wall suggestion: We could just ship with the necessary functions as part of gandalf - It's open source GPL code, we are GPL as long as we attribute it properly it should be fine...

2) The Ewald gravity is not included when stars are used. I've not changed this as it would also be needed to be added in the Nbody routines for consistency. I suspect that we should support this, however.

[rbooth200]

I said I'd look into what GIZMO is doing for its calculation of rdmdt. As far as I understand it, here is what I have found.

Like use, GIZMO is using a KDK scheme for the timestep update. However, it is doing the flux calculation at the end of the step, after the drift has been aplied and the gravitational accelerations have been calculated. This flux is then used again for the K at the start of the next step. This way rdmdt is up to date.

The GIZMO update is equivalent to what we are currently doing in SPH, but not with the Meshless, in which we compute the fluxes at the start of the KDK step. Essentially a single kick for dt/2 is done at the start of the simulation. From then on another dt is done for the double K (KDKKDKKDK) and a single kick for dt/2 is done at the end to finish the simulation.

We could experiment with a GIZMO style update, it would require an extra flux calculation for dt/2 in the post initial conditions setup, along with moving the density, gradients and flux calculation to the end of the time-step. There is also some subltely in how dQ should be applied doing it this way, especially with block time-steps.

Let me know your thoughts.

[dhubber]

Regarding GSL, yes I think this should be a simple exception thrown in the Ewald constructor if GSL is not compiled with the code. It currently runs through fine but actually not doing it correctly at all which is not as it should be. So, I'll add this on my next commit.

I don't think GSL should be shipped also with Gandalf. It is one of these standard libraries that almost all linux distros will have in their package managers (Mac also) so the way it works now is if you don't provide a path, it simply looks for the headers in the standard places (e.g. /usr/include) and if they are not there, it fails on compilation. The only addition I would make is perhaps there is some way we can check in the Makefile perhaps if those files are included, and if not, print an error that they cannot be found and cease compilation. All of this also applies really to the FFTW library.

[dhubber]

Changing the Meshless loop to mimic a LFKDK-style update does sound appealing. Maybe this would also make it easier to use the Runge-Kutta style integration (which is not currently working for gravity) for gravity also, if we decided to fully implement that approach (although we should concentrate just on MUSCL for now).

An extra flux-calculation in the set-up is nothing to worry about. Re-ordering the loop may cause issues. We'll have to think carefully, as you mentioned, about how the particles, W and Q variables are updated to make sure we don't screw something up, but should be do-able (obviously if it's already done in such a way in GIZMO). I'll have to look at the code a little more myself to think about how it should change before doing anything.

[dhubber]

Regarding the OpenMP-MPI compilation error, it is also done in SphSimulation.cpp at line 1311 in almost the exactly the same way. However, I can now see that the SphSimulation uses 'default(shared)' whereas the Meshless uses 'default(none)'. We've tried to use 'default(none)' everywhere but I guess this was changed (Giovanni?) to get the OpenMP and MPI to work together in harmony here. Nevertheless, I will change the Meshless to 'default(shared)' and see if it fixes the problem.

[rbooth200]

We'll have to think carefully, as you mentioned, about how the particles, W and Q variables are updated to make sure we don't screw something up, but should be do-able (obviously if it's already done in such a way in GIZMO). I'll have to look at the code a little more myself to think about how it should change before doing anything.

I will have a look into this next week. First with pure hydro to see whether it makes a difference on the gresho test -- I suspect it will. I'm going to start with a global time-step only implementation, which I already know how to do correctly. When I get that working I'll have a think about block timesteps. I've worked out how to do the conservation correctly with block timesteps, the main issue will be how to do the prediction of Qcons/Wprim properly for inactive particles. We should probably skype to chat about this before I get too far.

[giovanni-rosotti]

Sorry, just catching up now with the discussion.

1) About GSL, I think in theory we could ship it with GANDALF. The question is: how careful does one need to be in the compilation process? For example I wouldn't want to ship FFT with gandalf, as its optimization depends a lot on the machine and supercomputers typically have their own version (which is far more optimized than the one you would get simply by downloading the library and compiling it). In any case definitely the code that needs GSL should not run if compiled without GSL!!!!

2) I've double checked the level of threading we require in MPI. We require the minimum MPI_THREAD_FUNNELED (lower than this there is only MPI_THREAD_SINGLE, which means that you can't have multi-threaded code at all), where we promise that we will call MPI functions only from the master thread. That's what we are doing, so it's actually fine. Using default(shared) for this specific case is thus correct (although unpleasant).

3) can we move the discussion to reorganising the main loop to a separate thread? In any case I need to think more about this...

[rbooth200]

I don't think GSL does any special optimization. Personally I don't mind having it as a requirement for periodic gravity since it's so common and periodic gravity is not (cosmology aside). FFT not so much, many people will not make use of it - I don't think I've ever compiled with it enabled...

I'd make sure that the Ewald code only throws if we acutally try to use GSL functions. We often construct the Ewald object and never use it, so I'm not convinced we should throw on construction.

[dhubber]

I don't think there's any need at all to ship GSL with GANDALF. We've already provided a simple approach in the Makefile where the user only has to switch on the GSL (or FFTW for that library) flag to 1 and it will look in the usual/default place (which should cover all the Linux users and most of the Mac users) and if not, then they just add the library paths right there. I believe that the GSL functions are only used during set-up and therefore any kind of optimisation we could achieve through compiling statically with it would be negligible. But I should double-check this though!

Regarding throwing an exception, although the Ewald object is often constructed, the tables are only calculated (and therefore used) when at least one of the boundaries is periodic. So I was suggesting throwing an exception in this case, i.e. one of the boundaries is periodic AND GSL is not included.

[dhubber]

I pushed a one-line change with the 'default(shared)' edit but Travis didn't check it for some reason (there was another branch getting checked, but assumed it would be intelligent enough to queue them, but I guess not). Any way to force Travis to check the branch without actually needing to push another commit??

[giovanni-rosotti]

I closed the pull request and reopened it, and this was enough to force Travis (currently rebuilding).

Fine with throwing when the GSL stuff is actually used - please go ahead with it.

[rbooth200]

I think there are two things remaining to do on this branch:

• Add safety checks to the Ewald to make sure GSL is there when its needed (but not complain if its not there but not needed).
• Add Ewald summations for the stars.

I've already corrected the issue with the periodicity of the particles / cells. See commit 2c39d97

[giovanni-rosotti]

Unfortunately there are conflicts to be solved before we can merge...