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() |
29 |
gammaD = whereZero(x[0])+whereZero(x[1]) |
gammaD = whereZero(x[0])+whereZero(x[1]) |
30 |
# define PDE and get its solution u |
# define PDE and get its solution u |
31 |
mypde = Poisson(domain=mydomain) |
mypde = Poisson(domain=mydomain) |
32 |
mypde.setValue(f=1,q=gammaD) |
mypde.setValue(f=1,q=gammaD) |
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) |