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

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

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

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

Legend:
Removed from v.3773  
changed lines
  Added in v.3774

  ViewVC Help
Powered by ViewVC 1.1.26