--- trunk/doc/examples/usersguide/fluid.py 2009/12/01 06:10:10 2793 +++ trunk/doc/examples/usersguide/fluid.py 2014/04/09 00:26:32 4848 @@ -1,14 +1,18 @@ -######################################################## +from __future__ import print_function +############################################################################## # -# Copyright (c) 2008-2009 by University of Queensland -# Earth Systems Science Computational Center (ESSCC) -# http://www.uq.edu.au/esscc +# Copyright (c) 2008-2014 by University of Queensland +# http://www.uq.edu.au # # Primary Business: Queensland, Australia # Licensed under the Open Software License version 3.0 # http://www.opensource.org/licenses/osl-3.0.php # -######################################################## +# Development until 2012 by Earth Systems Science Computational Center (ESSCC) +# Development 2012-2013 by School of Earth Sciences +# Development from 2014 by Centre for Geoscience Computing (GeoComp) +# +############################################################################## ######## August 2008 ######## ########## Leon Graham ########## @@ -18,6 +22,7 @@ import esys.finley from esys.escript.linearPDEs import LinearPDE from esys.escript.models import StokesProblemCartesian +from esys.weipa import saveVTK #physical constants eta=1.0 @@ -37,7 +42,7 @@ H=2.0 L=1.0 W=1.0 -mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20) # use linear macro elements for pressure +mesh = esys.finley.Rectangle(l0=L, l1=H, order=-1, n0=20, n1=20, useElementsOnFace=0) # use linear macro elements for pressure coordinates = mesh.getX() #gravitational force @@ -60,22 +65,23 @@ while t <= t_end: - print " ----- Time step = %s -----"%( t ) - print "Time = %s seconds"%( time ) + print(" ----- Time step = %s -----"%( t )) + print("Time = %s seconds"%( time )) solution.initialize(fixed_u_mask=boundary_cond,eta=eta,f=Y) velocity,pressure=solution.solve(velocity,pressure,max_iter=max_iter,verbose=verbose,usePCG=True) - print "Max velocity =", Lsup(velocity), "m/s" + print("Max velocity =", Lsup(velocity), "m/s") #Courant condition dt=0.4*h/(Lsup(velocity)) - print "dt", dt + print("dt", dt) #displace the mesh displacement = velocity * dt coordinates = mesh.getX() - mesh.setX(coordinates + displacement) + newx=interpolate(coordinates + displacement, ContinuousFunction(mesh)) + mesh.setX(newx) time += dt