-
Notifications
You must be signed in to change notification settings - Fork 1
/
TestHapuaMod.py
131 lines (111 loc) · 6.03 KB
/
TestHapuaMod.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
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from hapuamod import loadmod, visualise, coast, riv, core, mor, out
#%% Test model load (particularly geometry processing)
ModelConfigFile = 'inputs\HurunuiModel.cnf'
(ModelName, FlowTs, WaveTs, SeaLevelTs, Origin, ShoreNormDir, ShoreX,
ShoreY, ShoreZ, RiverElev, OutletEndX, OutletEndWidth, OutletEndElev,
RiverWL, RiverVel, LagoonWL, LagoonVel, OutletDep, OutletVel, OutletEndDep,
OutletEndVel, Closed, TimePars, PhysicalPars, NumericalPars, OutputOpts) = \
loadmod.loadModel(ModelConfigFile)
plt.figure()
visualise.mapView(ShoreX, ShoreY, Origin, ShoreNormDir)
#%% Test longshore transport routine
EDir_h = WaveTs.EDir_h[0]
WavePower = WaveTs.WavePower[0]
WavePeriod = WaveTs.WavePeriod[0]
Wlen_h = WaveTs.Wlen_h[0]
plt.figure()
plt.subplot(2, 1, 1)
plt.plot(ShoreX, ShoreY[:,0])
LST = coast.longShoreTransport(ShoreY, NumericalPars['Dx'], WavePower,
WavePeriod, Wlen_h, EDir_h, PhysicalPars)
plt.subplot(2, 1, 2)
plt.plot((ShoreX[0:-1]+ShoreX[1:])/2, LST)
#%% Test river routines
# Join river and outlet through lagoon
(ChanDx, ChanElev, ChanWidth, LagArea, LagLen, ChanDep, ChanVel,
OnlineLagoon, OutletChanIx, ChanFlag, Closed) = \
riv.assembleChannel(ShoreX, ShoreY, ShoreZ,
OutletEndX, OutletEndWidth, OutletEndElev,
Closed, RiverElev,
np.zeros(RiverElev.size), np.zeros(RiverElev.size),
np.zeros(ShoreX.size), np.zeros(ShoreX.size),
np.zeros(ShoreX.size), np.zeros(ShoreX.size),
np.zeros(3), np.zeros(3), NumericalPars['Dx'],
PhysicalPars)
visualise.modelView(ShoreX, ShoreY, OutletEndX, OutletEndWidth, OutletChanIx, PhysicalPars['RiverWidth'], PhysicalPars['SpitWidth'])
# Steady state hydraulics
RivFlow = core.interpolate_at(FlowTs, pd.DatetimeIndex([TimePars['StartTime']])).values
SeaLevel = core.interpolate_at(SeaLevelTs, pd.DatetimeIndex([TimePars['StartTime']])).values
(ChanDep, ChanVel) = riv.solveSteady(ChanDx, ChanElev, ChanWidth,
PhysicalPars['Roughness'],
RivFlow[0], SeaLevel[0], NumericalPars)
# Store hydraulics and re-generate
(LagoonWL, LagoonVel, OutletDep, OutletVel, OutletEndDep, OutletEndVel) = \
riv.storeHydraulics(ChanDep, ChanVel, OnlineLagoon, OutletChanIx,
ChanFlag, ShoreZ[:,3], Closed)
(ChanDx, ChanElev, ChanWidth, LagArea, LagLen, ChanDep, ChanVel,
OnlineLagoon, OutletChanIx, ChanFlag, Closed) = \
riv.assembleChannel(ShoreX, ShoreY, ShoreZ,
OutletEndX, OutletEndWidth, OutletEndElev,
Closed, RiverElev,
ChanDep[ChanFlag==0], ChanVel[ChanFlag==0],
LagoonWL, LagoonVel, OutletDep, OutletVel,
OutletEndDep, OutletEndVel,
NumericalPars['Dx'], PhysicalPars)
# Bedload
Bedload = riv.calcBedload(ChanElev, ChanWidth, ChanDep, ChanVel, ChanDx,
PhysicalPars, NumericalPars['Psi'])
visualise.longSection(ChanDx, ChanElev, ChanWidth, ChanDep, ChanVel, Bedload)
ChanDist = np.insert(np.cumsum(ChanDx),0,0)
plt.figure()
plt.plot(ChanDist, ChanDep+ChanElev, 'b-')
# Unsteady hydraulics
riv.solveFullPreissmann(ChanElev, ChanWidth, LagArea, LagLen,
Closed, ChanDep, ChanVel, ChanDx,
TimePars['HydDt'], RivFlow, SeaLevel,
NumericalPars, PhysicalPars)
plt.plot(ChanDist, ChanDep+ChanElev, 'r:')
#%% Test transect plot
BeachSlope = PhysicalPars['BeachSlope']
BackshoreElev = PhysicalPars['BackshoreElev']
ClosureDepth = PhysicalPars['ClosureDepth']
BeachTopElev = PhysicalPars['BeachTopElev']
OutletWL = OutletDep + ShoreZ[:,1]
TransectFig = visualise.newTransectFig(ShoreX, ShoreY, ShoreZ, LagoonWL, OutletWL,
SeaLevel, BeachSlope, BackshoreElev,
ClosureDepth, BeachTopElev, -300)
#%% Runup/overtopping
WavePeriod = WaveTs.WavePeriod[0]
Hs_offshore = WaveTs.Hsig_Offshore[0]
Runup = coast.runup(WavePeriod, Hs_offshore, PhysicalPars['BeachSlope'])
(CST_tot, OverwashProp) = coast.overtopping(Runup, SeaLevel, ShoreY,
ShoreZ, PhysicalPars)
#%% Morphology updating
OldShoreY = ShoreY.copy()
MorDt = TimePars['MorDtMin']
(MorDt, Breach) = mor.updateMorphology(ShoreX, ShoreY, ShoreZ,
OutletEndX, OutletEndWidth, OutletEndElev,
RiverElev, OnlineLagoon,
OutletChanIx, LagoonWL, OutletDep,
ChanWidth, ChanDep, ChanDx, ChanFlag,
Closed, LST, Bedload, CST_tot, OverwashProp,
MorDt, PhysicalPars, TimePars, NumericalPars)
plt.plot(ShoreX, (ShoreY[:,0]-OldShoreY[:,0]))
#%% Create output netcdf file and write initial condition
out.newOutFile('test.nc', ModelName, TimePars['StartTime'],
ShoreX, NumericalPars['Dx'], RiverElev,
Origin, ShoreNormDir, PhysicalPars, False)
out.writeCurrent('test.nc', TimePars['StartTime'], SeaLevel[-1], RivFlow[-1],
Hs_offshore, EDir_h,
ShoreY, ShoreZ, LagoonWL, LagoonVel, np.zeros(ShoreX.size),
OutletDep, OutletVel, np.zeros(ShoreX.size),
np.zeros(ShoreX.size-1), np.zeros(ShoreX.size), np.zeros(ShoreX.size),
RiverElev, ChanDep[ChanFlag==0], ChanVel[ChanFlag==0], np.zeros(RiverElev.size),
OutletEndX, OutletEndElev, OutletEndWidth,
OutletEndDep, OutletEndVel, np.zeros(2), Closed)
#%% Test core timestepping
ModelConfigFile = 'inputs\HurunuiModel.cnf'
OutputTs = core.main(ModelConfigFile)