← Back to branch summary

~vcs-imports/escript-finley/trunk

« back to all changes in this revision

Viewing changes to dudley/test/python/PoissonSolverTest.py

  • Committer: jfenwick
  • Date: 2010-10-11 01:48:14 UTC
  • Revision ID: svn-v4:77569008-7704-0410-b7a0-a92fef0b09fd:trunk:3259
Merging dudley and scons updates from branches

Show diffs side-by-side

added added

removed removed

Lines of Context:
dudley/test/python/PoissonSolverTest.py
 
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()