/[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 1473 by gross, Thu Apr 3 09:07:07 2008 UTC revision 1474 by gross, Mon Apr 7 01:27:01 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, LoadMesh
4  from optparse import OptionParser  from optparse import OptionParser
5  from math import pi, ceil  from math import pi, ceil
6    def removeRestartDirectory(dir_name):
7       if os.path.isdir(dir_name):
8           for root, dirs, files in os.walk(dir_name, topdown=False):
9               for name in files: os.remove(os.path.join(root,name))
10               for name in dirs: os.remove(os.path.join(root,name))
11           os.rmdir(dir_name)
12           print "Restart files %s have been removed."%dir_name
13    
14  NE=20  NE=20
15  DIM=2  DIM=2
# Line 18  RA=1.e5 # Rayleigh number Line 25  RA=1.e5 # Rayleigh number
25  A=22.  # Arenious number  A=22.  # Arenious number
26  DI = 0.  # dissipation number  DI = 0.  # dissipation number
27  SUPG=False  SUPG=False
28    create_restartfiles_every_step=10
29    
30  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)  print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)
31    
32  # read options:  # read options:
33  parser = OptionParser(usage="%prog [-r [DIR]]")  parser = OptionParser(usage="%prog [-r [DIR]]")
34  parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory", default=False, action="store_true")  parser.add_option("-r", "--restart", dest="restart", help="restart from latest directory. It will be deleted after a new directory has been created.", default=False, action="store_true")
35  parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR",metavar="DIR", default=None)  parser.add_option("-d", "--dir", dest="restart_dir", help="restart from directory DIR. The directory will not be deleted but new restart directories are created.",metavar="DIR", default=None)
36  (options, args) = parser.parse_args()  (options, args) = parser.parse_args()
37  restart=options.restart or (options.restart_dir !=None)  restart=options.restart or (options.restart_dir !=None)
 print restart  
38  #  #
39  #   set up domain:  #   set up domain:
40  #  #
# Line 36  if restart: Line 44  if restart:
44        for f in os.listdir("."):        for f in os.listdir("."):
45            if f.startswith("restart"): restart_files.append(f)            if f.startswith("restart"): restart_files.append(f)
46        if len(restart_files)==0:        if len(restart_files)==0:
47            raise RunTimeError,"no restart files"            raise IOError,"no restart files"
48        restart_files.sort()        restart_files.sort()
49        f=restart_files[-1]        f=restart_files[-1]
50     else:     else:
51         f=options.restart_dir         f=options.restart_dir
52     print "restart from ",f     print "restart from directory ",f
53     dom=finley.load("mesh.nc")     dom=LoadMesh("mesh.nc")
54     velocity=load(os.path.join(f,"v.nc"),dom)     ff=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
    pressure=load(os.path.join(f,"p.nc"),dom)  
    ff=open(os.path.join(f,"s.%d"%dom.getMPIRank()),"r").read().split(";")  
55     t=float(ff[0])     t=float(ff[0])
56     t_out=float(ff[1])     t_out=float(ff[1])
57     n=int(ff[2])     n=int(ff[2])
58     n_out=int(ff[3])     n_out=int(ff[3])
59     dt=float(ff[4])     dt=float(ff[4])
60       v=load(os.path.join(f,"v.nc"),dom)
61       p=load(os.path.join(f,"p.nc"),dom)
62       T=load(os.path.join(f,"T.nc"),dom)
63       if n>1:
64          dt_a=float(ff[5])
65          a=load(os.path.join(f,"a.nc"),dom)
66       else:
67          dt_a=None
68          a=None
69     if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")     if dom.getMPIRank()==0: nusselt_file=open("nusselt.csv","a")
70  else:  else:
71    if DIM==2:    if DIM==2:
# Line 74  else: Line 89  else:
89    n=0    n=0
90    n_out=0    n_out=0
91    dt=None    dt=None
92      a=None
93      dt_a=None
94    
95  vol=integrate(1.,Function(dom))  vol=integrate(1.,Function(dom))
96  x=dom.getX()  x=dom.getX()
# Line 107  for d in range(DIM): Line 124  for d in range(DIM):
124  #  #
125  #  let the show begin:  #  let the show begin:
126  #  #
 dt_new=1.  
 a=None  
 a_alt=None  
 if dom.getMPIRank() ==0:  
      nusselt_file=open("nusselt.csv","w")  
127  while t<T_END:  while t<T_END:
128      v_last=v*1      v_last=v*1
129      print "============== solve for v ========================"      print "============== solve for v ========================"
# Line 131  while t<T_END: Line 143  while t<T_END:
143        t_out+=DT_OUT        t_out+=DT_OUT
144      Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)      Nu=1.+integrate(viscosity*length(grad(v))**2)/(RA*vol)
145      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))      if dom.getMPIRank() ==0: nusselt_file.write("%e %e\n"%(t,Nu))
     print "nusselt number = ",Nu  
     if a==None:  
        a=0  
     else:  
        a,a_alt = (v_last-v)/dt, a  
        dt_alt=dt  
        if a_alt!=None:  
           z=(a-a_alt)/((dt+dt_alt)/2)  
           f=Lsup(z)/Lsup(v)  
           print "estimated error ",f*dt**2  
           dt_new, dt_alt =min(2*dt,max(dt/2,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())      print "nusselt number = ",Nu,n
148        if n>0:
149            a,a_alt = (v_last-v)/dt, a
150            dt_a,dt_a_alt = dt, dt_a
151        if n>1:
152           z=(a-a_alt)/((dt_a+dt_a_alt)/2)
153           f=Lsup(z)/Lsup(v)
154           print "estimated error ",f*dt**2
155           dt_new=min(2*dt,max(dt/2,sqrt(0.05/f)))
156           dt=min(dt_new,heat.getSafeTimeStepSize())
157        else:
158           dt=heat.getSafeTimeStepSize()
159      print n,". time step t=",t," step size ",dt      print n,". time step t=",t," step size ",dt
160      print "============== solve for T ========================"      print "============== solve for T ========================"
161      T=heat.solve(dt, verbose=VERBOSE)        T=heat.solve(dt, verbose=VERBOSE)  
162      print "Temperature range ",inf(T),sup(T)      print "Temperature range ",inf(T),sup(T)
163      n+=1      n+=1
164      t+=dt      t+=dt
165        # =========================
166        #
167        #    create restart files:
168        #
169        if create_restartfiles_every_step>0:
170           if (n-1)%create_restartfiles_every_step == 0:
171             c=(n-1)/create_restartfiles_every_step
172             old_restart_dir="restart_%s_"%(c-1)
173             new_restart_dir="restart_%s_"%c
174    
175             print "Write restart files to ",new_restart_dir
176             if not os.path.isdir(new_restart_dir): os.mkdir(new_restart_dir)
177             v.dump(os.path.join(new_restart_dir,"v.nc"))
178             p.dump(os.path.join(new_restart_dir,"p.nc"))
179             T.dump(os.path.join(new_restart_dir,"T.nc"))
180             if n>1:
181                 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e; %e;\n"%(t, t_out, n, n_out, dt, dt_a))
182                 a.dump(os.path.join(new_restart_dir,"a.nc"))
183             else:
184                 file(os.path.join(new_restart_dir,"stamp.%d"%dom.getMPIRank()),"w").write("%e; %e; %s; %s; %e;\n"%(t, t_out, n, n_out, dt))
185             removeRestartDirectory(old_restart_dir)
186    
Legend:
Removed from v.1473  
changed lines
  Added in v.1474
  ViewVC Help
Powered by ViewVC 1.1.26