/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 8793 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

  ViewVC Help
Powered by ViewVC 1.1.26