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

Annotation of /branches/symbolic_from_3470/dudley/test/python/rayleigh_taylor_instabilty.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (hide annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 3 months ago) by caltinay
File MIME type: text/x-python
File size: 8619 byte(s)
Fast forward to latest trunk revision 3788.

1 ksteube 1810
2     ########################################################
3 ksteube 1312 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1810 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1810 # 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 ksteube 1312 #
12 ksteube 1810 ########################################################
13 ksteube 1312
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1810 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1810
22    
23 gross 868 ################################################
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 jfenwick 3087 import esys.dudley
36     from esys.dudley import dudley
37 gross 868 from esys.escript.linearPDEs import LinearPDE
38 gross 873 from esys.escript.pdetools import Projector, SaddlePointProblem
39 caltinay 3346 from esys.weipa import saveVTK
40 gross 868 import sys
41     import math
42    
43     ### DEFINITION OF THE DOMAIN ###
44     l0=1.
45     l1=1.
46 gross 877 n0=10 # IDEALLY 80...
47     n1=10 # IDEALLY 80...
48 jfenwick 3087 mesh=esys.dudley.Brick(l0=l0, l1=l1, l2=l0, order=2, n0=n0, n1=n1, n2=n0)
49 gross 868
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 gross 873 velocityPDE = LinearPDE(mesh, numEquations=numDim)
67    
68 gross 868 advectPDE = LinearPDE(mesh)
69     advectPDE.setReducedOrderOn()
70     advectPDE.setValue(D=1.0)
71     advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)
72 gross 873
73 gross 868 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 gross 873 velocity = Vector(0.0, ContinuousFunction(mesh))
93 gross 868
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 caltinay 3789 print("Advection step done")
115 gross 868 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 caltinay 3789 print("Reinitialisation started.")
127 gross 868 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 caltinay 3789 print("Reinitialisation iteration :", iter, " error:", error)
139 gross 868 previous = phi
140     iter +=1
141 caltinay 3789 print("Reinitialisation finalized.")
142 gross 868 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 gross 873 class StokesProblem(SaddlePointProblem):
160     """
161     simple example of saddle point problem
162     """
163 gross 877 def __init__(self,domain,debug=False):
164     super(StokesProblem, self).__init__(self,debug)
165 gross 873 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 gross 877 self.__pde_p.setValue(D=1/self.eta)
181 gross 873 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 gross 877 sol=StokesProblem(velocity.getDomain(),debug=True)
199 gross 873 def solve_vel_uszawa(rho, eta, velocity, pressure):
200 gross 868 ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
201 gross 873 Y = Vector(0.0,Function(mesh))
202     Y[1] -= rho*g
203     sol.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
204 gross 877 velocity,pressure=sol.solve(velocity,pressure,iter_max=100,tolerance=0.01) #,accepted_reduction=None)
205 gross 873 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 gross 868 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 gross 873 velocity = velocityPDE.getSolution()
231 gross 868 p_iter +=1
232     if p_iter >=500:
233 caltinay 3789 print("You're screwed...")
234 gross 868 sys.exit(1)
235    
236 gross 873 pressure -= penalty*eta*(trace(grad(velocity)))
237     error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))
238 caltinay 3789 print("\nPressure iteration number:", p_iter)
239     print("error", error)
240 gross 868 ref = pressure*1.0
241    
242 gross 873 return velocity, pressure
243 gross 868
244     ### MAIN LOOP, OVER TIME ###
245     while t_step <= t_step_end:
246 caltinay 3789 print("######################")
247     print("Time step:", t_step)
248     print("######################")
249 gross 868 rho = update_parameter(phi, rho1, rho2)
250     eta = update_parameter(phi, eta1, eta2)
251    
252 gross 873 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 gross 868
256     ### PSEUDO POST-PROCESSING ###
257 caltinay 3789 print("########## Saving image", t_step, " ###########")
258 gross 873 saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)
259 gross 868
260     t_step += 1
261 gross 873
262     # vim: expandtab shiftwidth=4:

  ViewVC Help
Powered by ViewVC 1.1.26