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

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

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

revision 3788 by caltinay, Tue Mar 15 04:23:54 2011 UTC revision 3789 by caltinay, Tue Jan 31 04:55:05 2012 UTC
# Line 111  def advect(phi, velocity, dt): Line 111  def advect(phi, velocity, dt):
111    advectPDE.setValue(Y=Y)        advectPDE.setValue(Y=Y)    
112    phi = advectPDE.getSolution()    phi = advectPDE.getSolution()
113    
114    print "Advection step done"    print("Advection step done")
115    return phi    return phi
116    
117  def reinitialise(phi):  def reinitialise(phi):
# Line 123  def reinitialise(phi): Line 123  def reinitialise(phi):
123    previous = 100.0    previous = 100.0
124    mask = whereNegative(abs(phi)-1.2*h)    mask = whereNegative(abs(phi)-1.2*h)
125    reinitPDE.setValue(q=mask, r=phi)    reinitPDE.setValue(q=mask, r=phi)
126    print "Reinitialisation started."    print("Reinitialisation started.")
127    while (iter<=reinit_max):    while (iter<=reinit_max):
128      prod_scal =0.0      prod_scal =0.0
129      for i in range(numDim):      for i in range(numDim):
# Line 135  def reinitialise(phi): Line 135  def reinitialise(phi):
135      reinitPDE.setValue(D=1.0, Y=phi+dtau*coeff-0.5*dtau**2*ps2)      reinitPDE.setValue(D=1.0, Y=phi+dtau*coeff-0.5*dtau**2*ps2)
136      phi = reinitPDE.getSolution()      phi = reinitPDE.getSolution()
137      error = Lsup((previous-phi)*whereNegative(abs(phi)-3.0*h))/h      error = Lsup((previous-phi)*whereNegative(abs(phi)-3.0*h))/h
138      print "Reinitialisation iteration :", iter, " error:", error      print("Reinitialisation iteration :", iter, " error:", error)
139      previous = phi      previous = phi
140      iter +=1      iter +=1
141    print "Reinitialisation finalized."    print("Reinitialisation finalized.")
142    return phi    return phi
143    
144  def update_phi(phi, velocity, dt, t_step):  def update_phi(phi, velocity, dt, t_step):
# Line 230  def solve_vel_penalty(rho, eta, velocity Line 230  def solve_vel_penalty(rho, eta, velocity
230      velocity = velocityPDE.getSolution()      velocity = velocityPDE.getSolution()
231      p_iter +=1      p_iter +=1
232      if p_iter >=500:      if p_iter >=500:
233        print "You're screwed..."        print("You're screwed...")
234        sys.exit(1)            sys.exit(1)    
235            
236      pressure -= penalty*eta*(trace(grad(velocity)))      pressure -= penalty*eta*(trace(grad(velocity)))
237      error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))      error = penalty*Lsup(trace(grad(velocity)))/Lsup(grad(velocity))
238      print "\nPressure iteration number:", p_iter      print("\nPressure iteration number:", p_iter)
239      print "error", error      print("error", error)
240      ref = pressure*1.0      ref = pressure*1.0
241            
242    return velocity, pressure    return velocity, pressure
243        
244  ### MAIN LOOP, OVER TIME ###  ### MAIN LOOP, OVER TIME ###
245  while t_step <= t_step_end:  while t_step <= t_step_end:
246    print "######################"    print("######################")
247    print "Time step:", t_step    print("Time step:", t_step)
248    print "######################"    print("######################")
249    rho = update_parameter(phi, rho1, rho2)    rho = update_parameter(phi, rho1, rho2)
250    eta = update_parameter(phi, eta1, eta2)    eta = update_parameter(phi, eta1, eta2)
251    
# Line 254  while t_step <= t_step_end: Line 254  while t_step <= t_step_end:
254    phi = update_phi(phi, velocity, dt, t_step)    phi = update_phi(phi, velocity, dt, t_step)
255    
256  ### PSEUDO POST-PROCESSING ###  ### PSEUDO POST-PROCESSING ###
257    print "##########  Saving image", t_step, " ###########"    print("##########  Saving image", t_step, " ###########")
258    saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)      saveVTK("phi3D.%2.2i.vtk"%t_step,layer=phi)  
259    
260    t_step += 1    t_step += 1

Legend:
Removed from v.3788  
changed lines
  Added in v.3789

  ViewVC Help
Powered by ViewVC 1.1.26