forked from respec/HSPsquared
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathADCALC.py
260 lines (220 loc) · 9.12 KB
/
ADCALC.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
''' Copyright (c) 2020 by RESPEC, INC.
Authors: Robert Heaphy, Ph.D. and Paul Duda
License: LGPL2
'''
from numpy import zeros, array
from numba import njit
from HSP2.utilities import make_numba_dict
# The clean way to get calculated data from adcalc() into advert() is to use a closure.
# This is not currently supported by Numba.
'''
split out the loop section of adcalc to allow numba caching
let omat be passed in advect to avoid recreating it each time
ADFLAG == 2 loop simplification as elif
'''
ERRMSG = []
def adcalc(store, siminfo, uci, ts):
'''Prepare to simulate advection of fully entrained constituents'''
errorsV = zeros(len(ERRMSG), dtype=int)
simlen = siminfo['steps']
ui = make_numba_dict(uci)
nexits = int(ui['NEXITS']) # table type GEN-INFO
ui['simlen'] = siminfo['steps']
ui['delts'] = siminfo['delt'] * 60.0 # delts is the simulation interval in seconds
ui['uunits'] = siminfo['units']
# calculated timeseries for advect()
if 'SROVOL' not in ts:
ts['SROVOL'] = zeros(simlen)
if 'EROVOL' not in ts:
ts['EROVOL'] = zeros(simlen)
############################################################################
errors = _adcalc_(ui, ts) # run ADCALC simulation code
############################################################################
if 'VOL' in ui:
vol = ui['VOL']
else:
vol = 0.0
VOL = ts['VOL']
SROVOL = ts['SROVOL']
EROVOL = ts['EROVOL']
SOVOL = zeros((simlen, nexits))
EOVOL = zeros((simlen, nexits))
for i in range(nexits):
SOVOL[:, i] = ts['SOVOL' + str(i + 1)]
EOVOL[:, i] = ts['EOVOL' + str(i + 1)]
uci['adcalcData'] = (nexits, vol, VOL, SROVOL, EROVOL, SOVOL, EOVOL)
return errorsV, ERRMSG
@njit(cache=True)
def _adcalc_(ui, ts):
''' Internal adcalc() loop for Numba'''
ADFG = ui['ADFG'] # table type ACTIVITY
simlen = int(ui['simlen'])
nexits = int(ui['NEXITS'])
delts = ui['delts']
uunits = ui['uunits']
# units conversion constants, 1 ACRE is 43560 sq ft. assumes input in acre-ft
VFACT = 43560.0
AFACT = 43560.0
if uunits == 2:
# si units conversion constants, 1 hectare is 10000 sq m, assumes area input in hectares, vol in Mm3
VFACT = 1.0e6
AFACT = 10000.0
# table ADCALC-DATA
if 'CRRAT' in ui:
crrat = ui['CRRAT']
else:
crrat = 1.5
if 'VOL' in ui:
vol = ui['VOL'] * VFACT
else:
vol = 0.0
ks = 0.0
if ADFG == 2:
ks = ui['KS'] # need to get from HYDR section
VOL = ts['VOL']
SROVOL = ts['SROVOL']
EROVOL = ts['EROVOL']
if 'SOVOL' in ts:
SOVOL = ts['SOVOL']
else:
SOVOL = zeros((simlen, nexits))
for index in range(nexits):
ts['SOVOL' + str(index+1)] = SOVOL[:,index]
if 'EOVOL' in ts:
EOVOL = ts['EOVOL']
else:
EOVOL = zeros((simlen, nexits))
for index in range(nexits):
ts['EOVOL' + str(index + 1)] = EOVOL[:, index]
ROS = ui['ROS']
# external time series
O = zeros((simlen, nexits))
if nexits > 1:
for i in range(nexits):
O[:, i] = ts['O' + str(i + 1)]
else:
O[:, 0] = ts['RO']
for loop in range(simlen):
vols = VOL[loop-1] * VFACT if loop > 0 else vol
o = O[loop]
os = O[loop-1] if loop > 0 else array([ROS])
ro = 0.0
ros= 0.0
for index in range(nexits):
ro += o[index]
ros += os[index]
# weighting factors to calculate mean outflow rate over ivl: constituents outflow rate at start (js) and end (cojs) of interval;
if ADFG >= 1: # first use standard method of computing advective weighting
if ros > 0.0: # calculate ratio of volume to outflow volume;
rat = vols / (ros * delts)
if rat < crrat:
# some outflow volume entered control volume as inflow during same interval; hence,
# concentration of inflowing material will affect outflow concentration => js will be < 1.0
js = rat / crrat
else:
# all water in outflow volume was contained in control volume at beginning of ivl;
# mean rate of outflow over ivl will be wholly dependent upon rate of outflow of constituents at start of ivl
js = 1.0
else: # reach/res has no outflow at start of ivl
js = 0.0
if ADFG == 2 and js > ks:
# use same weighting as used in flow routing i.e., JS = KS; this was added 6/2010 by brb & pbd
# based on suggestion of Sen Bai of Tetra Tech - ASCE J. Hydrol. Engr., Vol.15, No.3, March, 2010
js = ks
cojs = 1.0 - js # cojs is the complement of js
# calculate weighted volumes of outflow at start of ivl (srovol) and end of ivl (erovol)
SROVOL[loop] = js * ros * delts
EROVOL[loop] = cojs * ro * delts
# if nexits > 1: # determine weighted volume of outflow at start and end of ivl per exit
for index in range(nexits):
ts['SOVOL' + str(index + 1)][loop] = js * os[index] * delts
ts['EOVOL' + str(index + 1)][loop] = cojs * o[index] * delts
return
@njit(cache=True)
def advect(imat, conc, nexits, vols, vol, srovol, erovol, sovol, eovol):
''' Simulate advection of constituent totally entrained in water.
Originally designed to be called as: advect(loop, imat, conc, omat, *ui['adcalcData'])
but unit conversions in the calling routine make this impractical'''
# vols = VOL[loop-1] if loop > 0 else vol
# vol = VOL[loop]
# srovol = SROVOL[loop]
# erovol = EROVOL[loop]
# sovol = SOVOL[loop,:]
# eovol = EOVOL[loop,:]
omat = zeros((nexits))
if vol > 0.0: # reach/res contains water
concs = conc
conc = (imat + concs * (vols - srovol)) / (vol + erovol) # material entering during interval, weighted volume of outflow based on conditions at start of ivl (srovol), and weighted volume of outflow based on conditions at end of ivl (erovol)
romat = srovol * concs + erovol * conc # total material leaving reach/res in ivl
if nexits > 1: # material leaving through each exit gate
omat = sovol * concs + eovol * conc # qty.vol/l.ivl, array calculation
else: # reach/res has gone dry during the interval
romat = imat + (conc * vols) # total material leaving during interval = inflow + initial material
if nexits > 1 and srovol > 0: # calculate material leaving through each exit gate
omat = (sovol / srovol) * romat # array calculation
conc = -1.0e30
return conc, romat, omat
@njit(cache=True)
def oxrea(LKFG,wind,cforea,avvele,avdepe,tcginv,reamfg,reak,reakt,expred,exprev,len, delth,tw,delts,delt60,uunits):
''' Calculate oxygen reaeration coefficient'''
# DELTS - ???
# DELT60 - simulation time interval in hours
# UUNITS - system of units 1-english, 2-metric
if LKFG == 1: # this reach/res is a lake or reservoir
''' empirically reaeration coefficient based on windspeed, surface area, and volume;
windsp is windspeed in m/sec; wind is wind movement in m/ivl'''
windsp = wind / delts
windf = windsp * (-0.46 + 0.136 * windsp) if windsp > 6.0 else 2.0
korea = (0.032808 * windf * cforea / avdepe) * delt60
else:
# calculate reaeration coefficient for free-flowing reach
if reamfg == 1:
# calculate reaeration coefficient based on energy dissipation principles (tsivoglou method)
# convert length and drop in energy line along length of rchres to english units, if necessary
if uunits == 2:
lene = len * 3.28
delthe = delth * 3.28
else:
lene = len
delthe = delth
if abs(avvele) > 0.0:
flotim = lene / avvele
korea = reakt * (delthe / flotim) * (tcginv**(tw - 20.)) * delts
else:
korea = 0.0
else:
if reamfg == 2:
# calculate reaeration coefficient as a power function of average hydraulic
# depth and velocity; determine exponents to depth and velocity terms and assign value to reak
if avdepe <= 2.0: # use owen's formulation for reaeration
reak = 0.906
exprev = 0.67
expred = -1.85
else:
# calculate transition depth; transition depth determines which method
# of calculation is used given the current velocity
trandp = 0.0 if avvele < 1.7 else 0.4263 * (avvele**2.9135)
if avdepe - trandp <= 0.0: # use churchill's formulation for reaeration
reak = 0.484
exprev = 0.969
expred = -1.673
else: # use o'connor-dobbins formulation for reaeration
reak = 0.538
exprev = 0.5
expred = -1.5
korea = reak * avvele**exprev * avdepe**expred * tcginv**(tw - 20.0) * delt60 if tw < 66 else 0.999
if korea > 1.0:
korea = 0.999
return korea
@njit(cache=True)
def sink (vol, avdepe, kset, conc):
''' calculate quantity of material settling out of the control volume; determine the change in concentration as a result of sinking'''
if kset > 0.0 and avdepe > 0.17:
# calculate concentration change due to outgoing material; snkout is expressed in mass/liter/ivl; kset is expressed as ft/ivl and avdepe as feet
snkout = conc * (kset / avdepe) if kset < avdepe else conc # calculate portion of material which settles out of the control volume during time step; snkout is expressed as mass/liter.ivl; conc is the concentration of material in the control volume
conc -= snkout # calculate remaining concentration of material in the control volume
snkmat = snkout * vol # find quantity of material that sinks out; units are mass.ft3/l.ivl in english system, and mass.m3/l.ivl in metric system
else:
snkout = 0.0
snkmat = 0.0
return conc, snkmat