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

Revision 104 - (show annotations)
Fri Dec 17 07:43:12 2004 UTC (16 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 2964 byte(s)
```*** empty log message ***

```
 1 # \$Id\$ 2 3 from esys.escript import * 4 from esys.linearPDEs import Poisson 5 import esys.finley as finley 6 7 ne_list=[10,15,22,33,50,75] 8 height_list=[0.25,0.5,1.] 9 10 11 def getDomain(dim,ne,height): 12 13 if dim==2: 14 ne1=int(ne*height+0.5) 15 mydomain=finley.Rectangle(n0=ne,n1=ne1,l1=height,order=2) 16 totne=ne1*ne 17 else: 18 ne2=int(ne*height+0.5) 19 mydomain=finley.Brick(n0=ne,n1=ne,n2=ne2,l2=height,order=2) 20 totne=ne2*ne*ne 21 print "%d -dimensional domain generated."%dim 22 print "height of the domain is ",height 23 print "total number of elements is ",totne 24 return mydomain 25 26 27 def Solve1(mydomain,height): 28 print "Fully constraint solution" 29 l=[1.,1.,1.] 30 l[mydomain.getDim()-1]=height 31 cf=ContinuousFunction(mydomain) 32 x=cf.getX() 33 #construct exact solution: 34 u_ex=Scalar(1.,cf) 35 for i in range(mydomain.getDim()): 36 u_ex*=x[i]*(x[i]-l[i]) 37 #construct mask: 38 msk=Scalar(0.,cf) 39 for i in range(mydomain.getDim()): 40 msk+=x[i].whereZero()+(x[i]-l[i]).whereZero() 41 #construct right hand side 42 f=Scalar(0,cf) 43 for i in range(mydomain.getDim()): 44 f_p=Scalar(1,cf) 45 for j in range(mydomain.getDim()): 46 if i==j: 47 f_p*=-2. 48 else: 49 f_p*=x[j]*(x[j]-l[j]) 50 f+=f_p 51 52 mypde=Poisson(f=f,q=msk) 53 u=mypde.getSolution() 54 error=Lsup(u-u_ex)/Lsup(u_ex) 55 print "error = ",error 56 return error 57 58 def Solve2(mydomain,height): 59 print "Partially constraint solution" 60 l=[1.,1.,1.] 61 l[mydomain.getDim()-1]=height 62 print l 63 cf=ContinuousFunction(mydomain) 64 x=cf.getX() 65 #construct exact solution: 66 u_ex=Scalar(1.,cf) 67 for i in range(mydomain.getDim()): 68 u_ex*=x[i]*(2*l[i]-x[i]) 69 #construct mask: 70 msk=Scalar(0.,cf) 71 for i in range(mydomain.getDim()): 72 msk+=x[i].whereZero() 73 #construct right hand side 74 f=Scalar(0,cf) 75 for i in range(mydomain.getDim()): 76 f_p=Scalar(1,cf) 77 for j in range(mydomain.getDim()): 78 if i==j: 79 f_p*=2. 80 else: 81 f_p*=x[j]*(2*l[j]-x[j]) 82 f+=f_p 83 mypde=Poisson(f=f,q=msk) 84 u=mypde.getSolution() 85 error=Lsup(u-u_ex)/Lsup(u_ex) 86 print "error = ",error 87 return error 88 89 90 error=0 91 for ne in ne_list: 92 for dim in [2,3]: 93 for height in height_list: 94 print "***************************************************************" 95 mydomain= getDomain(dim,ne,height) 96 print "---------------------------------------------------------------" 97 error=max(error,Solve1(mydomain,height)) 98 print "---------------------------------------------------------------" 99 error=max(error,Solve2(mydomain,height)) 100 print "***************************************************************" 101 102 print "***************************************************************" 103 print "maximum error: ",error 104 print "***************************************************************"

Properties

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