# Contents of /trunk/finley/test/python/convection.py

Revision 1463 - (show annotations)
Tue Apr 1 05:20:57 2008 UTC (14 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 3143 byte(s)
```a better convection code
```
 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_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