/[escript]/trunk/finley/test/python/convection.py
ViewVC logotype

Annotation of /trunk/finley/test/python/convection.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26