# Diff of /trunk/doc/examples/usersguide/poisson.py

revision 2034 by artak, Thu Nov 13 03:03:42 2008 UTC revision 2035 by jfenwick, Thu Nov 13 03:16:09 2008 UTC
# Line 21  __url__="http://www.uq.edu.au/esscc/escr Line 21  __url__="http://www.uq.edu.au/esscc/escr
21
22  from esys.escript import *  from esys.escript import *
23  from esys.escript.linearPDEs import Poisson  from esys.escript.linearPDEs import Poisson
24  from esys.finley import Rectangle,Brick  from esys.finley import Rectangle
import time
from optparse import OptionParser

import sys

#sys.stderr=None
parser = OptionParser(usage="%prog [-e [NE]]")
parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=30)
(options, args) = parser.parse_args()
NE=int(options.NE)

25  # generate domain:  # generate domain:
26  mydomain = Rectangle(l0=1.,l1=1.,n0=NE, n1=NE)  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)
#mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
27  # define characteristic function of Gamma^D  # define characteristic function of Gamma^D
28  x = mydomain.getX()  x = mydomain.getX()
30  # define PDE and get its solution u  # define PDE and get its solution u
31  mypde = Poisson(domain=mydomain)  mypde = Poisson(domain=mydomain)
33    u = mypde.getSolution()
#mypde.setSolverMethod(mypde.PCG,mypde.JACOBI)
#t1 = time.time()
#u=mypde.getSolution(verbose=True)
#elapsed = time.time() - t1
#print "ESCRIPT Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
#print u

mypde.setSolverMethod(mypde.PCG,mypde.AMG)
t1 = time.time()
u=mypde.getSolution(verbose=True)
elapsed = time.time() - t1
print "AMG Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
print u

mypde.setSolverMethod(mypde.PCG,mypde.JACOBI)
t1 = time.time()
u=mypde.getSolution(verbose=True)
elapsed = time.time() - t1
print "ESCRIPT Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
print u
#u = mypde.getSolution()
34  # write u to an external file  # write u to an external file
35  #saveVTK("u.xml",sol=u)  saveVTK("u.xml",sol=u)

Legend:
 Removed from v.2034 changed lines Added in v.2035