MultirateGBR.py 39.2 KB
Newer Older
thomas's avatar
thomas committed
1
#"""
thomas's avatar
thomas committed
2
3
#   Simulations on the GBR with multirate explicit Runge-Kutta time-stepping.
#   All the parameters can be found and modified in the file "InpuGBR.py"
thomas's avatar
thomas committed
4
5
#"""

thomas's avatar
thomas committed
6
7
8
9
10
11
12
#   A) HOW TO LAUNCH A SIMULATION
#       Simulation on 1 processor  : rundgpy MultirateGBR.py
#       Simulation on n processors : mpirun -np n rundgpy MultirateGBR.py   ( !!! You need MPI !!! )
#
#
#   B) TO START SIMULATION FROM A PRE-EXISTING SOLUTION:
#     1. Set the following in Input.py:
thomas's avatar
thomas committed
13
14
15
16
17
#         -> set flag startFromExistingSolution = 1
#         -> set solutionFileName
#         -> set stepNbOfInitialSolution
#         -> set old_dt
#         -> then, Ti = old Ti + stepNbOfInitialSolution*old_dt
thomas's avatar
thomas committed
18
#     2. Set the following in MultirateGBR.py:
thomas's avatar
thomas committed
19
20
#         -> Nothing
#
thomas's avatar
thomas committed
21
22
#         Notes:  Must point to a 'merged' .idx file (use generateIdxFiles.py, and make sure you adjust nb_parts)
#                 Make sure idx file is complete (ie open it and check it looks right) - especially if previous job crashed out!
thomas's avatar
thomas committed
23
24
25
#
#   C) Warning: this script requires the presence of the proj_api library for C. 
#               --> If not installed, use MultirateGBR_oldProj.py instead
thomas's avatar
thomas committed
26
27
28
29
30
31
32
#"""

# Import for python
from dgpy import *
from gmshpy import *
from math import *
from InputGBR import *
thomas's avatar
thomas committed
33
34
from ProjectMesh import *
from DiffuseBathGBR import *
thomas's avatar
thomas committed
35
36
37
38
39
40
41
42
from termcolor import colored
import time, os, sys

if(Msg.GetCommRank() == 0):
  print''
  print(colored('SLIM: Multirate GBR simulation starting. Loading 1:Mesh, 2:Bathymetry 3:Reef Map 4:Conservation law.',"yellow"))
  print''

43
  print(colored('      ________ ___    ___  ___  ___         ',"red"))
44
45
46
  print(colored('     /  _____//  /   /  / /   \/   \        ',"red"))
  print(colored(' ____\____  //  /___/  /_/__/\__/\  \_____   ',"red"))
  print(colored('/__________/ \_____/__/___________\__\____\         [ Great Barrier Reef ]',"red"))
thomas's avatar
thomas committed
47
48
  print''
  
thomas's avatar
thomas committed
49
50
51
52
53
54
# Lauch Makefile to generate meshes and download forcings
if(Msg.GetCommRank() == 0):
  try : os.mkdir('./Meshes');
  except: 0;
  try : os.mkdir('./Data');
  except: 0;
thomas's avatar
thomas committed
55
56
  try : os.mkdir('./Bath');
  except: 0;
thomas's avatar
thomas committed
57
58
  try : os.mkdir('./GBR_Functions');
  except: 0;
thomas's avatar
thomas committed
59
60
61
62
  try : os.mkdir(outputDir+'/angle');
  except: 0;
  try : os.mkdir(outputDir+'/XYZ');
  except: 0;
thomas's avatar
thomas committed
63
64
65
  if(not os.path.exists("./Meshes/"+filename+".msh")):
    try : os.system("make "+filename+".msh");
    except: 0;
thomas's avatar
thomas committed
66
67
68
69
70
  if (gbrWind != 0):
    try : os.mkdir(outputDir+'/windUV');
    except: 0;
    try : os.mkdir(outputDir+'/windSourceAcc');
    except: 0;
thomas's avatar
thomas committed
71
72
73
  if (gbrWindFluxAtBoundaries == 1):
    try : os.mkdir(outputDir+'/fWind_depthFactor');
    except: 0;
thomas's avatar
thomas committed
74
75
76
77
78
79
80
81
82
  if (gbrWind == 2):
    try : os.system("make forcings_windNCEP")
    except: 0;
  elif (gbrWind == 3):
    try : os.system("make forcings_windCSFR")
    except: 0;
  if (gbrResidual == 2):
    try : os.system("make forcings_blink")
    except: 0;
thomas's avatar
thomas committed
83
84
    try : os.mkdir(outputDir+'/blinkSol');
    except: 0;
thomas's avatar
thomas committed
85
  elif (gbrResidual == 3):
thomas's avatar
thomas committed
86
    try : os.mkdir('./Data/currentsForcing_MYR_CCH')
thomas's avatar
thomas committed
87
88
89
    except: 0
    try : os.system("make forcings_GBROOS")
    except: 0;
thomas's avatar
thomas committed
90
91
92
93
    try : os.mkdir(outputDir+'/boundaryGBROOS_MYRSol');
    except: 0;
    try : os.mkdir(outputDir+'/boundaryGBROOS_CCHSol');
    except: 0;
thomas's avatar
thomas committed
94
95
96
97
98
99
  try : os.system("make forcings");
  except: 0;
  try : os.mkdir(outputDir);
  except: 0;
  try : os.mkdir(outputDir+'/Evaluator');
  except: 0;
thomas's avatar
thomas committed
100
101
  try : os.mkdir(outputDir+'/lonLat');
  except: 0;
thomas's avatar
thomas committed
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122

Msg.Barrier()

# Definition of Python functions
def printTime (t) :
  return "%dh%02d'%02d\"" % (t/3600,  (t%3600)/60,  t%60)

def printDate (t):
  date=time.gmtime(t)
  return "%02d/%02d/%04d %02dh%02d'%02d\""%(date[2], date[1], date[0], date[3], date[4], date[5])

# Compile libraries
if (Msg.GetCommRank() == 0):
  functionC.buildLibraryFromFile ("Libraries/GBR.cc", "Libraries/lib_gbr.so")
glib = "Libraries/lib_gbr.so"
Msg.Barrier()

# Load model
model = GModel()
if(Msg.GetCommSize() == 1):
  print(colored('Loading mesh for 1 processor (unpartitioned) ...',"red"))
thomas's avatar
thomas committed
123
  loadAndProjectMeshUTM(model, filename, projCodeLonLat, projCodeMesh)
thomas's avatar
thomas committed
124
else:
thomas's avatar
thomas committed
125
  if(os.path.exists("./Meshes/"+filename+"_utm_partitioned_"+str(Msg.GetCommSize())+".msh")):
thomas's avatar
thomas committed
126
    if(Msg.GetCommRank() == 0):
thomas's avatar
thomas committed
127
128
      print'Loading partitioned mesh: ./Meshes/'+filename+'_utm_partitioned_'+str(Msg.GetCommSize())+'.msh'
    model.load('./Meshes/'+filename+'_utm_partitioned_'+str(Msg.GetCommSize())+'.msh')
thomas's avatar
thomas committed
129
  else:
