run2d.py 15.9 KB
Newer Older
1
2
3
#!/usr/bin/env python
# -*- coding: utf-8 -*-

lambrechts's avatar
lambrechts committed
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
##
# @package run2d
# SLIM modelling of Congo River and estuarine.
#
# @author Thomas, Christopher
# @author Vallaeys, Valentin
# @author Le Bars, Yoann
# @author Lambrechts, Jonathan
# @author Delandmeter, Philippe
# @version 1.1
# @date 2013/07/18
# @date 2013/07/19
# @date 2013/07/24
# @date 2013/07/25
# @date 2013/07/27
# @date 2013/11/26
# @date 2014/03/19
# @date 2014/03/21
# @date 2014/03/24
# @date 2014/03/26
# @date 2014/03/27 
# @date 2014/04/03
# @date 2014/04/09
# @date 2014/04/10
# @date 2014/04/14
# @date 2014/04/15
# @date 2014/04/23
# @date 2014/04/24
# @date 2014/04/25
# @date 2014/04/28
# @date 2014/04/29
# @date 2014/05/30
# @date 2014/06/10
# @date 2014/06/30
# @date 2014/07/01
# @date 2014/07/17
40
41
42
# @date 2014/08/04
# @date 2014/08/06
# @date 2014/08/07
43
# @date 2014/09/16
44
# @date 2014/09/26
45
# @date 2014/12/10
46
# @date 2015/01/19
47
# 
lambrechts's avatar
lambrechts committed
48
49
# Congo River and neighbouring ocean modelling. Multirate explicit
# Runge-Kutta time-stepping.
50
# All the parameters can be found and modified in the file "param.py"
lambrechts's avatar
lambrechts committed
51
52
53
54
55
56

# Import for python
from dgpy import *
from math import *
import param
from termcolor import colored
vallaeys's avatar
vallaeys committed
57
import time, os
lambrechts's avatar
lambrechts committed
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
from functions import *
import util

if Msg.GetCommRank() == 0 :
  if not os.path.exists(param.outputDir) :
    os.mkdir(param.outputDir)
Msg.Barrier()

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])

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

# Current iteration time.
t = param.Ti
77
dt = 30
lambrechts's avatar
lambrechts committed
78
79
80
81
82
83
84
85
86

# Gmsh model (i.e. mesh)
model = GModel()
if Msg.GetCommSize() == 1 :
  model.load(param.meshDir + param.fileName + "_utm.msh")
else :
  model.load(param.meshDir + param.fileName + "_utm_part_" + str(Msg.GetCommSize()) + ".msh")

groups = dgGroupCollection(model)
87
groups.splitGroupsByPhysicalTag()
lambrechts's avatar
lambrechts committed
88
89
90

# Grouped parameters.
bath = dgDofContainer(groups, 1)
91
bath.importIdx("Bath/"+param.fileName+"_bath_smooth/"+param.fileName+"_bath_smooth.idx")
lambrechts's avatar
lambrechts committed
92
93
94
95
96
97
98
99
100
101
102
103
bath_PC = functionPrecomputed(groups, 3, 1)
bath_PC.compute(bath.getFunction())
bathGradient_PC = functionPrecomputed(groups, 3, 3)
bathGradient_PC.compute(bath.getFunctionGradient())

# Time format conversion.
timeFunction = functionConstant(t)
timeInitFunction = functionConstant(param.realTi)
julianTimeFunction = functionNumpy(1, epochToJulian, [timeFunction])

# Global conservation solution.
solution = dgDofContainer(groups, 3)
104
start = 0
105
106
if (param.realTi != param.Ti):
  start = int((t-param.realTi)/(dt*param.export))
107
  solution.importIdx(param.outputDir+"/idx/idx-"+str(start).zfill(6)+".idx")
108
109
else:
  solution.setAll(0.)
lambrechts's avatar
lambrechts committed
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

# Projection functions.
[angleProj, UTMtoLatLonDegrees, UTMtoLonLatDegrees] = util.setProj()
# Precomputed angle for projection.
angleProj_PC = functionPrecomputed(groups, 3, 2)
angleProj_PC.compute(angleProj)
# Precomputed conversion from UTM to latitude and longitude.
UTMtoLatLonDegrees_PC = functionPrecomputed(groups, 3, 3)
UTMtoLatLonDegrees_PC.compute(UTMtoLatLonDegrees)
# Precomputed conversion from UTM to longitude and latitude.
UTMtoLonLatDegrees_PC = functionPrecomputed(groups, 3, 3)
UTMtoLonLatDegrees_PC.compute(UTMtoLonLatDegrees)

