/[escript]/branches/stage3.1/doc/examples/usersguide/fluid.py
ViewVC logotype

Contents of /branches/stage3.1/doc/examples/usersguide/fluid.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2945 - (show annotations)
Wed Feb 24 00:17:46 2010 UTC (9 years ago) by jfenwick
File MIME type: text/x-python
File size: 2119 byte(s)
Bringing release stage up to trunk version 2944

1 ########################################################
2 #
3 # Copyright (c) 2008-2010 by University of Queensland
4 # Earth Systems Science Computational Center (ESSCC)
5 # http://www.uq.edu.au/esscc
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 ########################################################
12
13 ######## August 2008 ########
14 ########## Leon Graham ##########
15 ## Newtonian fluid using StokesProblemCartesian class##
16
17 from esys.escript import *
18 import esys.finley
19 from esys.escript.linearPDEs import LinearPDE
20 from esys.escript.models import StokesProblemCartesian
21
22 #physical constants
23 eta=1.0
24 rho=100.0
25 g=10.0
26
27 #solver settings
28 tolerance=1.0e-4
29 max_iter=200
30 t_end=50
31 t=0.0
32 time=0
33 verbose='TRUE'
34 useUzawa='TRUE'
35
36 #define mesh
37 H=2.0
38 L=1.0
39 W=1.0
40 mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20) # use linear macro elements for pressure
41 coordinates = mesh.getX()
42
43 #gravitational force
44 Y=Vector(0.0, Function(mesh))
45 Y[1]=-rho*g
46
47 #element spacing
48 h=Lsup(mesh.getSize())
49
50 #boundary conditions for slip at base
51 boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
52
53 #velocity and pressure vectors
54 velocity=Vector(0.0, Solution(mesh))
55 pressure=Scalar(0.0, ReducedSolution(mesh))
56
57 #Stokes Cartesian
58 solution=StokesProblemCartesian(mesh)
59 solution.setTolerance(tolerance)
60
61 while t <= t_end:
62
63 print " ----- Time step = %s -----"%( t )
64 print "Time = %s seconds"%( time )
65
66 solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)
67 velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
68
69 print "Max velocity =", Lsup(velocity), "m/s"
70
71 #Courant condition
72 dt=0.4*h/(Lsup(velocity))
73 print "dt", dt
74
75 #displace the mesh
76 displacement = velocity * dt
77 coordinates = mesh.getX()
78 mesh.setX(coordinates + displacement)
79
80 time += dt
81
82 vel_mag = length(velocity)
83
84 #save velocity and pressure output
85 saveVTK("vel.%2.2i.vtu"%(t),vel=vel_mag,vec=velocity,pressure=pressure)
86 t = t+1.0

Properties

Name Value
svn:mergeinfo /branches/lapack2681/doc/examples/usersguide/fluid.py:2682-2741 /branches/restext/doc/examples/usersguide/fluid.py:2610-2624 /trunk/doc/examples/usersguide/fluid.py:2872-2919,2924-2944

  ViewVC Help
Powered by ViewVC 1.1.26