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

  ViewVC Help
Powered by ViewVC 1.1.26