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

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_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