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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 # $Id$
2
3
4 __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 from esys.escript import *
10 from esys.escript.linearPDEs import Poisson
11 from esys import finley
12
13 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 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
22 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 msk+=whereZero(x[i])+whereZero(x[i]-l[i])
47 #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 mypde=Poisson(mydomain)
59 mypde.setTolerance(1.e-10)
60 mypde.setValue(f=f,q=msk)
61 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 msk+=whereZero(x[i])
81 #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 mypde=Poisson(mydomain)
92 mypde.setTolerance(1.e-10)
93 mypde.setValue(f=f,q=msk)
94 u=mypde.getSolution()
95 error=Lsup(u-u_ex)/Lsup(u_ex)
96 print "error = ",error
97 return error
98
99
100 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
114 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