/[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 155 - (show annotations)
Wed Nov 9 02:02:19 2005 UTC (13 years, 11 months ago) by jgs
File MIME type: text/x-python
File size: 3438 byte(s)
move all directories from trunk/esys2 into trunk and remove esys2

1 # $Id$
2
3
4 from esys.escript import *
5 from esys.escript.linearPDEs import Poisson
6 from esys import finley
7
8 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 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
17 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=ne2*ne*ne
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 mypde=Poisson(mydomain)
54 mypde.setTolerance(1.e-10)
55 mypde.setValue(f=f,q=msk)
56 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 mypde=Poisson(mydomain)
87 mypde.setTolerance(1.e-10)
88 mypde.setValue(f=f,q=msk)
89 u=mypde.getSolution()
90 error=Lsup(u-u_ex)/Lsup(u_ex)
91 print "error = ",error
92 return error
93
94
95 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
109 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()

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26