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

Revision 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 4184 byte(s)
```Change __url__ to launchpad site

```
 1 2 ######################################################## 3 # 4 # Copyright (c) 2003-2008 by University of Queensland 5 # Earth Systems Science Computational Center (ESSCC) 6 7 # 8 # Primary Business: Queensland, Australia 9 # Licensed under the Open Software License version 3.0 10 11 # 12 ######################################################## 13 14 __copyright__="""Copyright (c) 2003-2008 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 20 __url__= 21 22 from esys.escript import * 23 from esys.escript.linearPDEs import Poisson 24 from esys import finley 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=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1) 35 totne=ne1*ne 36 else: 37 ne2=int(ne*height+0.5) 38 mydomain=finley.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()

## Properties

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

 ViewVC Help Powered by ViewVC 1.1.26