# # # import numpy as np import os, sys import pickle 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 # # Gert the information from the configuration file. # import configparser config=configparser.ConfigParser({'SaveWeights':'true', "DiagLon":"0.0, 0.0", "DiagLat":"0.0, 0.0"}) config.read("run.def") saveweights=config.get("OverAll", "SaveWeights") nbasmax=int(config.getfloat("OverAll", "nbasmax")) OutGraphFile=config.get("OverAll","GraphFile") lonint=np.array(config.get("OverAll", "DiagLon").split(","),dtype=float) latint=np.array(config.get("OverAll", "DiagLat").split(","),dtype=float) # import ModelGrid as MG import Partition as PA import RPPtools as RPP import HydroGrid as HG import diagplots as DP import Interface as IF # # 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=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])) hydrogrid=HG.HydroGrid(modelgrid.box_land) icell = 0 sub_index = [] sub_area = [] sub_lon = [] sub_lat = [] # # Only compute the weights if the file does not exist. # if not os.path.exists(wdir+"/"+wfile) : INFO("Computing weights") for cell, area in zip(modelgrid.polylist, modelgrid.area) : overlap=np.zeros((hydrogrid.jjm,hydrogrid.iim), dtype=float) for hydrocell,index in zip(hydrogrid.polylist, hydrogrid.index) : if cell.intersects_poly(hydrocell) : overlap[index[0],index[1]] = cell.overlap(hydrocell) INFO(str(icell)+" Sum of overlap "+str(np.sum(overlap))+ " Nb of Hydro grid overlaping :"+str(np.count_nonzero(overlap))) sub_index.append(np.where(overlap > 0.0)) sub_area.append(np.extract(overlap > 0.0, overlap*area)) sub_lon.append(np.extract(overlap > 0.0, hydrogrid.lon)) sub_lat.append(np.extract(overlap > 0.0, hydrogrid.lat)) icell = icell + 1 # # Save the weights to a dedicated file for future use # if ( saveweights.lower() == 'true' ) : with open(wdir+"/"+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 (wdir+"/"+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() # # # nbpt = len(sub_index) sub_pts = np.array(list(len(sub_index[i][0]) for i in range(len(sub_index)))) nbvmax = max(sub_pts) print("nbpt : ", nbpt) print("nbvmax : ", nbvmax) print("nbasmax : ", nbasmax) # # Extract hydo data from file # INFO("hydrodata") hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, sub_index) INFO("nitiatmgrid") IF.initatmgrid(rank, nbcore, nbpt, modelgrid) INFO("hoverlap") hoverlap = IF.HydroOverlap(nbpt, nbvmax, sub_pts, sub_index, sub_area, sub_lon, sub_lat, modelgrid, hydrodata) # # Do some memory management and synchronize procs. # comm.Barrier() del sub_index del sub_area del sub_lon del sub_lat gc.collect() if rank ==0 : if lonint[0] != lonint[1] and latint[0] != latint[1] : DP.MAPPLOT("MapHydroGrid", lonint, latint, hoverlap, hoverlap.hierarchy_bx, modelgrid.polyll, title="Distance to ocean") hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap) # # Plot the hydrological supermesh # if rank == 0 : if lonint[0] != lonint[1] and latint[0] != latint[1] : DP.SUPERMESHPLOT("MapSuperGrid_Beforelinkup", lonint, latint, hoverlap, hsuper, modelgrid.polyll) print("=================== LINKUP ====================") hsuper.linkup(hydrodata) if rank ==0 : if lonint[0] != lonint[1] and latint[0] != latint[1] : DP.SUPERMESHPLOT("MapSuperGrid_Afterlinkup", lonint, latint, hoverlap, hsuper, modelgrid.polyll) # # Do some memory management and synchronize procs. # comm.Barrier() del hoverlap gc.collect() print("=================== Compute fetch ====================") hsuper.fetch(part) hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part) hgraph = IF.HydroGraph(nbasmax, hsuper, part) hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part) IF.closeinterface(comm)