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

Revision 617 - (hide annotations)
Wed Mar 22 02:58:17 2006 UTC (14 years ago) by elspeth
File MIME type: text/x-python
File size: 3721 byte(s)
```More copyright.

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

## Properties

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