examples/usersguide/fluid.py

########      August 2008      ########
##########    Leon Graham    ########## 
## Newtonian fluid using StokesProblemCartesian class##

from esys.escript import *
import esys.finley
from esys.escript.linearPDEs import LinearPDE
from esys.escript.models import StokesProblemCartesian

#physical constants
eta=1.0
rho=100.0
g=10.0 

#solver settings
tolerance=1.0e-4
max_iter=200
t_end=50
t=0.0
time=0
verbose='TRUE'
useUzawa='TRUE'

#define mesh 
H=2.0
L=1.0
W=1.0
mesh = esys.finley.Rectangle(l0=L, l1=H, order=2, n0=20, n1=20)
coordinates = mesh.getX()

#gravitational force
Y=Vector(0.0, Function(mesh))
Y[1]=-rho*g

#element spacing
h=Lsup(mesh.getSize())

#boundary conditions for slip at base
boundary_cond=whereZero(coordinates[1])*[0.0,1.0]

#velocity and pressure vectors
velocity=Vector(0.0, ContinuousFunction(mesh))
pressure=Scalar(0.0, ContinuousFunction(mesh))

#Stokes Cartesian
solution=StokesProblemCartesian(mesh)
solution.setTolerance(tolerance)

while t <= t_end:

  print " ----- Time step = %s -----"%( t )
  print "Time = %s seconds"%( time )  
 
  solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)
  velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
  
  print "Max velocity =", Lsup(velocity), "m/s"
  
  #Courant condition
  dt=0.4*h/(Lsup(velocity))
  print "dt", dt 
  
  #displace the mesh
  displacement = velocity * dt
  coordinates = mesh.getX()
  mesh.setX(coordinates + displacement)  
  
  time += dt
  
  vel_mag = length(velocity)

  #save velocity and pressure output
  saveVTK("vel.%2.2i.vtu"%(t),vel=vel_mag,vec=velocity,pressure=pressure)
  t = t+1.0
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