From 7ecc1533a546c8af86cd3647422a3fd2127a5631 Mon Sep 17 00:00:00 2001
From: Jan Polcher <jan.polcher@lmd.jussieu.fr>
Date: Wed, 26 Jun 2019 18:09:12 +0200
Subject: [PATCH] Version which now can copy core to halo with large data sets.
They are devided into chuncks.
---
DocumentationInterface | 5 +-
DomainDecompTests/DomainDecomp.py | 2 +-
F90subroutines/routing_interface.f90 | 49 ++++++++++---
F90subroutines/routing_reg.f90 | 30 ++------
Interface.py | 19 ++++-
Partition.py | 104 +++++++++++++++++++--------
6 files changed, 140 insertions(+), 69 deletions(-)
diff --git a/DocumentationInterface b/DocumentationInterface
index 58d11e5..499d46d 100644
--- a/DocumentationInterface
+++ b/DocumentationInterface
@@ -186,14 +186,14 @@ 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)
====================================================================
-basin_area_new,outflow_grid_new = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,fetch_basin,[nbpt,nwbas])
+outflow_grid_new = fetch(basin_count,basin_area,basin_id,basin_hierarchy,basin_fac,outflow_grid,outflow_basin,fetch_basin,[nbpt,nwbas])
Wrapper for ``fetch``.
Parameters
----------
basin_count : input rank-1 array('i') with bounds (nbpt)
-basin_area : input rank-2 array('f') with bounds (nbpt,nwbas)
+basin_area : in/output rank-2 array('f') with bounds (nbpt,nwbas)
basin_id : input rank-2 array('i') with bounds (nbpt,nwbas)
basin_hierarchy : input rank-2 array('f') with bounds (nbpt,nwbas)
basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
@@ -210,7 +210,6 @@ nwbas : input int, optional
Returns
-------
-basin_area_new : rank-2 array('f') with bounds (nbpt,nwbas)
outflow_grid_new : rank-2 array('i') with bounds (nbpt,nwbas)
====================================================================
routing_area,routing_cg,topo_resid,route_togrid,route_tobasin,route_nbintobas,global_basinid,route_outlet,route_type,origin_nbintobas,routing_fetch = 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])
diff --git a/DomainDecompTests/DomainDecomp.py b/DomainDecompTests/DomainDecomp.py
index 8281e5c..9049be0 100644
--- a/DomainDecompTests/DomainDecomp.py
+++ b/DomainDecompTests/DomainDecomp.py
@@ -132,7 +132,7 @@ print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landc
fx = part.zerocore(fx,order='F')
-if np.sum(part.sumcore(fx,order='F')) > 0 :
+if np.sum(part.sumcore(fx, order='F')) > 0 :
print "Error the core region was not set to zero on processor ", rank
sys.exit()
diff --git a/F90subroutines/routing_interface.f90 b/F90subroutines/routing_interface.f90
index 1d6a2e8..4d72ec3 100644
--- a/F90subroutines/routing_interface.f90
+++ b/F90subroutines/routing_interface.f90
@@ -286,8 +286,40 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
!
END SUBROUTINE linkup
+SUBROUTINE areanorm(nbpt, nwbas, basin_count, basin_area, basin_area_norm)
+ !
+ USE grid
+ !
+ !! INPUT VARIABLES
+ INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
+ INTEGER(i_std), INTENT (in) :: nwbas !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
+ !
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_norm !!
+ !
+ !! LOCAL VARIABLES
+ INTEGER(i_std) :: ib, ij
+ REAL(r_std) :: contarea !!
+ REAL(r_std) :: totbasins !!
+ !
+ ! Normalize the area of all basins
+ !
+ DO ib=1,nbpt
+ !
+ totbasins = SUM(basin_area(ib,1:basin_count(ib)))
+ contarea = area(ib)*contfrac(ib)
+ !
+ DO ij=1,basin_count(ib)
+ basin_area_norm(ib,ij) = basin_area(ib,ij)/totbasins*contarea
+ ENDDO
+ !
+ ENDDO
+ !
+END SUBROUTINE areanorm
+
SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
- & outflow_grid, outflow_basin, fetch_basin, basin_area_new, outflow_grid_new)
+ & outflow_grid, outflow_basin, fetch_basin, outflow_grid_new)
!
USE ioipsl
USE grid
@@ -299,30 +331,27 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nwbas !!
- INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_area !!
- INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac
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 !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
!
!! IN-OUTPUT VARIABLES
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
! Output
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: basin_area_new !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(out) :: outflow_grid_new !! Type of outflow on the grid box (unitless)
!
- basin_area_new(:,:) = basin_area(:,:)
outflow_grid_new(:,:) = outflow_grid(:,:)
!
- CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area_new, basin_id, basin_hierarchy, basin_fac,&
+ CALL routing_reg_fetch(nbpt, area, contfrac, nwbas, basin_count, basin_area, basin_id, basin_hierarchy, basin_fac,&
& outflow_grid_new, outflow_basin, fetch_basin)
!
END SUBROUTINE fetch
-
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, 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, routing_area, routing_cg, topo_resid, route_togrid, route_tobasin, route_nbintobas, &
@@ -337,7 +366,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
!
!! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
- INTEGER(i_std), INTENT (in) :: nbasmax
+ INTEGER(i_std), INTENT (in) :: nbxmax_in, nbasmax
INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
diff --git a/F90subroutines/routing_reg.f90 b/F90subroutines/routing_reg.f90
index c2f1c3f..6f4ad01 100644
--- a/F90subroutines/routing_reg.f90
+++ b/F90subroutines/routing_reg.f90
@@ -2703,41 +2703,25 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin
REAL(r_std), DIMENSION(nbpt), INTENT(in) :: contfrac !! Fraction of land in each grid box (unitless;0-1)
!
INTEGER(i_std) :: nwbas !!
- INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
- INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
+ INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac
- 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 !!
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: outflow_basin !!
!
-!! OUTPUT VARIABLES
- REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin !!
+ !! INPUT/OUTPUT VARIABLES
+ INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: outflow_grid !! Type of outflow on the grid box (unitless)
+ REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
!
!! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless)
- REAL(r_std) :: contarea !!
- REAL(r_std) :: totbasins !!
REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area !!
INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !!
!
INTEGER(i_std), PARAMETER :: maxiterations=10000
!_ ================================================================================================================================
- !
- !
- ! Normalize the area of all basins
- !
- DO ib=1,nbpt
- !
- totbasins = SUM(basin_area(ib,1:basin_count(ib)))
- contarea = gridarea(ib)*contfrac(ib)
- !
- DO ij=1,basin_count(ib)
- basin_area(ib,ij) = basin_area(ib,ij)/totbasins*contarea
- ENDDO
- !
- ENDDO
!
!
DO ib=1,nbpt
diff --git a/Interface.py b/Interface.py
index 32f1917..e264c61 100644
--- a/Interface.py
+++ b/Interface.py
@@ -162,9 +162,24 @@ class HydroSuper :
#
def fetch(self, part) :
#
+ fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
self.fetch_basin = np.zeros(self.basin_area.shape, dtype=np.float32, order='F')
- self.basin_area, self.outflow_grid = 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)
+ #
+ self.basin_area = routing_interface.areanorm(self.basin_count, self.basin_area)
+ partial_sum = np.copy(self.basin_area)
+ #
+ for i in range(part.size) :
+ fetch_basin[:,:] = 0.0
+ self.outflow_grid = routing_interface.fetch(self.basin_count, partial_sum, self.basin_id, self.basin_hierarchy, \
+ self.basin_fac, self.outflow_grid, self.outflow_basin, fetch_basin)
+ print("Rank : ", part.rank," iter =", i," Fetch_basin range :", np.min(fetch_basin), np.mean(fetch_basin), np.max(fetch_basin))
+ print("Rank : ", part.rank," Partial sum range :", np.min(partial_sum), np.mean(partial_sum), np.max(partial_sum))
+ partial_sum = np.copy(fetch_basin)
+ part.landsendtohalo(partial_sum, order='F')
+ partial_sum = part.copycore(partial_sum, self.basin_area, order='F')
+ print("Rank :", part.rank," iter =", i," Copy core done")
+
+ self.fetch_basin[:,:] = fetch_basin[:,:]
return
#
#
diff --git a/Partition.py b/Partition.py
index d7dc356..37adb4b 100644
--- a/Partition.py
+++ b/Partition.py
@@ -325,37 +325,52 @@ class partition :
#
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 :
+ # A simple vector
+ #
+ if len(x.shape) == 1 :
+ for i in range(max(self.nbreceive,self.nbsend)) :
+ if i < self.nbsend :
+ if len(self.landinnersend[i]) > 0 :
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()
- if i < self.nbreceive :
- if len(self.landhalosrc[i]) > 0 :
- if len(x.shape) == 1 :
+ if i < self.nbreceive :
+ if len(self.landhalosrc[i]) > 0 :
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()
+ #
+ # Working on matrices
+ #
+ elif len(x.shape) == 2 :
+ chksz = 100
+ if order == 'C' :
+ chkdim = x.shape[0]
+ elif order == 'F' :
+ chkdim = x.shape[1]
+ else :
+ ERROR("Unforessen order of the variable to be sent to halo of land points")
+ sys.exit()
+ chks = range(0, chkdim, chksz)
+ chks.append(chkdim)
+ #
+ for j in range(len(chks)-1) :
+ for i in range(max(self.nbreceive,self.nbsend)) :
+ if i < self.nbsend :
+ if len(self.landinnersend[i]) > 0 :
+ if order == 'C' :
+ self.comm.send(x[chks[j]:chks[j+1],self.landinnersend[i]], dest=self.sendto[i])
+ elif order == 'F' :
+ self.comm.send(x[self.landinnersend[i],chks[j]:chks[j+1]], dest=self.sendto[i])
+ if i < self.nbreceive :
+ if len(self.landhalosrc[i]) > 0 :
+ if order == 'C' :
+ x[chks[j]:chks[j+1],self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
+ elif order == 'F' :
+ x[self.landhalosrc[i],chks[j]:chks[j+1]] = self.comm.recv(source=self.receivefrom[i])
+ #
+ # Some higher ranked variable
+ #
+ else :
+ ERROR("Unforessen rank of the variable to be received in halo of land points")
+ sys.exit()
+
return
#
# Gather all fields partitioned in the 2D domain onto the root proc
@@ -446,7 +461,7 @@ class partition :
ERROR("Unforessen order for variable to be set to zero over core region.")
sys.exit()
#
- y = x
+ y = np.copy(x)
#
if len(x.shape) == 1 :
y[self.landcorelist] = 0.0
@@ -460,6 +475,35 @@ class partition :
sys.exit()
return y
#
+ # Copy all points in the core of y into x to produce z
+ #
+ def copycore(self, x, y, order='C') :
+ if order == 'C' :
+ if x.shape[-1] != self.nbland :
+ ERROR("Space dimension does not correspond to number of land points (case C)")
+ sys.exit()
+ elif order == 'F' :
+ if x.shape[0] != self.nbland :
+ ERROR("Space dimension does not correspond to number of land points (case F)")
+ sys.exit()
+ else :
+ ERROR("Unforessen order for variable to be set to zero over core region.")
+ sys.exit()
+ #
+ z = np.copy(x)
+ #
+ if len(x.shape) == 1 :
+ z[self.landcorelist] = y[self.landcorelist]
+ elif len(x.shape) == 2 :
+ if order == 'C' :
+ z[:,self.landcorelist] = y[:,self.landcorelist]
+ elif order == 'F' :
+ z[self.landcorelist,:] = y[self.landcorelist,:]
+ else :
+ ERROR("Unforessen rank for variable to summedibe set to zero over the core region.")
+ sys.exit()
+ return z
+ #
# Function to sum over the core regions
#
def sumcore(self, x, order='C') :
--
GitLab