Levelneib not updated with MPI

[giovanni-rosotti]

See for example here. The loop is over Nhydro, so only local particles have levelneib updated.

[dhubber]

So, I guess the solution would be :

• Also loop over the exported particles in these omp critical sections.
• Send levelneib back with all other exported info (e.g. accelerations) and do a max-reduce on levelneib in the import part of the code. Is it that simple or am I missing something?

[giovanni-rosotti]

I think the second part is already implemented, so we would need only the first one. Well the reason why I created an issue is to have a second opinion (and to carry on with what I was doing without forgetting to fix this), so if you also think it is that easy I'll just do that. But we should also worry about periodic ghosts, shouldn't we? Although we are not computing the force on them, if their levelneib has been modified this change should be propagated back to the original real particle, shouldn't it?

Other two things: -levelneib is not updated in ComputeSphGravForces. I don't see why, should I add it? -I think that critical section was the major bottleneck in the openMP scaling when I profiled the code, I'll think how to improve it. To start with: do we really need to loop over all particles? I guess so, levelneib might be modified also for inactive particles

[dhubber]

Okay, then hopefully it will be that simple then. You might be right about the periodic ghosts also. I was thinking maybe the original particles would do this but since the transmission of levelneib is effectively a scatter operation from live real particles, this will never happen for ghosts on real particles so I guess we also have to propagate this information from the ghosts to the real particles. Hopefully this will also be straightforward since we have iorig at hand. However, I'm not up-to-date will all ghost things after Richard's changes so maybe there is something more to it.

About the other things

• Yes, it should be added. The ComputeSphGravForces routine is only really used in tests (e.g. freefall collapse) but should be updated to be consistent with the AllForces routine.
• Okay, if this is true then it's definitely reason to optimise this. I guess I was hoping one simple critical section at the end (rather than many small ones through the loop) would not be a problem but guess that's not the case!

[rbooth200]

I think that periodic ghosts should be fine. The explanation follows (and hinges on the assumption that mirror ghosts are just MPI ghosts when MPI is used):

• Without MPI: The periodic ghosts are generated on the fly in the force calculation so we include them in the neighbour loop and thus include them in the local levelneib calculation for the active particles. (We don't need to worry about including the contribution for the periodic mirror back to the original particle since that contribution is done in the calculation for the periodic mirror itself in which the first particle is generated as a mirror).
• With MPI: The periodic mirrors are now just MPI ghosts, so any solution that works for general MPI will work for periodic mirrors when MPI is used,

[rbooth200]

-I think that critical section was the major bottleneck in the openMP scaling when I profiled the code, I'll think how to improve it. To start with: do we really need to loop over all particles? I guess so, levelneib might be modified also for inactive particles

There are two possible solutions:

1. Make levelneib an OMP shared variable, but I think we might run into problems with false sharing.
2. Stick with the per thread levelneib implementation, but allow all threads to access levelneib at the end. This way we can do the parallelization over the particles, rather than levelneib. I.e. turn the loop over particles into a parallel for and for each thread do a maximum over the level neib created for each thread.

Solution 2 will need a barrier and a flush at the end of the loop, so it's not completly obvious that it will be faster than 1, but I can't think of a 3rd solution.

[giovanni-rosotti]

This is now solved in #78. I think there is no need for a barrier since there is already an implicit barrier at the end of a "for" region