Commit 7ecc1533 authored by POLCHER Jan's avatar POLCHER Jan 🚴🏾
Browse files

Version which now can copy core to halo with large data sets. They are devided into chuncks.

parent f729ce16
...@@ -186,14 +186,14 @@ inflow_number : rank-2 array('i') with bounds (nbpt,8 * nbxmax_in) ...@@ -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_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) 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``. Wrapper for ``fetch``.
Parameters Parameters
---------- ----------
basin_count : input rank-1 array('i') with bounds (nbpt) 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_id : input rank-2 array('i') with bounds (nbpt,nwbas)
basin_hierarchy : input rank-2 array('f') 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) basin_fac : input rank-2 array('f') with bounds (nbpt,nwbas)
...@@ -210,7 +210,6 @@ nwbas : input int, optional ...@@ -210,7 +210,6 @@ nwbas : input int, optional
Returns Returns
------- -------
basin_area_new : rank-2 array('f') with bounds (nbpt,nwbas)
outflow_grid_new : rank-2 array('i') 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]) 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])
......
...@@ -132,7 +132,7 @@ print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landc ...@@ -132,7 +132,7 @@ print rank, "Sum over Core regions :", part.sumcore(fx,order='F')/len(part.landc
fx = part.zerocore(fx,order='F') 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 print "Error the core region was not set to zero on processor ", rank
sys.exit() sys.exit()
......
...@@ -286,8 +286,40 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas ...@@ -286,8 +286,40 @@ SUBROUTINE linkup(nbpt, nbxmax_in, nwbas, basin_count, basin_area, basin_id, bas
! !
END SUBROUTINE linkup 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,& 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 ioipsl
USE grid USE grid
...@@ -299,30 +331,27 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy ...@@ -299,30 +331,27 @@ SUBROUTINE fetch(nbpt, nwbas, basin_count, basin_area, basin_id, basin_hierarchy
!! INPUT VARIABLES !! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless)
INTEGER(i_std), INTENT (in) :: nwbas !! INTEGER(i_std), INTENT (in) :: nwbas !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! 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(inout) :: basin_area !!
INTEGER(i_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_id !! 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_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac 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(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 !! IN-OUTPUT VARIABLES
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !! REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: fetch_basin !!
! Output ! 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) 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(:,:) 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) & outflow_grid_new, outflow_basin, fetch_basin)
! !
END SUBROUTINE fetch END SUBROUTINE fetch
SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun, basin_area, basin_cg, basin_topoind,& 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,& & 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, & & 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, ...@@ -337,7 +366,7 @@ SUBROUTINE truncate(nbpt, nbxmax_in, nbasmax, nwbas, basin_count, basin_notrun,
! !
!! INPUT VARIABLES !! INPUT VARIABLES
INTEGER(i_std), INTENT (in) :: nbpt !! Domain size (unitless) 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), INTENT (in) :: nwbas !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt), INTENT(inout) :: basin_count !!
......
...@@ -2703,41 +2703,25 @@ SUBROUTINE routing_reg_fetch(nbpt, gridarea, contfrac, nwbas, basin_count, basin ...@@ -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) 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) :: nwbas !!
INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !! INTEGER(i_std), DIMENSION(nbpt), INTENT(in) :: basin_count !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(inout) :: basin_area !! 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,nwbas), INTENT(in) :: basin_id !!
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_hierarchy
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(in) :: basin_fac 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 !! INPUT/OUTPUT VARIABLES
REAL(r_std), DIMENSION(nbpt,nwbas), INTENT(out) :: fetch_basin !! 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 !! LOCAL VARIABLES
INTEGER(i_std) :: ib, ij, ff(1), it, itt, igrif, ibasf, nboutflow !! Indices (unitless) 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 !! REAL(r_std), DIMENSION(nbpt*nbvmax) :: tmp_area !!
INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !! INTEGER(i_std), DIMENSION(nbpt*nbvmax,2) :: tmpindex !!
! !
INTEGER(i_std), PARAMETER :: maxiterations=10000 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 DO ib=1,nbpt
......
...@@ -162,9 +162,24 @@ class HydroSuper : ...@@ -162,9 +162,24 @@ class HydroSuper :
# #
def fetch(self, part) : 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.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 return
# #
# #
......
...@@ -325,37 +325,52 @@ class partition : ...@@ -325,37 +325,52 @@ class partition :
# #
def landsendtohalo(self, x, order='C') : def landsendtohalo(self, x, order='C') :
# #
for i in range(max(self.nbreceive,self.nbsend)) : # A simple vector
if i < self.nbsend : #
if len(self.landinnersend[i]) > 0 : if len(x.shape) == 1 :
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]) self.comm.send(x[self.landinnersend[i]], dest=self.sendto[i])
elif len(x.shape) == 2 : if i < self.nbreceive :
if order == 'C' : if len(self.landhalosrc[i]) > 0 :
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 :
x[self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i]) x[self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i])
elif len(x.shape) == 2 : #
if order == 'C' : # Working on matrices
x[:,self.landhalosrc[i]] = self.comm.recv(source=self.receivefrom[i]) #
elif order == 'F' : elif len(x.shape) == 2 :
x[self.landhalosrc[i],:] = self.comm.recv(source=self.receivefrom[i]) chksz = 100
else : if order == 'C' :
ERROR("Unforessen order of the variable to be sent to halo of land points") chkdim = x.shape[0]
sys.exit() elif order == 'F' :
else : chkdim = x.shape[1]
ERROR("Unforessen rank of the variable to be received in halo of land points") else :
sys.exit() 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 return
# #
# Gather all fields partitioned in the 2D domain onto the root proc # Gather all fields partitioned in the 2D domain onto the root proc
...@@ -446,7 +461,7 @@ class partition : ...@@ -446,7 +461,7 @@ class partition :
ERROR("Unforessen order for variable to be set to zero over core region.") ERROR("Unforessen order for variable to be set to zero over core region.")
sys.exit() sys.exit()
# #
y = x y = np.copy(x)
# #
if len(x.shape) == 1 : if len(x.shape) == 1 :
y[self.landcorelist] = 0.0 y[self.landcorelist] = 0.0
...@@ -460,6 +475,35 @@ class partition : ...@@ -460,6 +475,35 @@ class partition :
sys.exit() sys.exit()
return y 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 # Function to sum over the core regions
# #
def sumcore(self, x, order='C') : def sumcore(self, x, order='C') :
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment