← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/convection.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/convection.py
 
1
########################################################
 
2
#
 
3
# Copyright (c) 2003-2010 by University of Queensland
 
4
# Earth Systems Science Computational Center (ESSCC)
 
5
# http://www.uq.edu.au/esscc
 
6
#
 
7
# Primary Business: Queensland, Australia
 
8
# Licensed under the Open Software License version 3.0
 
9
# http://www.opensource.org/licenses/osl-3.0.php
 
10
#
 
11
########################################################
 
12
"""
 
13
this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
 
14
 
 
15
It is solved in dimensionless form
 
16
 
 
17
"""
 
18
__copyright__="""Copyright (c) 2003-2010 by University of Queensland
 
19
Earth Systems Science Computational Center (ESSCC)
 
20
http://www.uq.edu.au/esscc
 
21
Primary Business: Queensland, Australia"""
 
22
__license__="""Licensed under the Open Software License version 3.0
 
23
http://www.opensource.org/licenses/osl-3.0.php"""
 
24
__url__="https://launchpad.net/escript-finley"
 
25
 
 
26
from esys.escript import *
 
27
from esys.escript.models import TemperatureCartesian, IncompressibleIsotropicFlowCartesian, Mountains, SubSteppingException
 
28
from esys.dudley import Rectangle, Brick, LoadMesh
 
29
from optparse import OptionParser
 
30
from math import pi, ceil
 
31
import sys
 
32
import time
 
33
 
 
34
# ======================= Default Values ==================================================
 
35
DIM=2                           # spatial dimension
 
36
H=1.                            # height
 
37
L=2*H                           # length
 
38
NE=30                           # number of elements in H-direction. 
 
39
PERT=0.15               # initial temperature perturbation
 
40
DT=1.e-4                        # initial time step size
 
41
CREATE_TOPOGRAPHY=False         # create topgraphy
 
42
DT_MIN=1.e-10                    # minumum time step size
 
43
T_END=10.                       # end time
 
44
 
 
45
RHO_0=100.                     # surface density  (lf ~ RHO_0**2) 
 
46
G=1.                            # gravitational constant
 
47
ALPHA_0=0.1                    # thermal expansion coefficient (Ra ~ RHO_0**2 * ALPHA_0 = lf * ALPHA_0)
 
48
T_0=0.                          # surface temperature
 
49
T_1=1.                          # bottom temperature
 
50
C_P=1                           # heat capacity 
 
51
K=1.                            # thermal conductivity
 
52
CHI=0.                          # Taylor Quinny coefficient
 
53
MUE=None                        # elastic shear modulus
 
54
TAU_Y=5*10**(2.5)               # Drucker-Prager cohesion factor
 
55
BETA=0                          # Drucker-Prager friction factor
 
56
TAU_0=2*10**(2.5)               # transition stress
 
57
N=3                             # power for power law
 
58
 
 
59
E=23*0                      # activation energy
 
60
V=18*0                        # activation volume 
 
61
T_OFFSET=1                      # temperature offset on surface (dimensionless formulation T_OFFSET=1 otherwise =0)
 
62
R=1                             # gas constant
 
63
ETA_N0=1.                       # viscosity at surface 
 
64
 
 
65
TOPO_SMOOTH=1e-5                # smoothing factor of extrapolation of surface velocity to interior
 
66
 
 
67
T_TOL=1.e-4                     # tolerance temperature transport
 
68
FLOW_TOL=1.e-3                  # tolerance for inconcompressible flow solver
 
69
TOPO_TOL=0.1                 # tolerance for update of topography
 
70
DIAGNOSTICS_FN="diagnostics.csv"
 
71
VERBOSE=True
 
72
DT_VIS=T_END/500                # time distane between two visulaization files
 
73
DN_VIS=5                        # maximum counter increment between two visulaization files
 
74
VIS_DIR="results"               # name of the director for vis files
 
75
DN_RESTART=1000000              # create a restart file every DN_RESTART step
 
76
PREFIX_RESTART="restart"        # name prefix for restart directories
 
77
TOPO_ITER_MAX=20                # maximum number of iteration steps to update topography
 
78
# =========================================================================================
 
79
 
 
80
def removeRestartDirectory(dir_name):
 
81
   if dom.onMasterProcessor() and os.path.isdir(dir_name):
 
82
       for root, dirs, files in os.walk(dir_name, topdown=False):
 
83
           for name in files: os.remove(os.path.join(root,name))
 
84
           for name in dirs: os.remove(os.path.join(root,name))
 
85
       os.rmdir(dir_name)
 
86
       print "Restart files %s have been removed."%dir_name
 
87
   dom.MPIBarrier()
 
88
 
 
89
# =========================================================================================
 
90
#
 
91
# read options:
 
92
#
 
93
parser = OptionParser(usage="%prog [Options]")
 
94
parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true")
 
95
parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None)
 
96
parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None)
 
97
(options, args) = parser.parse_args()
 
98
restart=options.restart or (options.restart_dir !=None)
 
99
#
 
100
#  overwrite the default options:
 
101
#
 
102
print "<%s> Execution started."%time.asctime()
 
103
if options.param !=None: 
 
104
     exec(open(options.param,'r'))
 
105
     print "Parameters are imported from file ",options.param
 
106
 
 
107
print "Input Parameters:"
 
108
print "\tDimension                     DIM\t\t=\t",DIM
 
109
print "\tHeight                        H\t\t\t=\t",H
 
110
print "\tLength                        L\t\t\t=\t",L
 
111
print "\telements in H                 NE\t\t=\t",NE
 
112
print "\ttemperature perturbation      PERT\t\t=\t",PERT
 
113
print "\tinitial time step size        DT\t\t=\t",DT
 
114
print "\tminimum time step size        DT_MIN\t\t=\t",DT_MIN
 
115
print "\tend time                      T_END\t\t=\t",T_END
 
116
print "\tcreate topgraphy              CREATE_TOPOGRAPHY\t=\t",CREATE_TOPOGRAPHY
 
117
print "\tsurface density               RHO_0\t\t=\t",RHO_0
 
118
print "\tgravitational constant        G\t\t\t=\t",G
 
119
print "\tthermal expansion coefficient ALPHA_0\t\t=\t",ALPHA_0
 
120
print "\tsurface temperature           T_0\t\t=\t",T_0
 
121
print "\tbottom temperature            T_1\t\t=\t",T_1
 
122
print "\theat capacity                 C_P\t\t=\t",C_P
 
123
print "\tthermal conductivity          K\t\t\t=\t",K
 
124
print "\tTaylor-Qinny coefficient      CHI\t\t=\t",CHI
 
125
print "\telastic shear modulus         MUE\t\t=\t",MUE
 
126
print "\tcohesion factor               TAU_Y\t\t=\t",TAU_Y
 
127
print "\tfriction factor               BETA\t\t=\t",BETA
 
128
print "\ttransition stress             TAU_0\t\t=\t",TAU_0
 
129
print "\tpower for power law           N\t\t\t=\t",N
 
130
print "\tviscosity at surface          ETA_N0\t\t=\t",ETA_N0
 
131
print "\tactivation energy             E\t\t\t=\t",E
 
132
print "\tactivation volume             V\t\t\t=\t",V
 
133
print "\ttemperature offset            T_OFFSET\t\t=\t",T_OFFSET
 
134
print "\tgas constant                  R\t\t\t=\t",R
 
135
print "\ttopography smoothing          TOPO_SMOOTH\t=\t",TOPO_SMOOTH
 
136
 
 
137
print "\ttolerance for topography      TOPO_TOL\t\t=\t",TOPO_TOL
 
138
print "\ttransport tolerance           T_TOL\t\t=\t",T_TOL
 
139
print "\tflow tolerance                FLOW_TOL\t\t=\t",FLOW_TOL
 
140
print "\tfile for diagnostics          DIAGNOSTICS_FN\t=\t",DIAGNOSTICS_FN
 
141
print "\tmin. time incr. for vis file  DT_VIS\t\t=\t",DT_VIS
 
142
print "\tmin. count incr. for vis file DN_VIS\t\t=\t",DN_VIS
 
143
print "\tvisualization dir             VIS_DIR\t\t=\t",VIS_DIR
 
144
print "\trestart counter incerement    DN_RESTART\t=\t",DN_RESTART
 
145
print "\tprefix for restart dirs       PREFIX_RESTART\t=\t",PREFIX_RESTART
 
146
print "\tverbosity                     VERBOSE\t\t=\t",VERBOSE
 
147
 
 
148
print "Control Parameters:"
 
149
t_REF=RHO_0*C_P*H**2/K
 
150
P_REF=ETA_N0/t_REF
 
151
Ar=E/R/(T_1-T_0)
 
152
if E>0 and V>0:
 
153
   V_REF=P_REF*V/E
 
154
else:
 
155
  V_REF=0
 
156
T_OFFSET_REF=T_OFFSET/(T_1-T_0)
 
157
Ra=RHO_0*G*H*(T_1-T_0)*ALPHA_0/P_REF
 
158
Di=ALPHA_0*G*H/C_P
 
159
CHI_REF=CHI*K*ETA_N0/(RHO_0**2*C_P**2*(T_1-T_0)*H**2)
 
160
if CREATE_TOPOGRAPHY:
 
161
   SURFACE_LOAD=RHO_0*G*H/P_REF
 
162
else:
 
163
   SURFACE_LOAD=0.
 
164
if MUE == None:
 
165
  De=None
 
166
else:
 
167
  De=ETA_N0/MUE/t_REF
 
168
ETA_BOT=exp(Ar*((1.+V_REF)/(T_OFFSET_REF+1)-1./T_OFFSET_REF))*ETA_N0
 
169
print "\ttotal #element                \t\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
 
170
print "\treference time                t_REF\t\t=\t%e"%t_REF
 
171
print "\treference pressure            P_REF\t\t=\t%e"%P_REF
 
172
print "\treference Taylor-Qinny        CHI_REF\t\t=\t%e"%CHI_REF
 
173
print "\tDissipation number            DI\t\t=\t%e"%Di
 
174
print "\tRayleigh number surface       Ra\t\t=\t%e"%Ra
 
175
if MUE == None:
 
176
   print "\tDebora number surface         De\t\t=\t",None
 
177
else:
 
178
   print "\tDebora number surface         De\t\t=\t%e"%De
 
179
     
 
180
print "\tBottom viscosity              \t\t\t=\t%e"%ETA_BOT
 
181
print "\tRayleigh number bottom        \t\t\t=\t%e"%(RHO_0*G*H*(T_1-T_0)*ALPHA_0*t_REF/ETA_BOT)
 
182
if MUE == None:
 
183
   print "\tDebora number bottom          \t\t\t=\t",None
 
184
else:
 
185
   print "\tDebora number bottom          \t\t\t=\t%d"%(ETA_BOT/MUE/t_REF)
 
186
print "\tArrhenius                     Ar\t\t=\t%e"%Ar
 
187
print "\tsurface load factor           SURFACE_LOAD\t=\t%e"%SURFACE_LOAD
 
188
print "\tscaled activation volume      V_REF\t\t=\t%e"%V_REF
 
189
#  some control variables (will be overwritten in  case of a restart:
 
190
#
 
191
t=0         # time stamp
 
192
n=0         # time step counter
 
193
dt=DT       # current time step size
 
194
t_vis=0     
 
195
n_vis=0
 
196
counter_vis=0
 
197
mkDir(VIS_DIR)
 
198
#=========================
 
199
#
 
200
#   set up domain or read restart file
 
201
#
 
202
if restart:
 
203
   if options.restart_dir ==None:
 
204
      restart_files=[]
 
205
      for f in os.listdir("."):
 
206
          if f.startswith(PREFIX_RESTART): restart_files.append(f)
 
207
      if len(restart_files)==0:
 
208
          raise IOError,"no restart files"
 
209
      restart_files.sort()
 
210
      f=restart_files[-1]
 
211
   else:
 
212
       f=options.restart_dir
 
213
   try:
 
214
      dom=LoadMesh("mesh.nc")
 
215
   except:
 
216
      pass
 
217
   FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
 
218
   t=float(FF[0])          # time stamp
 
219
   n=int(FF[3])            # time step counter
 