# ****************** Forcings ******************

# Tide value on each fronteer nodes.
[tideEta, tideUV] = util.setTide(UTMtoLonLatDegrees_PC, angleProj_PC,
                                 timeFunction, julianTimeFunction, bath_PC,
                                 solution.getFunction())
# Merge elevation and tide speed.
tide = functionNumpy(3, merge_ElevAndUV, [tideEta, tideUV])

# ..............................................

# ................OCEANIC CIRCULATION...................

if (param.oceanCirc == "hycom"):
  HYCOMInterpolatorU = slimStructDataInterpolatorMultilinear()
  HYCOMInterpolatorV = slimStructDataInterpolatorMultilinear()
  HYCOMInterpolatorEta = slimStructDataInterpolatorMultilinear()
  HYCOMUContainer = slimStructDataContainerNetcdf(param.hycomDir + "u_2011-04_depthAve_flux.nc","u", "lat", "lon")
  HYCOMVContainer = slimStructDataContainerNetcdf(param.hycomDir + "v_2011-04_depthAve_flux.nc","v", "lat", "lon")
  HYCOMEtaContainer = slimStructDataContainerNetcdf(param.hycomDir + "ssh_2011_fixed.nc","ssh", "lat", "lon")
  HYCOMU = slimFunctionStructData(HYCOMUContainer, HYCOMInterpolatorU, UTMtoLatLonDegrees_PC)
  HYCOMV = slimFunctionStructData(HYCOMVContainer, HYCOMInterpolatorV, UTMtoLatLonDegrees_PC)
  HYCOMEta = slimFunctionStructData(HYCOMEtaContainer, HYCOMInterpolatorEta, UTMtoLatLonDegrees_PC)
  HYCOMU.setTimeFunction(timeFunction)
  HYCOMV.setTimeFunction(timeFunction)
  HYCOMEta.setTimeFunction(timeFunction)

#test = dgDofContainer(groups, 1)
#test.L2Projection(HYCOMU)
#test.exportMsh("HYCOM_U_init")

  HYCOMUV = functionNumpy(3, rotateVector, [HYCOMU, HYCOMV, angleProj_PC])
  minFC = functionConstant(-3.0)
  HYCOMEta_cutOff = functionNumpy(1, cutOffMin, [HYCOMEta, minFC])

#test2 = dgDofContainer(groups, 3)
#test2.L2Projection(HYCOMUV)
#test2.exportMsh("HYCOM_UV_init")

  HYCOMuv = functionNumpy(2, transport2velocity, [HYCOMUV, bath_PC, solution.getFunction()])
  oceanic = functionNumpy(3, merge_ElevAndUV, [HYCOMEta_cutOff, HYCOMuv])

#test = dgDofContainer(groups, 3)
#test.L2Projection(HYCOMbnd)
#test.exportMsh("HYCOM_UVEta_init")
vallaeys's avatar
vallaeys committed
168
169
  if(Msg.GetCommRank() == 0):
    Msg.Info("Using Hycom for global circulation.")
lambrechts's avatar
lambrechts committed
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
else:
  
  oceanic = functionConstant([0., 0., 0.])

#bathyHycomInter = slimStructDataInterpolatorMultilinear()
#bathHycomContainer =  slimStructDataContainerNetcdf(param.hycomDir + "depth_GLBa0_08_09_Congo.nc", "bathymetry", "lat", "lon")
#bathHycom = slimFunctionStructData(bathHycomContainer, bathyHycomInter, UTMtoLatLonDegrees_PC)
#bathHycom.setTimeFunction(timeFunction)
#test.L2Projection(bathHycom)
#test.exportMsh("BATHYCOM")

#def difference(cmap, sol, f1, f2):
#  sol[:] = np.abs(f1 - f2);

#f3 = functionNumpy(1, difference, [bath_PC, bathHycom])
#test.L2Projection(f3)
#test.exportMsh("diffbath")

# .................. Wind     ..................

