Ghost corner cells

[fankiat]

Fixes #160

The code is a bit long, mainly because we are doing non-blocking communications for ghosting, and so we need to take care of the different parts of the ghost zones separately (i.e. corners). One can of course make blocking calls as shown here, which then would reduce the need to take care of different parts separately and rely on the order of communication to transfer the corners automatically.

In any case, this corner ghosting is needed when we have Eulerian-to-Lagrangian grid interpolation, since the Lagrangian grid might come close to the corner regions of the local domain, thus needed the ghosted corner cells. The implementation seems to work since (1) it passes all the update tests that includes checking corner cells and (2) Eulerian-to-Lagrangian interpolation now works (I will push the implementation in a different PR when addressing #157 )

[bhosale2]

@fankiat you requested a review but its still a draft PR, can you clarify?

[fankiat]

@fankiat you requested a review but its still a draft PR, can you clarify?

Right, this should be ready for review. Let me change the status.

[bhosale2]

The code seems to work and pass tests as seen, but I have comments regarding the large number of cases popping up and manual code repeat ensuing in that case. Please check how people do it in other MPI-based packages because maintaining or extending such code with 24 or so cases is error-prone and difficult. If a solution doesn't emerge soon please open an issue for this with medium to high priority.

Also, I'm not sure if the non-blocking is totally the cause of this issue. Can you see how they do it in PPM or some other well-established flow solver package?

[bhosale2]

After reviewing other packages if a non-blocking solution seems out of reach then maybe consider having a blocking option, which we can then consider if performance is not affected that much.

[fankiat]

The code seems to work and pass tests as seen, but I have comments regarding the large number of cases popping up and manual code repeat ensuing in that case. Please check how people do it in other MPI-based packages because maintaining or extending such code with 24 or so cases is error-prone and difficult. If a solution doesn't emerge soon please open an issue for this with medium to high priority.

Also, I'm not sure if the non-blocking is totally the cause of this issue. Can you see how they do it in PPM or some other well-established flow solver package?

Yes, I do realize that the code is cluttered and long, and it hurts my eyes as much as it hurts yours 😢 . Of course, I have also looked up how people have done it (including PPM, which IMHO is also cluttered, perhaps more confusing 😓 ), but as mentioned in the description of the PR, because we need both non-blocking MPI + corner cells (for Eulerian-to-Lagrangian interpolation kernels), we end up having to initialize sends and recvs for different parts of the ghost zones. This is not that common from what I see, given most approaches either (1) utilize kernels that rely on grid points along main axes of the structured grid and do not need the corner cells, or (2) use blocking communications to reduce the number of comm calls while achieving corner ghost cell exchange.

However, one thing I can try doing perhaps is to spend some time and carefully abstract the indices and in the data type init and ghosting calls, I have a loop to go over these preset indices. This should reduce clutter and make the code look cleaner. I am a little swamped these few weeks, perhaps I can open a separate PR for that as you suggested and set it at a lower priority for now. How does that sound?

[fankiat]

After reviewing other packages if a non-blocking solution seems out of reach then maybe consider having a blocking option, which we can then consider if performance is not affected that much.

@bhosale2 Yeah it seems that way, that a non-blocking solution without this many communication calls seems unavoidable. For example, in section 5.2 of this article where they talk about the corner cells, they mentioned (similar to what I meant by using blocking calls):
"A common way to solve this problem is to perform the border exchanges along each dimensional axis as independent
waves where each wave updates the halos in one direction"
This ordered, blocking communication pattern in 2D saves 4 exchanges (i.e. 4 sends + 4 recvs), which is also what we observe here, if we compare to the previous 2d ghosting code we had.

However, I think we should stick with asynchronous non-blocking ghost communication for the possible speedup it offers (depending on the decomposed subdomain size), as compared to no potential speedup in blocking ghost communication. We can remedy this code with the approach I suggested above. Please let me know what you think. Thanks!

[bhosale2]

However, one thing I can try doing perhaps is to spend some time and carefully abstract the indices and in the data type init and ghosting calls, I have a loop to go over these preset indices. This should reduce clutter and make the code look cleaner. I am a little swamped these few weeks, perhaps I can open a separate PR for that as you suggested and set it at a lower priority for now. How does that sound?

Can we do this? And if you want to delegate this to another issue let's keep it a mid-high priority, since I think this is a pretty important task that needs to be figured out to keep code maintenance and extension easier in the future. Also I am not clear which packages you looked but have a look at distributed immersed boundary packages if any like this one:
https://github.com/IBAMR/IBAMR

[bhosale2]

@fankiat I will approve this once the related issue tracking the comment above has been opened. I will then followup on other related PRs.

[fankiat]

However, one thing I can try doing perhaps is to spend some time and carefully abstract the indices and in the data type init and ghosting calls, I have a loop to go over these preset indices. This should reduce clutter and make the code look cleaner. I am a little swamped these few weeks, perhaps I can open a separate PR for that as you suggested and set it at a lower priority for now. How does that sound?

Can we do this? And if you want to delegate this to another issue let's keep it a mid-high priority, since I think this is a pretty important task that needs to be figured out to keep code maintenance and extension easier in the future. Also I am not clear which packages you looked but have a look at distributed immersed boundary packages if any like this one: https://github.com/IBAMR/IBAMR

I think I've seen this one before, but struggle to find where the ghosting code is implemented. I was also looking at a few others (most of which I need to dig and find again), but the overall idea I gather is we either have to sacrifice code length for non-blocking calls or sacrifice non-blocking for cleaner code. For example, here in their docstring they mentioned their approach, which is to use blocking calls to get corner cells across. The approach minimizes the number of MPI messages to send and maximizes the size of those messages, but serializes message traffic.

However, I have been thinking about this the past few days, and I think I have an idea worth trying to keep the current non-blocking implementation while making the code less lengthy and error-prone. The rough idea is to have a GhostField object, much like VectorField, and this GhostField will keep information about the ghost indices and so on. From the GhostField description, we automatically generate the required datatypes and send/recv in ghost communication. This would also be useful elsewhere when we try to refer to inner_cells and so on, and is also consistent with the way sopht use VectorField. I opened issue #165 to track the discussion here and will get working on that once other more important things are out of the way. I would like this functioning code to go through first for now (see other PRs for working illustrations), so that I can start developing cases needed for my defense (I spoke to Mattia yesterday on what needs to be done) and letting them run in the background on the cluster while I come back and work on this. Please let me know what you think.