/[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 2437 by jfenwick, Mon Mar 30 02:13:58 2009 UTC revision 2438 by gross, Tue May 26 02:26:52 2009 UTC
# Line 1  Line 1 
   
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2003-2008 by University of Queensland  # Copyright (c) 2003-2008 by University of Queensland
# Line 10  Line 9 
9  # http://www.opensource.org/licenses/osl-3.0.php  # http://www.opensource.org/licenses/osl-3.0.php
10  #  #
11  ########################################################  ########################################################
12    """
13    this is a convection simulation over a domain [0,L] X [0,L] x [0,H]
14    
15    """
16  __copyright__="""Copyright (c) 2003-2008 by University of Queensland  __copyright__="""Copyright (c) 2003-2008 by University of Queensland
17  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
18  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
# Line 24  from esys.escript.models import Temperat Line 26  from esys.escript.models import Temperat
26  from esys.finley import Rectangle, Brick, LoadMesh  from esys.finley import Rectangle, Brick, LoadMesh
27  from optparse import OptionParser  from optparse import OptionParser
28  from math import pi, ceil  from math import pi, ceil
   
   
 def removeRestartDirectory(dir_name):  
    if dom.onMasterProcessor() and os.path.isdir(dir_name):  
        for root, dirs, files in os.walk(dir_name, topdown=False):  
            for name in files: os.remove(os.path.join(root,name))  
            for name in dirs: os.remove(os.path.join(root,name))  
        os.rmdir(dir_name)  
        print "Restart files %s have been removed."%dir_name  
    dom.MPIBarrier()  
   
   
29  import sys  import sys
30  import time  import time
 t1 = time.time()  
31    
32  extratol=1  # ======================= Default Values ==================================================
33    DIM=2                           # spatial dimension
34    H=1.                            # height
35    L=4*H                           # length
36    NE=10                           # number of elements in H-direction.
37    PERTURBATION=0.1                # initial temperature perturbation
38    DT=1.e-4                        # initial time step size
39    USE_BACKWARD_EULER=False        # use backward Euler in time integartion of temperature otherwise Crank-Nicholson is used
40    CREATE_TOPGRAPHY=False          # create topgraphy
41    TOL=1.e-4                       # tolerance
42    DT_MIN=1.e-5                    # minumum time step size
43    T_END=0.1                       # end time
44    
 # read options:  
 parser = OptionParser(usage="%prog [-r [DIR]] [-e [NE]]")  
 parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=16)  
 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")  
 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)  
 (options, args) = parser.parse_args()  
 restart=options.restart or (options.restart_dir !=None)  
   
 NE=int(options.NE)  
   
 DIM=2  
 H=1.  
 L=4*H  
 USE_BACKWARD_EULER=False  
 TOL=1.e-4  
 PERTURBATION=0.1  
 DT=1.e-4  
 DT_MIN=1.e-5  
 T_END=0.1  
45  DT_OUT=T_END/500  DT_OUT=T_END/500
46  Dn_OUT=2  Dn_OUT=2
47  VERBOSE=True  VERBOSE=True
48  create_restartfiles_every_step=10  CREATE_RESTARTFILES_EVERY_STEP=10
49  if True:  if True:
50     # this is a simple linear Stokes model:     # this is a simple linear Stokes model:
51     RA=1.e6 # Rayleigh number     RA=1.e6 # Rayleigh number
# Line 103  else: Line 84  else:
84     ETAP0=ETA0     ETAP0=ETA0
85     useJAUMANNSTRESS=True     useJAUMANNSTRESS=True
86    
 print "total number of elements = ",NE**DIM*int(L/H)**(DIM-1)  
