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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 872 by gross, Tue Oct 10 22:32:13 2006 UTC revision 873 by gross, Mon Oct 16 04:07:33 2006 UTC
# Line 13  from esys.escript import * Line 13  from esys.escript import *
13  import esys.finley  import esys.finley
14  from esys.finley import finley  from esys.finley import finley
15  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
16  from esys.escript.pdetools import Projector  from esys.escript.pdetools import Projector, SaddlePointProblem
17  import sys  import sys
18  import math  import math
19    
# Line 40  numDim = mesh.getDim() Line 40  numDim = mesh.getDim()
40  smooth = h*2.0       # SMOOTHING PARAMETER FOR THE TRANSITION ACROSS THE INTERFACE  smooth = h*2.0       # SMOOTHING PARAMETER FOR THE TRANSITION ACROSS THE INTERFACE
41    
42  ### DEFINITION OF THE PDE ###  ### DEFINITION OF THE PDE ###
43  velocityPDE = LinearPDE(mesh, numEquations=3)  velocityPDE = LinearPDE(mesh, numEquations=numDim)
44  velocityPDE.setSolverMethod(solver=LinearPDE.DIRECT)  
45  advectPDE = LinearPDE(mesh)  advectPDE = LinearPDE(mesh)
46  advectPDE.setReducedOrderOn()  advectPDE.setReducedOrderOn()
47  advectPDE.setValue(D=1.0)  advectPDE.setValue(D=1.0)
48  advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)  advectPDE.setSolverMethod(solver=LinearPDE.DIRECT)
49    
50  reinitPDE = LinearPDE(mesh, numEquations=1)  reinitPDE = LinearPDE(mesh, numEquations=1)
51  reinitPDE.setReducedOrderOn()  reinitPDE.setReducedOrderOn()
52  reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)  reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING)
# Line 65  b_c = (bottom+top)*[1.0, 1.0, 1.0] + (le Line 66  b_c = (bottom+top)*[1.0, 1.0, 1.0] + (le
66  velocityPDE.setValue(q = b_c)  velocityPDE.setValue(q = b_c)
67    
68  pressure = Scalar(0.0, ContinuousFunction(mesh))  pressure = Scalar(0.0, ContinuousFunction(mesh))
69    velocity = Vector(0.0, ContinuousFunction(mesh))
70    
71  ### INITIALISATION OF THE INTERFACE ###  ### INITIALISATION OF THE INTERFACE ###
72  func = -(-0.1*cos(math.pi*xx/l0)*cos(math.pi*yy/l0)-zz+0.4)  func = -(-0.1*cos(math.pi*xx/l0)*cos(math.pi*yy/l0)-zz+0.4)
# Line 98  def reinitialise(phi): Line 100  def reinitialise(phi):
100    previous = 100.0    previous = 100.0
101    mask = whereNegative(abs(phi)-1.2*h)    mask = whereNegative(abs(phi)-1.2*h)
102    reinitPDE.setValue(q=mask, r=phi)    reinitPDE.setValue(q=mask, r=phi)
103      print "Reinitialisation started."
104    while (iter<=reinit_max):    while (iter<=reinit_max):
105      prod_scal =0.0      prod_scal =0.0
106      for i in range(numDim):      for i in range(numDim):
# Line 112  def reinitialise(phi): Line 115  def reinitialise(phi):
115      print "Reinitialisation iteration :", iter, " error:", error      print "Reinitialisation iteration :", iter, " error:", error
116      previous = phi      previous = phi
117      iter +=1      iter +=1
118      print "Reinitialisation finalized."
119    return phi    return phi
120    
121  def update_phi(phi, velocity, dt, t_step):  def update_phi(phi, velocity, dt, t_step):
# Line 129  def update_parameter(phi, param_neg, par Line 133  def update_parameter(phi, param_neg, par
133    param = param_pos*mask_pos + param_neg*mask_neg + ((param_pos+param_neg)/2 +(param_pos-param_neg)*phi/(2.*smooth))*mask_interface    param = param_pos*mask_pos + param_neg*mask_neg + ((param_pos+param_neg)/2 +(param_pos-param_neg)*phi/(2.*smooth))*mask_interface
134    return param    return param
135    
136  def solve_vel(rho, eta, pressure):  class StokesProblem(SaddlePointProblem):
137          """
138          simple example of saddle point problem
139          """
140          def __init__(self,domain):
141             super(StokesProblem, self).__init__(self)
142             self.domain=domain
143             self.__pde_u=LinearPDE(domain,numEquations=self.domain.getDim(),numSolutions=self.domain.getDim())
144             self.__pde_u.setSymmetryOn()
145    
146             self.__pde_p=LinearPDE(domain)
147             self.__pde_p.setReducedOrderOn()
148             self.__pde_p.setSymmetryOn()
149    
150          def initialize(self,f=Data(),fixed_u_mask=Data(),eta=1):
151             self.eta=eta
152             A =self.__pde_u.createCoefficientOfGeneralPDE("A")
153             for i in range(self.domain.getDim()):
154               for j in range(self.domain.getDim()):
155                 A[i,j,j,i] += self.eta
156                 A[i,j,i,j] += self.eta
157             self.__pde_p.setValue(D=1./self.eta)
158             self.__pde_u.setValue(A=A,q=fixed_u_mask,Y=f)
159    
160          def inner(self,p0,p1):
161             return integrate(p0*p1,Function(self.__pde_p.getDomain()))
162    
163          def solve_f(self,u,p,tol=1.e-8):
164             self.__pde_u.setTolerance(tol)
165             g=grad(u)
166             self.__pde_u.setValue(X=self.eta*symmetric(g)+p*kronecker(self.__pde_u.getDomain()))
167             return  self.__pde_u.getSolution()
168    
169          def solve_g(self,u,tol=1.e-8):
170             self.__pde_p.setTolerance(tol)
171             self.__pde_p.setValue(X=-u)
172             dp=self.__pde_p.getSolution()
173             return  dp
174    
175    def solve_vel_uszawa(rho, eta, velocity, pressure):
176  ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###  ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
177      sol=StokesProblem(velocity.getDomain())
178      Y = Vector(0.0,Function(mesh))
179      Y[1] -= rho*g
180      sol.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
181      velocity,pressure=sol.solve(velocity,pressure,relaxation=1.,iter_max=50,tolerance=0.01)
182      return velocity, pressure
183      
184    def solve_vel_penalty(rho, eta, velocity, pressure):
185    ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ###
186      velocityPDE.setSolverMethod(solver=LinearPDE.DIRECT)
187    error = 1.0    error = 1.0
188    ref = pressure*1.0    ref = pressure*1.0
189    p_iter=0    p_iter=0
# Line 151  def solve_vel(rho, eta, pressure): Line 204  def solve_vel(rho, eta, pressure):
204        X[i,i] += pressure        X[i,i] += pressure
205            
206      velocityPDE.setValue(A=A, X=X, Y=Y)      velocityPDE.setValue(A=A, X=X, Y=Y)
207      velocity_new = velocityPDE.getSolution()      velocity = velocityPDE.getSolution()
208      p_iter +=1      p_iter +=1
209      if p_iter >=500:      if p_iter >=500:
210        print "You're screwed..."        print "You're screwed..."
211        sys.exit(1)            sys.exit(1)    
212            
213      pressure -= penalty*eta*(trace(grad(velocity_new)))      pressure -= penalty*eta*(trace(grad(velocity)))
214      error = penalty*Lsup(trace(grad(velocity_new)))/Lsup(grad(velocity_new))      error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))
215      print "\nPressure iteration number:", p_iter      print "\nPressure iteration number:", p_iter
216      print "error", error      print "error", error
217      ref = pressure*1.0      ref = pressure*1.0
218            
219    return velocity_new, pressure    return velocity, pressure
220        
221  ### MAIN LOOP, OVER TIME ###  ### MAIN LOOP, OVER TIME ###
222  while t_step <= t_step_end:  while t_step <= t_step_end:
223      print "######################"
224      print "Time step:", t_step
225      print "######################"
226    rho = update_parameter(phi, rho1, rho2)    rho = update_parameter(phi, rho1, rho2)
227    eta = update_parameter(phi, eta1, eta2)    eta = update_parameter(phi, eta1, eta2)
228    
229    velocity_new, pressure = solve_vel(rho, eta, pressure)    velocity, pressure = solve_vel_uszawa(rho, eta,  velocity, pressure)
230    dt = 0.3*Lsup(mesh.getSize())/Lsup(velocity_new)    dt = 0.3*Lsup(mesh.getSize())/Lsup(velocity)
231    phi = update_phi(phi, velocity_new, dt, t_step)    phi = update_phi(phi, velocity, dt, t_step)
232    
233  ### PSEUDO POST-PROCESSING ###  ### PSEUDO POST-PROCESSING ###
234    print "##########  Saving image", t_step, " ###########"    print "##########  Saving image", t_step, " ###########"
235    phi.saveVTK("/home/laurent/results2006/instability/phi3D.%2.2i.vtk" % t_step)      saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)  
236    
   print "######################"  
   print "Time step:", t_step  
   print "######################"  
237    t_step += 1    t_step += 1
238    
239    # vim: expandtab shiftwidth=4:
Legend:
Removed from v.872  
changed lines
  Added in v.873
  ViewVC Help
Powered by ViewVC 1.1.26