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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2484 - (hide annotations)
Mon Jun 22 04:22:19 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 1651 byte(s)
numarray removed from docu; Locator revised.
1 lgraham 2191 ######## August 2008 ########
2     ########## Leon Graham ##########
3     ## Newtonian fluid using StokesProblemCartesian class##
4    
5     from esys.escript import *
6     import esys.finley
7     from esys.escript.linearPDEs import LinearPDE
8     from esys.escript.models import StokesProblemCartesian
9    
10     #physical constants
11     eta=1.0
12     rho=100.0
13     g=10.0
14    
15     #solver settings
16     tolerance=1.0e-4
17     max_iter=200
18     t_end=50
19     t=0.0
20     time=0
21     verbose='TRUE'
22     useUzawa='TRUE'
23    
24     #define mesh
25     H=2.0
26     L=1.0
27     W=1.0
28     mesh = esys.finley.Rectangle(l0=L, l1=H, order=2, n0=20, n1=20)
29     coordinates = mesh.getX()
30    
31     #gravitational force
32     Y=Vector(0.0, Function(mesh))
33     Y[1]=-rho*g
34    
35     #element spacing
36     h=Lsup(mesh.getSize())
37    
38     #boundary conditions for slip at base
39     boundary_cond=whereZero(coordinates[1])*[0.0,1.0]
40    
41     #velocity and pressure vectors
42     velocity=Vector(0.0, ContinuousFunction(mesh))
43     pressure=Scalar(0.0, ContinuousFunction(mesh))
44    
45     #Stokes Cartesian
46     solution=StokesProblemCartesian(mesh)
47     solution.setTolerance(tolerance)
48    
49     while t <= t_end:
50    
51     print " ----- Time step = %s -----"%( t )
52     print "Time = %s seconds"%( time )
53    
54     solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)
55 gross 2432 velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
56 lgraham 2191
57     print "Max velocity =", Lsup(velocity), "m/s"
58    
59     #Courant condition
60     dt=0.4*h/(Lsup(velocity))
61     print "dt", dt
62    
63     #displace the mesh
64     displacement = velocity * dt
65     coordinates = mesh.getX()
66     mesh.setX(coordinates + displacement)
67    
68     time += dt
69    
70     vel_mag = length(velocity)
71    
72     #save velocity and pressure output
73     saveVTK("vel.%2.2i.vtu"%(t),vel=vel_mag,vec=velocity,pressure=pressure)
74     t = t+1.0

Properties

Name Value
svn:mergeinfo

  ViewVC Help
Powered by ViewVC 1.1.26