87    
88    # =========================================================================================
89    
90    def removeRestartDirectory(dir_name):
91       if dom.onMasterProcessor() and os.path.isdir(dir_name):
92           for root, dirs, files in os.walk(dir_name, topdown=False):
93               for name in files: os.remove(os.path.join(root,name))
94               for name in dirs: os.remove(os.path.join(root,name))
95           os.rmdir(dir_name)
96           print "Restart files %s have been removed."%dir_name
97       dom.MPIBarrier()
98    # =========================================================================================
99  #  #
100  #   set up domain:  # read options:
101    #
102    parser = OptionParser(usage="%prog [Options]")
103    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")
104    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)
105    parser.add_option("-p", "--param", dest="param", help="name of file to be imported ",metavar="PARAM", default=None)
106    (options, args) = parser.parse_args()
107    restart=options.restart or (options.restart_dir !=None)
108    #
109    #  overwrite the default options:
110    #
111    print "Execution started ",time.asctime()
112    if options.param !=None:
113         exec(open(options.param,'r'))
114         print "Parameters are imported from file ",options.param
115    print "Parameters:"
116    print "\tDimension\tDIM=\t",DIM
117    print "\tHeight:\tH=\t",H
118    print "\tLength:\tL=\t",L
119    print "\telements in H:\tNE=\t",NE
120    print "\ttotal #element\t\t=\t",NE**DIM*int(L/H)**(DIM-1)
121    print "\ttolerance\tTOL=\t",TOL
122    print "\ttemperature perturbation\tPERTURBATION=\t",PERTURBATION
123    print "\tinitial time step size\tDT=\t",DT
124    print "\tminimum time step size\tDT_MIN=\t",tDT_MIN
125    print "\tend time\tT_END=\t",T_END
126    print "\tbackward Euler?\tUSE_BACKWARD_EULER=\t",USE_BACKWARD_EULER
127    print "\ttopgraphy?\tCREATE_TOPOGRAPHY=\t",CREATE_TOPOGRAPHY
128    #  some control variables (will be overwritten in  case of a restart:
129    #
130    t=0         # time stamp
131    n=0         # time step counter
132    dt=DT       # current time step size
133    
134    t_out=0     #
135    n_out=0
136    out_count=0
137    a=None
138    dt_a=None
139    #
140    #   set up domain or read restart file
141  #  #
142  if restart:  if restart:
143     if options.restart_dir ==None:     if options.restart_dir ==None:
# Line 119  if restart: Line 150  if restart:
150        f=restart_files[-1]        f=restart_files[-1]
151     else:     else:
152         f=options.restart_dir         f=options.restart_dir
153     print "restart from directory ",f     print ">>>Restart from directory ",f
154     try:     try:
155        dom=LoadMesh("mesh.nc")        dom=LoadMesh("mesh.nc")
156     except:     except:
157        pass        pass
158     FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")     FF=open(os.path.join(f,"stamp.%d"%dom.getMPIRank()),"r").read().split(";")
159     t=float(FF[0])     t=float(FF[0])          # time stamp
160     t_out=float(FF[1])     t_out=float(FF[1])      #
161     n_out=int(FF[2])     n_out=int(FF[2])
162     n=int(FF[3])     n=int(FF[3])
163     out_count=int(FF[4])     out_count=int(FF[4])
# Line 165  else: Line 196  else:
196    x2=ReducedSolution(dom).getX()    x2=ReducedSolution(dom).getX()
197    p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)    p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
198    if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")    if dom.onMasterProcessor(): nusselt_file=open("nusselt.csv","w")
   t=0  
   t_out=0  
   n_out=0  
   n=0  
   out_count=0  
   dt=DT  
   a=None  
   dt_a=None  
199    
200  vol=integrate(1.,Function(dom))  vol=integrate(1.,Function(dom))
201  p-=integrate(p)/vol  p-=integrate(p)/vol
# Line 181  x=dom.getX() Line 204  x=dom.getX()
204  #   set up heat problem:  #   set up heat problem:
205  #  #
206  heat=TemperatureCartesian(dom,useBackwardEuler=USE_BACKWARD_EULER)  heat=TemperatureCartesian(dom,useBackwardEuler=USE_BACKWARD_EULER)
207  heat.setTolerance(TOL*extratol)  heat.setTolerance(TOL)
   
208  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])  fixed_T_at=whereZero(x[DIM-1])+whereZero(H-x[DIM-1])
209  print "initial Temperature range ",inf(T),sup(T)  print "initial Temperature range ",inf(T),sup(T)
210  heat.setInitialTemperature(clip(T,0))  heat.setInitialTemperature(clip(T,0))
# Line 197  for d in range(DIM): Line 219  for d in range(DIM):
219         ll = H         ll = H
220      else:      else:
221         ll = L         ll = L
222      fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)      if CREATE_TOPOGRAPHY and d==DIM-1:
223           fixed_v_mask+=whereZero(x[d])*unitVector(d,DIM)
224        else:
225           fixed_v_mask+=(whereZero(x[d])+whereZero(x[d]-ll))*unitVector(d,DIM)
226    1/0
227  #  #
228  #   set up velovity problem  #   set up velovity problem
229  #  # ????????????????????
230  sp=StokesProblemCartesian(dom)  sp=StokesProblemCartesian(dom)
231  # ,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)  # ,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
232  # sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)  # sp=PlateMantelModel(dom,stress,v,p,t,useJaumannStress=useJAUMANNSTRESS)
# Line 210  sp.setTolerance(TOL*extratol) Line 236  sp.setTolerance(TOL*extratol)
236  #  #
237  #  let the show begin:  #  let the show begin:
238  #  #
239    t1 = time.time()
240  while t<T_END:  while t<T_END:
241      v_last=v*1      v_last=v*1
242      print "============== solve for v ========================"      print "============== solve for v ========================"

Legend:
Removed from v.2437  
changed lines
  Added in v.2438

  ViewVC Help
Powered by ViewVC 1.1.26