/[escript]/trunk/finley/test/python/rayleigh_taylor_instabilty.py
ViewVC logotype

Contents of /trunk/finley/test/python/rayleigh_taylor_instabilty.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1312 - (show annotations)
Mon Sep 24 06:18:44 2007 UTC (12 years, 6 months ago) by ksteube
File MIME type: text/x-python
File size: 8244 byte(s)
The MPI branch is hereby closed. All future work should be in trunk.

Previously in revision 1295 I merged the latest changes to trunk into trunk-mpi-branch.
In this revision I copied all files from trunk-mpi-branch over the corresponding
trunk files. I did not use 'svn merge', it was a copy.

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

  ViewVC Help
Powered by ViewVC 1.1.26