← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/rayleigh_taylor_instabilty.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/rayleigh_taylor_instabilty.py
 
1
 
 
2
########################################################
 
3
#
 
4
# Copyright (c) 2003-2010 by University of Queensland
 
5
# Earth Systems Science Computational Center (ESSCC)
 
6
# http://www.uq.edu.au/esscc
 
7
#
 
8
# Primary Business: Queensland, Australia
 
9
# Licensed under the Open Software License version 3.0
 
10
# http://www.opensource.org/licenses/osl-3.0.php
 
11
#
 
12
########################################################
 
13
 
 
14
__copyright__="""Copyright (c) 2003-2010 by University of Queensland
 
15
Earth Systems Science Computational Center (ESSCC)
 
16
http://www.uq.edu.au/esscc
 
17
Primary Business: Queensland, Australia"""
 
18
__license__="""Licensed under the Open Software License version 3.0
 
19
http://www.opensource.org/licenses/osl-3.0.php"""
 
20
__url__="https://launchpad.net/escript-finley"
 
21
 
 
22
 
 
23
################################################
 
24
##                                            ##
 
25
## October 2006                               ## 
 
26
##                                            ##
 
27
##  3D Rayleigh-Taylor instability benchmark  ##
 
28
##           by Laurent Bourgouin             ##
 
29
##                                            ##
 
30
################################################
 
31
 
 
32
 
 
33
### IMPORTS ###
 
34
from esys.escript import *
 
35
import esys.dudley
 
36
from esys.dudley import dudley
 
37
from esys.escript.linearPDEs import LinearPDE
 
38
from esys.escript.pdetools import Projector, SaddlePointProblem
 
39
import sys
 
40
import math
 
41
 
 
42
### DEFINITION OF THE DOMAIN ###
 
43
l0=1.
 
44
l1=1.
 
45
n0=10  # IDEALLY 80...
 
46
n1=10  # IDEALLY 80...
 
47
mesh=esys.dudley.Brick(l0=l0, l1=l1, l2=l0, order=2, n0=n0, n1=n1, n2=n0)
 
48
 
 
49
### PARAMETERS OF THE SIMULATION ###
 
50
rho1 = 1.0e3         # DENSITY OF THE FLUID AT THE BOTTOM
 
51
rho2 = 1.01e3        # DENSITY OF THE FLUID ON TOP
 
52
eta1 = 1.0e2         # VISCOSITY OF THE FLUID AT THE BOTTOM
 
53
eta2 = 1.0e2         # VISCOSITY OF THE FLUID ON TOP
 
54
penalty = 1.0e3      # PENALTY FACTOR FOT THE PENALTY METHOD
 
55
g=10.                # GRAVITY
 
56
t_step = 0
 
57
t_step_end = 2000
 
58
reinit_max = 30      # NUMBER OF ITERATIONS DURING THE REINITIALISATION PROCEDURE
 
59
reinit_each = 3      # NUMBER OF TIME STEPS BETWEEN TWO REINITIALISATIONS
 
60
h = Lsup(mesh.getSize())
 
61
numDim = mesh.getDim()
 
62
smooth = h*2.0       # SMOOTHING PARAMETER FOR THE TRANSITION ACROSS THE INTERFACE
 
63
 
 
64
### DEFINITION OF THE PDE ###
 
65
velocityPDE = LinearPDE(mesh, numEquations=numDim)
 
66
 
 
67
advectPDE = LinearPDE(mesh)
 
68
advectPDE.setReducedOrderOn()
 
69
advectPDE.setValue(D=1.0)
 
70
advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)
 
71
 
 
72
reinitPDE = LinearPDE(mesh, numEquations=1)
 
73
reinitPDE.setReducedOrderOn()
 
74
reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)
 
75
my_proj=Projector(mesh)
 
76
 
 
77
### BOUNDARY CONDITIONS ###
 
78
xx = mesh.getX()[0]
 
79
yy = mesh.getX()[1]
 
80
zz = mesh.getX()[2]
 
81
top = whereZero(zz-l1)
 
82
bottom = whereZero(zz)
 
83
left = whereZero(xx)
 
84
right = whereZero(xx-l0)
 
85
front = whereZero(yy)
 
86
back = whereZero(yy-l0)
 
