/[escript]/trunk-mpi-branch/finley/test/python/run_simplesolve.py
ViewVC logotype

Contents of /trunk-mpi-branch/finley/test/python/run_simplesolve.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1087 - (show annotations)
Thu Apr 12 10:01:47 2007 UTC (12 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 15253 byte(s)
the MPI version of PASO.PCG is running. 
There is a bug in the rectangular mesh generators but they need to be
revised in any case to clean up the code.


1 # $Id:$
2
3 """
4 Test suite for the linearPDE and pdetools test on finley
5
6 @remark:
7
8 @var __author__: name of author
9 @var __licence__: licence agreement
10 @var __url__: url entry point on documentation
11 @var __version__: version
12 @var __date__: date of the version
13 """
14
15 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
16 http://www.access.edu.au
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __author__="Lutz Gross, l.gross@uq.edu.au"
21 __url__="http://www.iservo.edu.au/esys/escript"
22 __version__="$Revision: 859 $"
23 __date__="$Date: 2006-09-26 12:19:18 +1000 (Tue, 26 Sep 2006) $"
24
25
26 import unittest
27
28 from esys.escript import *
29 from esys.finley import Rectangle,Brick
30 from esys.escript.linearPDEs import LinearPDE
31
32 try:
33 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
34 except KeyError:
35 FINLEY_TEST_DATA='.'
36
37 FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/"
38
39 # number of elements in the spatial directions
40 NE0=15
41 NE1=15
42 NE2=8
43
44 SOLVER_TOL=1.e-8
45 REL_TOL=1.e-6
46
47 FAC_DIAG=1.
48 FAC_OFFDIAG=-0.4
49
50 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
51 def test_solve(self):
52 domain=Rectangle(NE0,NE1,1)
53 x=Solution(domain).getX()
54 # --- set exact solution ----
55 u_ex=Scalar(0,Solution(domain))
56 u_ex=1.+2.*x[0]+3.*x[1]
57 # --- set exact gradient -----------
58 g_ex=Data(0.,(2,),Solution(domain))
59 g_ex[0]=2.
60 g_ex[1]=3.
61 # -------- test gradient --------------------------------
62 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
63 # -------- set-up PDE -----------------------------------
64 pde=LinearPDE(domain,numEquations=1)
65 mask=whereZero(x[0])
66 pde.setValue(r=u_ex,q=mask)
67 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
68 # -------- get the solution ---------------------------
69 pde.setTolerance(SOLVER_TOL)
70 pde.setSolverMethod(pde.PCG,pde.JACOBI)
71 pde.setSolverPackage(pde.PASO)
72 u=pde.getSolution()
73 # -------- test the solution ---------------------------
74 error=Lsup(u-u_ex)/Lsup(u_ex)
75 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
76 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
77 def test_solve(self):
78 domain=Rectangle(NE0,NE1,1)
79 x=Solution(domain).getX()
80 # --- set exact solution ----
81 u_ex=Vector(0,Solution(domain))
82 u_ex[0]=1.+2.*x[0]+3.*x[1]
83 u_ex[1]=-1.+3.*x[0]+2.*x[1]
84 # --- set exact gradient -----------
85 g_ex=Data(0.,(2,2),Solution(domain))
86 g_ex[0,0]=2.
87 g_ex[0,1]=3.
88 g_ex[1,0]=3.
89 g_ex[1,1]=2.
90 # -------- test gradient --------------------------------
91 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
92 # -------- set-up PDE -----------------------------------
93 pde=LinearPDE(domain,numEquations=2)
94 mask=whereZero(x[0])
95 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
96 A=Tensor4(0,Function(domain))
97 A[0,:,0,:]=kronecker(2)
98 A[1,:,1,:]=kronecker(2)
99 Y=Vector(0.,Function(domain))
100 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
101 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
102 pde.setValue(A=A,
103 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
104 Y=Y,
105 y=matrixmult(g_ex,domain.getNormal()))
106 # -------- get the solution ---------------------------
107 pde.setTolerance(SOLVER_TOL)
108 pde.setSolverMethod(pde.PCG,pde.JACOBI)
109 pde.setSolverPackage(pde.PASO)
110 u=pde.getSolution()
111 # -------- test the solution ---------------------------
112 error=Lsup(u-u_ex)/Lsup(u_ex)
113 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
114 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
115 def test_solve(self):
116 domain=Rectangle(NE0,NE1,2)
117 x=Solution(domain).getX()
118 # --- set exact solution ----
119 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
120 # --- set exact gradient -----------
121 g_ex=Data(0.,(2,),Solution(domain))
122 g_ex[0]=2.+8.*x[0]+5.*x[1]
123 g_ex[1]=3.+5.*x[0]+12.*x[1]
124 # -------- test gradient --------------------------------
125 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
126 # -------- set-up PDE -----------------------------------
127 pde=LinearPDE(domain,numEquations=1)
128 mask=whereZero(x[0])
129 pde.setValue(r=u_ex,q=mask)
130 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
131 # -------- get the solution ---------------------------
132 pde.setTolerance(SOLVER_TOL)
133 pde.setSolverMethod(pde.PCG,pde.JACOBI)
134 pde.setSolverPackage(pde.PASO)
135 u=pde.getSolution()
136 # -------- test the solution ---------------------------
137 error=Lsup(u-u_ex)/Lsup(u_ex)
138 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
139 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
140 def test_solve(self):
141 domain=Rectangle(NE0,NE1,2)
142 x=Solution(domain).getX()
143 # --- set exact solution ----
144 u_ex=Vector(0,Solution(domain))
145 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
146 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
147 # --- set exact gradient -----------
148 g_ex=Data(0.,(2,2),Solution(domain))
149 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
150 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
151 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
152 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
153 # -------- test gradient --------------------------------
154 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
155 # -------- set-up PDE -----------------------------------
156 pde=LinearPDE(domain,numEquations=2)
157 mask=whereZero(x[0])
158 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
159 A=Tensor4(0,Function(domain))
160 A[0,:,0,:]=kronecker(2)
161 A[1,:,1,:]=kronecker(2)
162 Y=Vector(0.,Function(domain))
163 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
164 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
165 pde.setValue(A=A,
166 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
167 Y=Y-[20.,10.],
168 y=matrixmult(g_ex,domain.getNormal()))
169 # -------- get the solution ---------------------------
170 pde.setTolerance(SOLVER_TOL)
171 pde.setSolverMethod(pde.PCG,pde.JACOBI)
172 pde.setSolverPackage(pde.PASO)
173 u=pde.getSolution()
174 # -------- test the solution ---------------------------
175 error=Lsup(u-u_ex)/Lsup(u_ex)
176 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
177 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
178 def test_solve(self):
179 domain=Brick(NE0,NE1,NE2,1)
180 x=Solution(domain).getX()
181 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
182 # --- set exact gradient -----------
183 g_ex=Data(0.,(3,),Solution(domain))
184 g_ex[0]=2.
185 g_ex[1]=3.
186 g_ex[2]=4.
187 # -------- test gradient --------------------------------
188 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
189 # -------- set-up PDE -----------------------------------
190 pde=LinearPDE(domain,numEquations=1)
191 mask=whereZero(x[0])
192 pde.setValue(r=u_ex,q=mask)
193 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
194 # -------- get the solution ---------------------------
195 pde.setTolerance(SOLVER_TOL)
196 pde.setSolverMethod(pde.PCG,pde.JACOBI)
197 pde.setSolverPackage(pde.PASO)
198 u=pde.getSolution()
199 # -------- test the solution ---------------------------
200 error=Lsup(u-u_ex)/Lsup(u_ex)
201 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
202 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
203 def test_solve(self):
204 domain=Brick(NE0,NE1,NE2,1)
205 x=Solution(domain).getX()
206 # --- set exact solution ----
207 u_ex=Vector(0,Solution(domain))
208 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
209 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
210 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
211 # --- set exact gradient -----------
212 g_ex=Data(0.,(3,3),Solution(domain))
213 g_ex[0,0]=2.
214 g_ex[0,1]=3.
215 g_ex[0,2]=4.
216 g_ex[1,0]=4.
217 g_ex[1,1]=1.
218 g_ex[1,2]=-2.
219 g_ex[2,0]=8.
220 g_ex[2,1]=4.
221 g_ex[2,2]=5.
222 # -------- test gradient --------------------------------
223 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
224 # -------- set-up PDE -----------------------------------
225 pde=LinearPDE(domain,numEquations=3)
226 mask=whereZero(x[0])
227 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
228 A=Tensor4(0,Function(domain))
229 A[0,:,0,:]=kronecker(3)
230 A[1,:,1,:]=kronecker(3)
231 A[2,:,2,:]=kronecker(3)
232 Y=Vector(0.,Function(domain))
233 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
234 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
235 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
236 pde.setValue(A=A,
237 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
238 Y=Y,
239 y=matrixmult(g_ex,domain.getNormal()))
240 # -------- get the solution ---------------------------
241 pde.setTolerance(SOLVER_TOL)
242 pde.setSolverMethod(pde.PCG,pde.JACOBI)
243 pde.setSolverPackage(pde.PASO)
244 u=pde.getSolution()
245 # -------- test the solution ---------------------------
246 error=Lsup(u-u_ex)/Lsup(u_ex)
247 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
248 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
249 def test_solve(self):
250 domain=Brick(NE0,NE1,NE2,2)
251 x=Solution(domain).getX()
252 # --- set exact solution ----
253 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
254 # --- set exact gradient -----------
255 g_ex=Data(0.,(3,),Solution(domain))
256 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
257 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
258 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
259 # -------- test gradient --------------------------------
260 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
261 # -------- set-up PDE -----------------------------------
262 pde=LinearPDE(domain,numEquations=1)
263 mask=whereZero(x[0])
264 pde.setValue(r=u_ex,q=mask)
265 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
266 # -------- get the solution ---------------------------
267 pde.setTolerance(SOLVER_TOL)
268 pde.setSolverMethod(pde.PCG,pde.JACOBI)
269 pde.setSolverPackage(pde.PASO)
270 u=pde.getSolution()
271 # -------- test the solution ---------------------------
272 error=Lsup(u-u_ex)/Lsup(u_ex)
273 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
274 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
275 def test_solve(self):
276 domain=Brick(NE0,NE1,NE2,2)
277 x=Solution(domain).getX()
278 # --- set exact solution ----
279 u_ex=Vector(0,Solution(domain))
280 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2
281 u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2
282 u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2
283 # --- set exact gradient -----------
284 g_ex=Data(0.,(3,3),Solution(domain))
285 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
286 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
287 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
288 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
289 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
290 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
291 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
292 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
293 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
294 # -------- test gradient --------------------------------
295 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
296 # -------- set-up PDE -----------------------------------
297 pde=LinearPDE(domain,numEquations=3)
298 mask=whereZero(x[0])
299 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
300 Y=Vector(0.,Function(domain))
301 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
302 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
303 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
304 A=Tensor4(0,Function(domain))
305 A[0,:,0,:]=kronecker(3)
306 A[1,:,1,:]=kronecker(3)
307 A[2,:,2,:]=kronecker(3)
308 pde.setValue(A=A,
309 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
310 Y=Y-numarray.array([60.,20.,22.]),
311 y=matrixmult(g_ex,domain.getNormal()))
312 # -------- get the solution ---------------------------
313 pde.setTolerance(SOLVER_TOL)
314 pde.setSolverMethod(pde.PCG,pde.JACOBI)
315 pde.setSolverPackage(pde.PASO)
316 u=pde.getSolution()
317 # -------- test the solution ---------------------------
318 error=Lsup(u-u_ex)/Lsup(u_ex)
319 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
320
321 if __name__ == '__main__':
322 suite = unittest.TestSuite()
323 # suite.addTest(unittest.makeSuite(SimpleSolve_Interval_Order1_SinglePDE_Paso_PCG_Jacobi))
324 # suite.addTest(unittest.makeSuite(SimpleSolve_Interval_Order1_SystemPDE_Paso_PCG_Jacobi))
325 # suite.addTest(unittest.makeSuite(SimpleSolve_Interval_Order2_SinglePDE_Paso_PCG_Jacobi))
326 # suite.addTest(unittest.makeSuite(SimpleSolve_Interval_Order2_SystemPDE_Paso_PCG_Jacobi))
327 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
328 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
329 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
330 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
331 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
332 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
333 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
334 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
335 s=unittest.TextTestRunner(verbosity=2).run(suite)

  ViewVC Help
Powered by ViewVC 1.1.26