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

  ViewVC Help
Powered by ViewVC 1.1.26