87
b_c = (bottom+top)*[1.0, 1.0, 1.0] + (left+right)*[1.0,0.0, 0.0] + (front+back)*[0.0, 1.0, 0.0]
 
88
velocityPDE.setValue(q = b_c)
 
89
 
 
90
pressure = Scalar(0.0, ContinuousFunction(mesh))
 
91
velocity = Vector(0.0, ContinuousFunction(mesh))
 
92
 
 
93
### INITIALISATION OF THE INTERFACE ###
 
94
func = -(-0.1*cos(math.pi*xx/l0)*cos(math.pi*yy/l0)-zz+0.4)
 
95
phi = func.interpolate(ReducedSolution(mesh))
 
96
 
 
97
 
 
98
def advect(phi, velocity, dt):
 
99
### SOLVES THE ADVECTION EQUATION ###
 
100
 
 
101
  Y = phi.interpolate(Function(mesh))
 
102
  for i in range(numDim):
 
103
    Y -= (dt/2.0)*velocity[i]*grad(phi)[i]
 
104
  advectPDE.setValue(Y=Y)    
 
105
  phi_half = advectPDE.getSolution()
 
106
 
 
107
  Y = phi
 
108
  for i in range(numDim):
 
109
    Y -= dt*velocity[i]*grad(phi_half)[i]
 
110
  advectPDE.setValue(Y=Y)    
 
111
  phi = advectPDE.getSolution()
 
112
 
 
113
  print "Advection step done"
 
114
  return phi
 
115
 
 
116
def reinitialise(phi):
 
117
### SOLVES THE REINITIALISATION EQUATION ###
 
118
  s = sign(phi.interpolate(Function(mesh)))
 
119
  w = s*grad(phi)/length(grad(phi))
 
120
  dtau = 0.3*h
 
121
  iter =0
 
122
  previous = 100.0
 
123
  mask = whereNegative(abs(phi)-1.2*h)
 
124
  reinitPDE.setValue(q=mask, r=phi)
 
125
  print "Reinitialisation started."
 
126
  while (iter<=reinit_max):
 
127
    prod_scal =0.0
 
128
    for i in range(numDim):
 
129
      prod_scal += w[i]*grad(phi)[i]
 
130
    coeff = s - prod_scal
 
131
    ps2=0
 
132
    for i in range(numDim):
 
133
      ps2 += w[i]*grad(my_proj(coeff))[i]
 
134
    reinitPDE.setValue(D=1.0, Y=phi+dtau*coeff-0.5*dtau**2*ps2)
 
135
    phi = reinitPDE.getSolution()
 
136
    error = Lsup((previous-phi)*whereNegative(abs(phi)-3.0*h))/h
 
137
    print "Reinitialisation iteration :", iter, " error:", error
 
138
    previous = phi
 
139
    iter +=1
 
140
  print "Reinitialisation finalized."
 
141
  return phi
 
142
 
 
143
def update_phi(phi, velocity, dt, t_step):
 
144
### CALLS THE ADVECTION PROCEDURE AND THE REINITIALISATION IF NECESSARY ###  
 
145
  phi=advect(phi, velocity, dt)
 
146
  if t_step%reinit_each ==0:
 
147
    phi = reinitialise(phi)
 
148
  return phi
 
149
  
 
150
def update_parameter(phi, param_neg, param_pos):
 
151
### UPDATES THE PARAMETERS TABLE USING THE SIGN OF PHI, A SMOOTH TRANSITION IS DONE ACROSS THE INTERFACE ###
 
152
  mask_neg = whereNonNegative(-phi-smooth)
 
153
  mask_pos = whereNonNegative(phi-smooth)
 
154
  mask_interface = whereNegative(abs(phi)-smooth)
 
155
  param = param_pos*mask_pos + param_neg*mask_neg + ((param_pos+param_neg)/2 +(param_pos-param_neg)*phi/(2.*smooth))*mask_interface
 
156
  return param
 
157
 
 
158
class StokesProblem(SaddlePointProblem):
 
159
      """
 
160
      simple example of saddle point problem
 
161
      """
 
162
      def __init__(self,domain,debug=False):
 
163
         super(StokesProblem, self).__init__(self,debug)
 
164
         self.domain=domain
 
165
         self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
 
