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

A first versino which write the weights into a netCDF files.

I verified that we can read the weight with another processor partition as the one use to compute them.
Small differences appear because of different rounding errors on the individual processors and the core proc who gives his values to the file.
parent 8b2082f5
......@@ -58,7 +58,9 @@ def corners(lon, lat) :
def gather(x, index) :
y=[]
for ia in index :
y.append(list(x[ia[0][i],ia[1][i]] for i in range(len(ia[0]))))
# print("XXXX shape of ia :", ia.shape, ia.dtype)
# print("XXXX ia :", ia)
y.append(list(x[ia[0,i],ia[1,i]] for i in range(ia.shape[1]) ))
return y
#
def getattrcontaining(nc, varname, substr) :
......@@ -72,9 +74,9 @@ def getattrcontaining(nc, varname, substr) :
class HydroGrid :
def __init__(self, lolacorners) :
#
filename=config.get("OverAll", "HydroFile")
INFO("Opening in HydroGrid :"+filename)
self.ncfile=Dataset(filename,'r')
self.source=config.get("OverAll", "HydroFile")
INFO("Opening in HydroGrid : "+self.source)
self.ncfile=Dataset(self.source,'r')
istr, iend, jstr, jend = getbox(self.ncfile, lolacorners)
self.box=[istr, iend, jstr, jend]
self.lon=np.copy(self.ncfile.variables["nav_lon"][jstr:jend,istr:iend])
......
......@@ -236,7 +236,7 @@ class HydroOverlap :
sub_lon[ib,0:sub_pts[ib]] = sub_lon_in[ib][:]
sub_lat[ib,0:sub_pts[ib]] = sub_lat_in[ib][:]
for ip in range(sub_pts[ib]) :
sub_index[ib,ip,:] = [sub_index_in[ib][0][ip],sub_index_in[ib][1][ip]]
sub_index[ib,ip,:] = sub_index_in[ib][:,ip]
#
part.landsendtohalo(np.array(sub_area), order='F')
#
......@@ -567,7 +567,7 @@ class HydroGraph :
if orig_type == "float":
var[np.isnan(var)] = NCFillValue
elif orig_type == "int":
var[var>=RPP.IntFillValue] = NCFillValue
var[var>=np.abs(RPP.IntFillValue)] = NCFillValue
if part.rank == 0:
ncvar = outnf.createVariable(name, vtyp, coord, fill_value=NCFillValue)
......
......@@ -233,8 +233,8 @@ class ModelGrid :
ERROR("Domain too small")
sys.exit()
#
filename=config.get("OverAll", "ModelGridFile")
geo=Dataset(filename,'r')
self.source = config.get("OverAll", "ModelGridFile")
geo=Dataset(self.source,'r')
#
# Get the coordinates from the grid file.
#
......@@ -413,9 +413,9 @@ class GlobalGrid :
lonrange=np.array(config.get("OverAll", "WEST_EAST").split(","),dtype=float)
latrange=np.array(config.get("OverAll", "SOUTH_NORTH").split(","),dtype=float)
filename=config.get("OverAll", "ModelGridFile")
INFO("Opening :"+filename)
geo=Dataset(filename,'r')
self.source=config.get("OverAll", "ModelGridFile")
INFO("Opening : "+self.source)
geo=Dataset(self.source,'r')
#
griddesc, lon_full, lat_full = getcoordinates(geo, 0, -1, 0, -1)
#
......
......@@ -201,14 +201,6 @@ def tolocal_index(x, istart, jstart) :
xout.append([c[0]-jstart,c[1]-istart])
return xout
#
# Convert indices to localglobal
#
def toglobal_index(x, istart, jstart) :
xout = []
for c in x :
xout.append([c[0]+jstart,c[1]+istart])
return xout
#
#
#
def toland_index(x,landmap) :
......@@ -238,6 +230,7 @@ class partition :
# Global info
self.nig = nig
self.njg = njg
self.nblandg = np.count_nonzero(land > 0)
self.allistart = []
self.alljstart = []
for i in range(self.size) :
......@@ -289,7 +282,7 @@ class partition :
self.sendto,self.innersend = coresendlist(innersend_map, self.istart, self.ihstart, self.ni, self.jstart, self.jhstart, self.nj)
self.nbsend = len(self.sendto)
# To global indexing
self.innersend_g = toglobal_index(self.innersend, self.ihstart, self.jhstart)
self.innersend_g = self.toglobal_index(self.innersend)
self.landinnersend = toland_index(self.innersend, landindmap)
if wunit != "None" :
wunit.write("+++++ Halo core send list : "+str(self.receivefrom)+"\n")
......@@ -302,6 +295,14 @@ class partition :
#
return
#
# Convert indices from local to global
#
def toglobal_index(self, x) :
xout = []
for c in x :
xout.append([c[0]+self.jhstart,c[1]+self.ihstart])
return xout
#
# Send halo to the other procs
#
def sendtohalo(self, x) :
......@@ -548,3 +549,22 @@ class partition :
#
def domainmax(self, x) :
return max(self.comm.allgather(x))
#
# Function to gather vectors of points on different procs on root
#
def gatherbypoint(self, x) :
maxlen = self.domainmax(len(x))
if self.rank == 0 :
xout = np.empty((maxlen,self.size))
# Receive from all procs
for lr in range(self.size) :
if lr == 0 :
xout[:,lr] = np.pad(x, [0,maxlen-len(x)], mode='constant', constant_values=np.nan)
else :
xout[:,lr] = self.comm.recv(source=lr)
else :
xout = None
self.comm.send(np.pad(x, [0,maxlen-len(x)], mode='constant', constant_values=np.nan), dest=0)
#
return xout
......@@ -4,8 +4,11 @@ from numba import jit
import os, sys
#
FillValue=np.nan
IntFillValue=999999
IntFillValue=-999999
vtyp=np.float64
vinttyp=np.int32
#
from netCDF4 import Dataset
import pickle
from spherical_geometry import vector
#
......@@ -77,9 +80,7 @@ def maxradius(cent, lon, lat) :
#
def dumpfield(x, lon, lat, filename, varname) :
#
from netCDF4 import Dataset
NCFillValue=1.0e20
vtyp=np.float64
#
print("Dumping overlap into file :", filename)
i=np.nonzero(np.mean(x,axis=0))
......@@ -113,72 +114,272 @@ def dumpfield(x, lon, lat, filename, varname) :
# Compute the weights of the overlap of modelgrif and hydrogrid
# Only if there is no file which already contains the information.
#
def compweights(wfile, modelgrid, hydrogrid) :
icell = 0
sub_index = []
sub_area = []
sub_lon = []
sub_lat = []
class compweights :
def __init__(self, wfile, part, modelgrid, hydrogrid) :
#
sub_index = []
sub_area = []
sub_lon = []
sub_lat = []
#
# Only compute the weights if the file does not exist.
#
if not os.path.exists(wfile) :
INFO("Computing weights")
for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
selected = hydrogrid.select(center, radius)
for isel in selected :
hydrocell = hydrogrid.polylist[isel]
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
# Another strange behaviour of SphericalHarmonics. There should be an overlap but the intersection polygone is empty.
if len(inter.polygons) > 0 :
inside = inter.polygons[0]._inside
if hydrocell.contains_point(inside) and cell.contains_point(inside):
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
##############
# Output of Overlap areas for each grid point
#
#ncfile="Area_"+str("%.2f"%llinside[0])+"_"+str("%.2f"%llinside[1])+"_"+str(rank+1)+"of"+str(nbcore)+".nc"
#RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
#
sub_index.append(np.array(np.where(area_in > 0.0)))
sub_area.append(np.extract(area_in > 0.0, area_in))
sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
#
# Save information in class variables.
#
self.index = sub_index
self.area = sub_area
self.lon = sub_lon
self.lat = sub_lat
self.hpts = [l.shape[0] for l in self.lon]
self.maxhpts = max(self.hpts)
for icell in range(len(self.lon)) :
INFO(str(icell)+" Sum of overlap "+str(np.sum(self.area[icell])/modelgrid.area[icell])+
" Nb of Hydro grid overlaping : "+str(self.hpts[icell]))
#
self.dumpweights(wfile, part, modelgrid, hydrogrid)
#
##self.printtest([15,19], part, modelgrid)
#
else :
#
INFO("Reading weights from "+wfile)
self.fetchweight(wfile, part, modelgrid, hydrogrid)
self.hpts = [l.shape[0] for l in self.lon]
self.maxhpts = max(self.hpts)
for icell in range(len(self.lon)) :
INFO(str(icell)+" Sum of overlap "+str(np.sum(self.area[icell])/modelgrid.area[icell])+
" Nb of Hydro grid overlaping : "+str(self.hpts[icell]))
#
##self.printtest([15,19], part, modelgrid)
#
return
#
# Only compute the weights if the file does not exist.
# Function to fetch the weights from weight file
#
if not os.path.exists(wfile) :
INFO("Computing weights")
for cell, center, radius, area in zip(modelgrid.polylist, modelgrid.centers, modelgrid.radius, modelgrid.area) :
llinside=vector.vector_to_lonlat(cell.polygons[0].inside[0],cell.polygons[0].inside[1],cell.polygons[0].inside[2])
area_in=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float)
selected = hydrogrid.select(center, radius)
INFO(str(icell)+" Number of HydroCells selected : "+str(len(selected))+ " out of "+str(len(hydrogrid.index)))
for isel in selected :
hydrocell = hydrogrid.polylist[isel]
index = hydrogrid.index[isel]
if cell.intersects_poly(hydrocell):
inter = cell.intersection(hydrocell)
# Another strange behaviour of SphericalHarmonics. There should be an overlap but the intersection polygone is empty.
if len(inter.polygons) > 0 :
inside = inter.polygons[0]._inside
if hydrocell.contains_point(inside) and cell.contains_point(inside):
area_in[index[0],index[1]] = inter.area()*(EarthRadius**2)
else:
area_in[index[0],index[1]] = (4*np.pi-inter.area())*(EarthRadius**2)
##############
# Output of Overlap areas for each grid point
def fetchweight(self, weightfile, part, modelgrid, hydrogrid) :
#
self.area = []
self.lon = []
self.lat = []
#
innf=Dataset(weightfile, 'r')
#
gindex=innf.variables["global_index"][:,:]
indexfill = innf.variables["global_index"]._FillValue
landind = part.toglobal_index(modelgrid.indP)
#
for pts in landind :
i = self.findinfile(pts, gindex)
maxhpts = np.array(np.where(innf.variables["hydro_index"][0,:,i] >= 0)).shape[1]
self.area.append(innf.variables["hydro_area"][0:maxhpts,i])
self.lon.append(innf.variables["hydro_lon"][0:maxhpts,i])
self.lat.append(innf.variables["hydro_lat"][0:maxhpts,i])
#
innf.close()
# Find the indicis of the points on the HydroGrid using the coordinates.
self.index = self.findhindex(self.lon, self.lat, hydrogrid)
#
return
#
# Function to find the indicis of the points on the hydrological grid for the current
# partition of the domain.
#
def findhindex(self, lon, lat, hydrogrid) :
ind=[]
for im in range(len(lon)) :
maxhpts = lon[im].shape[0]
hind=np.zeros((2,maxhpts))
for ih in range(maxhpts) :
dist=np.sqrt((hydrogrid.lon-lon[im][ih])**2+(hydrogrid.lat-lat[im][ih])**2)
hind[0,ih],hind[1,ih] = np.unravel_index(np.argmin(dist), dist.shape)
#
ind.append(np.array(hind, dtype=np.int32))
return ind
#
# TEST function to be deleted later
#
def printtest(self, printpts, part, modelgrid) :
landind = part.toglobal_index(modelgrid.indP)
ip = IntFillValue
for i in range(len(landind)) :
if printpts[0] == landind[i][0] and printpts[1] == landind[i][1] :
ip = i
if ip >= 0 :
print("H-Report on point ", landind[ip][:], " with ",self.index[ip].shape[1]," Hydro points.")
print("H- indices : ", self.index[ip].astype(int))
print("H- total area :", modelgrid.area[ip], np.sum(self.area[ip]), "Error =",modelgrid.area[ip]/np.sum(self.area[ip]))
print("H- areas :", self.area[ip])
print("H- lon :", self.lon[ip])
print("H- lat :", self.lat[ip])
return
#
# Function to find the index of the point in the file
#
def findinfile(self, pts, indlist) :
ib = IntFillValue
for i in range(indlist.shape[1]) :
if all(indlist[:,i] == pts) :
ib = i
return ib
#
# Function to dump all the weights into a netCDF file
#
def dumpweights(self, weightfile, part, modelgrid, hydrogrid) :
#
# This number needs to come from compweights
#
maxhydropts=part.domainmax(self.maxhpts)
landind = part.toglobal_index(modelgrid.indP)
#
# On the root proc set-up the netCDF file with all the variables needed.
#
if part.rank == 0 :
nf = self.generatenetcdffile(weightfile, part.nig, part.njg, part.nblandg, \
maxhydropts, modelgrid.source, hydrogrid.source)
#
#
#
maxlandpoint = part.domainmax(len(part.landcorelist))
print("Maximum points of Land points on any proc :", maxlandpoint, " Here : ", len(part.landcorelist), "Global Nbland = ", part.nblandg)
#
# Send each point on this proc to the roor_proc
#
ni = 0
for i in range(maxlandpoint) :
if i < len(part.landcorelist) :
# If this is a core point then send it.
glandind = part.gatherbypoint(landind[part.landcorelist[i]][:])
ghindj = part.gatherbypoint(self.index[part.landcorelist[i]][0,:])
ghindi = part.gatherbypoint(self.index[part.landcorelist[i]][1,:])
ghlon = part.gatherbypoint(self.lon[part.landcorelist[i]][:])
ghlat = part.gatherbypoint(self.lat[part.landcorelist[i]][:])
gharea = part.gatherbypoint(self.area[part.landcorelist[i]][:])
else :
# If it is not a core point then send nothing to the root_proc
glandind = part.gatherbypoint([])
ghindi = part.gatherbypoint([])
ghindj = part.gatherbypoint([])
ghlon = part.gatherbypoint([])
ghlat = part.gatherbypoint([])
gharea = part.gatherbypoint([])
#
#ncfile="Area_"+str("%.2f"%llinside[0])+"_"+str("%.2f"%llinside[1])+"_"+str(rank+1)+"of"+str(nbcore)+".nc"
#RPP.dumpfield(area_in, hydrogrid.lon, hydrogrid.lat, ncfile, "Intersect area")
# On the root proc which collected all the data we write data to the netCDF file
#
INFO(str(icell)+" Sum of overlap "+str(np.sum(area_in)/area)+
" Nb of Hydro grid overlaping : "+str(np.count_nonzero(area_in)))
sub_index.append(np.where(area_in > 0.0))
sub_area.append(np.extract(area_in > 0.0, area_in))
sub_lon.append(np.extract(area_in > 0.0, hydrogrid.lon))
sub_lat.append(np.extract(area_in > 0.0, hydrogrid.lat))
icell = icell + 1
#
# Save the weights to a dedicated file for future use
#
if ( saveweights.lower() == 'true' ) :
with open(wfile, 'wb') as fp:
pickle.dump(sub_index, fp)
pickle.dump(sub_area, fp)
pickle.dump(sub_lon, fp)
pickle.dump(sub_lat, fp)
fp.close()
else :
#
# Extract the weights from the file if it exists
#
INFO("Reading weights")
with open (wfile, 'rb') as fp:
sub_index = pickle.load(fp)
sub_area = pickle.load(fp)
sub_lon = pickle.load(fp)
sub_lat = pickle.load(fp)
fp.close()
#
#
#
return sub_lon, sub_lat, sub_index, sub_area
if part.rank == 0 :
ni = self.addtoncvar(ni, part.nig, part.njg, glandind, ghindi, ghindj, ghlon, ghlat, gharea)
#
if part.rank == 0 :
nf.close()
#
return
#
# Function to generate the netCDF file for the weight functions
#
def generatenetcdffile(self, weightfile, nig, njg, nbland, maxhydropts, mfile, hfile) :
#
outnf=Dataset(weightfile, 'w', format='NETCDF4_CLASSIC')
# Dimensions
print("Generating dimensions : ", nig, njg, nbland, maxhydropts)
outnf.createDimension('i', nig)
outnf.createDimension('j', njg)
outnf.createDimension('k', 2)
outnf.createDimension('nbland', nbland)
outnf.createDimension('nover', maxhydropts)
# Global attributes
outnf.HydroGrid=hfile
outnf.ModelGrid=mfile
#
self.mask = outnf.createVariable("mask", vtyp, ('j','i'), fill_value=FillValue)
self.mask.units="-"
self.mask.standard_name="Land sea mask"
self.mask.title="Land sea mask"
self.mask[:,:] = 0
#
self.proc = outnf.createVariable("proc_num", vtyp, ('j','i'), fill_value=FillValue)
self.proc.units="-"
self.proc.standard_name="Processor working on land point"
self.proc.title="Processor working on land point"
self.proc[:,:] = FillValue
#
self.gindex = outnf.createVariable("global_index", vinttyp, ('k','nbland'), fill_value=IntFillValue)
self.gindex.units="-"
self.gindex.standard_name="Index j,i on Global Grid"
self.gindex.title="Index j,i on Global Grid"
self.gindex[:,:] = IntFillValue
#
self.hindex = outnf.createVariable("hydro_index", vinttyp, ('k','nover','nbland'), fill_value=IntFillValue)
self.hindex.units="-"
self.hindex.standard_name="Index of HydroGrid point contributing to current land point"
self.hindex.title="Index of HydroGrid point contributing to current land point"
self.hindex[:,:,:] = IntFillValue
#
self.hlon = outnf.createVariable("hydro_lon", vtyp, ('nover','nbland'), fill_value=FillValue)
self.hlon.units="degrees_east"
self.hlon.standard_name="Longitude of HydroGrid point contributing to current land point"
self.hlon.title="Longitude of HydroGrid point contributing to current land point"
self.hlon[:,:] = FillValue
#
self.hlat = outnf.createVariable("hydro_lat", vtyp, ('nover','nbland'), fill_value=FillValue)
self.hlat.units="degrees_north"
self.hlat.standard_name="Latitude of HydroGrid point contributing to current land point"
self.hlat.title="Latitude of HydroGrid point contributing to current land point"
self.hlat[:,:] = FillValue
#
self.harea = outnf.createVariable("hydro_area", vtyp, ('nover','nbland'), fill_value=FillValue)
self.harea.units="m^2"
self.harea.standard_name="Area of HydroGrid point contributing to current land point"
self.harea.title="Area of HydroGrid point contributing to current land point"
self.harea[:,:] = FillValue
#
return outnf
#
# Add variables
#
def addtoncvar(self, nib, nig, njg, landind, hindi, hindj, hlon, hlat, harea) :
#
mxval,nbp=landind.shape
for ib in range(nbp) :
if ~np.isnan(landind[0,ib]) :
j,i=landind[:,ib]
self.mask[j,i] = 1.0
self.proc[j,i] = ib
self.gindex[0,nib] = j
self.gindex[1,nib] = i
self.hindex[0,:,nib] = hindj[:,ib]
self.hindex[1,:,nib] = hindi[:,ib]
self.hlon[:,nib] = hlon[:,ib]
self.hlat[:,nib] = hlat[:,ib]
self.harea[:,nib] = harea[:,ib]
nib = nib+1
return nib
......@@ -15,14 +15,15 @@ import time
# Gert the information from the configuration file.
#
import configparser
config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0", "numop":'100'})
config=configparser.ConfigParser({"DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"})
config.read("run.def")
nbasmax=config.getint("OverAll", "nbasmax")
numop=config.getint("OverAll", "numop")
numop=config.getint("OverAll", "numop", fallback=100)
OutGraphFile=config.get("OverAll","GraphFile")
DumpHydroSuper=config.getboolean("OverAll","DumpHydroSuper",fallback=False)
lonint=np.array(config.get("OverAll", "DiagLon").split(","),dtype=float)
latint=np.array(config.get("OverAll", "DiagLat").split(","),dtype=float)
wfile=config.get("OverAll","WeightFile",fallback="Weights.nc")
#
import ModelGrid as MG
import Partition as PA
......@@ -48,39 +49,24 @@ szhalo=1
nbcore=comm.Get_size()
rank=comm.Get_rank()
#
# Verify directories
#
wdir="Weights"
if rank == 0 :
if not os.path.exists(wdir) :
os.mkdir(wdir)
comm.Barrier()
#
# Region of grid to be treated
#
part = PA.partition(gg.ni, gg.nj, gg.land, comm, nbcore, szhalo, rank)
#
# Weight file to be created
#
wfile=os.path.basename(config.get("OverAll", "ModelGridFile").replace(".","_"))+"_i"+str(part.istart).zfill(4)+"di"+str(part.ni).zfill(4)+"j"+str(part.jstart).zfill(4)+"dj"+str(part.nj).zfill(4)
#
#
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]))
print("modelgrid.box_land",modelgrid.box_land)
#
hydrogrid=HG.HydroGrid(modelgrid.box_land)
#
# Computes weights of overlap of modelgrid and hydrogrid
#
sub_lon, sub_lat, sub_index, sub_area = RPP.compweights(wdir+"/"+wfile, modelgrid, hydrogrid)
w = RPP.compweights(wfile, part, modelgrid, hydrogrid)
#
#
#
nbpt = len(sub_index)
sub_pts = np.array(list(len(sub_index[i][0]) for i in range(len(sub_index))))
nbvmax = part.domainmax(max(sub_pts))
nbpt = len(w.index)
nbvmax = part.domainmax(max(w.hpts))
print("nbpt : ", nbpt)
print("nbvmax : ", nbvmax)
print("nbasmax : ", nbasmax)
......@@ -88,22 +74,19 @@ print("nbasmax : ", nbasmax)
# Extract hydo data from file
#
INFO("hydrodata")
hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, sub_index)
hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, w.index)
INFO("initiatmgrid")
IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
INFO("hoverlap")
hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, part, modelgrid, hydrodata)
hoverlap = IF.HydroOverlap(nbpt, nbvmax, w.hpts, w.index, w.area, w.lon, w.lat, part, modelgrid, hydrodata)
#
# Do some memory management and synchronize procs.
#
comm.Barrier()
del sub_index
del sub_area
del sub_lon
del sub_lat
del w
gc.collect()
if rank ==0 :
......
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