/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years, 1 month ago) by jfenwick
File MIME type: text/x-python
File size: 4366 byte(s)
Make everyone sad by touching all the files

Copyright dates update

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26