MultirateGBR.py 7.57 KB
Newer Older
bseny's avatar
bseny committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
"""
    Simulations on the GBR with multirate explicit Runge-Kutta time-stepping.
    All the parameters can be found and modified in the file "InpuGBR.py"
"""

"""
    *** How to launch the testcase ***

    Simulation on 1 processor  : rundgpy MultirateGBR.py
    Simulation on n processors : mpirun -np n rundgpy MultirateGBR.py   ( !!! You need MPI !!! )

    **********************************
"""

# Import for python
from dgpy import *
from math import *
from ProjectMesh import *
from InputGBR import *
from termcolor import colored
import time, os, sys

model = GModel()
bseny's avatar
update    
bseny committed
24
25
model.load(filename + '_partitioned_tan%d.msh'%(Msg.GetCommSize()))
#loadAndPartitionMesh(model,  filename,  mL, algo)
bseny's avatar
bseny committed
26
27
28
29
30
31
32
33
34
35
36

# Coordinate system
XYZ = function.getCoordinates();
lonLatDegrees = functionC(glib,"lonLatDegrees",3,[XYZ])
lonLat = functionC(glib,"lonLat",3,[XYZ])

latLonDegrees = functionC(glib,"latLonDegrees",3,[XYZ])
latLon = functionC(glib,"latLon",3,[XYZ])

groups = dgGroupCollection(model, dimension, order)

bseny's avatar
bseny committed
37
38
if(not os.path.exists(bathname)):
  Msg.Fatal('No valid bathymetry file :' + bathname + ' found,  please launch "rundgpy DiffuseBathymetryTsunami.py"')
bseny's avatar
bseny committed
39

bseny's avatar
bseny committed
40
bath = dgDofContainer(groups, 1);
41
bath.importIdx(bathname)
bseny's avatar
bseny committed
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

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

t = Ti
timeFunction = functionConstant(t)

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

# .................. Tides    ..................

# (a) Sinus tide
def tide (f):
  for i in range(0,f.size1()):
    f.set(i, 0, sin((t-Ti)*2*pi/3600/12))
    f.set(i, 1, 0)
    f.set(i, 2, 0)

# (b) Tpxo tide 
bseny's avatar
update    
bseny committed
60
61
62
63
64
65
66
tideEta = slimFunctionTpxo("h_tpxo7.2.nc","ha","hp","lon_z","lat_z")
tideEta.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideU = slimFunctionTpxo("u_tpxo7.2.nc","Ua","up","lon_u","lat_u")
tideU.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideV = slimFunctionTpxo("u_tpxo7.2.nc","Va","vp","lon_v","lat_v")
tideV.setCoordAndTimeFunctions(lonLatDegrees,timeFunction)
tideUV = functionC(glib,"latLonVector",3,[latLon,tideU,tideV])
bseny's avatar
bseny committed
67
68
69
70
71

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

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

bseny's avatar
update    
bseny committed
72
73
74
75
76
77
78
79
windInterpolator = slimStructDataInterpolatorMultilinear()
windUcontainer = slimStructDataContainerNetcdf("uwind-2000-2011.nc", "uwnd", "lat", "lon")
windVcontainer = slimStructDataContainerNetcdf("vwind-2000-2011.nc", "vwnd", "lat", "lon")
windU = slimFunctionStructData(windUcontainer, windInterpolator, latLonDegrees)
windV = slimFunctionStructData(windVcontainer, windInterpolator, latLonDegrees)
windU.setTimeFunction(timeFunction)
windV.setTimeFunction(timeFunction)
windUV = functionC(glib,"latLonVector",3,[latLon,windU,windV])
bseny's avatar
bseny committed
80
81
82
83
84
85
86
87
88
89
90
91

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

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


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

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

ld = functionConstant(0)
bseny's avatar
bseny committed
92
Cd = functionC(glib, "bottomDrag", 1, [function.getSolution(), bath.getFunction()])
bseny's avatar
bseny committed
93
94
95
96
97
98
99
100
if(gbrDiff == 0):
  Di = functionConstant(diff)
elif(gbrDiff == 1):
  Di = functionC(glib, "smagorinsky", 1, [function.getSolutionGradient()])
Coriolis = functionC(glib, "coriolis", 1, [XYZ])
if(gbrWind == 0):
  source = functionC(glib, "wind", 2, [XYZ])
elif(gbrWind == 1):
bseny's avatar
bseny committed
101
  source = functionC(glib, "windStress", 2, [windUV, function.getSolution(), bath.getFunction()])
bseny's avatar
bseny committed
102
103
104
105
106
107
108


claw.setCoriolisFactor(Coriolis)
claw.setQuadraticDissipation(Cd)
claw.setLinearDissipation(ld)
claw.setDiffusivity(Di)
claw.setSource(source)
bseny's avatar
bseny committed
109
110
claw.setBathymetry(bath.getFunction())
claw.setBathymetryGradient(bath.getFunctionGradient())
bseny's avatar
bseny committed
111
112

