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

Revision 302 - (show annotations)
Fri Dec 2 05:51:51 2005 UTC (16 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 3433 byte(s)
```runs again
```
 1 # \$Id\$ 2 3 4 from esys.escript import * 5 from esys.escript.linearPDEs import Poisson 6 from esys import finley 7 8 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 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=1) 17 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+=whereZero(x[i])+whereZero(x[i]-l[i]) 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 mypde=Poisson(mydomain) 54 mypde.setTolerance(1.e-10) 55 mypde.setValue(f=f,q=msk) 56 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+=whereZero(x[i]) 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 mypde=Poisson(mydomain) 87 mypde.setTolerance(1.e-10) 88 mypde.setValue(f=f,q=msk) 89 u=mypde.getSolution() 90 error=Lsup(u-u_ex)/Lsup(u_ex) 91 print "error = ",error 92 return error 93 94 95 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 109 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