/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 8756 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

  ViewVC Help
Powered by ViewVC 1.1.26