test/python/PoissonSolverTest.py

# $Id$


from esys.escript import *
from esys.escript.linearPDEs import Poisson
from esys import finley

ne_list=[10,15,22,33,50,75]
height_list=[0.25,0.5,1.]


def getDomain(dim,ne,height):

    if dim==2:
     ne1=int(ne*height+0.5)
     mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
     totne=ne1*ne
    else:
     ne2=int(ne*height+0.5)
     mydomain=finley.Brick(n0=ne,n1=ne,n2=ne2,l2=height,order=2)
     totne=ne2*ne*ne
    print "%d -dimensional domain generated."%dim
    print "height of the domain is ",height
    print "total number of elements is ",totne
    return mydomain


def Solve1(mydomain,height):
    print "Fully constraint solution"
    l=[1.,1.,1.]
    l[mydomain.getDim()-1]=height
    cf=ContinuousFunction(mydomain)
    x=cf.getX()
    #construct exact solution:
    u_ex=Scalar(1.,cf)
    for i in range(mydomain.getDim()):
      u_ex*=x[i]*(x[i]-l[i])
    #construct mask:
    msk=Scalar(0.,cf)
    for i in range(mydomain.getDim()):
      msk+=x[i].whereZero()+(x[i]-l[i]).whereZero()
    #construct right hand side 
    f=Scalar(0,cf)
    for i in range(mydomain.getDim()):
       f_p=Scalar(1,cf)
       for j in range(mydomain.getDim()):
          if i==j:
             f_p*=-2.
          else:
             f_p*=x[j]*(x[j]-l[j])
       f+=f_p

    mypde=Poisson(mydomain)
    mypde.setTolerance(1.e-10)
    mypde.setValue(f=f,q=msk)
    u=mypde.getSolution()
    error=Lsup(u-u_ex)/Lsup(u_ex)
    print "error = ",error
    return error

def Solve2(mydomain,height):
    print "Partially constraint solution"
    l=[1.,1.,1.]
    l[mydomain.getDim()-1]=height
    print l
    cf=ContinuousFunction(mydomain)
    x=cf.getX()
    #construct exact solution:
    u_ex=Scalar(1.,cf)
    for i in range(mydomain.getDim()):
      u_ex*=x[i]*(2*l[i]-x[i])
    #construct mask:
    msk=Scalar(0.,cf)
    for i in range(mydomain.getDim()):
      msk+=x[i].whereZero()
    #construct right hand side 
    f=Scalar(0,cf)
    for i in range(mydomain.getDim()):
       f_p=Scalar(1,cf)
       for j in range(mydomain.getDim()):
          if i==j:
             f_p*=2.
          else:
             f_p*=x[j]*(2*l[j]-x[j])
       f+=f_p
    mypde=Poisson(mydomain)
    mypde.setTolerance(1.e-10)
    mypde.setValue(f=f,q=msk)
    u=mypde.getSolution()
    error=Lsup(u-u_ex)/Lsup(u_ex)
    print "error = ",error
    return error


def main() :
    error=0
    for ne in ne_list:
       for dim in [2,3]:
       # for dim in [2]:
          for height in height_list:
             print "***************************************************************"
             mydomain= getDomain(dim,ne,height)
             print "---------------------------------------------------------------"
             error=max(error,Solve1(mydomain,height))
             print "---------------------------------------------------------------"
             error=max(error,Solve2(mydomain,height))
             print "***************************************************************"

    print "***************************************************************"
    print "maximum error: ",error
    print "***************************************************************"


import profile as Pr, pstats as Ps


if __name__ == "__main__":
    pr = Pr.Profile()
    pr.calibrate(10000)
    Pr.run('main()','eos_stats')
    stats = Ps.Stats('eos_stats')
    stats.strip_dirs()
    stats.sort_stats('time')
    stats.print_stats()
