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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3675 - (hide annotations)
Thu Nov 17 00:53:38 2011 UTC (7 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 2618 byte(s)
pasowrap joins the trunk.

1 lgraham 2129 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 caltinay 3346 from esys.weipa import saveVTK
7 lgraham 2129
8 gross 2502
9 lgraham 2129 #physical properties
10     rho1 = 1000 #fluid density on bottom
11     rho2 = 1010 #fluid density on top
12     eta1 = 100.0 #fluid viscosity on bottom
13     eta2 = 100.0 #fluid viscosity on top
14     g=10.0
15    
16     #solver settings
17     dt = 0.001
18     t_step = 0
19     t_step_end = 2000
20     TOL = 1.0e-5
21     max_iter=400
22     verbose=True
23     useUzawa=True
24    
25     #define mesh
26     l0=0.9142
27     l1=1.0
28     n0=200
29     n1=200
30    
31     mesh=Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1)
32     #get mesh dimensions
33     numDim = mesh.getDim()
34     #get element size
35     h = Lsup(mesh.getSize())
36    
37     #level set parameters
38     tolerance = 1.0e-6
39     reinit_max = 30
40     reinit_each = 3
41     alpha = 1
42     smooth = alpha*h
43    
44     #boundary conditions
45     x = mesh.getX()
46     #left + bottom + right + top
47     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]
48    
49     velocity = Vector(0.0, ContinuousFunction(mesh))
50     pressure = Scalar(0.0, ContinuousFunction(mesh))
51     Y = Vector(0.0,Function(mesh))
52    
53     #define initial interface between fluids
54     xx = mesh.getX()[0]
55     yy = mesh.getX()[1]
56     func = Scalar(0.0, ContinuousFunction(mesh))
57     h_interface = Scalar(0.0, ContinuousFunction(mesh))
58     h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2)
59     func = yy - h_interface
60     func_new = func.interpolate(ReducedSolution(mesh))
61    
62     #Stokes Cartesian
63     solution=StokesProblemCartesian(mesh,debug=True)
64     solution.setTolerance(TOL)
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 /branches/lapack2681/finley/test/python/RT2D.py:2682-2741 /branches/pasowrap/finley/test/python/RT2D.py:3661-3674 /branches/restext/finley/test/python/RT2D.py:2610-2624

  ViewVC Help
Powered by ViewVC 1.1.26