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

  ViewVC Help
Powered by ViewVC 1.1.26