RoutingPreProc.py 3.58 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12
#
#
#
import numpy as np
import os, sys
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
13
import time
14 15 16
#
# Gert the information from the configuration file.
#
17
ConfigFile="run.def"
18 19
import getargs
config = getargs.SetupConfig()
POLCHER Jan's avatar
POLCHER Jan committed
20
nbasmax=config.getint("OverAll", "nbasmax")
21
numop=config.getint("OverAll", "numop", fallback=100)
22
OutGraphFile=config.get("OverAll","GraphFile")
23
DumpHydroSuper=config.getboolean("OverAll","DumpHydroSuper",fallback=False)
24
wfile=config.get("OverAll","WeightFile",fallback=None)
25 26 27 28 29 30
#
import ModelGrid as MG
import Partition as PA
import RPPtools as RPP
import HydroGrid as HG
import Interface as IF
31
from Truncate import trunc as TR
32 33 34 35 36 37 38 39 40 41 42 43 44
#
# Logging in MPI : https://groups.google.com/forum/#!topic/mpi4py/SaNzc8bdj6U
#
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
45
nbcore=comm.Get_size()
46 47 48 49
rank=comm.Get_rank()
#
# Region of grid to be treated
#
50
part = PA.partition(gg.ni, gg.nj, gg.land, comm, nbcore, szhalo, rank)
51 52
#
modelgrid=MG.ModelGrid(part.ihstart+gg.igstart,part.nih,part.jhstart+gg.jgstart,part.njh)
53
#
54 55
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]))
56
#
57
hydrogrid=HG.HydroGrid(modelgrid.box_land, wfile)
58
#
59
# Computes weights of overlap of modelgrid and hydrogrid
60
w = RPP.compweights(wfile, part, modelgrid, hydrogrid)
61 62 63
#
#
#
64
nbpt = len(w.index)
65 66 67 68
nbvmax = part.domainmax(w.maxhpts)
#
nbasmax = min(nbasmax,nbvmax)
#
69 70 71
INFO("nbpt : {0}".format(nbpt))
INFO("nbvmax : {0}".format(nbvmax))
INFO("nbasmax : {0}".format(nbasmax))
72 73 74
#
# Extract hydo data from file
#
75
INFO("=================== HYDRODATA ====================")
76
hydrodata = HG.HydroData(hydrogrid.ncfile, hydrogrid.box, w.index, part, hydrogrid)
77
hydrotcst = HG.HydroParameter(hydrogrid)
78

79
INFO("=================== INITATMGRID ====================")
80
IF.initatmgrid(rank, nbcore, nbpt, modelgrid)
81

82
INFO("=================== HOVERLAP ====================")
83
hoverlap = IF.HydroOverlap(nbpt, nbvmax, w.hpts, w.index, w.area, w.lon, w.lat, part, modelgrid, hydrodata)
84 85 86
#
# Do some memory management and synchronize procs.
#
87

88
del w
89
gc.collect()
90
comm.Barrier()
91

92
hsuper = IF.HydroSuper(nbvmax, hydrodata, hoverlap, nbasmax, part)
93
#
94
INFO("=================== LINKUP ====================")
95
hsuper.linkup(hydrodata)
96

97 98 99
#
# Do some memory management and synchronize procs.
#
100

101
del hoverlap
102
gc.collect()
103
comm.Barrier()
104 105 106
#
hsuper.correct_outflows(part)
#
107
INFO("=================== Compute fetch ====================")
108
t = time.time()
POLCHER Jan's avatar
POLCHER Jan committed
109
hsuper.fetch(part)
110 111
comm.Barrier()
t1 = time.time()
112
print("Time for fetch: {:0.2f} s.".format(t1-t))
113 114
comm.Barrier()

115
if DumpHydroSuper :
116
    INFO("Dumping HydroSuper")
117
    hsuper.dumpnetcdf(OutGraphFile.replace(".nc","_HydroSuper.nc"), gg, modelgrid, part)
118

Anthony's avatar
Anthony committed
119
INFO("=================== Truncate ====================")
120
print("Rank : {0} - Basin_count Before Truncate : ".format(part.rank), hsuper.basin_count)
121
hs = TR(hsuper, part, comm, modelgrid, numop = numop)
122 123
print("Rank : {0} - Basin_count After Truncate : ".format(part.rank), hsuper.basin_count)

Anthony's avatar
Anthony committed
124
INFO("=================== HydroGraph ====================")
125
hgraph = IF.HydroGraph(nbasmax, hsuper, part, modelgrid)
126
hgraph.dumpnetcdf(OutGraphFile, gg, modelgrid, part, hydrotcst)
127 128

IF.closeinterface(comm)