thomas's avatar
thomas committed
130
    Msg.Fatal('No valid partitioned mesh file: ./Meshes/'+filename+'_utm_partitioned_[nbProcs].msh ; run PartitionGBR_Standalone.py first.')
thomas's avatar
thomas committed
131

thomas's avatar
thomas committed
132
# Group collection
thomas's avatar
thomas committed
133
134
groups = dgGroupCollection(model, dimension, order)

thomas's avatar
thomas committed
135
136
# Coordinate system
XYZ = function.getCoordinates()
thomas's avatar
thomas committed
137
138
XYZ_PC = functionPrecomputed(groups, integOrder, 3)
XYZ_PC.compute(XYZ)
thomas's avatar
thomas committed
139
XYZ_DC = dgDofContainer(groups, 3)
thomas's avatar
thomas committed
140
XYZ_DC.L2Projection(XYZ_PC)
thomas's avatar
thomas committed
141
142
143
144
145
146
147
148
exporterXYZ = dgIdxExporter(XYZ_DC, outputDir+'/XYZ')
exporterXYZ.exportIdx()

# Define functions to project UTM points --> lat/lon
UTMtoLatLonDegrees = functionC(glib,"UTMtoLatLonDegrees",3,[XYZ])
UTMtoLonLatDegrees = functionC(glib,"UTMtoLonLatDegrees",3,[XYZ])

# Place above functions into PrecomputedFunction objects
thomas's avatar
thomas committed
149
UTMtoLatLonDegrees_PC = functionPrecomputed(groups, integOrder, 3)
thomas's avatar
thomas committed
150
UTMtoLatLonDegrees_PC.compute(UTMtoLatLonDegrees)
thomas's avatar
thomas committed
151
UTMtoLonLatDegrees_PC = functionPrecomputed(groups, integOrder, 3)
thomas's avatar
thomas committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
UTMtoLonLatDegrees_PC.compute(UTMtoLonLatDegrees)

lonLat_DC = dgDofContainer(groups, 3)
lonLat_DC.L2Projection(UTMtoLonLatDegrees)
exporterLL = dgIdxExporter(lonLat_DC, outputDir+'/lonLat')
exporterLL.exportIdx()

angleProj = functionC(glib, "getAngleProjLonLat", 1, [XYZ])
angleProj_PC = functionPrecomputed(groups, 3, 1)
angleProj_PC.compute(angleProj)
angleProj_DC = dgDofContainer(groups, 1);
angleProj_DC.L2Projection(angleProj_PC)
exporterAngleProj = dgIdxExporter(angleProj_DC, outputDir+'/angle')
exporterAngleProj.exportIdx()
thomas's avatar
thomas committed
166

thomas's avatar
thomas committed
167
168
169
170
171
172
173
174
175
176
# ****************** Cons Law ******************
if (Msg.GetCommRank() == 0):
  print(colored('Setting up conservation law ...',"red"))

claw = dgConservationLawShallowWater2d()
solution = dgDofContainer(groups, claw.getNbFields())

# **********************************************

# ****************** Bathymetry ****************
thomas's avatar
thomas committed
177

thomas's avatar
thomas committed
178
bathDC = dgDofContainer(groups, 1)
thomas's avatar
thomas committed
179
180
if(os.path.exists("./Bath/"+filename+"_bath_smooth/"+filename+"_bath_smooth.idx")):
  Msg.Info("Smoothed bathymetry already exists.")
181
  bathDC.importIdx("./Bath/"+filename+"_bath_smooth/"+filename+"_bath_smooth.idx")
thomas's avatar
thomas committed
182
183
  Msg.Info("Smoothed bathymetry read.")
else:
thomas's avatar
thomas committed
184
  diffuseBathymetry(groups, UTMtoLonLatDegrees, bathDC, filename)
thomas's avatar
thomas committed
185
186
#if using full.bin: diffuseBathymetry(bathDC, groups, lonLat)

thomas's avatar
thomas committed
187
# **********************************************
thomas's avatar
thomas committed
188

thomas's avatar
thomas committed
189
# ****************** Reefs *********************
thomas's avatar
thomas committed
190
191
192
193
194
195
196
197
198
if (Msg.GetCommRank() == 0):
  print(colored('Loading reef map ...',"red"))

reefsDC = dgDofContainer(groups,1)

reefsDCfilename = "GBR_Functions/reefsFunction_"+filename+"/reefsFunction_"+filename+".idx"
if(os.path.exists(reefsDCfilename)):
  if (Msg.GetCommRank() == 0):
    print 'Reef Map file found: '+reefsDCfilename+' ; loading this.'
199
  reefsDC.importIdx(reefsDCfilename)
thomas's avatar
thomas committed
200
201
202
else:
  if (Msg.GetCommRank() == 0):
    print 'Reef Map file '+reefsDCfilename+' not found. Loading reef map ...'
thomas's avatar
thomas committed
203
204
205
206
207
208
209
210
211
    #1. Create a reefs object (initialise object and create reefs function)
    reefsObj = reefsSmooth(reefMapFile, UTMtoLonLatDegrees)
    #2. Interpolate the function onto DC
    reefsDC.interpolate(reefsObj)
    exporterReefsDC=dgIdxExporter(reefsDC, 'GBR_Functions/reefsFunction_'+filename)
    exporterReefsDC.exportIdx(0,0)
    #3. Delete reefs object data
    reefsObj.clearReefsData()
Msg.Barrier()
212
reefsDC.importIdx(reefsDCfilename)
thomas's avatar
thomas committed
213
214
215
216
217
218
219
220
221
222
223
224

# ************************************************

t = Ti
timeFunction = functionConstant(t)

# ****************** Forcings ******************
if (Msg.GetCommRank() == 0):
  print(colored('Setting up forcings ...',"red"))
# .................. Tides    ..................

# Tpxo tide 
225
226
227
tpxoTide = None
if (gbrTide == 1):
  tideEta = slimFunctionTpxo("./Data/h_tpxo7.2.nc","ha","hp","lon_z","lat_z")
thomas's avatar
thomas committed
228
  tideEta.setCoordAndTimeFunctions(UTMtoLonLatDegrees_PC,timeFunction)
229
  tideU = slimFunctionTpxo("./Data/u_tpxo7.2.nc","Ua","up","lon_u","lat_u")
thomas's avatar
thomas committed
230
  tideU.setCoordAndTimeFunctions(UTMtoLonLatDegrees_PC,timeFunction)
231
  tideV = slimFunctionTpxo("./Data/u_tpxo7.2.nc","Va","vp","lon_v","lat_v")
thomas's avatar
thomas committed
232
233
  tideV.setCoordAndTimeFunctions(UTMtoLonLatDegrees_PC,timeFunction)
  tideUV = functionC(glib, "lonLatVector_To_UTMVector", 3, [angleProj_PC, tideU, tideV])
234
235
236

  tideuv = functionC(glib,"transport2velocity",3,[tideUV, bathDC.getFunction(), solution.getFunction()])
  tpxoTide = functionC(glib, "merge", 3, [tideEta, tideuv])
thomas's avatar
thomas committed
237

thomas's avatar
thomas committed
238
239
# .................. Wind     ..................

