/[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 4657 - (show annotations)
Thu Feb 6 06:12:20 2014 UTC (5 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 4329 byte(s)
I changed some files.
Updated copyright notices, added GeoComp.



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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26