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

Contents of /trunk/finley/test/python/RT2D.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2484 - (show annotations)
Mon Jun 22 04:22:19 2009 UTC (9 years, 10 months ago) by gross
Original Path: trunk/doc/examples/usersguide/RT2D.py
File MIME type: text/x-python
File size: 2649 byte(s)
numarray removed from docu; Locator revised.
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

Properties

Name Value
svn:mergeinfo

  ViewVC Help
Powered by ViewVC 1.1.26