PICO: division by zero in some parallel runs about pism

[ckhroulev]

void PicoGeometry::get_basin_neighbors(...) is broken.

• It uses the same stencil everywhere in the domain (including domain boundaries). This means that a basin at the left edge of the domain would be considered a neighbor of a basin near the right edge of the domain. The variable basins in inputs.nc is constructed to avoid this issue, but we should not require such padding.
• The code attempting to combine results from sub-domains in a parallel run does not do what it is intended to do. (In short: GlobalSum() is not appropriate here.)
• It will fail if a basin has more than 2 neighbors. (We don't know if this is a possibility. We should add a check and stop with an error message saying that the current implementation does not support basins with more than two neighbors.)

diff good-1.log good-2.log produces this:

diff -u good-1.log good-2.log
--- good-1.log	2021-11-19 15:17:41.718276242 -0900
+++ good-2.log	2021-11-19 15:17:43.330264495 -0900
@@ -510,8 +510,8 @@
   done.
 PICO, get basin neighbors of b=1: b1=19 and b2=2 
 PICO, get basin neighbors of b=2: b1=1 and b2=3 
 PICO, get basin neighbors of b=3: b1=2 and b2=4 
-PICO, get basin neighbors of b=4: b1=5 and b2=3 
+PICO, get basin neighbors of b=4: b1=8 and b2=0 
 PICO, get basin neighbors of b=5: b1=6 and b2=4 
 PICO, get basin neighbors of b=6: b1=7 and b2=5 
 PICO, get basin neighbors of b=7: b1=8 and b2=6 
@@ -519,7 +519,7 @@
 PICO, get basin neighbors of b=9: b1=10 and b2=8 
 PICO, get basin neighbors of b=10: b1=12 and b2=9 
 PICO, get basin neighbors of b=11: b1=0 and b2=0 
-PICO, get basin neighbors of b=12: b1=13 and b2=10 
+PICO, get basin neighbors of b=12: b1=23 and b2=10 
 PICO, get basin neighbors of b=13: b1=12 and b2=14 
 PICO, get basin neighbors of b=14: b1=13 and b2=15 
 PICO, get basin neighbors of b=15: b1=14 and b2=16 

Note that 5 + 3 = 8 and 13 + 10 = 23. Also note that 23 is not a valid basin number.

This computation is cheap and it is performed only once, when PICO is initialized. It might be best to replace the buggy parallel implementation with a serial one. This almost certainly would not fix this issue, but it would remove one reason why results depend on the domain decomposition.

[talbrecht]

Oh yes, I see the problem...

The variable basins in inputs.nc is constructed to avoid this issue, but we should not require such padding.

Right, we always added a one cell-wide boundary basin with value 0 (in the example input data this boundary has been shifted due to remapping, but still avoiding the periodic boundary effect). I agree, that we just shouldn't have periodic boundary conditions for PICO.

(In short: GlobalSum() is not appropriate here.)

Right, a GlobalSum is not correct, maybe a GlobalMax would do a better job, assuming that there are exactly two neighbors?

It will fail if a basin has more than 2 neighbors. (We don't know if this is a possibility

The if-statement in get_basin_neighbors assumes basin neighbors on the ocean (with number larger than 0) as ocean-ward extension of the drainage basins (more or less following the longitude). One could add to the if-statement, that we only consider ocean cells next to the coast line?

It might be best to replace the buggy parallel implementation with a serial one

This would be worth a try...

PICO, get basin neighbors of b=11: b1=0 and b2=0

I merged basin 11 with 12, so it does not exist any more. But I also think it shouldn't cause any trouble.?!

[ckhroulev]

PICO, get basin neighbors of b=11: b1=0 and b2=0

I merged basin 11 with 12, so it does not exist any more. But I also think it shouldn't cause any trouble.?!

I don't think this is the problem.

This is, though:

mpiexec -n 1 gives this:
PICO, get basin neighbors of b=4: b1=5 and b2=3

mpiexec -n 2 gives this:
PICO, get basin neighbors of b=4: b1=8 and b2=0 

and

mpiexec -n 1:
PICO, get basin neighbors of b=12: b1=13 and b2=10

mpiexec -n 2:
PICO, get basin neighbors of b=12: b1=23 and b2=10 

where b1=23 is invalid. Results should not depend on the number of MPI processes and the parallel domain decomposition.

By the way: do we really have to assume that a basin has at most two neighbors?

[ckhroulev]

This issue appears to be resolved. The code in the current dev branch reports

PICO: basin 1 neighbors: 2, 19
PICO: basin 2 neighbors: 1, 3, 19
PICO: basin 3 neighbors: 2, 4
PICO: basin 4 neighbors: 3, 5
PICO: basin 5 neighbors: 4, 6
PICO: basin 6 neighbors: 5, 7
PICO: basin 7 neighbors: 6, 8
PICO: basin 8 neighbors: 7, 9
PICO: basin 9 neighbors: 8, 10
PICO: basin 10 neighbors: 9, 12
PICO: basin 12 neighbors: 10, 13
PICO: basin 13 neighbors: 12, 14
PICO: basin 14 neighbors: 13, 15
PICO: basin 15 neighbors: 14, 16
PICO: basin 16 neighbors: 15, 17
PICO: basin 17 neighbors: 16, 18
PICO: basin 18 neighbors: 17, 19
PICO: basin 19 neighbors: 1, 2, 18

in all three runs in the test case. Now PICO no longer attempts to divide by zero and no longer produces NaNs -- at least in this test case.

This also answers my question about the possible number of neighbors per basin. See the screenshot below. The cursor shows the location where basins 1, 2, and 19 come together so that basins 2 and 19 have three neighbors.

Results computed by PICO in this test case still depend on the number of MPI processes and the parallel domain decomposition. Since this test case is useful in tracking down that bug as well I'll wait to close this issue until that is taken care of.

[talbrecht]

Great, that the adjacency matrix works.

Within PICO only the ocean conditions next to the ice front (or on the continental shelf) should matter. If we consider only the coast line there should be exactly two neighbor basins left:

+++ b/src/coupler/ocean/PicoGeometry.cc
@@ -550,8 +550,10 @@ std::map<int,std::set<int> > PicoGeometry::basin_neighbors(const IceModelVec2Cel
     auto M = cell_type.star(i, j);
     auto B = basin_mask.star(i, j);
 
+    bool next_to_icefront = (ice_free_ocean(M.ij) and cell_type.next_to_ice(i,j));
+
     // skip the "dummy" basin and cells that are not in the "ice free ocean"
-    if (B.ij == 0 or not ice_free_ocean(M.ij)) {
+    if (B.ij == 0 or not ice_free_ocean(M.ij) or not next_to_icefront) {
       continue;
     }

I would assume, that in this case, this should be also true if the whole continental shelf area is considered instead of the coast line:

bool on_continental_shelf = (continental_shelf_mask.as_int(i, j) == 2);

But this would be tricky, as the basin_neighbors map is already calculated at init() and the continental_shelf_mask only at update().

[ckhroulev]

Results computed by PICO in this test case still depend on the number of MPI processes and the parallel domain decomposition.

Now these differences are small (small enough to explain them by differences in the order of summations used to compute shelf (box, basin) averages). Differences in summation orders are more-or-less unavoidable in a parallel implementation.

The second bug I was referring to is fixed by ef4964e.