if (param.windSource == "ncar"):
  # Longitudinal wind speed interpolation.
  windInterpolatorU = slimStructDataInterpolatorMultilinear()
  # Latitudinal wind speed interpolation.
  windInterpolatorV = slimStructDataInterpolatorMultilinear()
  # Raw longitudinal wind speed.
  windUcontainer = slimStructDataContainerNetcdf(param.windFile, "U_GRD_L103", "lat", "lon")
  # Raw latitudinal wind speed.
  windVcontainer = slimStructDataContainerNetcdf(param.windFile, "V_GRD_L103", "lat", "lon")
  # Longitudinal wind speed interpolated on model mesh.
  windU = slimFunctionStructData(windUcontainer, windInterpolatorU, UTMtoLatLonDegrees_PC)
  # Latitudinal wind speed interpolated on model mesh.
  windV = slimFunctionStructData(windVcontainer, windInterpolatorV, UTMtoLatLonDegrees_PC)
  windU.setTimeFunction(timeFunction)
  windV.setTimeFunction(timeFunction)
  # Wind speed vector.
  windUV = functionNumpy(3, rotateVector, [windU, windV, angleProj_PC])
else:
  windUV = functionConstant([0., 0., 0.])
# ..............................................

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


# ****************** Cons Law ******************

# Conservation law.
claw = dgConservationLawShallowWater2d()
218
claw.setIsLinear(False)
lambrechts's avatar
lambrechts committed
219
220
221

# Bottom friction.
Cd = functionNumpy(1, bottomDrag_const_n0235, [solution.getFunction(), bath_PC])
222
223
224
# Area of each element.
area = dgDofContainer(groups, 1)
area.computeArea()
lambrechts's avatar
lambrechts committed
225
# Diffusion (Smagorinsky)
226
Di = functionNumpy(1, smagorinskyTerm, [solution.getFunctionGradient(), area.getFunction()])
lambrechts's avatar
lambrechts committed
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
# Coriolis force.
Coriolis = functionNumpy(1, coriolisLonLatDegrees, [UTMtoLonLatDegrees_PC])
Coriolis_PC = functionPrecomputed(groups, 3, 1)
Coriolis_PC.compute(Coriolis)
# Wind stress from netcdf
windSource = functionNumpy(2, windStressSmithBanke, [windUV, solution.getFunction(), bath_PC])

claw.setCoriolisFactor(Coriolis_PC)
claw.setQuadraticDissipation(Cd)
claw.setDiffusivity(Di)
claw.setBathymetry(bath_PC)
claw.setBathymetryGradient(bathGradient_PC)

# Slope force.
ggradz = functionNumpy(2, ggradz, [UTMtoLonLatDegrees_PC])
ggradz_PC = functionPrecomputed(groups, 3, 2)
ggradz_PC.compute(ggradz)
source = functionNumpy(2, merge2D, [windSource, ggradz_PC])
claw.setSource(source)
# River discharge.
if (param.flowSource == "approx"):
  riverDischargeFunction = slimTemporalSerie(param.flowFile)
else:
  riverDischargeFunction = functionConstant(param.flowConstant)
riverDischargeFunctionSmooth = functionNumpy(1, rampFunctionSingle, [timeInitFunction, timeFunction, riverDischargeFunction])

# Flux integration over the section.
bathIntegralFM = fullMatrixDouble(1, 1)
dgFunctionIntegratorInterface(groups,  bath_PC).compute('source1',  bathIntegralFM)
# Flux integration type.
bathIntegralsFunc = functionConstant(bathIntegralFM(0, 0))
# Current normal speed vector.
currentUV = functionNumpy(3, fluxToVelocity, [function.getNormals(), bathIntegralsFunc,  riverDischargeFunctionSmooth])
260
boundarySmooth = functionNumpy(3, rampFunction, [timeInitFunction, timeFunction, tide])
lambrechts's avatar
lambrechts committed
261

262
if (param.nudging == True):
vallaeys's avatar
vallaeys committed
263
264
265
266
267
  coefNud = functionNumpy(2, coefAbs, [UTMtoLonLatDegrees])
  coefNud_PC = functionPrecomputed(groups, 3, 2)
  coefNud_PC.compute(coefNud)
  groups.exportFunctionMsh(coefNud,"t")
  claw.setAbsLayer('sea', coefNud_PC, oceanic)