1	jgs	102	# $Id$
2
3
4	jgs	149	from esys.escript import *
5			from esys.escript.linearPDEs import Poisson
6			from esys import finley
7	jgs	147
8	jgs	102	ne_list=[10,15,22,33,50,75]
9			height_list=[0.25,0.5,1.]
10
11
12			def getDomain(dim,ne,height):
13
14			if dim==2:
15			ne1=int(ne*height+0.5)
16	phornby	152	mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
17	jgs	102	totne=ne1*ne
18			else:
19			ne2=int(ne*height+0.5)
20			mydomain=finley.Brick(n0=ne,n1=ne,n2=ne2,l2=height,order=2)
21			totne=ne2nene
22			print "%d -dimensional domain generated."%dim
23			print "height of the domain is ",height
24			print "total number of elements is ",totne
25			return mydomain
26
27
28			def Solve1(mydomain,height):
29			print "Fully constraint solution"
30			l=[1.,1.,1.]
31			l[mydomain.getDim()-1]=height
32			cf=ContinuousFunction(mydomain)
33			x=cf.getX()
34			#construct exact solution:
35			u_ex=Scalar(1.,cf)
36			for i in range(mydomain.getDim()):
37			u_ex=x[i](x[i]-l[i])
38			#construct mask:
39			msk=Scalar(0.,cf)
40			for i in range(mydomain.getDim()):
41			msk+=x[i].whereZero()+(x[i]-l[i]).whereZero()
42			#construct right hand side
43			f=Scalar(0,cf)
44			for i in range(mydomain.getDim()):
45			f_p=Scalar(1,cf)
46			for j in range(mydomain.getDim()):
47			if i==j:
48			f_p*=-2.
49			else:
50			f_p=x[j](x[j]-l[j])
51			f+=f_p
52
53	jgs	108	mypde=Poisson(mydomain)
54	jgs	147	mypde.setTolerance(1.e-10)
55	jgs	108	mypde.setValue(f=f,q=msk)
56	jgs	102	u=mypde.getSolution()
57			error=Lsup(u-u_ex)/Lsup(u_ex)
58			print "error = ",error
59			return error
60
61			def Solve2(mydomain,height):
62			print "Partially constraint solution"
63			l=[1.,1.,1.]
64			l[mydomain.getDim()-1]=height
65			print l
66			cf=ContinuousFunction(mydomain)
67			x=cf.getX()
68			#construct exact solution:
69			u_ex=Scalar(1.,cf)
70			for i in range(mydomain.getDim()):
71			u_ex=x[i](2*l[i]-x[i])
72			#construct mask:
73			msk=Scalar(0.,cf)
74			for i in range(mydomain.getDim()):
75			msk+=x[i].whereZero()
76			#construct right hand side
77			f=Scalar(0,cf)
78			for i in range(mydomain.getDim()):
79			f_p=Scalar(1,cf)
80			for j in range(mydomain.getDim()):
81			if i==j:
82			f_p*=2.
83			else:
84			f_p=x[j](2*l[j]-x[j])
85			f+=f_p
86	jgs	108	mypde=Poisson(mydomain)
87	jgs	147	mypde.setTolerance(1.e-10)
88	jgs	108	mypde.setValue(f=f,q=msk)
89	jgs	102	u=mypde.getSolution()
90			error=Lsup(u-u_ex)/Lsup(u_ex)
91			print "error = ",error
92			return error
93
94
95	phornby	152	def main() :
96			error=0
97			for ne in ne_list:
98			for dim in [2,3]:
99			# for dim in [2]:
100			for height in height_list:
101			print "***************************************************************"
102			mydomain= getDomain(dim,ne,height)
103			print "---------------------------------------------------------------"
104			error=max(error,Solve1(mydomain,height))
105			print "---------------------------------------------------------------"
106			error=max(error,Solve2(mydomain,height))
107			print "***************************************************************"
108	jgs	102
109	phornby	152	print "***************************************************************"
110			print "maximum error: ",error
111			print "***************************************************************"
112
113
114
115			import profile as Pr, pstats as Ps
116
117
118			if __name__ == "__main__":
119			pr = Pr.Profile()
120			pr.calibrate(10000)
121			Pr.run('main()','eos_stats')
122			stats = Ps.Stats('eos_stats')
123			stats.strip_dirs()
124			stats.sort_stats('time')
125			stats.print_stats()
Name	Value
svn:eol-style	native
svn:keywords	Author Date Id Revision