/[escript]/trunk/esys2/finley/test/python/PoissonSolverTest.py
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 102 - (show annotations)
Wed Dec 15 07:08:39 2004 UTC (15 years, 8 months ago) by jgs
File MIME type: text/x-python
File size: 3152 byte(s)
*** empty log message ***

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26