Os = functionPython(3, tide,[])
bseny's avatar
update    
bseny committed
113
114
tideuv = functionC(glib,"transport2velocity",3,[tideUV,bath.getFunction(),function.getSolution()]);
tpxoTide = functionC(glib,"merge",3,[tideEta,tideuv]);
bseny's avatar
bseny committed
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131

if(gbrTide == 0):
  claw.addBoundaryCondition('Open Sea North',claw.newOutsideValueBoundary("Surface", Os))
  claw.addBoundaryCondition('Open Sea Central',claw.newOutsideValueBoundary("Surface", Os))
  claw.addBoundaryCondition('Open Sea South',claw.newOutsideValueBoundary("Surface", Os))
elif(gbrTide == 1):
  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))

claw.addBoundaryCondition('Coast',claw.newBoundaryWall())
claw.addBoundaryCondition('Islands',claw.newBoundaryWall())

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


# ****************** Time Int ******************
bseny's avatar
bseny committed
132
133
134
rk = dgMultirateERK(groups, claw, RK_TYPE)
dt = rk.splitGroupsForMultirate(mL, solution, [solution, bath])
rk.printMultirateInfo(0)
bseny's avatar
bseny committed
135
136
137
138

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

currentFunction = functionC(glib,"current",3,[function.getSolution()])
bseny's avatar
bseny committed
139
#tideCurrentFunction = functionC(glib,"current",3,[tpxoTide])
bseny's avatar
bseny committed
140
141
142
143
144

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

if(Msg.GetCommRank() == 0):
bseny's avatar
bseny committed
145
146
  if(os.path.exists(outputDir+'/mr-gbr')):
    try : os.system("rm -r "+outputDir+'/mr-gbr');
bseny's avatar
bseny committed
147
    except: 0;
bseny's avatar
bseny committed
148
exporterSol=dgIdxExporter(solution, outputDir+'/mr-gbr')
bseny's avatar
bseny committed
149
150
151

norm=solution.norm()

bseny's avatar
bseny committed
152
mass=functionPython(1,  massIntegrator,  [functionSolution.get(), bath.getFunction()])
bseny's avatar
bseny committed
153
154
155
156
157
158
159
160
161
162
integrator = dgFunctionIntegrator(groups,  mass,  solution)
totalMass=fullMatrixDouble(1, 1)
totalMassRef=fullMatrixDouble(1, 1)
integrator.compute(totalMassRef, "")

if (Msg.GetCommRank() == 0):
  print ''
  print(colored('******************** Initializing ********************', "blue"))
  print ''
  print '     - Number of processors:', Msg.GetCommSize()
bseny's avatar
bseny committed
163
  print '     - Multirate Runge Kutta Scheme:', RK_TYPE
bseny's avatar
bseny committed
164
165
166
167
  print '     - Initial Time (Ti):',  printDate(Ti)
  print '     - Final Time (Tf):',  printDate(Tf)
  print '     - Reference Time Step:', printTime(dt)
  print '     - Number of Iterartions:', nbSteps
168
  print '     - Theoretical Speedup=', rk.speedUp()
bseny's avatar
bseny committed
169
170
171
172
173
174
175
176
177
178
179
  print '     - Norm of Solution at Ti:',  norm
  print '     - Mass of Solution at Ti:',  totalMassRef(0, 0)
  print ''
  print(colored('******************************************************', "blue"))
  print ''


# ****************** Iterate  ******************


startcpu = time.clock()
bseny's avatar
bseny committed
180
181
for i in range(0, nbSteps):
  if(i == nbSteps - 1):
bseny's avatar
bseny committed
182
    dt=Tf-t
bseny's avatar
bseny committed
183
  norm = rk.iterate (solution, dt, t)
bseny's avatar
bseny committed
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
  t = t +dt
  timeFunction.set(t)
  if ( i%iter == 0 ):
     norm=solution.norm()
     integrator.compute(totalMass, "")
     relMass=totalMass(0, 0)
     if(Msg.GetCommRank() == 0):
       print(colored('|ITER| %d',  "red")%(i))
       print(colored('------------------------------------------------------', "red"))
       print'|TIME|', printDate(t)
       print'|NORM|',norm
       print'|MASS|',relMass 
       print'|CPUT|',printTime(time.clock() - startcpu)
       print(colored('------------------------------------------------------', "red"))
       print''
  if ( i%export  == 0 ):
bseny's avatar
update    
bseny committed
200
201
202
    #exporterSol.exportIdx(i, t)
    solution.exportMsh(outputScratch+"mr_gbr-%06d"% (i),t, nbExport)
    solution.exportFunctionMsh(currentFunction, outputScratch+'mr_gbr_current-%06d'%(i), t, nbExport, '')
bseny's avatar
bseny committed
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
    nbExport = nbExport + 1

endcpu=time.clock()
norm=solution.norm()
integrator.compute(totalMass, "")
relMass=totalMass(0, 0)
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 '     - Mass of Solution at Tf:', relMass
  print '     - Total CPU Time:', endcpu-startcpu
  print ''
  print(colored('******************************************************', "blue"))
  print ''
Msg.Exit(0)