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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (3 years 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 ksteube 1809
2 jfenwick 3981 ##############################################################################
3 ksteube 1312 #
4 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 ksteube 1312 #
7 ksteube 1809 # Primary Business: Queensland, Australia
8 jfenwick 6112 # Licensed under the Apache License, version 2.0
9     # http://www.apache.org/licenses/LICENSE-2.0
10 ksteube 1312 #
11 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 jfenwick 4657 # Development 2012-2013 by School of Earth Sciences
13     # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 jfenwick 3981 #
15     ##############################################################################
16 jgs 102
17 sshaw 5706 from __future__ import print_function, division
18    
19 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 ksteube 1809 Primary Business: Queensland, Australia"""
22 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
23     http://www.apache.org/licenses/LICENSE-2.0"""
24 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
25 ksteube 1809
26 jgs 149 from esys.escript import *
27     from esys.escript.linearPDEs import Poisson
28     from esys import finley
29 jgs 147
30 jgs 102 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 phornby 152 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1)
39 jgs 102 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 jfenwick 3772 print("%d -dimensional domain generated."%dim)
45     print("height of the domain is ",height)
46     print("total number of elements is ",totne)
47 jgs 102 return mydomain
48    
49    
50     def Solve1(mydomain,height):
51 jfenwick 3772 print("Fully constraint solution")
52 jgs 102 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 gross 302 msk+=whereZero(x[i])+whereZero(x[i]-l[i])
64 jgs 102 #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 jgs 108 mypde=Poisson(mydomain)
76 jgs 147 mypde.setTolerance(1.e-10)
77 jgs 108 mypde.setValue(f=f,q=msk)
78 jgs 102 u=mypde.getSolution()
79     error=Lsup(u-u_ex)/Lsup(u_ex)
80 jfenwick 3772 print("error = ",error)
81 jgs 102 return error
82    
83     def Solve2(mydomain,height):
84 jfenwick 3772 print("Partially constraint solution")
85 jgs 102 l=[1.,1.,1.]
86     l[mydomain.getDim()-1]=height
87 jfenwick 3772 print(l)
88 jgs 102 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 gross 302 msk+=whereZero(x[i])
98 jgs 102 #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 jgs 108 mypde=Poisson(mydomain)
109 jgs 147 mypde.setTolerance(1.e-10)
110 jgs 108 mypde.setValue(f=f,q=msk)
111 jgs 102 u=mypde.getSolution()
112     error=Lsup(u-u_ex)/Lsup(u_ex)
113 jfenwick 3772 print("error = ",error)
114 jgs 102 return error
115    
116    
117 phornby 152 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 jfenwick 3772 print("***************************************************************")
124 phornby 152 mydomain= getDomain(dim,ne,height)
125 jfenwick 3772 print("---------------------------------------------------------------")
126 phornby 152 error=max(error,Solve1(mydomain,height))
127 jfenwick 3772 print("---------------------------------------------------------------")
128 phornby 152 error=max(error,Solve2(mydomain,height))
129 jfenwick 3772 print("***************************************************************")
130 jgs 102
131 jfenwick 3772 print("***************************************************************")
132     print("maximum error: ",error)
133     print("***************************************************************")
134 phornby 152
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