Data mapping with MPI and duplicated vertices

Using preCICE
,
1

Hi, I’m working on developing a preCICE adapter for coupling two different parallel fluid codes in 3D. I’m using fortran with MPI for domain decomposition. I’d like to use the linear-cell-interpolation mapping method of preCICE, and as such I need to define tetrahedra between vertices.

Because of MPI + tetrahdra decomposition, I’m somewhat close to the case Use a single mesh and duplicate copied vertices in https://precice.org/couple-your-code-distributed-meshes.html

In other words, let’s say that I have a mesh defined by

x=0,\Delta x,2\Delta x,...,L_x=n_x\Delta x\\ y=\Delta y/2,3\Delta y/2,...,L_y-\Delta y/2\\ z=0,\Delta z,...,L_z

and the domain is split along x across multiple MPI tasks.

On the sender and receiver I have two different resolutions, i.e. different \Delta x and \Delta y, and the domain is periodic along y.

Because of periodicity (and more generally, if I’d like to use MPI also along y, which is the case), I thought creating two meshes, one on the sender, and one on the receiver. On the sender, I extend y to go from -\Delta y/2 to L_y+\Delta y/2 (for MPI decomposition along y, this would mean adding ghost cells at both ends). On the receiver, I keep y as given above. This way, I can put the periodic values var(x,-\Delta y/2)=var(x,L_y-\Delta y/2) and var(x,L_y+\Delta y/2)=var(x,\Delta y/2) into the send variable and, when performing the preCICE interpolation from one mesh to the other, ensure that the receive domain is fully contained into the domain of the sender (which is good for linear interpolations).

Now, I have the following problem: when I split x across MPI tasks, I will have the right boundary of rank-1 to be exactly the same as the left boundary of rank. For example, if L_x=1 and I use 2 MPI tasks, I will have x=0,dx,...,0.5 on rank 0 and x=0.5,0.5+dx,...,1.0 on rank 1. And therefore, I will have duplicated vertices among different MPI tasks. According to Use a single mesh and duplicate copied vertices, this doesn’t seem to be a major issue. However, when implementing it in a test code, it doesn’t work as expected.

I created an example, which I attach to this message (main.f90, inside the zip folder), together with the preCICE configure file I used. In short: I have one single program, which acts both as the sender and as the receiver. If it is invoked with the flag “0”, it acts as sender, if invoked with flag “1”, acts as receiver.
It then creates the 3D mesh (one for the sender and one for the receiver), as discussed above. Given that it is a structured cartesian grid with uniform \Delta x, \Delta y, and \Delta z, I use the standard Freudenthal method to decompose each hexaedron into thetrahedra.
I then create a 3D variable (varying only along z) and communicate via preCICE.

Everything runs fine when using 1 MPI task. On the other hand, when running with multiple processes, the program reports some filtering and missing thetraerda. As an example, I attach the output of the receiver for a run with

mpiexec -np 2 ../a.out 0
mpiexec -np 4 ../a.out 1

on two different terminals (the first command for the sender, the second for the receiver).

As can be seen in the out file inside files.zip, one of the receiver tasks report:

preCICE:WARNING:  1 of 2340 vertices of mesh meshSend have been filtered out since they have no influence on the mapping.

This is quite strange in my opinion, because if the problem are the duplicate vertices, I would expect more than one vertex to be filtered out. But, because of this filtering, I then see a bit below

preCICE: Linear Cell Interpolation is used, but some points from meshSend don't lie in the domain defined by the tetrahedra. These points have been projected on the domain boundary. This could come from non-matching discrete geometries or erroneous connectivity information. If distances seem too large, please check your mesh. The fallback statistics are: min:0.00349386 max:0.00349386 avg: 0.00349386 var: 0 cnt: 1

which is my main problem: the original send mesh fully contained the receive mesh, but because of the filtering of one vertex, I now have some missing tetrahedra.

This has further implications, as lower down I get the wrong result:

 #           2 : recv:   0.12000000000000002       0.20000000000000001

Here, I’d like to note that the wrong result 0.1200000, rather than 0.2000000, doesn’t appear at each run, but only sometimes (not always reproducible, even if I use preCICE compiled in debug mode, and I compile my code with the -O1 flag).

