/[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 2881 - (hide annotations)
Thu Jan 28 02:03:15 2010 UTC (9 years, 8 months ago) by jfenwick
File MIME type: text/x-python
File size: 2119 byte(s)
Don't panic.
Updating copyright stamps

1 jfenwick 2549 ########################################################
2     #
3 jfenwick 2881 # Copyright (c) 2008-2010 by University of Queensland
4 jfenwick 2549 # 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 lgraham 2191 ######## 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 gross 2793 mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20) # use linear macro elements for pressure
41 lgraham 2191 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 gross 2793 boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
52 lgraham 2191
53     #velocity and pressure vectors
54 gross 2793 velocity=Vector(0.0, Solution(mesh))
55     pressure=Scalar(0.0, ReducedSolution(mesh))
56 lgraham 2191
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 gross 2432 velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
68 lgraham 2191
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

  ViewVC Help
Powered by ViewVC 1.1.26