/[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 1472 - (show annotations)
Thu Apr 3 09:07:07 2008 UTC (11 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 4655 byte(s)
some first steps towards restart added
1 from esys.escript import *
2 from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3 from esys.finley import Rectangle, Brick
4 from optparse import OptionParser
5 from math import pi, ceil
6
7 NE=20
8 DIM=2
9 H=1.
10 L=4*H
11 THETA=0.5
12 TOL=1.e-3
13 PERTURBATION=0.1
14 T_END=0.3
15 DT_OUT=T_END/500
16 VERBOSE=False
17 RA=1.e5 # Rayleigh number
18 A=22. # Arenious number
19 DI = 0. # dissipation number
20 SUPG=False
21 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
22
23 # read options:
24 parser = OptionParser(usage="%prog [-r [DIR]]")
25 parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory", default=False, action="store_true")
26 parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR",metavar="DIR", default=None)
27 (options, args) = parser.parse_args()
28 restart=options.restart or (options.restart_dir !=None)
29 print restart
30 #
31 # set up domain:
32 #
33 if restart:
34 if options.restart_dir ==None:
35 restart_files=[]
36 for f in os.listdir("."):
37 if f.startswith("restart"): restart_files.append(f)
38 if len(restart_files)==0:
39 raise RunTimeError,"no restart files"
40 restart_files.sort()
41 f=restart_files[-1]
42 else:
43 f=options.restart_dir
44 print "restart from ",f
45 dom=finley.load("mesh.nc")
46 velocity=load(os.path.join(f,"v.nc"),dom)
47 pressure=load(os.path.join(f,"p.nc"),dom)
48 ff=open(os.path.join(f,"s.%d"%dom.getMPIRank()),"r").read().split(";")
49 t=float(ff[0])
50 t_out=float(ff[1])
51 n=int(ff[2])
52 n_out=int(ff[3])
53 dt=float(ff[4])
54 if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
55 else:
56 if DIM==2:
57 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
58 else:
59 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)
60 dom.dump("mesh.nc")
61 x=dom.getX()
62 T=Scalar(1,ReducedSolution(dom))
63 for d in range(DIM):
64 if d == DIM-1:
65 T*=sin(x[d]/H*pi)
66 else:
67 T*=cos(x[d]*((d+1)/L*pi))
68
69 T=1.-x[DIM-1]+PERTURBATION*T
70 v=Vector(0,Solution(dom))
71 if dom.getMPIRank() ==0: nusselt_file=open("nusselt.csv","w")
72 t=0
73 t_out=0
74 n=0
75 n_out=0
76 dt=None
77
78 vol=integrate(1.,Function(dom))
79 x=dom.getX()
80 #
81 # set up heat problem:
82 #
83 heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG)
84 heat.setTolerance(TOL)
85
86 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
87 heat.setInitialTemperature(T)
88 print "initial Temperature range ",inf(T),sup(T)
89 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
90 #
91 # set up velovity problem
92 #
93 sp=StokesProblemCartesian(dom)
94 sp.setTolerance(TOL)
95 sp.setToleranceReductionFactor(TOL)
96 x2=ReducedSolution(dom).getX()
97 p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
98
99 fixed_v_mask=Vector(0,Solution(dom))
100 for d in range(DIM):
101 if d == DIM-1:
102 ll = H
103 else:
104 ll = L
105 fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
106
107 #
108 # let the show begin:
109 #
110 dt_new=1.
111 a=None
112 a_alt=None
113 if dom.getMPIRank() ==0:
114 nusselt_file=open("nusselt.csv","w")
115 while t<T_END:
116 v_last=v*1
117 print "============== solve for v ========================"
118 viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))
119 print "viscosity range :", inf(viscosity), sup(viscosity)
120 sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)
121 v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')
122 # v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES')
123
124 for d in range(DIM):
125 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
126
127 if t>=t_out:
128 saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p)
129 print "visualization file %d for time step %e generated."%(n_out,t)
130 n_out+=1
131 t_out+=DT_OUT
132 Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)
133 if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
134 print "nusselt number = ",Nu
135 if a==None:
136 a=0
137 else:
138 a,a_alt = (v_last-v)/dt, a
139 dt_alt=dt
140 if a_alt!=None:
141 z=(a-a_alt)/((dt+dt_alt)/2)
142 f=Lsup(z)/Lsup(v)
143 print "estimated error ",f*dt**2
144 dt_new, dt_alt =min(2*dt,max(dt/2,sqrt(0.05/f))), dt
145 # dt_new, dt_alt =sqrt(0.05/f), dt
146 heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)
147 dt=min(dt_new,heat.getSafeTimeStepSize())
148 print n,". time step t=",t," step size ",dt
149 print "============== solve for T ========================"
150 T=heat.solve(dt, verbose=VERBOSE)
151 print "Temperature range ",inf(T),sup(T)
152 n+=1
153 t+=dt

  ViewVC Help
Powered by ViewVC 1.1.26