/[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 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years ago) by ksteube
File MIME type: text/x-python
File size: 4189 byte(s)
Copyright updated in all python files

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26