# Contents of /trunk/dudley/test/python/RT2D.py

Revision 2129 - (show annotations)
Thu Dec 4 03:53:02 2008 UTC (10 years, 10 months ago) by lgraham
Original Path: trunk/doc/examples/RT2D.py
File MIME type: text/x-python
File size: 2649 byte(s)
```Added RT2D.py python script for levelset tutorial to examples directory.

```
 1 from esys.escript import * 2 import esys.finley 3 from esys.escript.models import StokesProblemCartesian 4 from esys.finley import finley 5 from esys.finley import Rectangle 6 from LevelSet import * 7 8 #physical properties 9 rho1 = 1000 #fluid density on bottom 10 rho2 = 1010 #fluid density on top 11 eta1 = 100.0 #fluid viscosity on bottom 12 eta2 = 100.0 #fluid viscosity on top 13 g=10.0 14 15 #solver settings 16 dt = 0.001 17 t_step = 0 18 t_step_end = 2000 19 TOL = 1.0e-5 20 max_iter=400 21 verbose=True 22 useUzawa=True 23 24 #define mesh 25 l0=0.9142 26 l1=1.0 27 n0=200 28 n1=200 29 30 mesh=Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1) 31 #get mesh dimensions 32 numDim = mesh.getDim() 33 #get element size 34 h = Lsup(mesh.getSize()) 35 36 #level set parameters 37 tolerance = 1.0e-6 38 reinit_max = 30 39 reinit_each = 3 40 alpha = 1 41 smooth = alpha*h 42 43 #boundary conditions 44 x = mesh.getX() 45 #left + bottom + right + top 46 b_c = whereZero(x[0])*[1.0,0.0] + whereZero(x[1])*[1.0,1.0] + whereZero(x[0]-l0)*[1.0,0.0] + whereZero(x[1]-l1)*[1.0,1.0] 47 48 velocity = Vector(0.0, ContinuousFunction(mesh)) 49 pressure = Scalar(0.0, ContinuousFunction(mesh)) 50 Y = Vector(0.0,Function(mesh)) 51 52 #define initial interface between fluids 53 xx = mesh.getX()[0] 54 yy = mesh.getX()[1] 55 func = Scalar(0.0, ContinuousFunction(mesh)) 56 h_interface = Scalar(0.0, ContinuousFunction(mesh)) 57 h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2) 58 func = yy - h_interface 59 func_new = func.interpolate(ReducedSolution(mesh)) 60 61 #Stokes Cartesian 62 solution=StokesProblemCartesian(mesh,debug=True) 63 solution.setTolerance(TOL) 64 solution.setSubProblemTolerance(TOL**2) 65 66 #level set 67 levelset = LevelSet(mesh, func_new, reinit_max, reinit_each, tolerance, smooth) 68 69 while t_step <= t_step_end: 70 #update density and viscosity 71 rho = levelset.update_parameter(rho1, rho2) 72 eta = levelset.update_parameter(eta1, eta2) 73 74 #get velocity and pressure of fluid 75 Y[1] = -rho*g 76 solution.initialize(fixed_u_mask=b_c,eta=eta,f=Y) 77 velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,useUzawa=useUzawa) 78 79 #update the interface 80 func = levelset.update_phi(velocity, dt, t_step) 81 82 print "##########################################################" 83 print "time step:", t_step, " completed with dt:", dt 84 print "Velocity: min =", inf(velocity), "max =", Lsup(velocity) 85 print "##########################################################" 86 87 #save interface, velocity and pressure 88 saveVTK("phi2D.%2.4i.vtu"%t_step,interface=func,velocity=velocity,pressure=pressure) 89 #courant condition 90 dt = 0.4*Lsup(mesh.getSize())/Lsup(velocity) 91 t_step += 1