I don’t understand if I’m doing something wrong, or if this setup is simply not supported by preCICE. Any help on this issue would be greatly appreciated.

Best,
Fabio
precice-config.xml (1.3 KB)
files.zip (3.6 KB)

2

Hi everyone, me again. I’m still investigating my issue, but maybe an additional, simpler question:
Does preCICE handle domain decomposition with MPI + mesh connectivity + different spatial meshes between sender and receiver? If yes, how to handle mesh connectivity between different MPI domains? Is the only possible way to either avoid mesh connectivity at MPI boundaries, or to have the global mesh on the master MPI process (so no vertex duplication), and then scatter the result to everyone else?

Best,
Fabio

3

Does preCICE handle domain decomposition with MPI + mesh connectivity + different spatial meshes between sender and receiver?

I have not done what you are trying to do with volume coupling, but I have done it with surface coupling and in that context the answer to this question is 100% yes. So I would expect it also to work with volume coupling.

The message about vertices being filtered out has nothing to do with vertices on multiple ranks but with the same coordinates. PreCICE constructs an interpolation from data on your source mesh to your target mesh, and depending on where you are asking for data on your target mesh, it may be that some vertices on the source mesh do not contribute to the interpolation. PreCICE “filters out” those source vertices, meaning that their data is ignored when deciding what values the receiver will receive. Seeing a warning about vertices being filtered out is not inherently cause for concern. However, if your two solvers use the exact same domain (same extents, same partitioning into hexahedra, same division into tetrahedra), no interpolation should be necessary so no warnings should appear. You can test this by setting the same Δ⁢𝑥 and Δ⁢𝑦 for both solvers.

I think the second warning:

some points from meshSend don’t lie in the domain defined by the tetrahedra

is not related to filtering out a vertex. This message suggests to me that you have a send-receive going the wrong direction in your configuration somewhere, because there should be no problem with a sending volume enclosing a receiving volume and that is what you intended.

4

Thank you for your feedback.
I am aware that the message about vertices being filtered out is not, by itself, an issue. What is suspicious to me is that just one single vertex (and not a full line of vertices, that are exactly duplicated) is filtered out in my example. I’ll come back to this point later.

Before that, I’d like to report on some additional tests. You said that, if the domain is exactly the same between the sender and the receiver, I should get no warning on filtering out. While I can confirm that this is the case without MPI, if I add MPI, and therefore duplication of vertices, I get a warning on filtering out.
For example, with a mesh that is (256+2) x (8 + 2) x 3 (the +2 here means ghost cells), and using 2 MPI ranks both for the sender and the receiver, I get a message that 30 vertices are filtered out. This makes perfectly sense, I’m duplicating exactly 1 x (8+2) x 3 vertices between the two MPI ranks, and preCICE is filtering out 30 vertices. In this case, no errors about missing tetrahedra is reported.

I then repeated the test with same resolution, but where the receiver is a subset of the domain of the sender (so, again, no interpolation). Also in this case, with MPI, I get the warning about filtering out vertices (in this case, more vertices are filtered out, which is expected), but no errors about missing tetrahedra (I tested this with different number of MPI tasks, and there are no problems).

Going back to my particular example above: what is strange is that only one vertex is filtered out, and after the filtering, one tetrahedra is missing (cnt: 1). Not several tetrahedra, but just one (and it is just on one MPI rank, not for all of them). And this tetrahedra is present if I don’t use MPI (without MPI, everything works fine).

Given your feedback, I also created a second example, with a 2D mesh, where I use the set_quads routine to define the connectivity of my mesh. I can confirm that, if I just use the 2D mapping, everything works even with duplicated vertices + MPI decomposition + different resolution (depending on the MPI settings, I sometimes get a warning about filtering out vertices, which is expected, but never saw the problem of missing triangles).

I’m now puzzled. I will double-check again the tetrahedra definition. But if that is correct, I think I’m left either with the rbf mapping (which does not require mesh connectivity, but is a bit more problematic for me to implement with proper periodic boundary conditions), or with duplicating many different 2D mesh (anyhow, I’m not interpolating in the 3D dimension, because the two solvers use the same \Delta z).