[mpas-developers] Fix for isolated cells

[Mark Petersen]

mpas developers,

We have had trouble with our global domains with land, and I have a
bug fix that solves the problem.  The error occurs on runs with more
than one processor, where a processor has an isolated ocean cell.  For
example:

1 1 1 L L L 2 2 2
1 1 1 L L 1 2 2 2
1 1 1 L L L 2 2 2

These are cells owned by processor 1 and 2, and land cells L are shown
but are not really cells.

In module block_decomp, subroutine block_decomp_add_halo, the section

       do i=1,local_graph_info % nVertices
          do j=1,local_graph_info % nAdjacent(i)
             if (local_graph_info % adjacencyList(j,i) /= 0) then
                if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then
                   call hash_insert(h, local_graph_info % adjacencyList(j,i))
                end if
             end if
          end do
       end do

only adds cells to the hash table if they have neighbors owned by that
processor.  Thus isolated cells are not added,

       local_graph_with_halo % nVerticesTotal = hash_size(h)

is too small, the warning

block_decomp_add_halo: Somehow we don't have the right number of total 
cells.

appears, and we end up reading outside of memory bounds on arrays like
local_graph_with_halo % vertexID.

My fix is to add all cells to the hash table first, and then add
neighbors of cells:

lo1-fe> pwd
/usr/projects/climate/mpeterse/coy/mpas_100603/trunk/mpas/src/framework
lo1-fe> svn diff module_block_decomp.F
Index: module_block_decomp.F
===================================================================
--- module_block_decomp.F       (revision 528)
+++ module_block_decomp.F       (working copy)
@@ -236,6 +236,10 @@
        call hash_init(h)

        do i=1,local_graph_info % nVertices
+         call hash_insert(h, local_graph_info % vertexID(i))
+      end do
+
+      do i=1,local_graph_info % nVertices
           do j=1,local_graph_info % nAdjacent(i)
              if (local_graph_info % adjacencyList(j,i) /= 0) then
                 if (.not. hash_search(h, local_graph_info % adjacencyList(j,i))) then

A more elegant solution would be for the domain decomposition to not
have isolated cells, but kmetis handles this and I don't know how to
change it.  Also, users should be able to choose any grid and
decomposition they like without unnecessary errors.  A warning would be 
more appropriate.

I will make the above change tomorrow, unless I hear other feedback.

Mark