Skip to content

Commit

Permalink
bug fix numerical imprecisions warm start
Browse files Browse the repository at this point in the history
  • Loading branch information
@StefaniaGrimaldi
StefaniaGrimaldi committed Dec 7, 2022
1 parent 579cdc0 commit 6043971
Show file tree
Hide file tree
Showing 3 changed files with 11 additions and 24 deletions.
15 changes: 1 addition & 14 deletions src/lisflood/Lisflood_dynamic.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -196,21 +196,8 @@ def splitlanduse(array1, array2=None, array3=None):
# sum of both lines
# CrossSection2Area = pcraster.max(scalar(0.0), (self.Chan2M3Kin - self.Chan2M3Start) / self.ChanLength)

# temporary subtract routing ChanQ here, to be added after checks for small values
self.sumDisDay -= self.ChanQ

self.TotalCrossSectionArea = self.ChanM3 * self.InvChanLength

if option['SplitRouting']:
self.TotalCrossSectionAreaSAFE = self.TotalCrossSectionArea
self.CrossSection2AreaSAFE = self.CrossSection2Area
self.TotalCrossSectionArea = np.where(np.abs(self.TotalCrossSectionAreaSAFE - self.CrossSection2AreaSAFE)<0.0000001,0.0,self.TotalCrossSectionArea)
self.CrossSection2Area = np.where(np.abs(self.TotalCrossSectionAreaSAFE - self.CrossSection2AreaSAFE)<0.0000001,0.0,self.CrossSection2Area)
self.ChanQ = np.where(np.abs(self.TotalCrossSectionAreaSAFE - self.CrossSection2AreaSAFE)<0.0000001,0.0,self.ChanQ)
self.ChanM3 = np.where(np.abs(self.TotalCrossSectionAreaSAFE - self.CrossSection2AreaSAFE)<0.0000001,0.0,self.ChanM3)

# add ChanQ again after small values check
self.sumDisDay += self.ChanQ
self.sumDis += self.sumDisDay
self.ChanQAvg = self.sumDisDay/self.NoRoutSteps