thomas's avatar
thomas committed
240
241
242
243
244
245
246
247
248
249
250
251
252
windInterpolatorU = slimStructDataInterpolatorMultilinear()
windInterpolatorV = slimStructDataInterpolatorMultilinear()
windUcontainer = None
windVcontainer = None
if gbrWind == 1:
  windUcontainer = slimStructDataContainerTemporalSerie("./Data/AIMS_wind_Davies_U.txt")
  windVcontainer = slimStructDataContainerTemporalSerie("./Data/AIMS_wind_Davies_V.txt")
elif gbrWind == 2:
  windUcontainer = slimStructDataContainerNetcdf("./Data/uwind-2000-2011.nc", "uwnd", "lat", "lon")
  windVcontainer = slimStructDataContainerNetcdf("./Data/vwind-2000-2011.nc", "vwnd", "lat", "lon")
elif gbrWind == 3:
  windUcontainer = slimStructDataContainerNetcdf("./Data/wnd10m_csfr_2008-2010.grb2.nc", "U_GRD_L103", "lat", "lon")
  windVcontainer = slimStructDataContainerNetcdf("./Data/wnd10m_csfr_2008-2010.grb2.nc", "V_GRD_L103", "lat", "lon")
thomas's avatar
thomas committed
253
254
windU = slimFunctionStructData(windUcontainer, windInterpolatorU, UTMtoLatLonDegrees_PC)
windV = slimFunctionStructData(windVcontainer, windInterpolatorV, UTMtoLatLonDegrees_PC)
thomas's avatar
thomas committed
255
256
257
windU.setTimeFunction(timeFunction)
windV.setTimeFunction(timeFunction)

thomas's avatar
thomas committed
258
#windU_PostSetTimeFn_DC = dgDofContainer(groups, 1)
thomas's avatar
thomas committed
259
#windU_PostSetTimeFn_DC.L2Projection(windU)
thomas's avatar
thomas committed
260
#windV_PostSetTimeFn_DC = dgDofContainer(groups, 1)
thomas's avatar
thomas committed
261
#windV_PostSetTimeFn_DC.L2Projection(windV)
thomas's avatar
thomas committed
262
263
264
#windU_PostSetTimeFn_DC.exportMsh('windU_PostSetTimeFn')
#windV_PostSetTimeFn_DC.exportMsh('windV_PostSetTimeFn')

thomas's avatar
thomas committed
265
windUV = functionC(glib, "lonLatVector_To_UTMVector", 3, [angleProj_PC, windU, windV])
thomas's avatar
thomas committed
266

thomas's avatar
thomas committed
267
268
# .... Define wind-induced current functions
currentFromWind_North = None
thomas's avatar
thomas committed
269
270
currentFromWind_Central = None
currentFromWind_South = None
thomas's avatar
thomas committed
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
if gbrWindFluxAtBoundaries == 1 and gbrWind != 0:
  #currentFromWind_North = functionC(glib, "functionMultiply", 2, [windUV, northMultFactor])
  #currentFromWind_Central = functionC(glib, "functionMultiply", 2, [windUV, centralMultFactor])
  #currentFromWind_South = functionC(glib, "functionMultiply", 2, [windUV, southMultFactor])
  
  functionCurrentAsFWind_depthFactor = functionC(glib, "functionCurrentAsFWind_calcDepthFactor", 1, [bathDC.getFunction()])
  functionCurrentAsFWind_depthFactor_PC = functionPrecomputed(groups, integOrder, 1)
  functionCurrentAsFWind_depthFactor_PC.compute(functionCurrentAsFWind_depthFactor)

  depthFactorDC = dgDofContainer(groups, 1)
  depthFactorDC.L2Projection(functionCurrentAsFWind_depthFactor_PC)
  depthFactorExporter = dgIdxExporter(depthFactorDC, outputDir+'/fWind_depthFactor')
  depthFactorExporter.exportIdx(0,0)
  
  currentFromWind_North = functionC(glib, "functionCurrentAsFWind", 2, [windUV, northMultFactor, functionCurrentAsFWind_depthFactor_PC])
  currentFromWind_Central = functionC(glib, "functionCurrentAsFWind", 2, [windUV, centralMultFactor, functionCurrentAsFWind_depthFactor_PC])
  currentFromWind_South = functionC(glib, "functionCurrentAsFWind", 2, [windUV, southMultFactor, functionCurrentAsFWind_depthFactor_PC])
  

# ............ Boundary current functions    .............
boundaryUVElev_North = None
boundaryUVElev_Central = None
boundaryUVElev_South = None
boundaryUVElev_North_WIP = None
boundaryUVElev_Central_WIP = None
boundaryUVElev_South_WIP = None
boundaryUVElevAndTide_North = None
boundaryUVElevAndTide_Central = None
boundaryUVElevAndTide_South = None
boundaryUVElevAndTide_North_WIP = None
boundaryUVElevAndTide_Central_WIP = None
boundaryUVElevAndTide_South_WIP = None
thomas's avatar
thomas committed
303
304
305
306

# .... Constant NCJ
if gbrResidual == 1:
  #Put data into functions
thomas's avatar
thomas committed
307
308
309
  data_A = functionConstant([F_A, vx_A_bulbSouth, vy_A_bulbSouth])
  data_B = functionConstant([F_B, vx_B_bulbSouth, vy_B_bulbSouth])
  data_C = functionConstant([F_C, vx_C_bulbSouth, vy_C_bulbSouth])
thomas's avatar
thomas committed
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325

  # NCJ Current forcing
  bathTimesVector_A = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_A])
  bathTimesVector_B = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_B])
  bathTimesVector_C = functionC(glib, "bathVec", 3, [bathDC.getFunction(), function.getNormals(), data_C])
  bathIntegralNorthFM = fullMatrixDouble(1,1)
  bathIntegralCentralFM = fullMatrixDouble(1,1)
  bathIntegralSouthFM = fullMatrixDouble(1,1)
  dgFunctionIntegratorInterface(groups, bathTimesVector_A).compute('Open Sea North', bathIntegralNorthFM)
  dgFunctionIntegratorInterface(groups, bathTimesVector_B).compute('Open Sea Central', bathIntegralCentralFM)
  dgFunctionIntegratorInterface(groups, bathTimesVector_C).compute('Open Sea South', bathIntegralSouthFM)
  bathIntegralsFuncNorth = functionConstant(bathIntegralNorthFM(0,0))
  bathIntegralsFuncCentral = functionConstant(bathIntegralCentralFM(0,0))
  bathIntegralsFuncSouth = functionConstant(bathIntegralSouthFM(0,0))

  # Generate SEC UV functions
thomas's avatar
thomas committed
326
327
328
  boundaryUVElev_North_WIP = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncNorth, function.getNormals(), data_A])
  boundaryUVElev_Central_WIP = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncCentral, function.getNormals(), data_B])
  boundaryUVElev_South_WIP = functionC(glib, "SEC_UV", 3, [bathIntegralsFuncSouth, function.getNormals(), data_C])
