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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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

Legend:
Removed from v.3675  
changed lines
  Added in v.3772

  ViewVC Help
Powered by ViewVC 1.1.26