/[escript]/trunk/doc/examples/usersguide/fluid.py
ViewVC logotype

Diff of /trunk/doc/examples/usersguide/fluid.py

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

revision 2793 by gross, Tue Dec 1 06:10:10 2009 UTC revision 3914 by jfenwick, Wed Jun 20 04:51:47 2012 UTC
# Line 1  Line 1 
1  ########################################################  ########################################################
2  #  #
3  # Copyright (c) 2008-2009 by University of Queensland  # Copyright (c) 2008-2012 by University of Queensland
4  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
5  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
6  #  #
# Line 18  from esys.escript import * Line 18  from esys.escript import *
18  import esys.finley  import esys.finley
19  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
20  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian
21    from esys.weipa import saveVTK
22    
23  #physical constants  #physical constants
24  eta=1.0  eta=1.0
# Line 37  useUzawa='TRUE' Line 38  useUzawa='TRUE'
38  H=2.0  H=2.0
39  L=1.0  L=1.0
40  W=1.0  W=1.0
41  mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20) # use linear macro elements for pressure  mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20, useElementsOnFace=0) # use linear macro elements for pressure
42  coordinates = mesh.getX()  coordinates = mesh.getX()
43    
44  #gravitational force  #gravitational force
# Line 60  solution.setTolerance(tolerance) Line 61  solution.setTolerance(tolerance)
61    
62  while t <= t_end:  while t <= t_end:
63    
64    print " ----- Time step = %s -----"%( t )    print(" ----- Time step = %s -----"%( t ))
65    print "Time = %s seconds"%( time )      print("Time = %s seconds"%( time ))  
66    
67    solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)    solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)
68    velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)    velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
69        
70    print "Max velocity =", Lsup(velocity), "m/s"    print("Max velocity =", Lsup(velocity), "m/s")
71        
72    #Courant condition    #Courant condition
73    dt=0.4*h/(Lsup(velocity))    dt=0.4*h/(Lsup(velocity))
74    print "dt", dt    print("dt", dt)
75        
76    #displace the mesh    #displace the mesh
77    displacement = velocity * dt    displacement = velocity * dt
78    coordinates = mesh.getX()    coordinates = mesh.getX()
79    mesh.setX(coordinates + displacement)      newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))
80      mesh.setX(newx)  
81        
82    time += dt    time += dt
83        

Legend:
Removed from v.2793  
changed lines
  Added in v.3914

  ViewVC Help
Powered by ViewVC 1.1.26