/[escript]/trunk/doc/user/levelset.tex
ViewVC logotype

Diff of /trunk/doc/user/levelset.tex

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

revision 2115 by lgraham, Tue Dec 2 01:21:07 2008 UTC revision 2120 by lgraham, Tue Dec 2 06:21:49 2008 UTC
# Line 240  The Rayleigh-Taylor instability problem Line 240  The Rayleigh-Taylor instability problem
240  \begin{figure}  \begin{figure}
241  \center  \center
242  \scalebox{0.7}{\includegraphics{figures/RT2Dsetup.eps}}  \scalebox{0.7}{\includegraphics{figures/RT2Dsetup.eps}}
243  \caption{Parameters, initial interface and boundary conditions for the Rayleigh-Taylor instability problem. The interface is defined as $\phi=0.02cos(\frac{\pi x}{\lambda}) + 0.2$. The fluids have been assigned different densities and equal viscosity (isoviscous).}  \caption{Parameters, initial interface and boundary conditions for the Rayleigh-Taylor instability problem. The interface is defined as $\phi=0.02cos(\frac{\pi x}{\lambda}) + 0.2$. The fluids have been assigned different densities and equal viscosity (isoviscous) \cite{BOURGOUIN2006}.}
244  \label{RT2DSETUP}  \label{RT2DSETUP}
245  \end{figure}  \end{figure}
246    
# Line 248  The Rayleigh-Taylor instability problem Line 248  The Rayleigh-Taylor instability problem
248    
249  \begin{python}  \begin{python}
250    
251    from esys.escript import *
252    import esys.finley
253    from esys.escript.models import StokesProblemCartesian
254    from esys.finley import finley
255    from LevelSet import *
256    
257    #physical properties
258    rho1 = 1000     #fluid density on bottom
259    rho2 = 1010     #fluid density on top
260    eta1 = 100.0        #fluid viscosity on bottom
261    eta2 = 100.0        #fluid viscosity on top
262    penalty = 100.0
263    g=10.0
264    
265    #solver settings
266    dt = 0.001
267    t_step = 0
268    t_step_end = 2000
269    TOL = 1.0e-5
270    max_iter=400
271    verbose=True
272    useUzawa=True
273    
274    #define mesh
275    l0=0.9142
276    l1=1.0
277    n0=100      
278    n1=100
279    
280    mesh=esys.finley.Rectangle(l0=l0, l1=l1, order=2, n0=n0, n1=n1)
281    #get mesh dimensions
282    numDim = mesh.getDim()
283    #get element size
284    h = Lsup(mesh.getSize())
285    print "element size",h
286    
287    #level set parameters
288    tolerance = 1.0e-6
289    reinit_max = 30
290    reinit_each = 3
291    alpha = 1
292    smooth = alpha*h
293    
294    #boundary conditions
295    x = mesh.getX()
296    #left + bottom + right + top
297    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]
298    
299    velocity = Vector(0.0, ContinuousFunction(mesh))
300    pressure = Scalar(0.0, ContinuousFunction(mesh))
301    Y = Vector(0.0,Function(mesh))
302    
303    #define initial interface between fluids
304    xx = mesh.getX()[0]
305    yy = mesh.getX()[1]
306    func = Scalar(0.0, ContinuousFunction(mesh))
307    h_interface = Scalar(0.0, ContinuousFunction(mesh))
308    h_interface = h_interface + (0.02*cos(math.pi*xx/l0) + 0.2)
309    func = yy - h_interface
310    func_new = func.interpolate(ReducedSolution(mesh))
311    
312    #Stokes cartesian
313    solution=StokesProblemCartesian(mesh,debug=True)
314    solution.setTolerance(TOL)
315    solution.setSubProblemTolerance(TOL**2)
316    
317    #level set
318    levelset = LevelSet(mesh, func_new, reinit_max, reinit_each, tolerance, smooth)    
319    
320    while t_step <= t_step_end:
321      #update density and viscosity
322      rho = levelset.update_parameter(rho1, rho2)
323      eta = levelset.update_parameter(eta1, eta2)
324    
325      #get velocity and pressue of fluid
326      Y[1] = -rho*g
327      solution.initialize(fixed_u_mask=b_c,eta=eta,f=Y)
328      velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,useUzawa=useUzawa)
329      
330      #update the interface
331      func = levelset.update_phi(velocity, dt, t_step)  
332    
333      print "##########################################################"
334      print "time step:", t_step, " completed with dt:", dt
335      print "Velocity: min =", inf(velocity), "max =", Lsup(velocity)
336      print "##########################################################"
337    
338      #save interface, velocity and pressure
339      saveVTK("phi2D.%2.4i.vtu"%t_step,interface=func,velocity=velocity,pressure=pressure)
340      #courant condition
341      dt = 0.4*Lsup(mesh.getSize())/Lsup(velocity)
342      t_step += 1
343    
344  \end{python}  \end{python}

Legend:
Removed from v.2115  
changed lines
  Added in v.2120

  ViewVC Help
Powered by ViewVC 1.1.26