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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26