Prepares and tests functions to echange halo and make some operations on core...

Prepares and tests functions to echange halo and make some operations on core points only (zero and sum).

DocumentationInterface

@@ -186,7 +186,7 @@ inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in)
 inflow_grid : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
 inflow_basin : rank-3 array('i') with bounds (nbpt,8 * nbxmax_in,8 * nbxmax_in)
 ====================================================================
-fetch_basin = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,[nbpt,nwbas])
+fetch_basin_min,fetch_basin_max = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,fetch_basin,[nbpt,nwbas])
 
 Wrapper for ``fetch``.
 
@@ -199,6 +199,7 @@ basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
 basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
 outflow_grid : in/output rank-2 array('i') with bounds (nbpt,nwbas)
 outflow_basin : input rank-2 array('i') with bounds (nbpt,nwbas)
+fetch_basin : in/output rank-2 array('f') with bounds (nbpt,nwbas)
 
 Other Parameters
 ----------------
@@ -209,7 +210,8 @@ nwbas : input int, optional
 
 Returns
 -------
-fetch_basin : rank-2 array('f') with bounds (nbpt,nwbas)
+fetch_basin_min : float
+fetch_basin_max : float
 ====================================================================
 routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas = truncate(nbasmax,basin_count,basin_notrun,basin_area,basin_cg,basin_topoind,fetch_basin,basin_id,basin_coor,basin_type,basin_flowdir,outflow_grid,outflow_basin,inflow_number,inflow_grid,inflow_basin,[nbpt,nbxmax_in,nwbas])
 

DomainDecompTests/DomainDecomp.py

@@ -109,12 +109,38 @@ x[jhoff:jhoff+part.nj,ihoff:ihoff+part.ni]=rank
 lx = modelgrid.landgather(x)
 part.landsendtohalo(lx)
 xland = part.landscatgat(modelgrid, lx)
- 
+
 if rank == 0 :
     xland[~np.isnan(xland)] = 1
     xland[np.isnan(xland)] = 0
     print("Difference between global and gathered version of land only :",np.sum(xland-gg.land))
 
+#
+# Test a 2D landsendtohalo with a FORTRAN ordering
+#
+fx=np.zeros((modelgrid.nbland,nbc), dtype=np.float32, order='F')
+fx[:,:] = rank
+
+part.landsendtohalo(fx, order='F')
+
+if rank == 3 :
+    print "Before Zerored core on proc 1"
+    yy=modelgrid.landscatter(fx[:,3])
+    print yy
+
+print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landcorelist)
+    
+part.zerocore(fx,order='F')
+
+if np.sum(part.sumcore(fx,order='F')) > 0 :
+    print "Error the core region was not set to zero on processor ", rank
+    sys.exit()
+
+if rank == 3 :
+    print "Zerored core on proc 3"
+    yy=modelgrid.landscatter(fx[:,3])
+    print yy
+    
 print("len areas : ", len(modelgrid.area))
 areas = part.landscatgat(modelgrid, np.array(modelgrid.area))
 if rank == 0 :
@@ -125,4 +151,4 @@ y=np.zeros((nbc,modelgrid.nbland))
 print("Test : ", y.shape, modelgrid.nbland)
 y[rank,:]=rank
 
-savefile.dumpnetcdf("TestDomainDecomp.nc", gg, modelgrid, part, rland, x, y)
+savefile.dumpnetcdf("TestDomainDecomp.nc", gg, modelgrid, part, rland, x, y, fx)

DomainDecompTests/savefile.py

@@ -4,7 +4,7 @@ from netCDF4 import Dataset
 
 FillValue=1.0e20
 
-def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland) :
+def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland, fx) :
         #
         vtyp=np.float64
         cornerind=[0,2,4,6]
@@ -87,7 +87,17 @@ def dumpnetcdf(filename, globalgrid, procgrid, part, x1dland, x1dxy, x2dland) :
                 x2dl.standard_name="Some 2D land variable"
                 x2dl.title="2D Land"
                 x2dl[:,:,:] = fullX2D[:,:,:]
-        
+        #
+        # Test a Fortran gathered field
+        #
+        fg = procgrid.landscatter(fx[:,:], order='F')
+        if part.rank == 0 :                           
+            l2id = outnf.createVariable("land2d", vtyp, ('htu','lat','lon'), fill_value=FillValue)
+            l2id.title = "Land 2d in F order"
+            l2id.units = "-"
+        else :
+            l2id = np.zeros((1,1,1))
+        l2id[:,:,:] = part.gather(fg)
         #
         # Proc by land
         if part.rank == 0 :

F90subroutines/routing_interface.f90

@@ -287,7 +287,7 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
 END SUBROUTINE linkup
 
 SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
-     & outflow_grid, outflow_basin, fetch_basin)
+     & outflow_grid, outflow_basin, fetch_basin, fetch_basin_min, fetch_basin_max)
   !
   USE ioipsl
   USE grid
@@ -307,11 +307,16 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
   INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
   !
-  !! OUTPUT VARIABLES
-  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin !!
+  !! IN-OUTPUT VARIABLES
+  REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
+  ! Output
+  REAL(r_std), INTENT(out) :: fetch_basin_min, fetch_basin_max
   
   CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
        & outflow_grid, outflow_basin, fetch_basin)
+  fetch_basin_min = MINVAL(fetch_basin)
+  fetch_basin_max = MAXVAL(fetch_basin)
+  
   !
 END SUBROUTINE fetch
 

F90subroutines/routing_reg.f90

