/[escript]/branches/symbolic_from_3470/dudley/test/python/PoissonSolverTest.py
ViewVC logotype

Contents of /branches/symbolic_from_3470/dudley/test/python/PoissonSolverTest.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3789 - (show annotations)
Tue Jan 31 04:55:05 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 4199 byte(s)
Fast forward to latest trunk revision 3788.

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # 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 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="https://launchpad.net/escript-finley"
21
22 from esys.escript import *
23 from esys.escript.linearPDEs import Poisson
24 from esys import dudley
25
26 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 mydomain=dudley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
35 totne=ne1*ne
36 else:
37 ne2=int(ne*height+0.5)
38 mydomain=dudley.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 msk+=whereZero(x[i])+whereZero(x[i]-l[i])
60 #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 mypde=Poisson(mydomain)
72 mypde.setTolerance(1.e-10)
73 mypde.setValue(f=f,q=msk)
74 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 msk+=whereZero(x[i])
94 #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 mypde=Poisson(mydomain)
105 mypde.setTolerance(1.e-10)
106 mypde.setValue(f=f,q=msk)
107 u=mypde.getSolution()
108 error=Lsup(u-u_ex)/Lsup(u_ex)
109 print("error = ",error)
110 return error
111
112
113 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
127 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