220
   t_vis=float(FF[1])      # 
 
221
   n_vis=int(FF[2])
 
222
   counter_vis=int(FF[4])
 
223
   dt=float(FF[5])
 
224
   stress=load(os.path.join(f,"stress.nc"),dom)
 
225
   v=load(os.path.join(f,"v.nc"),dom)
 
226
   p=load(os.path.join(f,"p.nc"),dom)
 
227
   T=load(os.path.join(f,"T.nc"),dom)
 
228
   if CREATE_TOPOGRAPHY:
 
229
       topography=load(os.path.join(f,"topo.nc"),dom)
 
230
   
 
231
   diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=True)
 
232
   print "<%s> Restart from file %s at time step %s (t=%s) completed."%(time.asctime(),f,t,n)
 
233
else:
 
234
  if DIM==2:
 
235
    dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L/H,l1=1,order=-1,optimize=True)
 
236
  else:
 
237
    dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L/H,l1=L/H,l2=1,order=-1,optimize=True)
 
238
  try:
 
239
     dom.dump("mesh.nc")
 
240
  except:
 
241
     pass
 
242
  x=dom.getX()
 
243
  T=Scalar(1,Solution(dom))
 
244
  for d in range(DIM):
 
245
      if d == DIM-1: 
 
246
         T*=sin(x[d]/H*pi)
 
247
      else:
 
248
         T*=cos(x[d]/L*pi)
 
249
 
 
250
  T=(1.-x[DIM-1])+PERT*T
 
251
  v=Vector(0,Solution(dom))
 
252
  stress=Tensor(0,Function(dom))
 
253
  x2=ReducedSolution(dom).getX()
 
254
  p=Ra*(x2[DIM-1]-0.5*x2[DIM-1]**2-0.5)
 
255
 
 
256
  if CREATE_TOPOGRAPHY:
 
257
      topography=Scalar(0.,Solution(dom))
 
258
  diagnostics_file=FileWriter(DIAGNOSTICS_FN,append=False)
 
259
  diagnostics_file.write("Ra = %e Lambda= %e\n"%(Ra, SURFACE_LOAD))
 
260
 
 
261
p_last=p
 
262
x=dom.getX()
 
263
#
 
264
#   set up heat problem:
 
265
#
 
266
heat=TemperatureCartesian(dom,useBackwardEuler=False)
 
267
print "<%s> Temperature transport has been set up."%time.asctime()
 
268
heat.getSolverOptions().setTolerance(T_TOL)
 
269
heat.getSolverOptions().setVerbosity(VERBOSE)
 
270
fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
 
271
heat.setInitialTemperature(clip(T,T_0))
 
272
heat.setValue(rhocp=1,k=1,given_T_mask=fixed_T_at)
 
273
#
 
274
#   velocity constraints:
 
275
#
 
276
fixed_v_mask=Vector(0,Solution(dom))
 
277
faces=Scalar(0.,Solution(dom))
 
278
for d in range(DIM):
 
279
    if d == DIM-1: 
 
280
       ll = H
 
281
    else:
 
282
       ll = L
 
283
    if CREATE_TOPOGRAPHY and d==DIM-1:
 
284
       fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
 
285
    else:
 
286
       s=whereZero(x[d])+whereZero(x[d]-ll)
 
287
       faces+=s
 
288
       fixed_v_mask+=s*unitVector(d,DIM)
 
289
#
 
290
#   set up velovity problem
 
291
#
 
292
flow=IncompressibleIsotropicFlowCartesian(dom, stress=stress, v=v, p=p, t=t, numMaterials=2, verbose=VERBOSE)
 
293
flow.setDruckerPragerLaw(tau_Y=TAU_Y/P_REF+BETA*(1.-Function(dom).getX()[DIM-1]))
 
294
 
 
295
flow.setElasticShearModulus(MUE)
 
296
flow.setTolerance(FLOW_TOL)
 
297
flow.setEtaTolerance(FLOW_TOL)
 
298
flow.setExternals(fixed_v_mask=fixed_v_mask)
 
299
print "<%s> Flow solver has been set up."%time.asctime()
 
300
#
 
301
# topography set-up
 
