/[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 1516 - (show annotations)
Fri Apr 18 01:02:25 2008 UTC (11 years, 3 months ago) by ksteube
File MIME type: text/x-python
File size: 6739 byte(s)
Reduced L and changed Arenious from 22.0 to 0.0

1 from esys.escript import *
2 from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3 from esys.finley import Rectangle, Brick, LoadMesh
4 from optparse import OptionParser
5 from math import pi, ceil
6 def removeRestartDirectory(dir_name):
7 if os.path.isdir(dir_name):
8 for root, dirs, files in os.walk(dir_name, topdown=False):
9 for name in files: os.remove(os.path.join(root,name))
10 for name in dirs: os.remove(os.path.join(root,name))
11 os.rmdir(dir_name)
12 print "Restart files %s have been removed."%dir_name
13
14
15 import sys
16 import time
17 t1 = time.time()
18
19 if (len(sys.argv)>=3):
20 NE=int(sys.argv[2])
21 else:
22 NE=20
23
24 if (len(sys.argv)>=2):
25 solver=sys.argv[1]
26 else:
27 solver='PCG'
28
29 if solver!='PCG':
30 extratol=0.001
31 else:
32 extratol=1
33
34 DIM=3
35 H=1.
36 L=2*H
37 THETA=0.5
38 TOL=1.e-3
39 PERTURBATION=0.1
40 T_END=0.3
41 DT_OUT=T_END/500
42 VERBOSE=False
43 RA=1.e5 # Rayleigh number
44 A=0. # Arenious number
45 DI = 0. # dissipation number
46 SUPG=False
47 create_restartfiles_every_step=10
48
49 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
50
51 # read options:
52 parser = OptionParser(usage="%prog [-r [DIR]]")
53 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")
54 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)
55 (options, args) = parser.parse_args()
56 restart=options.restart or (options.restart_dir !=None)
57 #
58 # set up domain:
59 #
60 if restart:
61 if options.restart_dir ==None:
62 restart_files=[]
63 for f in os.listdir("."):
64 if f.startswith("restart"): restart_files.append(f)
65 if len(restart_files)==0:
66 raise IOError,"no restart files"
67 restart_files.sort()
68 f=restart_files[-1]
69 else:
70 f=options.restart_dir
71 print "restart from directory ",f
72 try:
73 dom=LoadMesh("mesh.nc")
74 except:
75 pass
76 ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
77 t=float(ff[0])
78 t_out=float(ff[1])
79 n=int(ff[2])
80 n_out=int(ff[3])
81 dt=float(ff[4])
82 v=load(os.path.join(f,"v.nc"),dom)
83 p=load(os.path.join(f,"p.nc"),dom)
84 T=load(os.path.join(f,"T.nc"),dom)
85 if n>1:
86 dt_a=float(ff[5])
87 a=load(os.path.join(f,"a.nc"),dom)
88 else:
89 dt_a=None
90 a=None
91 if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
92 else:
93 if DIM==2:
94 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
95 else:
96 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)
97 try:
98 dom.dump("mesh.nc")
99 except:
100 pass
101 x=dom.getX()
102 T=Scalar(1,ReducedSolution(dom))
103 for d in range(DIM):
104 if d == DIM-1:
105 T*=sin(x[d]/H*pi)
106 else:
107 T*=cos(x[d]*((d+1)/L*pi))
108
109 T=1.-x[DIM-1]+PERTURBATION*T
110 v=Vector(0,Solution(dom))
111 if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")
112 t=0
113 t_out=0
114 n=0
115 n_out=0
116 dt=None
117 a=None
118 dt_a=None
119
120 vol=integrate(1.,Function(dom))
121 x=dom.getX()
122 #
123 # set up heat problem:
124 #
125 heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG)
126 heat.setTolerance(TOL)
127
128 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
129 heat.setInitialTemperature(T)
130 print "initial Temperature range ",inf(T),sup(T)
131 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
132 #
133 # set up velovity problem
134 #
135 sp=StokesProblemCartesian(dom)
136 sp.setTolerance(TOL*extratol)
137 sp.setToleranceReductionFactor(TOL)
138 x2=ReducedSolution(dom).getX()
139 p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
140
141 fixed_v_mask=Vector(0,Solution(dom))
142 for d in range(DIM):
143 if d == DIM-1:
144 ll = H
145 else:
146 ll = L
147 fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
148
149 #
150 # let the show begin:
151 #
152 while t<T_END:
153 v_last=v*1
154 print "============== solve for v ========================"
155 viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-1./2.))
156 print "viscosity range :", inf(viscosity), sup(viscosity)
157 sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)
158 #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')
159 #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES')
160 #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='MINRES')
161 v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver=solver)
162
163 for d in range(DIM):
164 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
165
166 if t>=t_out:
167 saveVTK("state.%d.vtu"%n_out,T=T,v=v)
168 print "visualization file %d for time step %e generated."%(n_out,t)
169 n_out+=1
170 t_out+=DT_OUT
171 Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)
172 if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
173 heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)
174 print "nusselt number = ",Nu,n
175 if n>0:
176 a,a_alt = (v_last-v)/dt, a
177 dt_a,dt_a_alt = dt, dt_a
178 if n>1:
179 z=(a-a_alt)/((dt_a+dt_a_alt)/2)
180 f=Lsup(z)/Lsup(v)
181 print "estimated error ",f*dt**2
182 dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
183 dt=min(dt_new,heat.getSafeTimeStepSize())
184 else:
185 dt=heat.getSafeTimeStepSize()
186 print n,". time step t=",t," step size ",dt
187 print "============== solve for T ========================"
188 T=heat.solve(dt, verbose=VERBOSE)
189 print "Temperature range ",inf(T),sup(T)
190 n+=1
191 t+=dt
192 # =========================
193 #
194 # create restart files:
195 #
196 if create_restartfiles_every_step>0:
197 if (n-1)%create_restartfiles_every_step == 0:
198 c=(n-1)/create_restartfiles_every_step
199 old_restart_dir="restart_%s_"%(c-1)
200 new_restart_dir="restart_%s_"%c
201
202 print "Write restart files to ",new_restart_dir
203 if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
204 v.dump(os.path.join(new_restart_dir,"v.nc"))
205 p.dump(os.path.join(new_restart_dir,"p.nc"))
206 T.dump(os.path.join(new_restart_dir,"T.nc"))
207 if n>1:
208 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e; %e;\n"%(t, t_out, n, n_out, dt, dt_a))
209 a.dump(os.path.join(new_restart_dir,"a.nc"))
210 else:
211 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e;\n"%(t, t_out, n, n_out, dt))
212 removeRestartDirectory(old_restart_dir)
213 elapsed = time.time() - t1
214 print "plot","\t",NE,"\t",elapsed
215

  ViewVC Help
Powered by ViewVC 1.1.26