thomas's avatar
thomas committed
329
330
331
332
333
334
335
336
337
338
339
340
  
  if (Msg.GetCommRank() == 0):
    outputVec = open('GBR_Functions/CurrentUnitVectors.pos', 'w')
    outputVec.write('View \"CurrentUnitVectors\" {\n')
    outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_A + ',' + '%.5f'%vy_A + ',0};\n')
    outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_B + ',' + '%.5f'%vy_B + ',0};\n')
    outputVec.write('VP (' + '%.5f'%0 + ',' + '%.5f'%0 + ',0) {' + '%.5f'%vx_C + ',' + '%.5f'%vy_C + ',0};\n')
    outputVec.write('};')
    outputVec.close()

# .... Bluelink
elif gbrResidual == 2:
thomas's avatar
thomas committed
341
342
343
  blinkInterpolatorU = slimStructDataInterpolatorMultilinear()
  blinkInterpolatorV = slimStructDataInterpolatorMultilinear()
  blinkInterpolatorEta = slimStructDataInterpolatorMultilinear()
thomas's avatar
thomas committed
344
345
346
  blinkUContainer = slimStructDataContainerNetcdf("./Data/bluelink_transport_27nov2007_35days.nc", "u", "lat", "lon")
  blinkVContainer = slimStructDataContainerNetcdf("./Data/bluelink_transport_27nov2007_35days.nc", "v", "lat", "lon")
  blinkEtaContainer = slimStructDataContainerNetcdf("./Data/bluelink_eta_27nov2007_35days.nc", "eta", "lat", "lon") 
thomas's avatar
thomas committed
347
348
349
  blinkU = slimFunctionStructData(blinkUContainer, blinkInterpolatorU, UTMtoLatLonDegrees_PC)
  blinkV = slimFunctionStructData(blinkVContainer, blinkInterpolatorV, UTMtoLatLonDegrees_PC)
  blinkEta = slimFunctionStructData(blinkEtaContainer, blinkInterpolatorEta, UTMtoLatLonDegrees_PC)
thomas's avatar
thomas committed
350
351
352
353
354
  blinkU.setTimeFunction(timeFunction)
  blinkV.setTimeFunction(timeFunction)
  blinkEta.setTimeFunction(timeFunction)
  #blinkEta.setMinValue(-3.0)
  #blinkEta.setMaxValue(3.0)
thomas's avatar
thomas committed
355
  blinkUV = functionC(glib, "lonLatVector_To_UTMVector", 3, [angleProj_PC, blinkU, blinkV])
thomas's avatar
thomas committed
356
357
358
  minFC = functionConstant(-3.0)
  blinkEta_cutOff = functionC(glib, "cutOffMin", 1, [blinkEta, minFC])

thomas's avatar
thomas committed
359
  blinkuv = functionC(glib,"transport2velocity",3,[blinkUV, bathDC.getFunction(), solution.getFunction()])
thomas's avatar
thomas committed
360
361
362
  boundaryUVElev_North_WIP = functionC(glib, "merge", 3, [blinkEta_cutOff, blinkuv])
  boundaryUVElev_Central_WIP = functionC(glib, "merge", 3, [blinkEta_cutOff, blinkuv])
  boundaryUVElev_South_WIP = functionC(glib, "merge", 3, [blinkEta_cutOff, blinkuv])
thomas's avatar
thomas committed
363
364

# .... GBROOS Boundary Currents
thomas's avatar
thomas committed
365
#NB: -9 to GBROOS input files if necessary (GMT time offset)
thomas's avatar
thomas committed
366
elif gbrResidual == 3:
thomas's avatar
thomas committed
367
368
369
370
  #Put data into functions
  data_Central = functionConstant([vx_B, vy_B])
  data_South = functionConstant([vx_C, vy_C])

thomas's avatar
thomas committed
371
372
373
  boundaryInterpolatorU_MYR = slimStructDataInterpolatorMultilinear()
  boundaryInterpolatorV_MYR = slimStructDataInterpolatorMultilinear()
  boundaryInterpolatorEta_MYR = slimStructDataInterpolatorMultilinear()
thomas's avatar
thomas committed
374
375
376
  boundaryElevContainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_elev_novdec07janfeb08.txt")
  boundaryUcontainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_ucur_novdec07janfeb08.txt")
  boundaryVcontainer_MYR = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/MYR_vcur_novdec07janfeb08.txt")
thomas's avatar
thomas committed
377
378
379
  boundaryElev_MYR = slimFunctionStructData(boundaryElevContainer_MYR, boundaryInterpolatorEta_MYR, UTMtoLatLonDegrees_PC)
  boundaryU_MYR = slimFunctionStructData(boundaryUcontainer_MYR, boundaryInterpolatorU_MYR, UTMtoLatLonDegrees_PC)
  boundaryV_MYR = slimFunctionStructData(boundaryVcontainer_MYR, boundaryInterpolatorV_MYR, UTMtoLatLonDegrees_PC)
thomas's avatar
thomas committed
380
381
382
  boundaryElev_MYR.setTimeFunction(timeFunction)
  boundaryU_MYR.setTimeFunction(timeFunction)
  boundaryV_MYR.setTimeFunction(timeFunction)
thomas's avatar
thomas committed
383
  boundaryUV_MYR = functionC(glib, "lonLatVector_To_UTMVector", 3, [angleProj_PC, boundaryU_MYR, boundaryV_MYR])
thomas's avatar
thomas committed
384
  boundaryUV_modulated_MYR = functionC(glib,"currentsAlongVecDotNormals", 3, [function.getNormals(), data_Central, boundaryUV_MYR])
thomas's avatar
thomas committed
385
  boundaryUVElev_Central_WIP = functionC(glib, "merge", 3, [boundaryElev_MYR, boundaryUV_modulated_MYR])
thomas's avatar
thomas committed
386

thomas's avatar
thomas committed
387
388
389
  boundaryInterpolatorU_CCH = slimStructDataInterpolatorMultilinear()
  boundaryInterpolatorV_CCH = slimStructDataInterpolatorMultilinear()
  boundaryInterpolatorEta_CCH = slimStructDataInterpolatorMultilinear()
thomas's avatar
thomas committed
390
391
392
  boundaryElevContainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_elev_novdec07janfeb08.txt")
  boundaryUcontainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_ucur_novdec07janfeb08.txt")
  boundaryVcontainer_CCH = slimStructDataContainerTemporalSerie("./Data/currentsForcing_MYR_CCH/CCH_vcur_novdec07janfeb08.txt")
thomas's avatar
thomas committed
393
394
395
  boundaryElev_CCH = slimFunctionStructData(boundaryElevContainer_CCH, boundaryInterpolatorEta_CCH, UTMtoLatLonDegrees_PC)
  boundaryU_CCH = slimFunctionStructData(boundaryUcontainer_CCH, boundaryInterpolatorU_CCH, UTMtoLatLonDegrees_PC)
  boundaryV_CCH = slimFunctionStructData(boundaryVcontainer_CCH, boundaryInterpolatorV_CCH, UTMtoLatLonDegrees_PC)
thomas's avatar
thomas committed
396
397
398
  boundaryElev_CCH.setTimeFunction(timeFunction)
  boundaryU_CCH.setTimeFunction(timeFunction)
  boundaryV_CCH.setTimeFunction(timeFunction)
thomas's avatar
thomas committed
399
  boundaryUV_CCH = functionC(glib, "lonLatVector_To_UTMVector", 3, [angleProj_PC, boundaryU_CCH, boundaryV_CCH])
