/[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 5705 - (hide annotations)
Mon Jun 29 01:07:07 2015 UTC (4 years, 3 months ago) by sshaw
File MIME type: text/x-python
File size: 2523 byte(s)
no more tab characters in python files
1 jfenwick 3774 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 sshaw 5705 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 jfenwick 3774 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/4.0fordebian/dudley/test/python/RT2D.py:5567-5588 /branches/lapack2681/finley/test/python/RT2D.py:2682-2741 /branches/pasowrap/dudley/test/python/RT2D.py:3661-3674 /branches/py3_attempt2/dudley/test/python/RT2D.py:3871-3891 /branches/restext/finley/test/python/RT2D.py:2610-2624 /branches/ripleygmg_from_3668/dudley/test/python/RT2D.py:3669-3791 /branches/symbolic_from_3470/dudley/test/python/RT2D.py:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/test/python/RT2D.py:3517-3974 /release/4.0/dudley/test/python/RT2D.py:5380-5406 /trunk/ripley/test/python/dudley/test/python/RT2D.py:3480-3515

  ViewVC Help
Powered by ViewVC 1.1.26