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

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
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