Expand Down Expand Up @@ -271,4 +258,4 @@ def splitlanduse(array1, array2=None, array3=None):
self.stateVar_module.dynamic()
if option['wateruse'] and option['indicator'] and self.monthend:
self.indicatorcalc_module.dynamic_setzero()
# setting monthly and yearly dindicator to zero at the end of the month (year)
# setting monthly and yearly dindicator to zero at the end of the month (year)
18 changes: 9 additions & 9 deletions src/lisflood/global_modules/default_options.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -127,12 +127,12 @@
'wateruseRegion',
'indicator'],
monthly=True, yearly=False),
'ChSideEnd': ReportedMap(name='ChSideEnd', output_var='Sideflow1Chan', unit='m',
'ChSideEnd': ReportedMap(name='ChSideEnd', output_var='Sideflow1Chan', unit='m3/s',
end=['repEndMaps'], steps=[], all=[],
restrictoption=['SplitRouting'],
monthly=False, yearly=False),
'ChSideState': ReportedMap(name='ChSideState', output_var='Sideflow1Chan',
unit='m', end=[], steps=['repStateMaps'],
unit='m3/s', end=[], steps=['repStateMaps'],
all=[], restrictoption=['nonInit', 'SplitRouting'],
monthly=False, yearly=False),
'ChanCrossSectionEnd': ReportedMap(name='ChanCrossSectionEnd',
Expand Down Expand Up @@ -208,12 +208,12 @@
'indicator'],
monthly=True, yearly=False),
'CrossSection2End': ReportedMap(name='CrossSection2End',
output_var='CrossSection2Area', unit='m',
output_var='CrossSection2Area', unit='m2',
end=['repEndMaps'], steps=[], all=[],
restrictoption=['SplitRouting'],
monthly=False, yearly=False),
'CrossSection2State': ReportedMap(name='CrossSection2State',
output_var='CrossSection2Area', unit='m',
output_var='CrossSection2Area', unit='m2',
end=[], steps=['repStateMaps'], all=[],
restrictoption=['nonInit', 'SplitRouting'],
monthly=False, yearly=False),
Expand All @@ -240,7 +240,7 @@
yearly=False),
'CumInterceptionForestMaps': ReportedMap(name='CumInterceptionForestMaps',
output_var='CumInterception[1]',
unit='m', end=[], steps=[],
unit='mm', end=[], steps=[],
all=['repCumInterCeptionMaps'],
restrictoption=['nonInit'],
monthly=False, yearly=False),
Expand All @@ -261,7 +261,7 @@
unit='mm', end=[], steps=['repStateMaps'], all=[],
restrictoption=['nonInit'], monthly=False, yearly=False),
'CumInterceptionMaps': ReportedMap(name='CumInterceptionMaps',
output_var='CumInterception[0]', unit='m',
output_var='CumInterception[0]', unit='mm',
end=[], steps=[],
all=['repCumInterCeptionMaps'],
restrictoption=['nonInit'], monthly=False,
Expand Down Expand Up @@ -791,7 +791,7 @@
restrictoption=[],
monthly=False, yearly=False),
'ReservoirFillState': ReportedMap(name='ReservoirFillState',
output_var='ReservoirFill', unit='m',
output_var='ReservoirFill', unit='-',
end=[], steps=['repStateMaps'], all=[],
restrictoption=['simulateReservoirs'],
monthly=False, yearly=False),
Expand Down Expand Up @@ -1087,7 +1087,7 @@
end=['repEndMaps'], steps=[], all=[],
restrictoption=[], monthly=False,
yearly=False),
'UZForestMaps': ReportedMap(name='UZForestMaps', output_var='UZ[1]', unit='m',
'UZForestMaps': ReportedMap(name='UZForestMaps', output_var='UZ[1]', unit='mm',
end=[], steps=[], all=['repUZMaps'],
restrictoption=['nonInit'], monthly=False,
yearly=False),
Expand All @@ -1103,7 +1103,7 @@
unit='mm', end=[], steps=['repStateMaps'],
all=[], restrictoption=['nonInit'],
monthly=False, yearly=False),
'UZMaps': ReportedMap(name='UZMaps', output_var='UZ[0]', unit='m', end=[],
'UZMaps': ReportedMap(name='UZMaps', output_var='UZ[0]', unit='mm', end=[],
steps=[], all=['repUZMaps'], restrictoption=['nonInit'],
monthly=False, yearly=False),
'UZOutflowForestMaps': ReportedMap(name='UZOutflowForestMaps',
Expand Down
2 changes: 1 addition & 1 deletion src/lisflood/hydrological_modules/routing.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -390,7 +390,7 @@ def initialSecond(self):
self.var.Chan2M3Kin = self.var.CrossSection2Area * self.var.ChanLength + self.var.Chan2M3Start
self.var.ChanM3Kin = self.var.ChanM3 - self.var.Chan2M3Kin + self.var.Chan2M3Start

self.var.ChanM3Kin = np.where((self.var.ChanM3Kin < 0.0) & (self.var.ChanM3Kin > -0.0000001),0.0,self.var.ChanM3Kin)
self.var.ChanM3Kin = np.where((self.var.ChanM3Kin < 0.0) & (self.var.ChanM3Kin > -0.0000001),0.0,self.var.ChanM3Kin) # this line prevents the warmstart from failing in case of small numerical imprecisions when writing and reading the maps

self.var.Chan2QKin = (self.var.Chan2M3Kin * self.var.InvChanLength * self.var.InvChannelAlpha2) ** (self.var.InvBeta)
self.var.ChanQKin = (self.var.ChanM3Kin * self.var.InvChanLength * self.var.InvChannelAlpha) ** (self.var.InvBeta)
Expand Down

0 comments on commit 6043971

Please sign in to comment.