Commit 67878875 authored by Anthony's avatar Anthony
Browse files

Correction / Improvement of Truncate

parent c2b62560
......@@ -17,6 +17,10 @@ localdir=os.path.dirname(getframeinfo(currentframe()).filename)
sys.path.append(localdir+'/F90subroutines')
import routing_interface
import getargs
config = getargs.SetupConfig()
frac_lim=config.getfloat("OverAll", "frac_lim", fallback=0.3)
limit_step5=config.getfloat("OverAll","limit_step5",fallback=5)
log_master, log_world = getargs.getLogger(__name__)
INFO, DEBUG, ERROR = log_master.info, log_master.debug, log_world.error
......@@ -41,7 +45,12 @@ class trunc:
"""
#
# Initialize some variables
#
#
# Fraction of grid limit to operate the different step
self.frac_lim = frac_lim
# Limit of % of grid for HTUs to kill during step 5
self.limit_step5 = limit_step5
self.nbpt = int(hydrosuper.nbpt)
self.nbasmax = hydrosuper.nbasmax
self.gridarea = modelgrid.area
......@@ -59,6 +68,8 @@ class trunc:
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 1 iter. {0}, maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num, max_bas))
INFO("Step 1 iter. {0}, maxops: {1}->max_bas : {2}".format(i, self.max_ops_num, max_bas))
i+=1
hydrosuper.fetch(part)
......@@ -75,6 +86,8 @@ class trunc:
hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 2 iter. {0} (core) maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num_core, max_bas))
INFO("Step 2 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_core, max_bas))
comm.Barrier()
hydrosuper.fetch(part)
......@@ -87,6 +100,8 @@ class trunc:
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 2 iter. {0} (halo) maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num_halo, max_bas))
INFO("Step 2 iter. {0} (halo) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_halo, max_bas))
comm.Barrier()
hydrosuper.fetch(part)
......@@ -105,6 +120,8 @@ class trunc:
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 3 iter. {0} (core) maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num_core, max_bas))
INFO("Step 3 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_core, max_bas))
hydrosuper.fetch(part)
comm.Barrier()
......@@ -115,6 +132,8 @@ class trunc:
hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 3 iter. {0} (halo) maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num_halo, max_bas))
INFO("Step 3 iter. {0} (halo) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num_halo, max_bas))
hydrosuper.fetch(part)
comm.Barrier()
......@@ -129,17 +148,24 @@ class trunc:
hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
INFO("Step 4 iter. {0} (core) maxops: {1}->max_bas : {2}".format(i, self.max_ops_num, max_bas))
with open("check.out","a+") as foo:
foo.write("Step 4 iter. {0} maxops: {1}->max_bas : {2}\n".format(i, self.max_ops_num, max_bas))
INFO("Step 4 iter. {0} maxops: {1}->max_bas : {2}".format(i, self.max_ops_num, max_bas))
hydrosuper.fetch(part)
comm.Barrier()
i += 1
# Step 5 : Brutal method
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step5\n")
self.step5(hydrosuper, part)
hydrosuper.killbas(self.tokill, self.totakeover, self.numops)
max_bas = part.domainmax(np.max(self.basin_count))
if part.rank == 0:
with open("check.out","a+") as foo:
foo.write("Step 5, maxops: {0}->max_bas : {1}".format(self.numopmax, max_bas))
INFO("Step 5, maxops: {0}->max_bas : {1}".format(self.numopmax, max_bas))
hydrosuper.fetch(part)
......@@ -187,7 +213,7 @@ class trunc:
for k in range(len(self.landcorelist)):
changes = self.get_changes_step2(k, neigh_domain)
if len(changes)>0:
if len(changes)>1:
for ch in changes:
orig_ops = self.numops[self.landcorelist[k]]
if orig_ops < self.numop:
......@@ -217,7 +243,7 @@ class trunc:
for k in range(len(self.landcorelist)):
changes = self.get_changes_step3(k, neigh_domain)
if len(changes)>0:
if len(changes)>1:
for ch in changes:
orig_ops = self.numops[self.landcorelist[k]]
if orig_ops < self.numop:
......@@ -247,7 +273,7 @@ class trunc:
for k in range(len(self.landcorelist)):
tokk, totakk = self.get_changes_step4(k)
if len(tokk)>0:
if (len(tokk)>0 and len(totakk)>0):
ops = min(self.numop,len(tokk))
self.tokill[self.landcorelist[k],:ops] = tokk[:ops]
......@@ -268,7 +294,6 @@ class trunc:
# Step 5 : htus inside point
"""
# Limit of the area of basin to delete by step 5 (in %)
limit_step5 = 5
self.load_data(hydrosuper)
# Number of basin to delete in each point of the core
......@@ -280,7 +305,7 @@ class trunc:
if self.numopmax>self.numop:
print("Error - cannot finish to resolve truncate by step 5")
print("Try with an higher nbasmax")
print("Try with an higher numops, here numopsmax = {0}".format(self.numopmax))
sys.exit()
......@@ -290,13 +315,15 @@ class trunc:
for k, k_land in enumerate(self.landcorelist):
if num2del[k] > 0:
B = np.argsort(self.fetch[k])[:num2del[k]]
tokk = [k+1 for k in B]
totakk = [np.argsort(self.fetch[k])[-1] + 1] * num2del[k]
area_basin_to_truncate = np.sum([self.basin_area[k][l] for l in B])
if (area_basin_to_truncate/np.sum(self.basin_area[k]))*100 < limit_step5:
B = np.argsort(self.basin_area[k])[:num2del[k]+1]
f = [self.fetch[k][m] for m in B]
f = np.argsort(f)
tokk = [B[m]+1 for m in f[:-1]]
totakk = [B[f[-1]] + 1] * num2del[k]
area_basin_to_truncate = np.sum([self.basin_area[k][B[m]] for m in f[:-1]])
AREA = (area_basin_to_truncate/np.sum(self.basin_area[k]))*100
if AREA < self.limit_step5:
ops = len(tokk)
self.tokill[k_land,:ops] = tokk
self.totakeover[k_land,:ops] = totakk
......@@ -306,7 +333,7 @@ class trunc:
else:
print("Error - cannot finish to resolve truncate by step 5")
print("Try with an higher nbasmax")
print("Try with an higher nbasmax, total area : {0}".format(AREA))
sys.exit()
part.landsendtohalo(self.numops, order='F')
......@@ -347,7 +374,7 @@ class trunc:
"""
# List of neighbour which are outflow
outgrid = self.outflow_grid[npt]
outgrid = ma.masked_where(outgrid == self.landcorelist[npt], outgrid)
outgrid = ma.masked_where(outgrid == self.landcorelist_f[npt], outgrid)
outgrid = ma.masked_where(outgrid <= 0, outgrid)
outgrid = ma.unique(outgrid)
......@@ -370,7 +397,7 @@ class trunc:
if len(A)>1:
B = self.fetch[npt][A]
C = self.basin_area[npt][A]
changes.append([A[k]+1 for k in np.argsort(B) if C[k]/area < 0.1])
changes.append([A[k]+1 for k in np.argsort(B) if C[k]/area < self.frac_lim])
elif neigh_domain == "halo":
for k in halo:
......@@ -378,7 +405,7 @@ class trunc:
if len(A)>1:
B = self.fetch[npt][A]
C = self.basin_area[npt][A]
changes.append([A[k]+1 for k in np.argsort(B) if C[k]/area < 0.1])
changes.append([A[k]+1 for k in np.argsort(B) if C[k]/area < self.frac_lim])
return changes
......@@ -389,42 +416,49 @@ class trunc:
"""
# List of basin ID included in the grid point
basid = self.basin_id[npt]
basid = ma.unique(basid)
basid_un = ma.unique(basid)
core = []; halo = []
# For each basin id, list of point with outflow in the core / halo
for bid in basid:
for bid in basid_un:
htus = np.where(basid == bid)[0]
outgrid = [self.outflow_grid[npt][htu] for htu in htus]
#outgrid = [self.outflow_grid[npt][htu] for htu in htus]
c = []; h = []
for g,htu in zip(outgrid,htus):
if g in self.landhalolist_f:
h.append(htu)
elif g in self.landcorelist_f:
c.append(htu)
if len(h)>1:
halo.append(h)
if len(c)>1:
core.append(c)
#for g,htu in zip(outgrid,htus):
# if g in self.landhalolist:
# h.append(htu)
# elif g in self.landcorelist:
# c.append(htu)
#if len(h)>1:
# halo.append(h)
#if len(c)>1:
# core.append(c)
if len(htus)>1: core.append(htus)
# Establishing the list of changes
changes = []
area = np.sum(self.basin_area[npt])
for c in core:
FETCH = [self.fetch[npt][m] for m in c]
AREA = [self.basin_area[npt][m] for m in c]
changes.append([c[k]+1 for k in np.argsort(FETCH) if AREA[k]/area < self.frac_lim])
"""
if neigh_domain == "core":
for c in core:
B = self.fetch[npt][c]
C = self.basin_area[npt][A]
changes.append([c[k]+1 for k in np.argsort(B) if C[k]/area < 0.1])
B = [self.fetch[npt][m] for m in c]
C = [self.basin_area[npt][m] for m in c]
changes.append([c[k]+1 for k in np.argsort(B) if C[k]/area < 0.2])
elif neigh_domain == "halo":
for h in halo:
B = self.fetch[npt][h]
C = self.basin_area[npt][h]
changes.append([h[k]+1 for k in np.argsort(B) if C[k]/area < 0.1])
changes.append([h[k]+1 for k in np.argsort(B) if C[k]/area < 0.2])
"""
return changes
def get_changes_step4(self, npt):
......@@ -433,13 +467,13 @@ class trunc:
# List of HTUs flowing in the same grid
outgrid = self.outflow_grid[npt]
htus = np.where(outgrid == self.landcorelist[npt])[0]
htus = np.where(outgrid == self.landcorelist_f[npt])[0]
# Establishing the lists of changes
if len(htus)>0:
if len(htus)>1:
B = self.fetch[npt][htus]
# np.flip -> to reduce the HTUs nummber over the main river
tok = [htus[k]+1 for k in np.flip(np.argsort(B))]
tok = [htus[k]+1 for k in np.argsort(B)]
totak = [self.outflow_basin[npt][htus[k]] for k in np.argsort(B)]
return tok, totak
......@@ -467,7 +501,7 @@ class trunc:
HTUs_to_truncate : List of HTUs to truncate, ordered by fetch
"""
orig_ops = self.numops[self.landcorelist[k]]
poss_ops = len(HTUs_to_truncate)-1
poss_ops = max(len(HTUs_to_truncate)-1,0)
ops = min(poss_ops,self.numop-orig_ops)
self.tokill[self.landcorelist[k],orig_ops:orig_ops+ops] = HTUs_to_truncate[:ops]
......
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