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

revision 1471 by artak, Wed Apr 2 03:33:43 2008 UTC revision 1472 by gross, Thu Apr 3 09:07:07 2008 UTC
# Line 1  Line 1
1  from esys.escript import *  from esys.escript import *
2  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian  from esys.escript.models import TemperatureCartesian, StokesProblemCartesian
3  from esys.finley import Rectangle, Brick  from esys.finley import Rectangle, Brick
4    from optparse import OptionParser
5  from math import pi, ceil  from math import pi, ceil
6  NE=40
7  DIM=3  NE=20
8    DIM=2
9  H=1.  H=1.
10  L=4*H  L=4*H
11  THETA=0.5  THETA=0.5
# Line 17  A=22.  # Arenious number Line 19  A=22.  # Arenious number
19  DI = 0.  # dissipation number  DI = 0.  # dissipation number
20  SUPG=False  SUPG=False
21  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
22
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  #  #
31  #   set up domain:  #   set up domain:
32  #  #
33  if DIM==2:  if restart:
34    dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True)     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
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  else:  else:
56    dom=Brick(int(ceil(L*NE/H)),int(ceil(L*NE/H)),NE,l0=L,l1=L,l2=H,order=2, useFullElementOrder=True)    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
69      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  vol=integrate(1.,Function(dom))  vol=integrate(1.,Function(dom))
79  x=dom.getX()  x=dom.getX()
# Line 34  heat=TemperatureCartesian(dom,theta=THET Line 84  heat=TemperatureCartesian(dom,theta=THET
84  heat.setTolerance(TOL)  heat.setTolerance(TOL)
85
86  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
T=Scalar(1,ReducedSolution(dom))
for d in range(DIM):
if d == DIM-1:
T*=sin(x[d]/H*pi)
else:
T*=cos(x[d]/H*pi)
T=1.-x[DIM-1]+PERTURBATION*T
87  heat.setInitialTemperature(T)  heat.setInitialTemperature(T)
88  print "initial Temperature range ",inf(T),sup(T)  print "initial Temperature range ",inf(T),sup(T)
# Line 52  sp.setTolerance(TOL) Line 95  sp.setTolerance(TOL)
95  sp.setToleranceReductionFactor(TOL)  sp.setToleranceReductionFactor(TOL)
96  x2=ReducedSolution(dom).getX()  x2=ReducedSolution(dom).getX()
97  p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)  p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
v=Vector(0,Solution(dom))
98
100  for d in range(DIM):  for d in range(DIM):
# Line 65  for d in range(DIM): Line 107  for d in range(DIM):
107  #  #
108  #  let the show begin:  #  let the show begin:
109  #  #
t=0
t_out=0
n=0
n_out=0
dt=None
110  dt_new=1.  dt_new=1.
111  a=0  a=None
112    a_alt=None
113  if dom.getMPIRank() ==0:  if dom.getMPIRank() ==0:
114       nusselt_file=open("nusselt.csv","w")       nusselt_file=open("nusselt.csv","w")
115  while t<T_END:  while t<T_END:
116      v_last=v*1.      v_last=v*1
117      print "============== solve for v ========================"      print "============== solve for v ========================"
118      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-0.5))      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))
119      print "viscosity range :", inf(viscosity), sup(viscosity)      print "viscosity range :", inf(viscosity), sup(viscosity)
121      v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='GMRES')      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
124      for d in range(DIM):      for d in range(DIM):
125           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
126
127      if t>=t_out:      if t>=t_out:
128        saveVTK("state.%d.vtu"%n_out,T=T,v=v)        saveVTK("state.%d.vtu"%n_out,T=T,v=v,p=p)
129        print "visualization file %d for time step %e generated."%(n_out,t)        print "visualization file %d for time step %e generated."%(n_out,t)
130        n_out+=1        n_out+=1
131        t_out+=DT_OUT        t_out+=DT_OUT
133      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
134      print "nusselt number = ",Nu      print "nusselt number = ",Nu
135      if n>0:      if a==None:
136           a=0
137        else:
138         a,a_alt = (v_last-v)/dt, a         a,a_alt = (v_last-v)/dt, a
139         dt_alt=dt         dt_alt=dt
140      if n>1:         if a_alt!=None:
141         z=(a-a_alt)/((dt+dt_alt)/2)            z=(a-a_alt)/((dt+dt_alt)/2)
142         f=Lsup(z)/Lsup(v)            f=Lsup(z)/Lsup(v)
143         print "estimated error ",f*dt**2            print "estimated error ",f*dt**2
144         dt_new, dt_alt =min(2*dt,max(dt/2,sqrt(0.05/f))), dt            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            # dt_new, dt_alt =sqrt(0.05/f), dt