/[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 1482 by gross, Mon Apr 7 01:27:01 2008 UTC revision 1483 by artak, Wed Apr 9 03:20:00 2008 UTC
# Line 11  def removeRestartDirectory(dir_name): Line 11  def removeRestartDirectory(dir_name):
11         os.rmdir(dir_name)         os.rmdir(dir_name)
12         print "Restart files %s have been removed."%dir_name         print "Restart files %s have been removed."%dir_name
13    
14  NE=20  
15    import sys
16    import time
17    t1 = time.time()
18    
19    if (len(sys.argv)>=3):
20     NE=int(sys.argv[2])
21    else:
22     NE=20
23    
24    if (len(sys.argv)>=2):
25     solver=sys.argv[1]
26    else:
27     solver='PCG'
28    
29    if solver!='PCG':
30     extratol=0.001
31    else:
32     extratol=1
33    
34  DIM=2  DIM=2
35  H=1.  H=1.
36  L=4*H  L=4*H
# Line 50  if restart: Line 69  if restart:
69     else:     else:
70         f=options.restart_dir         f=options.restart_dir
71     print "restart from directory ",f     print "restart from directory ",f
72     dom=LoadMesh("mesh.nc")     try:
73          dom=LoadMesh("mesh.nc")
74       except:
75          pass
76     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(";")
77     t=float(ff[0])     t=float(ff[0])
78     t_out=float(ff[1])     t_out=float(ff[1])
# Line 72  else: Line 94  else:
94      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)      dom=Rectangle(int(ceil(L*NE/H)),NE,l0=L,l1=H,order=2, useFullElementOrder=True,optimize=True)
95    else:    else:
96      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)      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)
97    dom.dump("mesh.nc")    try:
98         dom.dump("mesh.nc")
99      except:
100         pass
101    x=dom.getX()    x=dom.getX()
102    T=Scalar(1,ReducedSolution(dom))    T=Scalar(1,ReducedSolution(dom))
103    for d in range(DIM):    for d in range(DIM):
# Line 108  heat.setValue(rhocp=Scalar(1.,Function(d Line 133  heat.setValue(rhocp=Scalar(1.,Function(d
133  #   set up velovity problem  #   set up velovity problem
134  #  #
135  sp=StokesProblemCartesian(dom)  sp=StokesProblemCartesian(dom)
136  sp.setTolerance(TOL)  sp.setTolerance(TOL*extratol)
137  sp.setToleranceReductionFactor(TOL)  sp.setToleranceReductionFactor(TOL)
138  x2=ReducedSolution(dom).getX()  x2=ReducedSolution(dom).getX()
139  p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)  p=-RA*(x2[DIM-1]-0.5*x2[DIM-1]**2)
# Line 130  while t<T_END: Line 155  while t<T_END:
155      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))      viscosity=exp(A*(1./(1+T.interpolate(Function(dom)))-2./3.))
156      print "viscosity range :", inf(viscosity), sup(viscosity)      print "viscosity range :", inf(viscosity), sup(viscosity)
157      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)
158      v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')      #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='PCG')
159      # 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='GMRES')
160        #v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver='MINRES')
161        v,p=sp.solve(v,p,show_details=VERBOSE, verbose=True,max_iter=500,solver=solver)
162    
163      for d in range(DIM):      for d in range(DIM):
164           print "range %d-velocity"%d,inf(v[d]),sup(v[d])           print "range %d-velocity"%d,inf(v[d]),sup(v[d])
# Line 183  while t<T_END: Line 210  while t<T_END:
210           else:           else:
211               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))               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))
212           removeRestartDirectory(old_restart_dir)           removeRestartDirectory(old_restart_dir)
213    elapsed = time.time() - t1
214    print "plot","\t",NE,"\t",elapsed
215    
Legend:
Removed from v.1482  
changed lines
  Added in v.1483
  ViewVC Help
Powered by ViewVC 1.1.26