lambrechts's avatar
lambrechts committed
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
claw.addBoundaryCondition('ocean',claw.newOutsideValueBoundary("Surface", boundarySmooth))
claw.addBoundaryCondition('source1', claw.newOutsideValueBoundary("Surface", currentUV))
claw.addBoundaryCondition('coast',claw.newBoundaryWall())
claw.addBoundaryCondition('removedElements',claw.newBoundaryWall())

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

# ****************** Time Int ******************

#temporal iterator
rk = 0
# Time step in seconds.
dt = 0

if param.temporalScheme == "implicit" :
  sys = linearSystemPETScBlockDouble()
  sys.setParameter("petscOptions", "-ksp_rtol 1e-3")
  dof = dgDofManager.newDGBlock(groups, 3, sys)
  rk = dgDIRK(claw, dof, 2)
  #rk.getNewton().setVerb(10)
  rk.getNewton().setRtol(1e-5)
  dt = 900
elif param.temporalScheme == "explicit" :
  rk = dgERK(claw, None, DG_ERK_44)
  dt = 2
elif param.temporalScheme == "multirate" :
  rk = dgMultirateERK(groups, claw, param.multirateScheme)
295
  dt = rk.splitGroupsForMultirate(param.nMRLevels, solution, [solution, bath, area]) 
lambrechts's avatar
lambrechts committed
296
297
298
299
300
301
302
303
304
305
306
307
308
  dt = 30
  rk.printMultirateInfo(0)
else :
  Msg.Fatal("unkown temporal scheme : %s" % param.temporalScheme)

# groups splitted, need to recompute...
angleProj_PC.compute(angleProj)
UTMtoLatLonDegrees_PC.compute(UTMtoLatLonDegrees)
UTMtoLonLatDegrees_PC.compute(UTMtoLonLatDegrees)
ggradz_PC.compute(ggradz)
bath_PC.compute(bath.getFunction())
bathGradient_PC.compute(bath.getFunctionGradient())
Coriolis_PC.compute(Coriolis)
309
310
groups.splitGroupsByPhysicalTag()
if (param.nudging == True):
vallaeys's avatar
vallaeys committed
311
  coefNud_PC.compute(coefNud)
312

lambrechts's avatar
lambrechts committed
313
314
315
316
317
318
319
320

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

# Number of time steps.
nbSteps = int(ceil((param.Tf-param.Ti)/dt))
# Number of export made.
nbExport = 0

vallaeys's avatar
vallaeys committed
321
if (param.outputIdx):
322
  if(Msg.GetCommRank() == 0 and param.realTi == param.Ti):
vallaeys's avatar
vallaeys committed
323
324
    if(os.path.exists(param.outputDir+'/idx')):
      try : os.system("rm -r "+param.outputDir+'/idx');
lambrechts's avatar
lambrechts committed
325
      except: 0;
vallaeys's avatar
vallaeys committed
326
327
  # Exporting object to export solution in idx.
  exporterIdx = dgIdxExporter(solution, param.outputDir+'/idx')
328
329
330
  if (param.realTi == param.Ti):
    exporterIdx.exportIdx(0, t)
  if(Msg.GetCommRank() == 0 and param.realTi == param.Ti):
vallaeys's avatar
vallaeys committed
331
332
    Msg.Info("Export initial solution in idx format... done.")
if (param.outputNetCDF):
333
  if(Msg.GetCommSize() == 1):
vallaeys's avatar
vallaeys committed
334
335
    if(os.path.exists(param.outputDir+'/netCDF')):
      try : os.system("rm "+param.outputDir+'/netCDF/*');
lambrechts's avatar
lambrechts committed
336
337
      except: 0;
    else:
vallaeys's avatar
vallaeys committed
338
339
340
341
342
      os.mkdir(param.outputDir+'/netCDF')
    # Exporting object to export solution in netCDF.
    exporterNetCDF = dgExportNetCDF(groups, param.proj, param.latMin, param.latMax, param.lonMin, param.lonMax, param.dLat, param.dLon)
    exporterNetCDF.exportNetCDF(solution, param.outputDir+'/netCDF/out000000.nc', t/3600.)
    Msg.Info("Export initial solution in netCDF format... done.")
343
  else:
344
    Msg.Info("Exports in netCDF do not work with multi-processor runs!") 
vallaeys's avatar
vallaeys committed
345
346
if (not param.outputNetCDF and not param.outputIdx):
  Msg.Fatal("Unknown output format, no output created !")
