test/python/PoissonSolverTest.py

# $Id$


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