fix: SegFault when a MPI rank has no data of a sub-region

[MelReyCG]

Closing #3528

My goal here is to solve the crash, and to propose a way to reduce (Min/Max) pairs lexicographically over MPI ranks. That can typically be useful to get the min/max value (pressure, temperature...) along with its globalIndex in the mesh (while ensuring that the globalIndex will be stable).
It could be generalized to tuples, let me know if that could be useful.

[corbett5]

I'm sure there's a reason I haven't thought about, but why not use std::pair? And those mpiPairType::get functions can be more succinctly written as templated variables.

Example: https://godbolt.org/z/3sxTeGrzs

[MelReyCG]

Hi @corbett5!

1. Yes, I indeed started with std::pair, but I did not continue with it because is not trivially copiable. As it seem for the least that it could be used in performance critical places, it seems like a reasonable requirement.
   I wonder if this is why the original code here did not go for std::pair.
   https://godbolt.org/z/KE7TEMx7d
   https://stackoverflow.com/questions/64978278/what-is-the-performance-problem-with-stdpair-assignment-operator

2. IIUC, you suggest the following?

/* no default get() implementation, please add a template specialization and add it in the "testMpiWrapper" unit test. */
template< typename FIRST, typename SECOND >
MPI_Datatype const mpiPairType;

template<> MPI_Datatype const mpiPairType< float, int > = MPI_FLOAT_INT;
template<> MPI_Datatype const mpiPairType< double, int > = MPI_DOUBLE_INT;
template<> MPI_Datatype const mpiPairType< int, int > = MPI_2INT;
template<> MPI_Datatype const mpiPairType< long int, int > = MPI_LONG_INT;
template<> MPI_Datatype const mpiPairType< long int, long int > = getMpiCustomPairType< long int, long int >();
template<> MPI_Datatype const mpiPairType< long long int, long long int > = getMpiCustomPairType< long long int, long long int >();
template<> MPI_Datatype const mpiPairType< double, long int > = getMpiCustomPairType< double, long int >();
template<> MPI_Datatype const mpiPairType< double, long long int > = getMpiCustomPairType< double, long long int >();
template<> MPI_Datatype const mpiPairType< double, double > = getMpiCustomPairType< double, double >();

[corbett5]

@MelReyCG

1. Gotcha, yeah I think using std::pair would just work, but I think using the custom struct is not a bad idea.
2. I forgot about the getMpiCustomPairType function, on second thought would it be cleaner if we eliminated the mpiPairType variable completely and just put all the logic in getMpiCustomPairType? You can't template specialize functions, but it could be done with if constexpr checks of the types I believe.

[MelReyCG]

@corbett5 I added a new proposal in the "proposal #1" commit.
For now, I did not go for the if constexpr solution, but it would easily be feasible if you guys think it would be cleaner than this proposal.

[CusiniM]

@MelReyCG

1. Gotcha, yeah I think using `std::pair` would just work, but I think using the custom struct is not a bad idea.

I remember an old discussion about the std::pair layout not being guaranteed to be contiguous in memory which I think is the reason why it was not used in the first place. I don't recall the details of how it was being used though.

2. I forgot about the `getMpiCustomPairType` function, on second thought would it be cleaner if we eliminated the `mpiPairType` variable completely and just put all the logic in `getMpiCustomPairType`? You can't template specialize functions, but it could be done with `if constexpr` checks of the types I believe.

Personally I don't have a strong preference between using if constexpr in a single function or using template specializations as @MelReyCG has done. They way it's written now seems pretty tidy and easy to read but something equally good could probably be obtained using if constexpr. Maybe this second option is a bit more modern?

[corbett5]

Personally I don't have a strong preference between using if constexpr in a single function or using template specializations as @MelReyCG has done. They way it's written now seems pretty tidy and easy to read but something equally good could probably be obtained using if constexpr. Maybe this second option is a bit more modern?

I think the single function is cleaner, but this is in my opinion a small improvement to a minor part of this PR and I think the current implementation is fine.

[MelReyCG]

@CusiniM

I remember an old discussion about the std::pair layout not being guaranteed to be contiguous in memory which I think is the reason why it was not used in the first place. I don't recall the details of how it was being used though.

This is exactly why I chose the struct, it seems that there is no guarantee that the compilator will be allowed to optimize std::pair<> copies as much as structures. Another source: https://www.reddit.com/r/cpp/comments/ar4ghs/stdpair_disappointing_performance/?rdt=48069

About the if constexpr, you can see that I used it in getMpiPairReductionOp(), I tend to use it more when it concerns implementation details, or when it can simplify overcomplicated std::enable_if.
For such "switch on type T" usage, I prefer the good old template. I think that's a matter of style / tastes. If the code readability is not affected, I think we can keep template specializations.

[MelReyCG]

I finished to write the last comments I think, we may be good to go?

[CusiniM]

LGTM. Thanks @MelReyCG

[corbett5]

Why does this need to be static? mpiOp is returned by value, so it shouldn't matter if it gets destroyed after the function returns.

[MelReyCG]

The static storage here does matter as MPI_Op (or MPI_Datatype for custom types) represents a persistent MPI resource that must:

• be initialized only once to avoid mem-leaks (MPI_X_create() functions allocates internal resources),
• remain valid for the entire MPI lifetime (MPI_Init->MPI_Finalize),
• be shared across all calls to this function.

The static construction I used here is a variation of the "Meyer's Singleton", which lazily constructs a unique instance of an object when the resource is first requested (by calling the lambda in this case).

At first I did not mind calling MPI_Op_free() since we already call MPI_Finalize() but I found the following in the MPI documentation:

The call to MPI_FINALIZE does not free objects created by MPI calls; these objects are freed using MPI_ XXX_FREE, MPI_COMM_DISCONNECT, or MPI_FILE_CLOSE calls.

To completely manage these resources lifetime I will work a bit more on this PR.

[corbett5]

Oh god, yeah I get it. If performance isn't a concern (and I doubt it is) I would simply createa new op every time and free it every time as well. You could do this pretty easily with a RAII object that had a move constructor and deleted copy constructor.

[MelReyCG]

@corbett5 Indeed, a RAII object would do the job.
Yesterday I finished with another approach in which I kept using the MPI operation / data-type caching and stored their trace to release them on MPI_Finalize(). Let me know if my approach is right (commit).

[corbett5]

I prefer my approach (in which these objects are create and destroyed on the fly) since it would eliminate the need for static variables, but this is still quite the improvement and I'm not sure it's worth your time to rewrite it.