166
         self.__pde_u.setSymmetryOn()
 
167
 
 
168
         self.__pde_p=LinearPDE(domain)
 
169
         self.__pde_p.setReducedOrderOn()
 
170
         self.__pde_p.setSymmetryOn()
 
171
 
 
172
      def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
 
173
         self.eta=eta
 
174
         A =self.__pde_u.createCoefficientOfGeneralPDE("A")
 
175
         for i in range(self.domain.getDim()):
 
176
           for j in range(self.domain.getDim()):
 
177
             A[i,j,j,i] += self.eta
 
178
             A[i,j,i,j] += self.eta
 
179
         self.__pde_p.setValue(D=1/self.eta)
 
180
         self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
 
181
 
 
182
      def inner(self,p0,p1):
 
183
         return integrate(p0*p1,Function(self.__pde_p.getDomain()))
 
184
 
 
185
      def solve_f(self,u,p,tol=1.e-8):
 
186
         self.__pde_u.setTolerance(tol)
 
187
         g=grad(u)
 
188
         self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
 
189
         return  self.__pde_u.getSolution()
 
190
 
 
191
      def solve_g(self,u,tol=1.e-8):
 
192
         self.__pde_p.setTolerance(tol)
 
193
         self.__pde_p.setValue(X=-u) 
 
194
         dp=self.__pde_p.getSolution()
 
195
         return  dp
 
196
 
 
197
sol=StokesProblem(velocity.getDomain(),debug=True)
 
198
def solve_vel_uszawa(rho, eta, velocity, pressure):
 
199
### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
 
200
  Y = Vector(0.0,Function(mesh))
 
201
  Y[1] -= rho*g
 
202
  sol.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
 
203
  velocity,pressure=sol.solve(velocity,pressure,iter_max=100,tolerance=0.01) #,accepted_reduction=None)
 
204
  return velocity, pressure
 
205
  
 
206
def solve_vel_penalty(rho, eta, velocity, pressure):
 
207
### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
 
208
  velocityPDE.setSolverMethod(solver=LinearPDE.DIRECT)
 
209
  error = 1.0
 
210
  ref = pressure*1.0
 
211
  p_iter=0
 
212
  while (error >= 1.0e-2):
 
213
  
 
214
    A=Tensor4(0.0, Function(mesh))
 
215
    for i in  range(numDim):
 
216
      for j in range(numDim):
 
217
        A[i,j,i,j] += eta
 
218
        A[i,j,j,i] += eta
 
219
        A[i,i,j,j] += penalty*eta
 
220
 
 
221
    Y = Vector(0.0,Function(mesh))
 
222
    Y[1] -= rho*g
 
223
    
 
224
    X = Tensor(0.0, Function(mesh))
 
225
    for i in range(numDim):
 
226
      X[i,i] += pressure
 
227
    
 
228
    velocityPDE.setValue(A=A, X=X, Y=Y)
 
229
    velocity = velocityPDE.getSolution()
 
230
    p_iter +=1
 
231
    if p_iter >=500:
 
232
      print "You're screwed..."
 
233
      sys.exit(1)    
 
234
    
 
235
    pressure -= penalty*eta*(trace(grad(velocity)))
 
236
    error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))
 
237
    print "\nPressure iteration number:", p_iter
 
238
    print "error", error
 
239
    ref = pressure*1.0
 
240
    
 
241
  return velocity, pressure
 
242
  
 
243
### MAIN LOOP, OVER TIME ###
 
244
while t_step <= t_step_end:
 
245
  print "######################"
 
246
  print "Time step:", t_step
 
247
  print "######################"
 
248
  rho = update_parameter(phi, rho1, rho2)
 
249
  eta = update_parameter(phi, eta1, eta2)
 
250
 
 
251
  velocity, pressure = solve_vel_uszawa(rho, eta,  velocity, pressure)
 
252
  dt = 0.3*Lsup(mesh.getSize())/Lsup(velocity)
 
253
  phi = update_phi(phi, velocity, dt, t_step)
 
254
 
 
255
### PSEUDO POST-PROCESSING ###
 
256
  print "##########  Saving image", t_step, " ###########" 
 
257
  saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)  
 
258
 
 
259
  t_step += 1
 
260
 
 
261
# vim: expandtab shiftwidth=4: