Diff of /trunk/doc/examples/usersguide/fluid.py

revision 4853 by jfenwick, Wed Apr 9 05:41:57 2014 UTC revision 5288 by sshaw, Tue Dec 2 23:18:40 2014 UTC
# Line 20  from __future__ import print_function Line 20  from __future__ import print_function
20  ## Newtonian fluid using StokesProblemCartesian class##  ## Newtonian fluid using StokesProblemCartesian class##
21
22  from esys.escript import *  from esys.escript import *
23  import esys.finley  try:
24        from esys.finley import Rectangle
25        HAVE_FINLEY = True
26    except ImportError:
27        HAVE_FINLEY = False
28  from esys.escript.linearPDEs import LinearPDE  from esys.escript.linearPDEs import LinearPDE
29  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian
30  from esys.weipa import saveVTK  from esys.weipa import saveVTK
31    if not HAVE_FINLEY:
32  #physical constants      print("Finley module not available")
33  eta=1.0  else:
34  rho=100.0      #physical constants
35  g=10.0      eta=1.0
36        rho=100.0
37  #solver settings      g=10.0
38  tolerance=1.0e-4
39  max_iter=200      #solver settings
40  t_end=50      tolerance=1.0e-4
41  t=0.0      max_iter=200
42  time=0      t_end=50
43  verbose='TRUE'      t=0.0
44  useUzawa='TRUE'      time=0
45        verbose='TRUE'
46  #define mesh      useUzawa='TRUE'
47  H=2.0
48  L=1.0      #define mesh
49  W=1.0      H=2.0
50  mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20, useElementsOnFace=0) # use linear macro elements for pressure      L=1.0
51  coordinates = mesh.getX()      W=1.0
52        mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20, useElementsOnFace=0) # use linear macro elements for pressure
53  #gravitational force      coordinates = mesh.getX()
54  Y=Vector(0.0, Function(mesh))
55  Y[1]=-rho*g      #gravitational force
56        Y=Vector(0.0, Function(mesh))
57  #element spacing      Y[1]=-rho*g
58  h=Lsup(mesh.getSize())
59        #element spacing
60  #boundary conditions for slip at base      h=Lsup(mesh.getSize())
61  boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
62        #boundary conditions for slip at base
63  #velocity and pressure vectors      boundary_cond=whereZero(coordinates[1])*[0.0,1.0]+whereZero(coordinates[0])*[1.0,0.0]
64  velocity=Vector(0.0, Solution(mesh))
65  pressure=Scalar(0.0, ReducedSolution(mesh))      #velocity and pressure vectors
66        velocity=Vector(0.0, Solution(mesh))
67  #Stokes Cartesian      pressure=Scalar(0.0, ReducedSolution(mesh))
68  solution=StokesProblemCartesian(mesh)
69  solution.setTolerance(tolerance)      #Stokes Cartesian
70        solution=StokesProblemCartesian(mesh)
71  while t <= t_end:      solution.setTolerance(tolerance)
72
73    print(" ----- Time step = %s -----"%( t ))      while t <= t_end:
74    print("Time = %s seconds"%( time ))
75          print(" ----- Time step = %s -----"%( t ))
76    solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y)        print("Time = %s seconds"%( time ))
77    velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
79    print("Max velocity =", Lsup(velocity), "m/s")        velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True)
80
81    #Courant condition        print("Max velocity =", Lsup(velocity), "m/s")
82    dt=0.4*h/(Lsup(velocity))
83    print("dt", dt)        #Courant condition
84            dt=0.4*h/(Lsup(velocity))
85    #displace the mesh        print("dt", dt)
86    displacement = velocity * dt
87    coordinates = mesh.getX()        #displace the mesh
88    newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))        displacement = velocity * dt
89    mesh.setX(newx)          coordinates = mesh.getX()
90            newx=interpolate(coordinates + displacement, ContinuousFunction(mesh))
91    time += dt        mesh.setX(newx)
92
93    vel_mag = length(velocity)        time += dt
94
95    #save velocity and pressure output        vel_mag = length(velocity)
96    saveVTK("vel.%2.2i.vtu"%(t),vel=vel_mag,vec=velocity,pressure=pressure)
97    t = t+1.0        #save velocity and pressure output
98          saveVTK("vel.%2.2i.vtu"%(t),vel=vel_mag,vec=velocity,pressure=pressure)
99          t = t+1.0

Legend:
 Removed from v.4853 changed lines Added in v.5288