/[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 1466 - (show annotations)
Wed Apr 2 03:33:43 2008 UTC (11 years, 7 months ago) by artak
File MIME type: text/x-python
File size: 3158 byte(s)


1 from esys.escript import *
2 from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3 from esys.finley import Rectangle, Brick
4 from math import pi, ceil
5 NE=40
6 DIM=3
7 H=1.
8 L=4*H
9 THETA=0.5
10 TOL=1.e-3
11 PERTURBATION=0.1
12 T_END=0.3
13 DT_OUT=T_END/500
14 VERBOSE=False
15 RA=1.e5 # Rayleigh number
16 A=22. # Arenious number
17 DI = 0. # dissipation number
18 SUPG=False
19 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
20 #
21 # set up domain:
22 #
23 if DIM==2:
24 dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True)
25 else:
26 dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True)
27
28 vol=integrate(1.,Function(dom))
29 x=dom.getX()
30 #
31 # set up heat problem:
32 #
33 heat=TemperatureCartesian(dom,theta=THETA,useSUPG=SUPG)
34 heat.setTolerance(TOL)
35
36 fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
37 T=Scalar(1,ReducedSolution(dom))
38 for d in range(DIM):
39 if d == DIM-1:
40 T*=sin(x[d]/H*pi)
41 else:
42 T*=cos(x[d]/H*pi)
43 T=1.-x[DIM-1]+PERTURBATION*T
44 heat.setInitialTemperature(T)
45 print "initial Temperature range ",inf(T),sup(T)
46 heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
47 #
48 # set up velovity problem
49 #
50 sp=StokesProblemCartesian(dom)
51 sp.setTolerance(TOL)
52 sp.setToleranceReductionFactor(TOL)
53 x2=ReducedSolution(dom).getX()
54 p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
55 v=Vector(0,Solution(dom))
56
57 fixed_v_mask=Vector(0,Solution(dom))
58 for d in range(DIM):
59 if d == DIM-1:
60 ll = H
61 else:
62 ll = L
63 fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
64
65 #
66 # let the show begin:
67 #
68 t=0
69 t_out=0
70 n=0
71 n_out=0
72 dt=None
73 dt_new=1.
74 a=0
75 if dom.getMPIRank() ==0:
76 nusselt_file=open("nusselt.csv","w")
77 while t<T_END:
78 v_last=v*1.
79 print "============== solve for v ========================"
80 viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-0.5))
81 print "viscosity range :", inf(viscosity), sup(viscosity)
82 sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)
83 v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES')
84
85 for d in range(DIM):
86 print "range %d-velocity"%d,inf(v[d]),sup(v[d])
87
88 if t>=t_out:
89 saveVTK("state.%d.vtu"%n_out,T=T,v=v)
90 print "visualization file %d for time step %e generated."%(n_out,t)
91 n_out+=1
92 t_out+=DT_OUT
93 Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)
94 if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
95 print "nusselt number = ",Nu
96 if n>0:
97 a,a_alt = (v_last-v)/dt, a
98 dt_alt=dt
99 if n>1:
100 z=(a-a_alt)/((dt+dt_alt)/2)
101 f=Lsup(z)/Lsup(v)
102 print "estimated error ",f*dt**2
103 dt_new, dt_alt =min(2*dt,max(dt/2,sqrt(0.05/f))), dt
104 # dt_new, dt_alt =sqrt(0.05/f), dt
105 heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)
106 dt=min(dt_new,heat.getSafeTimeStepSize())
107 print n,". time step t=",t," step size ",dt
108 print "============== solve for T ========================"
109 T=heat.solve(dt, verbose=VERBOSE)
110 print "Temperature range ",inf(T),sup(T)
111 n+=1
112 t+=dt

  ViewVC Help
Powered by ViewVC 1.1.26