Added a __main__ to only compute the weight. Usefull for dealing with big problems.

Partition.py

@@ -62,8 +62,6 @@ def adjustpart(partin, land, nbcore) :
         addnbland(partin, land)
         
     nbptproc=np.array([dom['nbi']*dom['nbj'] for dom in partin])
-    INFO("Nb of land points per proc : Mean = "+str(np.mean(nbptproc))+" Min = "+str(np.min(nbptproc))+\
-         " Max = "+str(np.max(nbptproc))+" Var :"+str((np.max(nbptproc)-np.min(nbptproc))/np.mean(nbptproc)*100)+"%")
     if np.min(nbptproc) < 3 :
         ERROR("The smallest domain in the partition is of a size less than 5. This will probably not work.")
         sys.exit()        
@@ -71,7 +69,8 @@ def adjustpart(partin, land, nbcore) :
 #
 # Partition domain in two steps : 1) halving, 2) adjusting by nb land points
 #
-def partitiondom(nig, njg, land, nbcore) :
+def partitiondom(land, nbcore, rank) :
+    njg,nig = land.shape
     partout=[{"nbi":nig,"nbj":njg,"istart":0,"jstart":0, "nbland":int(np.sum(land))}]
     while len(partout) <= nbcore/2 :
         partout = halfpartition(partout)
@@ -79,14 +78,96 @@ def partitiondom(nig, njg, land, nbcore) :
     #
     adjustpart(partout, land, nbcore)
     #
+    gindland=np.zeros((njg,nig), dtype=np.int32)
+    gindland[:,:]=-1
+    n=0
+    for i in range(nig) :
+        for j in range(njg) :
+            if (land[j,i] > 0 ) :
+                gindland[j,i] = n
+                n += 1
+    #
+    for ic in range(len(partout)) :
+        istart = partout[ic]["istart"]
+        nbi = partout[ic]["nbi"]
+        jstart = partout[ic]["jstart"]
+        nbj = partout[ic]["nbj"]
+        locgind=np.copy(gindland[jstart:jstart+nbj,istart:istart+nbi])
+        loccore=np.copy(gindland[jstart:jstart+nbj,istart:istart+nbi])
+        loccore[loccore >= 0] = ic
+        partout[ic]["glandindex"] = locgind
+        partout[ic]["corenb"] = loccore
+        partout[ic]["nblocland"] = np.count_nonzero(locgind >= 0)
+    #
+    return partout
+#
+# Build global map of processors
+#
+def buildprocmap(nig, njg, part) :
     procmap=np.zeros((njg,nig), dtype=np.int8)
     procmap[:,:] = -1
-    for i in range(len(partout)) :
-        for ij in range(partout[i]["nbj"]) :
-            for ii in range(partout[i]["nbi"]) :
-                procmap[ij+partout[i]["jstart"],ii+partout[i]["istart"]] = i
-
-    return partout, procmap
+    for i in range(len(part)) :
+        for ij in range(part[i]["nbj"]) :
+            for ii in range(part[i]["nbi"]) :
+                procmap[ij+part[i]["jstart"],ii+part[i]["istart"]] = i
+    #
+    return procmap
+#
+# Partion of only land points in the case where we do not need halos
+#
+def landpartition(land, nbcore, rank) :
+    nbland =  np.count_nonzero(land > 0)
+    optland = int(nbland/nbcore)+(nbland%nbcore > 0)
+    print("RANK",rank,"Nbland = ", np.count_nonzero(land > 0), " nbcore = ", nbcore, " optland = ",optland)
+    #
+    distppc=np.zeros(nbcore)
+    distppc[:]=0
+    for i in range(nbcore-1) :
+        ptsleft = int(nbland-np.sum(distppc[0:i]))
+        crsleft = nbcore-i
+        if  ptsleft > optland and ptsleft/2 > crsleft :
+            distppc[i] = optland
+        else :
+            print("On proc ", i," Still have ", ptsleft, " points to distribute and ", crsleft, " cores")
+            distppc[i] = int(ptsleft/crsleft)
+        #
+    distppc[-1]=nbland-np.sum(distppc[0:-1])
+    print("Distpcc :", distppc)
+    #
+    njg,nig=land.shape
+    gindland=np.zeros((njg,nig), dtype=np.int32)
+    gindland[:,:]=-1
+    glandcore=np.zeros((njg,nig), dtype=np.int32)
+    glandcore[:,:]=-1
+    n=0
+    for i in range(nig) :
+        for j in range(njg) :
+            if (land[j,i] > 0 ) :
+                gindland[j,i] = n
+                cc = np.argmax(distppc > 0)
+                glandcore[j,i] = cc
+                distppc[cc] = distppc[cc]-1
+                n += 1
+    #
+    # Decompose domain according to land point positions
+    #
+    partout=[]
+    for ic in range(nbcore) :
+        ff = np.nonzero(glandcore == ic)
+        # Region is at least 3x3
+        nbi=max(max(ff[1])-min(ff[1])+1, 3)
+        nbj=max(max(ff[0])-min(ff[0])+1, 3)
+        istart=min(min(ff[1]),nig-nbi)
+        jstart=min(min(ff[0]),njg-nbj)
+        #
+        loccore=np.copy(glandcore[jstart:jstart+nbj,istart:istart+nbi])
+        loccore[loccore != ic] = -1
+        locgind=np.copy(gindland[jstart:jstart+nbj,istart:istart+nbi])
+        locgind[loccore < 0] = -1
+        partout.append({"nbi":nbi, "nbj":nbj, "istart":istart, "jstart":jstart, "nbland":ff[0].shape[0], \
+                        "corenb":loccore, "glandindex":locgind, "nblocland":np.count_nonzero(locgind >= 0)})
+    #
+    return partout
 #
 # Add halo to the partition
 #