302
#
 
303
boundary=FunctionOnBoundary(dom).getX()[DIM-1]
 
304
top_boundary_mask=whereZero(boundary-sup(boundary))
 
305
surface_area=integrate(top_boundary_mask)
 
306
if CREATE_TOPOGRAPHY:
 
307
    mts=Mountains(dom,eps=TOPO_SMOOTH)
 
308
    mts.setTopography(topography)
 
309
    print "<%s> topography has been set up."%time.asctime()
 
310
#
 
311
#  let the show begin:
 
312
#
 
313
t1 = time.time()
 
314
print "<%s> Start time step %s (t=%s)."%(time.asctime(),n,t)
 
315
while t<T_END:
 
316
    if CREATE_TOPOGRAPHY: topography_old=topography
 
317
    v_old, p_old, stress_old=v, p, stress
 
318
    T_old=T
 
319
    #======= solve for velovity ====================================================================
 
320
    eta_N=exp(Ar*((1.+V_REF*(1-Function(dom).getX()[DIM-1]))/(T_OFFSET_REF+interpolate(T,Function(dom)))-1./T_OFFSET_REF))
 
321
    print "viscosity range :", inf(eta_N), sup(eta_N)
 
322
    flow.setPowerLaws([eta_N, eta_N ], [ 1., TAU_0],  [1,N])
 
323
    flow.setExternals(F=Ra*T*unitVector(DIM-1,DIM))
 
324
    # if dt<=0 or not CREATE_TOPOGRAPHY:
 
325
    if not CREATE_TOPOGRAPHY:
 
326
            flow.update(dt, iter_max=100, verbose=False)
 
327
    else:
 
328
        topography_last=topography
 
329
        Topo_norm, error_Topo=1,1
 
330
        i=0
 
331
        print "DDDDD : ====",dt
 
332
        while error_Topo > TOPO_TOL * Topo_norm:
 
333
            flow.setStatus(t, v_old, p_old, stress_old)
 
334
            flow.setExternals(f=-SURFACE_LOAD*(topography-dt*v)*unitVector(DIM-1,DIM)*top_boundary_mask, restoration_factor=SURFACE_LOAD*dt*top_boundary_mask) 
 
335
            flow.update(dt, iter_max=100, verbose=False)
 
336
            v=flow.getVelocity()
 
337
            mts.setTopography(topography_old)
 
338
            mts.setVelocity(v)
 
339
            topography_last, topography=topography, mts.update(dt, allow_substeps=True)
 
340
            error_Topo=sqrt(integrate(((topography-topography_last)*top_boundary_mask)**2))
 
341
            Topo_norm=sqrt(integrate((topography*top_boundary_mask)**2))
 
342
            print "DDDDD :", "input=",integrate(v*top_boundary_mask)[DIM-1],"output=", integrate(topography*top_boundary_mask)/Lsup(topography), error_Topo, Topo_norm
 
343
            print "topography update step %s error = %e, norm = %e."%(i, error_Topo, Topo_norm), Lsup(v)
 
344
            i+=1
 
345
            if i > TOPO_ITER_MAX: 
 
346
               raise RuntimeError,"topography did not converge after %s steps."%TOPO_ITER_MAX
 
347
    v=flow.getVelocity()
 
348
    for d in range(DIM):
 
349
         print "range %d-velocity"%d,inf(v[d]),sup(v[d])
 
350
    print "Courant = ",inf(dom.getSize()/(length(v)+1e-19)), inf(dom.getSize()**2)
 
351
    print "<%s> flow solver completed."%time.asctime()
 
352
    n+=1
 
353
    t+=dt
 
354
    # print "influx= ",integrate(inner(v,dom.getNormal())), sqrt(integrate(length(v)**2,FunctionOnBoundary(dom))), integrate(1., FunctionOnBoundary(dom))
 
355
    print "<%s> Time step %s (t=%s) completed."%(time.asctime(),n,t)
 
356
    #======= setup Temperature problem ====================================================================
 
357
    #
 
358
    heat.setValue(v=v,Q=CHI_REF*flow.getTau()**2/flow.getCurrentEtaEff())
 
359
    dt=heat.getSafeTimeStepSize()
 
