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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (2 years, 1 month 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 ksteube 1810
2 jfenwick 3981 ##############################################################################
3 ksteube 1312 #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1312 #
7 ksteube 1810 # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1312 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 ksteube 1312
17 sshaw 5706 from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 ksteube 1810 Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1810
26    
27 gross 868 ################################################
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 caltinay 3346 from esys.weipa import saveVTK
42 gross 868 from esys.escript.linearPDEs import LinearPDE
43 gross 873 from esys.escript.pdetools import Projector, SaddlePointProblem
44 gross 868 import sys
45     import math
46    
47     ### DEFINITION OF THE DOMAIN ###
48     l0=1.
49     l1=1.
50 gross 877 n0=10 # IDEALLY 80...
51     n1=10 # IDEALLY 80...
52 gross 868 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 gross 873 velocityPDE = LinearPDE(mesh, numEquations=numDim)
71    
72 gross 868 advectPDE = LinearPDE(mesh)
73     advectPDE.setReducedOrderOn()
74     advectPDE.setValue(D=1.0)
75     advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)
76 gross 873
77 gross 868 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 gross 873 velocity = Vector(0.0, ContinuousFunction(mesh))
97 gross 868
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 jfenwick 3772 print("Advection step done")
119 gross 868 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 jfenwick 3772 print("Reinitialisation started.")
131 gross 868 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 jfenwick 3772 print("Reinitialisation iteration :", iter, " error:", error)
143 gross 868 previous = phi
144     iter +=1
145 jfenwick 3772 print("Reinitialisation finalized.")
146 gross 868 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 gross 873 class StokesProblem(SaddlePointProblem):
164     """
165     simple example of saddle point problem
166     """
167 gross 877 def __init__(self,domain,debug=False):
168     super(StokesProblem, self).__init__(self,debug)
169 gross 873 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 gross 877 self.__pde_p.setValue(D=1/self.eta)
185 gross 873 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 gross 877 sol=StokesProblem(velocity.getDomain(),debug=True)
203 gross 873 def solve_vel_uszawa(rho, eta, velocity, pressure):
204 gross 868 ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
205 gross 873 Y = Vector(0.0,Function(mesh))
206     Y[1] -= rho*g
207     sol.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
208 gross 877 velocity,pressure=sol.solve(velocity,pressure,iter_max=100,tolerance=0.01) #,accepted_reduction=None)
209 gross 873 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 gross 868 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 sshaw 4576 A[i,i,j,j] += penalty*eta
225 gross 868
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 gross 873 velocity = velocityPDE.getSolution()
235 gross 868 p_iter +=1
236     if p_iter >=500:
237 jfenwick 3772 print("You're screwed...")
238 gross 868 sys.exit(1)
239    
240 gross 873 pressure -= penalty*eta*(trace(grad(velocity)))
241     error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))
242 jfenwick 3772 print("\nPressure iteration number:", p_iter)
243     print("error", error)
244 gross 868 ref = pressure*1.0
245    
246 gross 873 return velocity, pressure
247 gross 868
248     ### MAIN LOOP, OVER TIME ###
249     while t_step <= t_step_end:
250 jfenwick 3772 print("######################")
251     print("Time step:", t_step)
252     print("######################")
253 gross 868 rho = update_parameter(phi, rho1, rho2)
254     eta = update_parameter(phi, eta1, eta2)
255    
256 gross 873 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 gross 868
260     ### PSEUDO POST-PROCESSING ###
261 jfenwick 3772 print("########## Saving image", t_step, " ###########")
262 gross 873 saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)
263 gross 868
264     t_step += 1
265 gross 873
266     # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26