lambrechts's avatar
lambrechts committed
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372

# SLIM function which computes total water mass.
mass=functionNumpy(1,  massIntegrator,  [function.getSolution(), bath_PC])
# Mass integration type.
integrator = dgFunctionIntegrator(groups, mass, solution)
# Current total water mass.
totalMass=fullMatrixDouble(1, 1)
# Initial total water mass.
totalMassRef=fullMatrixDouble(1, 1)

if (Msg.GetCommRank() == 0):
  print( '')
  print((colored('******************** Initializing ********************', "blue")))
  print( '')
  print( '     - Number of processors:', Msg.GetCommSize())
  print( '     - Multirate Runge Kutta Scheme:', param.multirateScheme)
 # print( '     - Multirate Speedup:', rk.speedUp())
  print( '     - Initial Time (Ti):',  printDate(param.Ti))
  print( '     - Final Time (Tf):',  printDate(param.Tf))
  print( '     - Reference Time Step:', printTime(dt))
  print( '     - Number of Iterations:', nbSteps)
  print( '     - Mass of Solution at Ti:',  totalMassRef(0, 0))
  print( '')
  print((colored('******************************************************', "blue")))
  print( '')

373
374
375
376
377
378
379
#Write info file for LPTracker
if (Msg.GetCommRank() == 0):
  simInfo = [str(dt)+'\n', str(nbSteps)+'\n', str(param.export)+'\n', str(param.Ti)+'\n', str(param.outputDir)+'\n', str(param.fileName + "_utm.msh")]
  simInfoFile = open(param.outputDir+'/simulationInfo.txt','w')
  simInfoFile.writelines(simInfo)
  simInfoFile.close() 

lambrechts's avatar
lambrechts committed
380
# ****************** Iterate  ******************
vallaeys's avatar
vallaeys committed
381
relMass = 0.
lambrechts's avatar
lambrechts committed
382
383
# System time when starting computations.
startcpu = time.clock()
384
385
startI = start * param.export
for i in range(startI, startI + nbSteps):
lambrechts's avatar
lambrechts committed
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
  if(i == nbSteps - 1):
    dt = param.Tf - t
  rk.iterate (solution, dt, t)
  t = t + dt
  timeFunction.set(t)
  if ((i+1) % param.display == 0):
    integrator.compute(totalMass, "")
    # Current total mass as a scalar value.
    relMass=totalMass(0, 0)
    if(Msg.GetCommRank() == 0):
      print(colored('|ITER| %d',  "red")%(i+1))
      print(colored('------------------------------------------------------', "red"))
      print('|TIME|', printDate(t))
      print('|MASS|',relMass)
      print('|CPUT|',printTime(time.clock() - startcpu))
      print(colored('------------------------------------------------------', "red"))
      print('')
  if ( (i+1)%param.export  == 0 ):
vallaeys's avatar
vallaeys committed
404
405
406
407
408
    if (param.outputIdx):
      exporterIdx.exportIdx((i+1) // param.export, t)
      if(Msg.GetCommRank() == 0):
        Msg.Info("Export solution in idx format... done.")
    if (param.outputNetCDF):
409
410
411
412
413
      if(Msg.GetCommSize() == 1):
        # Number of hours since 1970/01/01 00:00
        hours = t/3600.
        exporterNetCDF.exportNetCDF(solution, param.outputDir+'/netCDF/out'+str((i+1) // param.export).zfill(6)+".nc", hours)
        Msg.Info("Export solution in netCDF format... done.")
lambrechts's avatar
lambrechts committed
414
    nbExport = nbExport + 1
415
  if (isnan(relMass)):
416
    Msg.Fatal("NaN detected !")  
lambrechts's avatar
lambrechts committed
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
# System time when stopping computations.
endcpu=time.clock()
integrator.compute(totalMass, "")
relMass=totalMass(0, 0)
Msg.Barrier()
if (Msg.GetCommRank() == 0):
  print( '')
  print(colored('********************     End      ********************', "blue"))
  print( '')
  print( '     - Final Time (param.Tf):',  printDate(t))
  print( '     - Mass of Solution at Tf:', relMass)
  print( '     - Total CPU Time:', endcpu-startcpu)
  print( '')
  print(colored('******************************************************', "blue"))
  print( '')
Msg.Exit(0)