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

Revision 2742 - (show annotations)
Thu Nov 12 06:03:37 2009 UTC (9 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 2054 byte(s)
```Merging changes from the lapack branch.

The inverse() operation has been moved into c++. [No lazy support for this operation yet.]
Optional Lapack support has been added for matrices larger than 3x3.
service0 is set to use mkl_lapack.

```
 1 ######################################################## 2 # 3 # Copyright (c) 2008-2009 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=2, n0=20, n1=20) 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] 52 53 #velocity and pressure vectors 54 velocity=Vector(0.0, ContinuousFunction(mesh)) 55 pressure=Scalar(0.0, ContinuousFunction(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