@@ -2738,10 +2738,6 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
        !
     ENDDO
     !
-    ! Compute the area upstream of each basin
-    !
-    fetch_basin(:,:) = zero
-    !
     !
     DO ib=1,nbpt
        !

Interface.py

@@ -161,8 +161,11 @@ class HydroSuper :
         return
     #
     def fetch(self) :
-        self.fetch_basin = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
-                                                   self.basin_fac, self.outflow_grid, self.outflow_basin)
+        #
+        self.fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
+        self.fetch_basin_min, self.fetch_basin_max = routing_interface.fetch(self.basin_count, self.basin_area, self.basin_id, self.basin_hierarchy, \
+                                                                             self.basin_fac, self.outflow_grid, self.outflow_basin, self.fetch_basin)
+        print("## Extrema of basin fetch :", self.fetch_basin_min, self.fetch_basin_max/1.0e6)
         return
 #
 #

ModelGrid.py

@@ -297,7 +297,7 @@ class ModelGrid :
 #
     def landscatter(self, var, order='C') :
         #
-        # Some arrays can be in FORTRAN convention and thus the dimension to be scattered is not the last but the first.
+        # Some arrays can be in FORTRAN convention and thus the dimension to be scattered is not the last but the first.
         #
         dims = var.shape
         transpose=False

Partition.py

@@ -147,6 +147,18 @@ def coresendlist(innersend_map, istart, ihstart, ni, jstart, jhstart, nj) :
         innersend[ic][1][:] = innersend[ic][1][:]+ihoff
     return sendto, innersend
 #
+# All points belonging to the core of the domain
+#
+def listcoreland(istart, ihstart, ni, jstart, jhstart, nj, landind) :
+    corelandlist=[]
+    ihoff=istart-ihstart
+    jhoff=jstart-jhstart
+    for i in range(ni) :
+        for j in range(nj) :
+            if landind[jhoff+j,ihoff+i] >= 0 :
+                corelandlist.append(landind[jhoff+j,ihoff+i])
+    return corelandlist
+#
 # Get list of halo points which need to receive data
 #
 def haloreceivelist(halosource_map, rank) :
@@ -285,6 +297,10 @@ class partition :
             for ic in range(self.nbsend) :
                 wunit.write("To "+str(self.sendto[ic])+" send innersend :"+str(self.innersend[ic])+"\n")
                 wunit.write("To "+str(self.sendto[ic])+" send landinnersend :"+str(self.landinnersend[ic])+"\n")
+        #
+        self.landcorelist = listcoreland(self.istart, self.ihstart, self.ni, self.jstart, self.jhstart, self.nj, landindmap)
+        #
+        return
     #
     # Send halo to the other procs
     #
@@ -307,13 +323,21 @@ class partition :
     #
     # For field gathered in land point send halo to the other procs
     #
-    def landsendtohalo(self, x) :
+    def landsendtohalo(self, x, order='C') :
         #
         for i in range(max(self.nbreceive,self.nbsend)) :
             if i < self.nbsend :
                 if len(self.landinnersend[i]) > 0 :
                     if len(x.shape) == 1 :
                         self.comm.send(x[self.landinnersend[i]], dest=self.sendto[i])
+                    elif len(x.shape) == 2 :
+                        if order == 'C' :
+                            self.comm.send(x[:,self.landinnersend[i]], dest=self.sendto[i])
+                        elif order == 'F' :
+                            self.comm.send(x[self.landinnersend[i],:], dest=self.sendto[i])
+                        else :
+                            ERROR("Unforessen order of the variable to be sent to halo of land points")
+                            sys.exit()
                     else :
                         ERROR("Unforessen rank of the variable to be sent to halo of land points")
                         sys.exit()
@@ -321,6 +345,14 @@ class partition :
                 if len(self.landhalosrc[i]) > 0 :
                     if len(x.shape) == 1 :
                         x[self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
+                    elif len(x.shape) == 2 :
+                        if order == 'C' :
+                            x[:,self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
+                        elif order == 'F' :
+                            x[self.landhalosrc[i],:] = self.comm.recv(source=self.receivefrom[i])
+                        else :
+                            ERROR("Unforessen order of the variable to be sent to halo of land points")
+                            sys.exit()
                     else :
                         ERROR("Unforessen rank of the variable to be received in halo of land points")
                         sys.exit()
@@ -398,3 +430,34 @@ class partition :
             ERROR("The first dimension does not have the length of the number of land points")
             sys.exit()
         return y
+    #
+    # Set to zero all points in the core
+    #
+    def zerocore(self, x, order='C') :
+        if len(x.shape) == 1 :
+            x[self.landcorelist] = 0.0
+        elif len(x.shape) == 2 :
+            if order == 'C' :
+                x[:,self.landcorelist] = 0.0
+            elif order == 'F' :
+                x[self.landcorelist,:] = 0.0
+            else :
+                ERROR("Unforessen order for variable to zero")
+                sys.exit()
+        
+        return
+    #
+    # Function to sum over the core regions
+    #
+    def sumcore(self, x, order='C') :
+        if len(x.shape) == 1 :
+            y = np.nansum(x[self.landcorelist])
+        elif len(x.shape) == 2 :
+            if order == 'C' :
+               y =  np.nansum(x[:,self.landcorelist], axis=1)
+            elif order == 'F' :
+               y = np.nansum(x[self.landcorelist,:], axis=0)
+            else :
+                ERROR("Unforessen order for variable to zero")
+                sys.exit()
+        return y