360
    print "<%s> New time step size is %e"%(time.asctime(),dt)
 
361
    if n == 10: 1/0
 
362
    #======= set-up topography ==================================================================================
 
363
    if CREATE_TOPOGRAPHY:
 
364
        dt=min(mts.getSafeTimeStepSize()*0.5,dt)
 
365
        print "<%s> New time step size is %e"%(time.asctime(),dt)
 
366
    print "<%s> Start time step %s (t=%s)."%(time.asctime(),n+1,t+dt)
 
367
    #
 
368
    #  solve temperature:
 
369
    #
 
370
    T=heat.getSolution(dt)
 
371
    print "Temperature range ",inf(T),sup(T)
 
372
    print "<%s> temperature update completed."%time.asctime()
 
373
    #======= anaysis ==================================================================================
 
374
    #
 
375
    #  .... Nusselt number
 
376
    #
 
377
    dTdz=grad(T)[DIM-1]
 
378
    Nu=1.-integrate(v[DIM-1]*T)/integrate(dTdz)
 
379
    eta_bar=integrate(flow.getTau())/integrate(flow.getTau()/flow.getCurrentEtaEff())
 
380
    Ra_eff= (t_REF*RHO_0*G*H*(T_1-T_0)*ALPHA_0)/eta_bar
 
381
    print "nusselt number = %e"%Nu
 
382
    print "av. eta = %e"%eta_bar
 
383
    print "effective Rayleigh number = %e"%Ra_eff
 
384
    if CREATE_TOPOGRAPHY:
 
385
       topo_level=integrate(topography*top_boundary_mask)/surface_area
 
386
       valleys_deep=inf(topography)
 
387
       mountains_heigh=sup(topography)
 
388
       print "topography level = ",topo_level
 
389
       print "valleys deep = ", valleys_deep
 
390
       print "mountains_heigh = ", mountains_heigh
 
391
       diagnostics_file.write("%e %e %e %e %e %e %e\n"%(t,Nu, topo_level, valleys_deep, mountains_heigh, eta_bar, Ra_eff))
 
392
    else:
 
393
       diagnostics_file.write("%e %e %e %e\n"%(t,Nu, eta_bar, Ra_eff))
 
394
    #
 
395
    #  .... visualization
 
396
    #
 
397
    if t>=t_vis or n>n_vis:
 
398
      if CREATE_TOPOGRAPHY:
 
399
         saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff(), topography=topography_old, vex=mts.getVelocity())
 
400
      else:
 
401
         saveVTK(os.path.join(VIS_DIR,"state.%d.vtu"%counter_vis),T=T,v=v,eta=flow.getCurrentEtaEff())
 
402
      print "<%s> Visualization file %d for time step %e generated."%(time.asctime(),counter_vis,t)
 
403
      counter_vis+=1
 
404
      t_vis+=DT_VIS
 
405
      n_vis+=DN_VIS
 
406
    # =========================
 
407
    #
 
408
    #    create restart files:
 
409
    #
 
410
    if DN_RESTART>0:
 
411
       if (n-1)%DN_RESTART == 0:
 
412
         c=(n-1)/DN_RESTART
 
413
         old_restart_dir="%s_%s_"%(PREFIX_RESTART,c-1)
 
414
         new_restart_dir="%s_%s_"%(PREFIX_RESTART,c)
 
415
         mkDir(new_restart_dir)
 
416
         dom.MPIBarrier()
 
417
         file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_vis, n_vis, n, counter_vis, dt))
 
418
         v.dump(os.path.join(new_restart_dir,"v.nc"))
 
419
         p.dump(os.path.join(new_restart_dir,"p.nc"))
 
420
         T.dump(os.path.join(new_restart_dir,"T.nc"))
 
421
         if CREATE_TOPOGRAPHY: topography.dump(os.path.join(new_restart_dir,"topo.nc"))
 
422
         removeRestartDirectory(old_restart_dir)
 
423
         print "<%s> Restart files written to %s."%(time.asctime(),new_restart_dir)
 
424
print "<%s> Calculation finalized after %s sec."%(time.asctime(),time.time() - t1)
 
425