thomas's avatar
thomas committed
400
  boundaryUV_modulated_CCH = functionC(glib,"currentsAlongVecDotNormals", 3, [function.getNormals(), data_South, boundaryUV_CCH])
thomas's avatar
thomas committed
401
402
403
  boundaryUVElev_South_WIP = functionC(glib, "merge", 3, [boundaryElev_CCH, boundaryUV_modulated_CCH])
  
  boundaryUVElev_North_WIP = functionConstant([0.0,0.0,0.0])
thomas's avatar
thomas committed
404

thomas's avatar
thomas committed
405
# .... Constant boundary currents
thomas's avatar
thomas committed
406
elif gbrResidual == 4:
thomas's avatar
thomas committed
407
408
409
410
  elevFactor = functionConstant([0.0])
  northAddFactor = functionConstant([0.0, 0.0])
  #centralAddFactor = functionConstant([0.0, -0.10])
  centralAddFactor = functionC(glib, "getRes3", 2, [UTMtoLonLatDegrees_PC])
thomas's avatar
thomas committed
411
412
  southAddFactor = functionConstant([0.08, -0.08])
  
thomas's avatar
thomas committed
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
  #Res3DC = dgDofContainer(groups, 2)
  #Res3DC.L2Projection(centralAddFactor)
  #Res3DC.exportMsh("Res3DC")
  
  boundaryUVElev_North_WIP = functionC(glib, "merge", 3, [elevFactor, northAddFactor])
  boundaryUVElev_Central_WIP = functionC(glib, "merge", 3, [elevFactor, centralAddFactor])
  boundaryUVElev_South_WIP = functionC(glib, "merge", 3, [elevFactor, southAddFactor])


# .... Add tides to boundary conditions
if gbrTide == 1 and gbrResidual != 0:
  boundaryUVElevAndTide_North_WIP = functionC(glib,"merge_3DFunctions",3,[tpxoTide, boundaryUVElev_North_WIP])
  boundaryUVElevAndTide_Central_WIP = functionC(glib, "merge_3DFunctions", 3, [tpxoTide, boundaryUVElev_Central_WIP])
  boundaryUVElevAndTide_South_WIP = functionC(glib, "merge_3DFunctions", 3, [tpxoTide, boundaryUVElev_South_WIP])

    
# .... Add wind-induced boundary currents: UV = UV_existing + f(windUV)
if gbrWindFluxAtBoundaries == 1 and gbrWind != 0:
  if gbrResidual == 0:
    if gbrTide == 0:
      eta = functionConstant([0.0])
      boundaryUVElev_North = functionC(glib, "merge", 3, [eta, currentFromWind_North])
      boundaryUVElev_Central = functionC(glib, "merge", 3, [eta, currentFromWind_Central])
      boundaryUVElev_South = functionC(glib, "merge", 3, [eta, currentFromWind_South])
    elif gbrTide == 1:
      boundaryUVElevAndTide_North = functionC(glib, "merge3D_with_UV", 3, [tpxoTide, currentFromWind_North])
      boundaryUVElevAndTide_Central = functionC(glib, "merge3D_with_UV", 3, [tpxoTide, currentFromWind_Central])
      boundaryUVElevAndTide_South = functionC(glib, "merge3D_with_UV", 3, [tpxoTide, currentFromWind_South])
  else:
    if gbrTide == 0:
      boundaryUVElev_North = functionC(glib, "merge3D_with_UV", 3, [boundaryUVElev_North_WIP, currentFromWind_North])
      boundaryUVElev_Central = functionC(glib, "merge3D_with_UV", 3, [boundaryUVElev_Central_WIP, currentFromWind_Central])
      boundaryUVElev_South = functionC(glib, "merge3D_with_UV", 3, [boundaryUVElev_South_WIP, currentFromWind_South])
    elif gbrTide == 1:
      boundaryUVElevAndTide_North = functionC(glib,"merge3D_with_UV",3,[boundaryUVElevAndTide_North_WIP, currentFromWind_North])
      boundaryUVElevAndTide_Central = functionC(glib,"merge3D_with_UV",3,[boundaryUVElevAndTide_Central_WIP, currentFromWind_Central])
      boundaryUVElevAndTide_South = functionC(glib,"merge3D_with_UV",3,[boundaryUVElevAndTide_South_WIP, currentFromWind_South])
else:
  if gbrTide == 0:
    boundaryUVElev_North = boundaryUVElev_North_WIP
    boundaryUVElev_Central = boundaryUVElev_Central_WIP
    boundaryUVElev_South = boundaryUVElev_South_WIP
  elif gbrTide == 1:
    boundaryUVElevAndTide_North = boundaryUVElevAndTide_North_WIP
    boundaryUVElevAndTide_Central = boundaryUVElevAndTide_Central_WIP
    boundaryUVElevAndTide_South = boundaryUVElevAndTide_South_WIP

thomas's avatar
thomas committed
460
  
thomas's avatar
thomas committed
461
462
463
464
465
# **********************************************

# ** OPTIONAL 1/3: Load an existing DC as the initial solution DC

if (startFromExistingSolution == 1):
466
  solution.importIdx( solutionFileName )
thomas's avatar
thomas committed
467
468
469
470
471
472

# **********************************************

ld = functionConstant(0)

#Cd = functionC(glib, "bottomDrag_old", 1, [solution.getFunction(), bathDC.getFunction()])
thomas's avatar
thomas committed
473
#Cd = functionC(glib, "bottomDrag", 1, [ solution.getFunction(), bathDC.getFunction(), reefsDC.getFunction() ])
thomas's avatar
thomas committed
474
Cd = functionC(glib, "bottomDrag_const_n016", 1, [ solution.getFunction(), bathDC.getFunction() ])
thomas's avatar
thomas committed
475
476
477
478
if(gbrDiff == 0):
  Di = functionConstant(diff)
elif(gbrDiff == 1):
  Di = functionC(glib, "smagorinsky", 1, [solution.getFunctionGradient()]) #***
thomas's avatar
thomas committed
479
Coriolis = functionC(glib, "coriolisLonLatDegrees", 1, [UTMtoLonLatDegrees])
thomas's avatar
thomas committed
480
481
if(gbrWind == 0):
  source = functionConstant([0.0, 0.0])
thomas's avatar
thomas committed
482
else:
thomas's avatar
thomas committed
483
  source = functionC(glib, "windStressSmithBanke", 2, [windUV, solution.getFunction(), bathDC.getFunction()])
thomas's avatar
thomas committed
484
485
486
487
488
489
490
491
492
493

claw.setCoriolisFactor(Coriolis)
claw.setQuadraticDissipation(Cd)
claw.setLinearDissipation(ld)
claw.setDiffusivity(Di)
claw.setSource(source)
claw.setBathymetry(bathDC.getFunction())
claw.setBathymetryGradient(bathDC.getFunctionGradient())

if(gbrResidual == 0):
thomas's avatar
thomas committed
494
  if(gbrTide == 0 and gbrWindFluxAtBoundaries == 0):
