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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1877 - (show annotations)
Tue Oct 14 02:58:39 2008 UTC (11 years, 8 months 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
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 from esys.escript import *
23 from esys.escript.models import TemperatureCartesian, PlateMantelModel
24 from esys.finley import Rectangle, Brick, LoadMesh
25 from optparse import OptionParser
26 from math import pi, ceil
27
28
29 def removeRestartDirectory(dir_name):
30 if dom.onMasterProcessor() and os.path.isdir(dir_name):
31 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 dom.MPIBarrier()
37
38
39 import sys
40 import time
41 t1 = time.time()
42
43 extratol=1
44
45 # 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
50 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
55 solver=options.solver
56 NE=int(options.NE)
57
58 DIM=2
59 H=1.
60 L=4*H
61 THETA=0.5 # time stepping THETA=0.5 cranck nicolson
62 TOL=1.e-4
63 PERTURBATION=0.1
64 DT=1.e-4
65 DT_MIN=1.e-5
66 T_END=0.1
67 DT_OUT=T_END/500
68 Dn_OUT=2
69 VERBOSE=True
70 create_restartfiles_every_step=10
71 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 # 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 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
109 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
110
111 #
112 # set up domain:
113 #
114 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 raise IOError,"no restart files"
121 restart_files.sort()
122 f=restart_files[-1]
123 else:
124 f=options.restart_dir
125 print "restart from directory ",f
126 try:
127 dom=LoadMesh("mesh.nc")
128 except:
129 pass
130 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 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 dt_a=float(FF[6])
143 a=load(os.path.join(f,"a.nc"),dom)
144 else:
145 dt_a=None
146 a=None
147 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","a")
148 else:
149 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 try:
154 dom.dump("mesh.nc")
155 except:
156 pass
157 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 T*=cos(x[d]/L*pi)
164
165 T=1.-x[DIM-1]+PERTURBATION*T
166 v=Vector(0,Solution(dom))
167 stress=Tensor(0,Function(dom))
168 x2=ReducedSolution(dom).getX()
169 p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
170 if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")
171 t=0
172 t_out=0
173 n_out=0
174 n=0
175 out_count=0
176 dt=DT
177 a=None
178 dt_a=None
179
180 vol=integrate(1.,Function(dom))
181 p-=integrate(p)/vol
182 x=dom.getX()
183 #
184 # set up heat problem:
185 #
186 heat=TemperatureCartesian(dom,theta=THETA)
187 heat.setTolerance(TOL*extratol)
188
189 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
190 print "initial Temperature range ",inf(T),sup(T)
191 heat.setInitialTemperature(clip(T,0))
192 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
193 #
194 # velocity constraints:
195 #
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 # 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 # let the show begin:
213 #
214 while t<T_END:
215 v_last=v*1
216 print "============== solve for v ========================"
217 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 sp.update(dt,max_inner_iter=20, verbose=VERBOSE, show_details=False, tol=10., solver=solver)
221 v=sp.getVelocity()
222
223 for d in range(DIM):
224 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
225 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 t_out+=DT_OUT
230 n_out+=Dn_OUT
231 # calculation of nusselt number:
232 se=sp.getMechanicalPower()
233 print "Xse:",inf(se),sup(se)
234 Nu=1.+integrate(se)/(RA*vol)
235 if dom.onMasterProcessor(): nusselt_file.write("%e %e\n"%(t,Nu))
236 heat.setValue(v=interpolate(v,ReducedSolution(dom)),Q=DI/RA*se)
237 print "Xnusselt number = ",Nu, "dt =",dt
238 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 # dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
246 dt_new=sqrt(0.05/f)
247 dt=min(dt_new,heat.getSafeTimeStepSize())
248 else:
249 dt=heat.getSafeTimeStepSize()
250 dt=max(DT_MIN,dt)
251 print n,". time step t=",t," step size ",dt
252 print "============== solve for T ========================"
253 T=heat.getSolution(dt, verbose=VERBOSE)
254 print "Temperature range ",inf(T),sup(T)
255 n+=1
256 t+=dt
257 # =========================
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 if dom.onMasterProcessor() and not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
269 dom.MPIBarrier()
270 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 T.dump(os.path.join(new_restart_dir,"T.nc"))
274 if n>1:
275 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 a.dump(os.path.join(new_restart_dir,"a.nc"))
277 else:
278 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 removeRestartDirectory(old_restart_dir)
280 elapsed = time.time() - t1
281 print "plot","\t",NE,"\t",elapsed
282

  ViewVC Help
Powered by ViewVC 1.1.26