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

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 6 # 7 # Primary Business: Queensland, Australia 8 # Licensed under the Open Software License version 3.0 9 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