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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 152 - (hide annotations)
Fri Oct 21 08:32:21 2005 UTC (15 years ago) by phornby
File MIME type: text/x-python
File size: 3438 byte(s)
DataConstantTestCase.cpp, DataExpandedTestCase.cpp:

Replace

DataExpanded testData1(FunctionSpace(),....)

with

FunctionSpace tmp_fns;
DataExpanded testData1(tmp_fns,shape,data);

for GCC's sake.




Added some profiling to PoissonSolverTest.py.

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