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

Revision 3346 - (show annotations)
Fri Nov 12 01:19:02 2010 UTC (8 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 2566 byte(s)
```Replaced usage of esys.escript.util.saveVTK by weipa.saveVTK in all python
scripts.

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

## Properties

Name Value
svn:mergeinfo /branches/lapack2681/finley/test/python/RT2D.py:2682-2741 /branches/restext/finley/test/python/RT2D.py:2610-2624