thomas's avatar
thomas committed
495
496
497
    claw.addBoundaryCondition('Open Sea North',claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea Central',claw.newBoundaryWall())
    claw.addBoundaryCondition('Open Sea South',claw.newBoundaryWall())
thomas's avatar
thomas committed
498
499
    claw.addBoundaryCondition('Open Sea Neutral',claw.newBoundaryWall())
  elif(gbrTide == 1 and gbrWindFluxAtBoundaries == 0):
thomas's avatar
thomas committed
500
501
502
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", tpxoTide))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", tpxoTide))
thomas's avatar
thomas committed
503
504
505
506
507
508
509
510
511
512
513
514
    claw.addBoundaryCondition('Open Sea Neutral',claw.newOutsideValueBoundary("Surface", tpxoTide))
  if(gbrTide == 0 and gbrWindFluxAtBoundaries == 1):
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", boundaryUVElev_North))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", boundaryUVElev_Central))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", boundaryUVElev_South))
    claw.addBoundaryCondition('Open Sea Neutral',claw.newOutsideValueBoundary("Surface", boundaryUVElev_Central))
  elif(gbrTide == 1 and gbrWindFluxAtBoundaries == 1):
    claw.addBoundaryCondition('Open Sea North', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_North))
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_Central))
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_South))
    claw.addBoundaryCondition('Open Sea Neutral',claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_Central))
elif(gbrResidual == 1 or gbrResidual == 2 or gbrResidual == 3 or gbrResidual == 4):
thomas's avatar
thomas committed
515
  if(gbrTide == 0):
thomas's avatar
thomas committed
516
517
518
519
    claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", boundaryUVElev_North))
    claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", boundaryUVElev_Central))
    claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", boundaryUVElev_South))
    claw.addBoundaryCondition('Open Sea Neutral',claw.newOutsideValueBoundary("Surface", boundaryUVElev_Central))
thomas's avatar
thomas committed
520
  elif(gbrTide == 1):
thomas's avatar
thomas committed
521
522
523
524
    claw.addBoundaryCondition('Open Sea North', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_North))
    claw.addBoundaryCondition('Open Sea Central', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_Central))
    claw.addBoundaryCondition('Open Sea South', claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_South))
    claw.addBoundaryCondition('Open Sea Neutral',claw.newOutsideValueBoundary("Surface", boundaryUVElevAndTide_Central))
thomas's avatar
thomas committed
525

thomas's avatar
thomas committed
526
claw.addBoundaryCondition('Coast',claw.newBoundaryWall())
thomas's avatar
thomas committed
527
claw.addBoundaryCondition('Island',claw.newBoundaryWall())
thomas's avatar
thomas committed
528
529
530
531
532
533
claw.addBoundaryCondition('Islands',claw.newBoundaryWall())

# **********************************************
#NB: Must make 1st argument of Cd solution.getFunction() to use the following:
bottomDragDC = dgDofContainer(groups, 1)
bottomDragDC.interpolate(Cd)
thomas's avatar
thomas committed
534
exporterQuadDissip = dgIdxExporter(bottomDragDC, 'GBR_Functions/quadDissipation_'+filename)
thomas's avatar
thomas committed
535
536
537
538
539
540
541
542
543
544
545
546
exporterQuadDissip.exportIdx(0,0)

# ****************** Time Int ******************
#WARNING: in MultiRate, all DofContainers are split into groups
#       ---> Must add every DC used by claw to splitForMultiRate!!!
if (Msg.GetCommRank() == 0):
  print(colored('Setting up time integrator ...',"red"))

rk = dgMultirateERK(groups, claw, RK_TYPE)
dt = rk.splitGroupsForMultirate(mL, solution, [solution, bathDC, reefsDC])
rk.printMultirateInfo(0)

thomas's avatar
thomas committed
547
548
549
550
# Recompute Precomputed functions after splitting of element groups
UTMtoLonLatDegrees_PC.compute(UTMtoLonLatDegrees)
UTMtoLatLonDegrees_PC.compute(UTMtoLatLonDegrees)
angleProj_PC.compute(angleProj)
thomas's avatar
thomas committed
551
functionCurrentAsFWind_depthFactor_PC.compute(functionCurrentAsFWind_depthFactor)
thomas's avatar
thomas committed
552

thomas's avatar
thomas committed
553
554
555
556
557
#OPTIONAL 2/3: if starting from existing solution, set dt = old_dt
if (startFromExistingSolution == 1):
  dt = old_dt
# **********************************************

thomas's avatar
thomas committed
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
##Write dissipation term
#dissipationTermFC = functionC(glib, "dissipationTerm", 2, [solution.getFunction(), solution.getFunctionGradient(), Di, bathDC.getFunction(), bathDC.getFunctionGradient()])
#dissipationTermDC = dgDofContainer(groups, 2)
#dissipationTermDC.interpolate(dissipationTermFC)
#exporterDissipationTerm = dgIdxExporter(dissipationTermDC, outputDir+'/DissipationTerm/dissipationTerm')
#exporterDissipationTerm.exportIdx(0,0)

#Write windUV term
windUVDC = None
if (gbrWind != 0):
  windUVDC = dgDofContainer(groups, 3)
  windUVDC.L2Projection(windUV)
  exporterWindUV = dgIdxExporter(windUVDC, outputDir+'/windUV')
  exporterWindUV.exportIdx(0, t)
  
#Write windStress term
windSourceDC = None
if (gbrWind != 0):
  windSourceDC = dgDofContainer(groups, 2)
thomas's avatar
thomas committed
577
  windSourceDC.L2Projection(source)
thomas's avatar
thomas committed
578
579
580
581
582
583
584
  exporterWindSourceAcc = dgIdxExporter(windSourceDC, outputDir+'/windSourceAcc')
  exporterWindSourceAcc.exportIdx(0, t)
  
#Write blink term
blinkSolDC = None
if(gbrResidual == 2):
  blinkSolDC = dgDofContainer(groups, 3)
thomas's avatar
thomas committed
585
586
  blinkSolDC.L2Projection(boundaryUVElev_Central)
  exporterBlinkSol=dgIdxExporter(blinkSolDC, outputDir+'/blinkSol_boundaryTerm')
thomas's avatar
thomas committed
587
588
  exporterBlinkSol.exportIdx(0, t)
  
thomas's avatar
thomas committed
589
590
591
592
593
594
595
596
597
598
599
600
##Write GBROOS boundary currents/elev term
#boundaryGBROOS_MYRSol_DC = None
#boundaryGBROOS_CCHSol_DC = None
#if(gbrResidual == 3):
  #boundaryGBROOS_MYRSol_DC = dgDofContainer(groups, 3)
  #boundaryGBROOS_MYRSol_DC.L2Projection(boundaryUVElev_MYR)
  #exporterGBROOS_MYRSol=dgIdxExporter(boundaryGBROOS_MYRSol_DC, outputDir+'/boundaryGBROOS_MYRSol')
  #exporterGBROOS_MYRSol.exportIdx(0, t)
  #boundaryGBROOS_CCHSol_DC = dgDofContainer(groups, 3)
  #boundaryGBROOS_CCHSol_DC.L2Projection(boundaryUVElev_CCH)
  #exporterGBROOS_CCHSol=dgIdxExporter(boundaryGBROOS_CCHSol_DC, outputDir+'/boundaryGBROOS_CCHSol')
  #exporterGBROOS_CCHSol.exportIdx(0, t)