@@ -217,9 +298,6 @@ def toland_index(x,landmap) :
 #
 class partition :
     def __init__ (self, nig, njg, land, mpicomm, nbcore, halosz, rank, wunit="None") :
-        #
-        part, procmap = partitiondom(nig, njg, land, nbcore)
-        halosource_map, innersend_map = addhalo(nig, njg, part, procmap, halosz)
         #
         # Self varibales
         #
@@ -231,12 +309,27 @@ class partition :
         self.nig = nig
         self.njg = njg
         self.nblandg = np.count_nonzero(land > 0)
+        #
+        part = partitiondom(land, nbcore, rank)
+        procmap = buildprocmap(nig, njg, part)
+        # To be activated when halosz = 0 and we can deal with single point domain grids.
+        #
+        #     part = landpartition(land, nbcore, rank)
+        #     procmap = buildprocmap(nig, njg, part)
+        #
+        halosource_map, innersend_map = addhalo(nig, njg, part, procmap, halosz)
+        #
+        INFO(" Nb of land points on proc "+str(rank)+" = "+str(part[rank]["nblocland"])+" from "+str(self.nblandg)+\
+         " -- "+"%.2f"%(part[rank]["nblocland"]/self.nblandg*100)+"% of total.")
+        #
         self.allistart = []
         self.alljstart = []
         for i in range(self.size) :
             self.allistart.append(part[i]["istart"])
             self.alljstart.append(part[i]["jstart"])
-        # Local info
+        #
+        # Information of domain on local processor
+        #
         self.ni = part[rank]["nbi"]
         self.nj = part[rank]["nbj"]
         self.istart = part[rank]["istart"]

RoutingPreProc.py

@@ -27,7 +27,6 @@ import ModelGrid as MG
 import Partition as PA
 import RPPtools as RPP
 import HydroGrid as HG
-import diagplots as DP
 import Interface as IF
 from Truncate import trunc as TR
 #

WeightsOnly.py

+#
+#
+#
+import numpy as np
+import os, sys
+from mpi4py import MPI
+import gc
+#
+import matplotlib as mpl
+mpl.use('Agg')
+from matplotlib.backends.backend_pdf import PdfPages
+import matplotlib.pyplot as plt
+import time
+#
+# Gert the information from the configuration file.
+#
+import configparser
+config=configparser.ConfigParser()
+config.read("run.def")
+nbasmax=config.getint("OverAll", "nbasmax")
+numop=config.getint("OverAll", "numop", fallback=100)
+OutGraphFile=config.get("OverAll","GraphFile")
+DumpHydroSuper=config.getboolean("OverAll","DumpHydroSuper",fallback=False)
+wfile=config.get("OverAll","WeightFile",fallback="Weights.nc")
+#
+import ModelGrid as MG
+import Partition as PA
+import RPPtools as RPP
+import HydroGrid as HG
+#
+# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
+#
+import getargs
+log_master, log_world = getargs.getLogger()
+INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
+INFO_ALL, DEBUG_ALL = log_world.info, log_world.debug
+#
+# Read full grid and partition the domain.
+#
+gg = MG.GlobalGrid()
+#
+comm = MPI.COMM_WORLD
+szhalo=0
+nbcore=comm.Get_size()
+rank=comm.Get_rank()
+#
+# Region of grid to be treated
+#
+part = PA.partition(gg.ni, gg.nj, gg.land, comm, nbcore, szhalo, rank)
+#
+modelgrid=MG.ModelGrid(part.ihstart+gg.igstart,part.nih,part.jhstart+gg.jgstart,part.njh)
+INFO("Longitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[0][0])+" : "+str(modelgrid.box_land[0][1]))
+INFO("Latitude interval on proc "+str(rank)+" = "+str(modelgrid.box_land[1][0])+" : "+str(modelgrid.box_land[1][1]))
+#
+hydrogrid=HG.HydroGrid(modelgrid.box_land)
+#
+# Computes weights of overlap of modelgrid and hydrogrid
+#
+w = RPP.compweights(wfile, part, modelgrid, hydrogrid)
+#
+#
+