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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3774 - (show annotations)
Wed Jan 18 06:29:34 2012 UTC (7 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 8619 byte(s)
dudley, pasowrap, pycad

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

  ViewVC Help
Powered by ViewVC 1.1.26