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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 lgraham 2129 from esys.escript import *
2     from esys.escript.models import StokesProblemCartesian
3 jfenwick 3087 from esys.dudley import Rectangle
4 lgraham 2129
5 gross 2502
6 lgraham 2129 #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