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

revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC revision 2034 by artak, Thu Nov 13 03:03:42 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  from esys.finley import Rectangle,Brick
25    import time
26    from optparse import OptionParser
27
28    import sys
29
30    #sys.stderr=None
31    parser = OptionParser(usage="%prog [-e [NE]]")
32    parser.add_option("-e", "--elements", dest="NE", help="number of elements in one direction.",metavar="NE", default=30)
33    (options, args) = parser.parse_args()
34    NE=int(options.NE)
35
36  # generate domain:  # generate domain:
37  mydomain = Rectangle(l0=1.,l1=1.,n0=40, n1=20)  mydomain = Rectangle(l0=1.,l1=1.,n0=NE, n1=NE)
38    #mydomain = Brick(l0=1.,l1=1., l2=1.,n0=10, n1=10, n2=10)
39  # define characteristic function of Gamma^D  # define characteristic function of Gamma^D
40  x = mydomain.getX()  x = mydomain.getX()
41  gammaD = whereZero(x[0])+whereZero(x[1])  gammaD = whereZero(x[0])+whereZero(x[1])
42  # define PDE and get its solution u  # define PDE and get its solution u
43  mypde = Poisson(domain=mydomain)  mypde = Poisson(domain=mydomain)
45  u = mypde.getSolution()
46    #mypde.setSolverMethod(mypde.PCG,mypde.JACOBI)
47    #t1 = time.time()
48    #u=mypde.getSolution(verbose=True)
49    #elapsed = time.time() - t1
50    #print "ESCRIPT Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
51    #print u
52
53    mypde.setSolverMethod(mypde.PCG,mypde.AMG)
54    t1 = time.time()
55    u=mypde.getSolution(verbose=True)
56    elapsed = time.time() - t1
57    print "AMG Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
58    print u
59
60    mypde.setSolverMethod(mypde.PCG,mypde.JACOBI)
61    t1 = time.time()
62    u=mypde.getSolution(verbose=True)
63    elapsed = time.time() - t1
64    print "ESCRIPT Job completed","\t",NE,"\tin\t",elapsed,"\t(seconds)\t",elapsed/3600.0, " (hours)"
65    print u
66    #u = mypde.getSolution()
67  # write u to an external file  # write u to an external file
68  saveVTK("u.xml",sol=u)  #saveVTK("u.xml",sol=u)

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