/[escript]/trunk/finley/test/python/convection.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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    
23    # read options:
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
45       dom=finley.load("mesh.nc")
46       velocity=load(os.path.join(f,"v.nc"),dom)
47       pressure=load(os.path.join(f,"p.nc"),dom)
48       ff=open(os.path.join(f,"s.%d"%dom.getMPIRank()),"r").read().split(";")
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)
89  heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)  heat.setValue(rhocp=Scalar(1.,Function(dom)),k=Scalar(1.,Function(dom)),given_T_mask=fixed_T_at)
# 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    
99  fixed_v_mask=Vector(0,Solution(dom))  fixed_v_mask=Vector(0,Solution(dom))
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)
120      sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)      sp.initialize(f=T*(RA*unitVector(DIM-1,DIM)),eta=viscosity,fixed_u_mask=fixed_v_mask)
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
132      Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)      Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)
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
146      heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)      heat.setValue(v=v,Q=DI/RA*viscosity*length(symmetric(grad(v)))**2)
147      dt=min(dt_new,heat.getSafeTimeStepSize())      dt=min(dt_new,heat.getSafeTimeStepSize())
148      print n,". time step t=",t," step size ",dt      print n,". time step t=",t," step size ",dt

Legend:
Removed from v.1471  
changed lines
  Added in v.1472

  ViewVC Help
Powered by ViewVC 1.1.26