/[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 1483 - (hide annotations)
Wed Apr 9 03:20:00 2008 UTC (11 years, 3 months ago) by artak
File MIME type: text/x-python
File size: 6744 byte(s)
command line input parameters manipulation and execution time calculation are added
1 gross 1417 from esys.escript import *
2     from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3 gross 1474 from esys.finley import Rectangle, Brick, LoadMesh
4 gross 1472 from optparse import OptionParser
5 gross 1417 from math import pi, ceil
6 gross 1474 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 gross 1472
14 artak 1483
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 gross 1472 DIM=2
35 gross 1417 H=1.
36 gross 1463 L=4*H
37 gross 1417 THETA=0.5
38     TOL=1.e-3
39     PERTURBATION=0.1
40 gross 1463 T_END=0.3
41     DT_OUT=T_END/500
42 gross 1417 VERBOSE=False
43 gross 1463 RA=1.e5 # Rayleigh number
44     A=22. # Arenious number
45 gross 1417 DI = 0. # dissipation number
46 gross 1463 SUPG=False
47 gross 1474 create_restartfiles_every_step=10
48    
49 gross 1463 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
50 gross 1472
51     # read options:
52     parser = OptionParser(usage="%prog [-r [DIR]]")
53 gross 1474 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 gross 1472 (options, args) = parser.parse_args()
56     restart=options.restart or (options.restart_dir !=None)
57 gross 1417 #
58     # set up domain:
59     #
60 gross 1472 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 gross 1474 raise IOError,"no restart files"
67 gross 1472 restart_files.sort()
68     f=restart_files[-1]
69     else:
70     f=options.restart_dir
71 gross 1474 print "restart from directory ",f
72 artak 1483 try:
73     dom=LoadMesh("mesh.nc")
74     except:
75     pass
76 gross 1474 ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
77 gross 1472 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 gross 1474 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 gross 1472 if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
92 gross 1417 else:
93 gross 1472 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 artak 1483 try:
98     dom.dump("mesh.nc")
99     except:
100     pass
101 gross 1472 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 gross 1417
109 gross 1472 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 gross 1474 a=None
118     dt_a=None
119 gross 1472
120 gross 1417 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 artak 1483 sp.setTolerance(TOL*extratol)
137 gross 1417 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 gross 1472 v_last=v*1
154 gross 1417 print "============== solve for v ========================"
155 gross 1472 viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))
156 gross 1417 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 artak 1483 #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 gross 1417
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 gross 1472 saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p)
168 gross 1417 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 gross 1474 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 gross 1472 else:
185 gross 1474 dt=heat.getSafeTimeStepSize()
186 gross 1417 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 gross 1474 # =========================
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 artak 1483 elapsed = time.time() - t1
214     print "plot","\t",NE,"\t",elapsed
215 gross 1474

  ViewVC Help
Powered by ViewVC 1.1.26