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

Contents of /trunk/doc/examples/fluid.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2432 - (show annotations)
Wed May 20 06:06:20 2009 UTC (9 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 1651 byte(s)
power law and incompressible flow is running now
1 ######## 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 velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
56
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

  ViewVC Help
Powered by ViewVC 1.1.26