thomas's avatar
thomas committed
601

thomas's avatar
thomas committed
602
currentFunction = functionC(glib,"current",3,[solution.getFunction()]) #***
thomas's avatar
thomas committed
603
604
if gbrTide != 0:
  tideCurrentFunction = functionC(glib,"current",3,[tpxoTide])
thomas's avatar
thomas committed
605
606
607
608
609
610
611
612
613
614
615

nbSteps = int(ceil((Tf-Ti)/dt))
nbExport = 0

if(Msg.GetCommRank() == 0):
  if(os.path.exists(outputDir+'/mr-gbr')):
    try : os.system("rm -r "+outputDir+'/mr-gbr');
    except: 0;
exporterSol=dgIdxExporter(solution, outputDir+'/mr-gbr')
#exporterSol=dgIdxExporter(currentFunction, outputDir+'/mr-gbr-func')

616
617
618
619
620
621
622
norm = solution.norm()
waterDepth = functionC(glib, "getWaterDepth", 1, [solution.getFunction(), bathDC.getFunction()])
integrator = dgFunctionIntegrator(groups,  waterDepth,  solution)
totalMass=fullMatrixDouble(1, 1)
totalMassRef=fullMatrixDouble(1, 1)
integrator.compute(totalMassRef, "")

thomas's avatar
thomas committed
623
624
625
626
627
628
629
630
631
632
633

if (Msg.GetCommRank() == 0):
  print ''
  print(colored('******************** Initialising ********************', "blue"))
  print ''
  print '     - Number of processors:', Msg.GetCommSize()
  print '     - Multirate Runge Kutta Scheme:', RK_TYPE
  print '     - Initial Time (Ti):',  printDate(Ti)
  print '     - Final Time (Tf):',  printDate(Tf)
  print '     - Simulation Length (days):',  (Tf-Ti)/3600.0/24.0  
  print '     - Reference Time Step:', printTime(dt), ', which is', str(dt), 's'
thomas's avatar
thomas committed
634
  print '     - Number of Iterations:', nbSteps
thomas's avatar
thomas committed
635
  print '     - Theoretical Speedup:', rk.speedUp()
thomas's avatar
thomas committed
636
  print '     - Norm of Solution at Ti:',  norm
637
  print '     - Mass of Solution at Ti:',  totalMassRef(0, 0)
thomas's avatar
thomas committed
638
  print ''
thomas's avatar
thomas committed
639
  print '     - Tides:', gbrTideLabels[gbrTide], ' || Wind:', gbrWindLabels[gbrWind], ' || Wind-induced boundary UV:', gbrWindFluxAtBoundariesLabels[gbrWindFluxAtBoundaries], ' || Residual current:', gbrResidualLabels[gbrResidual]
thomas's avatar
thomas committed
640
  print ''
thomas's avatar
thomas committed
641
642
  print '     - Output Directory:', outputDir
  print ''
thomas's avatar
thomas committed
643
644
645
646
647
648
649
650
651
  print(colored('******************************************************', "blue"))
  print ''
  
  #Write info file for LPTracker
  simInfo = [str(dt)+'\n', str(nbSteps)+'\n', str(export)+'\n', str(Ti)+'\n', str(outputDir)+'\n', str(filename)]
  simInfoFile = open(outputDir+'/simulationInfo.txt','w')
  simInfoFile.writelines(simInfo)
  simInfoFile.close()

652
#Initialise Function Evaluators
thomas's avatar
thomas committed
653
654
655
FunctionEval = dgFunctionEvaluator(groups, currentFunction)
BathEval = dgFunctionEvaluator(groups, bathDC.getFunction())
ElevationEval = dgFunctionEvaluator(groups, solution.getFunction())
thomas's avatar
thomas committed
656
657
658
WindUVEval = None
if (gbrWind != 0):
  WindUVEval = dgFunctionEvaluator(groups, windUV)
thomas's avatar
thomas committed
659
660
res = fullMatrixDouble()

thomas's avatar
thomas committed
661
#Evaluate Bathymetry Evaluator
thomas's avatar
thomas committed
662
663
outBathEval = open(outputDir+'/Evaluator/evalBathymetry_StdLocations.dat',"w")
outBathGBROOSEval = open(outputDir+'/Evaluator/evalBathymetry_GBROOSLocations.dat',"w")
thomas's avatar
thomas committed
664
665
666
667
668
for point in range(0,len(evalPointsLatLon)):
  #Only write to file if we are on correct cpu:
  BathEval.compute(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2], res)
  outBathEval.write(str(res(0,0))+'\n')
outBathEval.close()
thomas's avatar
thomas committed
669
670
671
672
673
for point in range(0,len(evalPointsGBROOSLatLon)):
  #Only write to file if we are on correct cpu:
  BathEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
  outBathGBROOSEval.write(str(res(0,0))+'\n')
outBathGBROOSEval.close()
thomas's avatar
thomas committed
674

675
676
677
678
679
680
#Create total mass file
if (Msg.GetCommRank() == 0):
  outMass = open(outputDir+'/totalMass.dat',"w")
  outMass.write(str(t/3600-Ti/3600)+'\t'+str(totalMassRef(0,0))+'\n')
  outMass.close()

thomas's avatar
thomas committed
681
682
683
684
685
686
# **********************************************
exporterSol.exportIdx(0, t)
# ****************** Iterate  ******************
t_exportStart = Ti + t_exportStart
startcpu = time.clock()

687
688
689
690
691
##TEST: OUTPUT VALUE OF WIND (PART 1 OF 2)
#if (Msg.GetCommRank() == 0):
  #for siteId in range(0,len(evalPointsGBROOSLatLon)):
    #outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[siteId][3])+'.dat',"w")
    #outWind.close()
thomas's avatar
thomas committed
692

thomas's avatar
thomas committed
693
694
695
696
697
698
699
700
701
# OPTIONAL 3/3: if starting from an initial solution, adjust start & end steps accordingly
firstStep = 1
lastStep = 1
if (startFromExistingSolution == 1):
  firstStep = stepNbOfInitialSolution + 1
  lastStep = stepNbOfInitialSolution + nbSteps + 1
else:
  firstStep = 1
  lastStep = nbSteps + 1
thomas's avatar
thomas committed
702

thomas's avatar
thomas committed
703
704
705
if (Msg.GetCommRank() == 0):
  print 'Starting simulation from iteration #',firstStep,' ...'
  print ''
thomas's avatar
thomas committed
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
for i in range(firstStep, lastStep):
  if(i == lastStep-1):
    dt=Tf-t
  norm = rk.iterate (solution, dt, t)
  t = t +dt
  timeFunction.set(t)
  if ( i%iter == 0 ):
     norm=solution.norm()
     if(Msg.GetCommRank() == 0):
       print''
       print(colored('|ITER| %d of %d',  "red")%(i,lastStep))
       print(colored('------------------------------------------------------', "red"))
       print'|TIME|', printDate(t)
       print'|TIME ELAPSED|', printTime((t - Ti))
       print'| DT |','%.2f'%(dt)
       print'|NORM|','%.4f'%(norm)
       print'|CPUT|',printTime(time.clock() - startcpu)
       print(colored('------------------------------------------------------', "red"))
       print''
  if ( i%export  == 0 and t > t_exportStart ):
    exporterSol.exportIdx(i, t)
