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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/x-python
File size: 2534 byte(s)
Merging dudley and scons updates from branches

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

  ViewVC Help
Powered by ViewVC 1.1.26