/[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 1877 - (hide annotations)
Tue Oct 14 02:58:39 2008 UTC (12 years ago) by ksteube
File MIME type: text/x-python
File size: 8913 byte(s)
convection.py checkpointing uses mkdir/rmdir, and under MPI there
was a race condition.

mkdir needs to be run on only one CPU and then a barrier to prevent
working processors from using the directory before it exists.

Added methods domain.MPIBarrier and domain.onMasterProcessor() to
implement this technique.

A more general solution might be possible in the future.

1 ksteube 1811
2     ########################################################
3     #
4     # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7     #
8     # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11     #
12     ########################################################
13    
14     __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 gross 1417 from esys.escript import *
23 gross 1639 from esys.escript.models import TemperatureCartesian, PlateMantelModel
24 gross 1474 from esys.finley import Rectangle, Brick, LoadMesh
25 gross 1472 from optparse import OptionParser
26 gross 1417 from math import pi, ceil
27 gross 1639
28    
29 gross 1474 def removeRestartDirectory(dir_name):
30 ksteube 1877 if dom.onMasterProcessor() and os.path.isdir(dir_name):
31 gross 1474 for root, dirs, files in os.walk(dir_name, topdown=False):
32     for name in files: os.remove(os.path.join(root,name))
33     for name in dirs: os.remove(os.path.join(root,name))
34     os.rmdir(dir_name)
35     print "Restart files %s have been removed."%dir_name
36 ksteube 1877 dom.MPIBarrier()
37 gross 1472
38 artak 1483
39     import sys
40     import time
41     t1 = time.time()
42    
43 artak 1674 extratol=1
44 artak 1483
45 artak 1674 # read options:
46     parser = OptionParser(usage="%prog [-r [DIR]] [-e [NE]] [-s [solver]]")
47     parser.add_option("-s", "--solver", dest="solver", help="solver to be used for saddle point problem. The possible options are PCG, GMRES, NewtonGMRES, MINRES and TFQMR.", metavar="solver", default="PCG")
48     parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=16)
49 artak 1483
50 artak 1674 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")
51     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)
52     (options, args) = parser.parse_args()
53     restart=options.restart or (options.restart_dir !=None)
54 artak 1483
55 artak 1674 solver=options.solver
56     NE=int(options.NE)
57 gross 1639
58     DIM=2
59 gross 1417 H=1.
60 gross 1639 L=4*H
61     THETA=0.5 # time stepping THETA=0.5 cranck nicolson
62     TOL=1.e-4
63 gross 1417 PERTURBATION=0.1
64 gross 1639 DT=1.e-4
65     DT_MIN=1.e-5
66 gross 1552 T_END=0.1
67 gross 1463 DT_OUT=T_END/500
68 gross 1552 Dn_OUT=2
69 gross 1639 VERBOSE=True
70 gross 1474 create_restartfiles_every_step=10
71 gross 1639 if True:
72     # this is a simple linear Stokes model:
73     RA=1.e6 # Rayleigh number
74     A=22 # Arenious number
75     DI = 0. # dissipation number
76     MU=None
77     ETA0=1.
78     TAU0=None
79     N=None
80     NPL=None
81     ETAP0=ETA0
82     TAUY=None
83     useJAUMANNSTRESS=False
84 gross 1659 # this is a simple linear Stokes model:
85     RA=1.e5 # Rayleigh number
86     A=0 # Arenious number
87     DI = 0. # dissipation number
88     MU=None
89     ETA0=1.
90     TAU0=250.
91     N=None
92     NPL=None
93     ETAP0=ETA0
94     TAUY=TAU0
95     useJAUMANNSTRESS=False
96 gross 1639 else:
97     RA=1.e4 # Rayleigh number
98     A=22 # Arenious number
99     DI = 0. # dissipation number
100     MU=1.e4
101     ETA0=1.
102     TAU0=0.866*10**2.5
103     N=3
104     NPL=14
105     TAUY=TAU0
106     ETAP0=ETA0
107     useJAUMANNSTRESS=True
108 gross 1474
109 gross 1463 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
110 gross 1472
111 gross 1417 #
112     # set up domain:
113     #
114 gross 1472 if restart:
115     if options.restart_dir ==None:
116     restart_files=[]
117     for f in os.listdir("."):
118     if f.startswith("restart"): restart_files.append(f)
119     if len(restart_files)==0:
120 gross 1474 raise IOError,"no restart files"
121 gross 1472 restart_files.sort()
122     f=restart_files[-1]
123     else:
124     f=options.restart_dir
125 gross 1474 print "restart from directory ",f
126 artak 1483 try:
127     dom=LoadMesh("mesh.nc")
128     except:
129     pass
130 gross 1639 FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
131     t=float(FF[0])
132     t_out=float(FF[1])
133     n_out=int(FF[2])
134     n=int(FF[3])
135     out_count=int(FF[4])
136     dt=float(FF[5])
137     stress=load(os.path.join(f,"stress.nc"),dom)
138 gross 1474 v=load(os.path.join(f,"v.nc"),dom)
139     p=load(os.path.join(f,"p.nc"),dom)
140     T=load(os.path.join(f,"T.nc"),dom)
141     if n>1:
142 gross 1639 dt_a=float(FF[6])
143 gross 1474 a=load(os.path.join(f,"a.nc"),dom)
144     else:
145     dt_a=None
146     a=None
147 ksteube 1877 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","a")
148 gross 1417 else:
149 gross 1472 if DIM==2:
150     dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
151     else:
152     dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True,optimize=True)
153 artak 1483 try:
154     dom.dump("mesh.nc")
155     except:
156     pass
157 gross 1472 x=dom.getX()
158     T=Scalar(1,ReducedSolution(dom))
159     for d in range(DIM):
160     if d == DIM-1:
161     T*=sin(x[d]/H*pi)
162     else:
163 gross 1552 T*=cos(x[d]/L*pi)
164 gross 1417
165 gross 1472 T=1.-x[DIM-1]+PERTURBATION*T
166     v=Vector(0,Solution(dom))
167 gross 1639 stress=Tensor(0,Function(dom))
168     x2=ReducedSolution(dom).getX()
169     p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
170 ksteube 1877 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")
171 gross 1472 t=0
172     t_out=0
173 gross 1552 n_out=0
174 gross 1472 n=0
175 gross 1552 out_count=0
176 gross 1639 dt=DT
177 gross 1474 a=None
178     dt_a=None
179 gross 1472
180 gross 1417 vol=integrate(1.,Function(dom))
181 gross 1639 p-=integrate(p)/vol
182 gross 1417 x=dom.getX()
183     #
184     # set up heat problem:
185     #
186 gross 1841 heat=TemperatureCartesian(dom,theta=THETA)
187 gross 1552 heat.setTolerance(TOL*extratol)
188 gross 1417
189     fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
190     print "initial Temperature range ",inf(T),sup(T)
191 gross 1841 heat.setInitialTemperature(clip(T,0))
192 gross 1417 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
193     #
194 gross 1639 # velocity constraints:
195 gross 1417 #
196     x2=ReducedSolution(dom).getX()
197     fixed_v_mask=Vector(0,Solution(dom))
198     for d in range(DIM):
199     if d == DIM-1:
200     ll = H
201     else:
202     ll = L
203     fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
204     #
205 gross 1639 # set up velovity problem
206     #
207     sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
208     sp.initialize(mu=MU, tau_0=TAU0, n=N, eta_Y=ETAP0, tau_Y=TAU0, n_Y=NPL, q=fixed_v_mask)
209     sp.setTolerance(TOL*extratol)
210     sp.setToleranceReductionFactor(TOL)
211     #
212 gross 1417 # let the show begin:
213     #
214     while t<T_END:
215 gross 1472 v_last=v*1
216 gross 1417 print "============== solve for v ========================"
217 gross 1639 FF=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
218     print "viscosity range :", inf(FF)*ETA0, sup(FF)*ETA0
219     sp.initialize(eta_N=ETA0*FF, eta_0=ETA0*FF, F=T*(RA*unitVector(DIM-1,DIM)))
220 artak 1674 sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver=solver)
221 gross 1639 v=sp.getVelocity()
222 gross 1417
223     for d in range(DIM):
224     print "range %d-velocity"%d,inf(v[d]),sup(v[d])
225 gross 1552 if t>=t_out or n>n_out:
226     saveVTK("state.%d.vtu"%out_count,T=T,v=v)
227     print "visualization file %d for time step %e generated."%(out_count,t)
228     out_count+=1
229 gross 1417 t_out+=DT_OUT
230 gross 1552 n_out+=Dn_OUT
231 gross 1639 # calculation of nusselt number:
232 gross 1659 se=sp.getMechanicalPower()
233     print "Xse:",inf(se),sup(se)
234 gross 1639 Nu=1.+integrate(se)/(RA*vol)
235 ksteube 1877 if dom.onMasterProcessor(): nusselt_file.write("%e %e\n"%(t,Nu))
236 gross 1639 heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
237 gross 1659 print "Xnusselt number = ",Nu, "dt =",dt
238 gross 1474 if n>0:
239     a,a_alt = (v_last-v)/dt, a
240     dt_a,dt_a_alt = dt, dt_a
241     if n>1:
242     z=(a-a_alt)/((dt_a+dt_a_alt)/2)
243     f=Lsup(z)/Lsup(v)
244     print "estimated error ",f*dt**2
245 gross 1552 # dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
246     dt_new=sqrt(0.05/f)
247 gross 1474 dt=min(dt_new,heat.getSafeTimeStepSize())
248 gross 1472 else:
249 gross 1474 dt=heat.getSafeTimeStepSize()
250 gross 1639 dt=max(DT_MIN,dt)
251 gross 1417 print n,". time step t=",t," step size ",dt
252     print "============== solve for T ========================"
253 gross 1841 T=heat.getSolution(dt, verbose=VERBOSE)
254 gross 1417 print "Temperature range ",inf(T),sup(T)
255     n+=1
256     t+=dt
257 gross 1474 # =========================
258     #
259     # create restart files:
260     #
261     if create_restartfiles_every_step>0:
262     if (n-1)%create_restartfiles_every_step == 0:
263     c=(n-1)/create_restartfiles_every_step
264     old_restart_dir="restart_%s_"%(c-1)
265     new_restart_dir="restart_%s_"%c
266    
267     print "Write restart files to ",new_restart_dir
268 ksteube 1877 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
269     dom.MPIBarrier()
270 gross 1639 sp.getStress().dump(os.path.join(new_restart_dir,"stress.nc"))
271     sp.getVelocity().dump(os.path.join(new_restart_dir,"v.nc"))
272     sp.getPressure().dump(os.path.join(new_restart_dir,"p.nc"))
273 gross 1474 T.dump(os.path.join(new_restart_dir,"T.nc"))
274     if n>1:
275 gross 1552 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e; %e;\n"%(t, t_out, n_out, n, out_count, dt, dt_a))
276 gross 1474 a.dump(os.path.join(new_restart_dir,"a.nc"))
277     else:
278 gross 1552 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %s; %e;\n"%(t, t_out, n_out, n, out_count, dt))
279 gross 1474 removeRestartDirectory(old_restart_dir)
280 artak 1483 elapsed = time.time() - t1
281     print "plot","\t",NE,"\t",elapsed
282 gross 1474

  ViewVC Help
Powered by ViewVC 1.1.26