thomas's avatar
thomas committed
727
728
729
  if ( i%export_extraTerms == 0 ):
    #exporterDissipationTerm.exportIdx(i, t)
    if(gbrWind!=0):
thomas's avatar
thomas committed
730
      windUVDC.L2Projection(windUV)
thomas's avatar
thomas committed
731
      exporterWindUV.exportIdx(i, t)
thomas's avatar
thomas committed
732
      windSourceDC.L2Projection(source)
thomas's avatar
thomas committed
733
734
      exporterWindSourceAcc.exportIdx(i, t)
    if(gbrResidual == 2):
thomas's avatar
thomas committed
735
      blinkSolDC.L2Projection(blink_VelAndElev)
thomas's avatar
thomas committed
736
737
      exporterBlinkSol.exportIdx(i, t)
    if(gbrResidual == 3):
thomas's avatar
thomas committed
738
      boundaryGBROOS_MYRSol_DC.L2Projection(boundaryUVElev_MYR)
thomas's avatar
thomas committed
739
      exporterGBROOS_MYRSol.exportIdx(i, t)
thomas's avatar
thomas committed
740
      boundaryGBROOS_CCHSol_DC.L2Projection(boundaryUVElev_CCH)
thomas's avatar
thomas committed
741
      exporterGBROOS_CCHSol.exportIdx(i, t)
thomas's avatar
thomas committed
742
743
  if ( i%export_fnEval  == 0 ):    
    #Export Function Evaluator results:
744
745
746
747
748
749
    integrator.compute(totalMass, "")
    relMass=totalMass(0, 0)
    if (Msg.GetCommRank() == 0):
      outMass = open(outputDir+'/totalMass.dat',"a")
      outMass.write(str(t/3600-Ti/3600)+'\t'+str(relMass)+'\n')
      outMass.close()
thomas's avatar
thomas committed
750
751
752
753
754
755
756
757
758
    #*** Current evaluator: ***
    for point in range(0,len(evalPointsCart)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      FunctionEval.compute(evalPointsCart[point][0],evalPointsCart[point][1],evalPointsCart[point][2], res)
      #Project velocities onto lonLat space
      res_uv = fullMatrixDouble(1,2)
thomas's avatar
thomas committed
759
      #xyVector(res_uv, evalPointsCart[point][0], evalPointsCart[point][1], res)
thomas's avatar
thomas committed
760
      if (Msg.GetCommRank()+1 == whichPartition):
thomas's avatar
thomas committed
761
762
763
        outEval = open(outputDir+'/Evaluator/evalCurrent_%02d'%(point)+'.dat',"a")
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\t'+str(res.get(0,1))+'\n')
        outEval.close()
thomas's avatar
thomas committed
764
765
    #*** Elevation evaluator - GBROOS: ***
    for point in range(0,len(evalPointsGBROOSLatLon)):
thomas's avatar
thomas committed
766
      #Only write to file if we are on correct cpu:
thomas's avatar
thomas committed
767
      evalPointSP = SPoint3(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2])
thomas's avatar
thomas committed
768
769
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
thomas's avatar
thomas committed
770
771
772
      ElevationEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res)
      res2 = fullMatrixDouble()
      FunctionEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res2)
thomas's avatar
thomas committed
773
      if (Msg.GetCommRank()+1 == whichPartition):
thomas's avatar
thomas committed
774
        outElevEval = open(outputDir+'/Evaluator/evalElevGBROOS_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
thomas's avatar
thomas committed
775
776
        outElevEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\n')
        outElevEval.close()
thomas's avatar
thomas committed
777
        outEval = open(outputDir+'/Evaluator/evalCurrentGBROOS_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
thomas's avatar
thomas committed
778
779
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res2.get(0,0))+'\t'+str(res2.get(0,1))+'\n')
        outEval.close()
thomas's avatar
thomas committed
780
        
thomas's avatar
thomas committed
781
782
783
784
      ##TEST: OUTPUT VALUE OF WIND (PART 2 OF 2) (Doesn't work if using Precomputed functions in wind)
      #if gbrWind != 0:
        #res3 = fullMatrixDouble()
        #WindUVEval.compute(evalPointsGBROOSCart[point][0],evalPointsGBROOSCart[point][1],evalPointsGBROOSCart[point][2], res3)
thomas's avatar
thomas committed
785
        #if (Msg.GetCommRank()+1 == whichPartition):
thomas's avatar
thomas committed
786
787
788
          #outWind = open(outputDir+'/Evaluator/WINDUV_%02d'%(evalPointsGBROOSLatLon[point][3])+'.dat',"a")
          #outWind.write(str(t/3600-Ti/3600)+'\t'+str(res3.get(0,0))+'\t'+str(res3.get(0,1))+'\n')
          #outWind.close()
thomas's avatar
thomas committed
789
        
thomas's avatar
thomas committed
790
791
792
793
794
795
796
    #*** Elevation evaluator - Seafarer: ***
    for point in range(0,len(tidalPointsSeafarerLatLon)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(tidalPointsSeafarerCart[point][0],tidalPointsSeafarerCart[point][1],tidalPointsSeafarerCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      ElevationEval.compute(tidalPointsSeafarerCart[point][0],tidalPointsSeafarerCart[point][1],tidalPointsSeafarerCart[point][2], res)
thomas's avatar
thomas committed
797
      if (Msg.GetCommRank()+1 == whichPartition):
thomas's avatar
thomas committed
798
799
800
801
802
803
804
805
806
807
        outElevEval = open(outputDir+'/Evaluator/evalElevSeafarer_%02d'%(point)+'.dat',"a")
        outElevEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\n')
        outElevEval.close()
    #*** Reef Points: current & elevation evaluator: ***
    for point in range(0,len(reefPointsLatLon)):
      #Only write to file if we are on correct cpu:
      evalPointSP = SPoint3(reefPointsCart[point][0],reefPointsCart[point][1],reefPointsCart[point][2])
      evalPointElem = model.getMeshElementByCoord(evalPointSP)
      whichPartition = evalPointElem.getPartition()
      ElevationEval.compute(reefPointsCart[point][0],reefPointsCart[point][1],reefPointsCart[point][2], res)
thomas's avatar
thomas committed
808
      if (Msg.GetCommRank()+1 == whichPartition):
thomas's avatar
thomas committed
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
        outEval = open(outputDir+'/Evaluator/evalReefPoint_%02d'%(point)+'.dat',"a")
        outEval.write(str(t/3600-Ti/3600)+'\t'+str(res.get(0,0))+'\t'+str(res.get(0,1))+'\t'+str(res.get(0,2))+'\n')
        outEval.close()

endcpu=time.clock()
norm=solution.norm()
Msg.Barrier()
if (Msg.GetCommRank() == 0):
  print ''
  print(colored('********************     End      ********************', "blue"))
  print ''
  print '     - Final Time (Tf):',  printDate(t)
  print '     - Norm of Solution at Tf:',  norm
  print '     - Total CPU Time:', endcpu-startcpu
  print ''
  print(colored('******************************************************', "blue"))
  print ''
Msg.Exit(0)