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

revision 3675 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
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.3675 changed lines Added in v.3774