/[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 1472 - (hide annotations)
Thu Apr 3 09:07:07 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 4655 byte(s)
some first steps towards restart added
1 gross 1417 from esys.escript import *
2     from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3     from esys.finley import Rectangle, Brick
4 gross 1472 from optparse import OptionParser
5 gross 1417 from math import pi, ceil
6 gross 1472
7     NE=20
8     DIM=2
9 gross 1417 H=1.
10 gross 1463 L=4*H
11 gross 1417 THETA=0.5
12     TOL=1.e-3
13     PERTURBATION=0.1
14 gross 1463 T_END=0.3
15     DT_OUT=T_END/500
16 gross 1417 VERBOSE=False
17 gross 1463 RA=1.e5 # Rayleigh number
18     A=22. # Arenious number
19 gross 1417 DI = 0. # dissipation number
20 gross 1463 SUPG=False
21     print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
22 gross 1472
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 gross 1417 #
31     # set up domain:
32     #
33 gross 1472 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 gross 1417 else:
56 gross 1472 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 gross 1417
69 gross 1472 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 gross 1417 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 gross 1472 a=None
112     a_alt=None
113 gross 1417 if dom.getMPIRank() ==0:
114 gross 1463 nusselt_file=open("nusselt.csv","w")
115 gross 1417 while t<T_END:
116 gross 1472 v_last=v*1
117 gross 1417 print "============== solve for v ========================"
118 gross 1472 viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))
119 gross 1417 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 gross 1472 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 gross 1417
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 gross 1472 saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p)
129 gross 1417 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 gross 1472 if a==None:
136     a=0
137     else:
138 gross 1417 a,a_alt = (v_last-v)/dt, a
139     dt_alt=dt
140 gross 1472 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 gross 1417 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