/[escript]/trunk/finley/test/python/PoissonSolverTest.py
ViewVC logotype

Annotation of /trunk/finley/test/python/PoissonSolverTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years ago) by elspeth
File MIME type: text/x-python
File size: 3721 byte(s)
More copyright.

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=ne2*ne*ne
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()

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26