# Contents of /trunk/finley/test/python/rayleigh_taylor_instabilty.py

Revision 868 - (show annotations)
Tue Oct 10 22:32:13 2006 UTC (14 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 5816 byte(s)
```new benchmark problem
```
 1 ################################################ 2 ## ## 3 ## October 2006 ## 4 ## ## 5 ## 3D Rayleigh-Taylor instability benchmark ## 6 ## by Laurent Bourgouin ## 7 ## ## 8 ################################################ 9 10 11 ### IMPORTS ### 12 from esys.escript import * 13 import esys.finley 14 from esys.finley import finley 15 from esys.escript.linearPDEs import LinearPDE 16 from esys.escript.pdetools import Projector 17 import sys 18 import math 19 20 ### DEFINITION OF THE DOMAIN ### 21 l0=1. 22 l1=1. 23 n0=20 # IDEALLY 80... 24 n1=20 # IDEALLY 80... 25 mesh=esys.finley.Brick(l0=l0, l1=l1, l2=l0, order=2, n0=n0, n1=n1, n2=n0) 26 27 ### PARAMETERS OF THE SIMULATION ### 28 rho1 = 1.0e3 # DENSITY OF THE FLUID AT THE BOTTOM 29 rho2 = 1.01e3 # DENSITY OF THE FLUID ON TOP 30 eta1 = 1.0e2 # VISCOSITY OF THE FLUID AT THE BOTTOM 31 eta2 = 1.0e2 # VISCOSITY OF THE FLUID ON TOP 32 penalty = 1.0e3 # PENALTY FACTOR FOT THE PENALTY METHOD 33 g=10. # GRAVITY 34 t_step = 0 35 t_step_end = 2000 36 reinit_max = 30 # NUMBER OF ITERATIONS DURING THE REINITIALISATION PROCEDURE 37 reinit_each = 3 # NUMBER OF TIME STEPS BETWEEN TWO REINITIALISATIONS 38 h = Lsup(mesh.getSize()) 39 numDim = mesh.getDim() 40 smooth = h*2.0 # SMOOTHING PARAMETER FOR THE TRANSITION ACROSS THE INTERFACE 41 42 ### DEFINITION OF THE PDE ### 43 velocityPDE = LinearPDE(mesh, numEquations=3) 44 velocityPDE.setSolverMethod(solver=LinearPDE.DIRECT) 45 advectPDE = LinearPDE(mesh) 46 advectPDE.setReducedOrderOn() 47 advectPDE.setValue(D=1.0) 48 advectPDE.setSolverMethod(solver=LinearPDE.DIRECT) 49 reinitPDE = LinearPDE(mesh, numEquations=1) 50 reinitPDE.setReducedOrderOn() 51 reinitPDE.setSolverMethod(solver=LinearPDE.LUMPING) 52 my_proj=Projector(mesh) 53 54 ### BOUNDARY CONDITIONS ### 55 xx = mesh.getX()[0] 56 yy = mesh.getX()[1] 57 zz = mesh.getX()[2] 58 top = whereZero(zz-l1) 59 bottom = whereZero(zz) 60 left = whereZero(xx) 61 right = whereZero(xx-l0) 62 front = whereZero(yy) 63 back = whereZero(yy-l0) 64 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] 65 velocityPDE.setValue(q = b_c) 66 67 pressure = Scalar(0.0, ContinuousFunction(mesh)) 68 69 ### INITIALISATION OF THE INTERFACE ### 70 func = -(-0.1*cos(math.pi*xx/l0)*cos(math.pi*yy/l0)-zz+0.4) 71 phi = func.interpolate(ReducedSolution(mesh)) 72 73 74 def advect(phi, velocity, dt): 75 ### SOLVES THE ADVECTION EQUATION ### 76 77 Y = phi.interpolate(Function(mesh)) 78 for i in range(numDim): 79 Y -= (dt/2.0)*velocity[i]*grad(phi)[i] 80 advectPDE.setValue(Y=Y) 81 phi_half = advectPDE.getSolution() 82 83 Y = phi 84 for i in range(numDim): 85 Y -= dt*velocity[i]*grad(phi_half)[i] 86 advectPDE.setValue(Y=Y) 87 phi = advectPDE.getSolution() 88 89 print "Advection step done" 90 return phi 91 92 def reinitialise(phi): 93 ### SOLVES THE REINITIALISATION EQUATION ### 94 s = sign(phi.interpolate(Function(mesh))) 95 w = s*grad(phi)/length(grad(phi)) 96 dtau = 0.3*h 97 iter =0 98 previous = 100.0 99 mask = whereNegative(abs(phi)-1.2*h) 100 reinitPDE.setValue(q=mask, r=phi) 101 while (iter<=reinit_max): 102 prod_scal =0.0 103 for i in range(numDim): 104 prod_scal += w[i]*grad(phi)[i] 105 coeff = s - prod_scal 106 ps2=0 107 for i in range(numDim): 108 ps2 += w[i]*grad(my_proj(coeff))[i] 109 reinitPDE.setValue(D=1.0, Y=phi+dtau*coeff-0.5*dtau**2*ps2) 110 phi = reinitPDE.getSolution() 111 error = Lsup((previous-phi)*whereNegative(abs(phi)-3.0*h))/h 112 print "Reinitialisation iteration :", iter, " error:", error 113 previous = phi 114 iter +=1 115 return phi 116 117 def update_phi(phi, velocity, dt, t_step): 118 ### CALLS THE ADVECTION PROCEDURE AND THE REINITIALISATION IF NECESSARY ### 119 phi=advect(phi, velocity, dt) 120 if t_step%reinit_each ==0: 121 phi = reinitialise(phi) 122 return phi 123 124 def update_parameter(phi, param_neg, param_pos): 125 ### UPDATES THE PARAMETERS TABLE USING THE SIGN OF PHI, A SMOOTH TRANSITION IS DONE ACROSS THE INTERFACE ### 126 mask_neg = whereNonNegative(-phi-smooth) 127 mask_pos = whereNonNegative(phi-smooth) 128 mask_interface = whereNegative(abs(phi)-smooth) 129 param = param_pos*mask_pos + param_neg*mask_neg + ((param_pos+param_neg)/2 +(param_pos-param_neg)*phi/(2.*smooth))*mask_interface 130 return param 131 132 def solve_vel(rho, eta, pressure): 133 ### SOLVES THE VELOCITY PROBLEM USING A PENALTY METHOD FOR THE INCOMPRESSIBILITY ### 134 error = 1.0 135 ref = pressure*1.0 136 p_iter=0 137 while (error >= 1.0e-2): 138 139 A=Tensor4(0.0, Function(mesh)) 140 for i in range(numDim): 141 for j in range(numDim): 142 A[i,j,i,j] += eta 143 A[i,j,j,i] += eta 144 A[i,i,j,j] += penalty*eta 145 146 Y = Vector(0.0,Function(mesh)) 147 Y[1] -= rho*g 148 149 X = Tensor(0.0, Function(mesh)) 150 for i in range(numDim): 151 X[i,i] += pressure 152 153 velocityPDE.setValue(A=A, X=X, Y=Y) 154 velocity_new = velocityPDE.getSolution() 155 p_iter +=1 156 if p_iter >=500: 157 print "You're screwed..." 158 sys.exit(1) 159 160 pressure -= penalty*eta*(trace(grad(velocity_new))) 161 error = penalty*Lsup(trace(grad(velocity_new)))/Lsup(grad(velocity_new)) 162 print "\nPressure iteration number:", p_iter 163 print "error", error 164 ref = pressure*1.0 165 166 return velocity_new, pressure 167 168 ### MAIN LOOP, OVER TIME ### 169 while t_step <= t_step_end: 170 171 rho = update_parameter(phi, rho1, rho2) 172 eta = update_parameter(phi, eta1, eta2) 173 174 velocity_new, pressure = solve_vel(rho, eta, pressure) 175 dt = 0.3*Lsup(mesh.getSize())/Lsup(velocity_new) 176 phi = update_phi(phi, velocity_new, dt, t_step) 177 178 ### PSEUDO POST-PROCESSING ### 179 print "########## Saving image", t_step, " ###########" 180 phi.saveVTK("/home/laurent/results2006/instability/phi3D.%2.2i.vtk" % t_step) 181 182 print "######################" 183 print "Time step:", t_step 184 print "######################" 185 t_step += 1

 ViewVC Help Powered by ViewVC 1.1.26