/[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 1556 - (show annotations)
Mon May 12 00:54:58 2008 UTC (12 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 17131 byte(s)
Modification to allow mixed mode execution. 
In order to keep the code portable accross platform all MPI calls within
parallel regions have been moved. 


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 # setNumberOfThreads(2)
47
48 try:
49 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
50 except KeyError:
51 FINLEY_TEST_DATA='.'
52
53 FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/"
54
55 # number of elements in the spatial directions
56 NE0=8
57 NE1=10
58 NE2=12
59
60 NE0=12
61 NE1=12
62 NE2=8
63
64 SOLVER_TOL=1.e-8
65 REL_TOL=1.e-6
66
67 FAC_DIAG=1.
68 FAC_OFFDIAG=-0.4
69
70
71 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
72 def test_solve(self):
73 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
74 x=Solution(domain).getX()
75 # --- set exact solution ----
76 u_ex=Scalar(0,Solution(domain))
77 u_ex=1.+2.*x[0]+3.*x[1]
78 # --- set exact gradient -----------
79 g_ex=Data(0.,(2,),Solution(domain))
80 g_ex[0]=2.
81 g_ex[1]=3.
82 # -------- test gradient --------------------------------
83 g=grad(u_ex)
84 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
85 # -------- set-up PDE -----------------------------------
86 pde=LinearPDE(domain,numEquations=1)
87 mask=whereZero(x[0])
88 pde.setValue(r=u_ex,q=mask)
89 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
90 # -------- get the solution ---------------------------
91 pde.setTolerance(SOLVER_TOL)
92 pde.setSolverMethod(pde.BICGSTAB,pde.JACOBI)
93 pde.setSolverPackage(pde.PASO)
94 u=pde.getSolution(verbose=SOLVER_VERBOSE)
95 # -------- test the solution ---------------------------
96 error=Lsup(u-u_ex)/Lsup(u_ex)
97 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
98 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
99 def test_solve(self):
100 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
101 x=Solution(domain).getX()
102 # --- set exact solution ----
103 u_ex=Scalar(0,Solution(domain))
104 u_ex=1.+2.*x[0]+3.*x[1]
105 # --- set exact gradient -----------
106 g_ex=Data(0.,(2,),Solution(domain))
107 g_ex[0]=2.
108 g_ex[1]=3.
109 # -------- test gradient --------------------------------
110 g=grad(u_ex)
111 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
112 # -------- set-up PDE -----------------------------------
113 pde=LinearPDE(domain,numEquations=1)
114 mask=whereZero(x[0])
115 pde.setValue(r=u_ex,q=mask)
116 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
117 # -------- get the solution ---------------------------
118 pde.setTolerance(SOLVER_TOL)
119 pde.setSolverMethod(pde.PCG,pde.JACOBI)
120 pde.setSolverPackage(pde.PASO)
121 u=pde.getSolution(verbose=SOLVER_VERBOSE)
122 # -------- test the solution ---------------------------
123 error=Lsup(u-u_ex)/Lsup(u_ex)
124 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
125 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
126 def test_solve(self):
127 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
128 x=Solution(domain).getX()
129 # --- set exact solution ----
130 u_ex=Vector(0,Solution(domain))
131 u_ex[0]=1.+2.*x[0]+3.*x[1]
132 u_ex[1]=-1.+3.*x[0]+2.*x[1]
133 # --- set exact gradient -----------
134 g_ex=Data(0.,(2,2),Solution(domain))
135 g_ex[0,0]=2.
136 g_ex[0,1]=3.
137 g_ex[1,0]=3.
138 g_ex[1,1]=2.
139 # -------- test gradient --------------------------------
140 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
141 # -------- set-up PDE -----------------------------------
142 pde=LinearPDE(domain,numEquations=2)
143 mask=whereZero(x[0])
144 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
145 A=Tensor4(0,Function(domain))
146 A[0,:,0,:]=kronecker(2)
147 A[1,:,1,:]=kronecker(2)
148 Y=Vector(0.,Function(domain))
149 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
150 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
151 pde.setValue(A=A,
152 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
153 Y=Y,
154 y=matrixmult(g_ex,domain.getNormal()))
155 # -------- get the solution ---------------------------
156 pde.setTolerance(SOLVER_TOL)
157 pde.setSolverMethod(pde.PCG,pde.JACOBI)
158 pde.setSolverPackage(pde.PASO)
159 u=pde.getSolution(verbose=SOLVER_VERBOSE)
160 # -------- test the solution ---------------------------
161 error=Lsup(u-u_ex)/Lsup(u_ex)
162 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
163 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
164 def test_solve(self):
165 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
166 x=Solution(domain).getX()
167 # --- set exact solution ----
168 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
169 # --- set exact gradient -----------
170 g_ex=Data(0.,(2,),Solution(domain))
171 g_ex[0]=2.+8.*x[0]+5.*x[1]
172 g_ex[1]=3.+5.*x[0]+12.*x[1]
173 # -------- test gradient --------------------------------
174 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
175 # -------- set-up PDE -----------------------------------
176 pde=LinearPDE(domain,numEquations=1)
177 mask=whereZero(x[0])
178 pde.setValue(r=u_ex,q=mask)
179 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
180 # -------- get the solution ---------------------------
181 pde.setTolerance(SOLVER_TOL)
182 pde.setSolverMethod(pde.PCG,pde.JACOBI)
183 pde.setSolverPackage(pde.PASO)
184 u=pde.getSolution(verbose=SOLVER_VERBOSE)
185 # -------- test the solution ---------------------------
186 error=Lsup(u-u_ex)/Lsup(u_ex)
187 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
188 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
189 def test_solve(self):
190 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
191 x=Solution(domain).getX()
192 # --- set exact solution ----
193 u_ex=Vector(0,Solution(domain))
194 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
195 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
196 # --- set exact gradient -----------
197 g_ex=Data(0.,(2,2),Solution(domain))
198 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
199 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
200 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
201 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
202 # -------- test gradient --------------------------------
203 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
204 # -------- set-up PDE -----------------------------------
205 pde=LinearPDE(domain,numEquations=2)
206 mask=whereZero(x[0])
207 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
208 A=Tensor4(0,Function(domain))
209 A[0,:,0,:]=kronecker(2)
210 A[1,:,1,:]=kronecker(2)
211 Y=Vector(0.,Function(domain))
212 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
213 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
214 pde.setValue(A=A,
215 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
216 Y=Y-[20.,10.],
217 y=matrixmult(g_ex,domain.getNormal()))
218 # -------- get the solution ---------------------------
219 pde.setTolerance(SOLVER_TOL)
220 pde.setSolverMethod(pde.PCG,pde.JACOBI)
221 pde.setSolverPackage(pde.PASO)
222 u=pde.getSolution(verbose=SOLVER_VERBOSE)
223 # -------- test the solution ---------------------------
224 error=Lsup(u-u_ex)/Lsup(u_ex)
225 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
226 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
227 def test_solve(self):
228 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
229 x=Solution(domain).getX()
230 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
231 # --- set exact gradient -----------
232 g_ex=Data(0.,(3,),Solution(domain))
233 g_ex[0]=2.
234 g_ex[1]=3.
235 g_ex[2]=4.
236 # -------- test gradient --------------------------------
237 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
238 # -------- set-up PDE -----------------------------------
239 pde=LinearPDE(domain,numEquations=1)
240 mask=whereZero(x[0])
241 pde.setValue(r=u_ex,q=mask)
242 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
243 # -------- get the solution ---------------------------
244 pde.setTolerance(SOLVER_TOL)
245 pde.setSolverMethod(pde.PCG,pde.JACOBI)
246 pde.setSolverPackage(pde.PASO)
247 u=pde.getSolution(verbose=SOLVER_VERBOSE)
248 # -------- test the solution ---------------------------
249 error=Lsup(u-u_ex)/Lsup(u_ex)
250 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
251 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
252 def test_solve(self):
253 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
254 x=Solution(domain).getX()
255 # --- set exact solution ----
256 u_ex=Vector(0,Solution(domain))
257 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
258 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
259 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
260 # --- set exact gradient -----------
261 g_ex=Data(0.,(3,3),Solution(domain))
262 g_ex[0,0]=2.
263 g_ex[0,1]=3.
264 g_ex[0,2]=4.
265 g_ex[1,0]=4.
266 g_ex[1,1]=1.
267 g_ex[1,2]=-2.
268 g_ex[2,0]=8.
269 g_ex[2,1]=4.
270 g_ex[2,2]=5.
271 # -------- test gradient --------------------------------
272 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
273 # -------- set-up PDE -----------------------------------
274 pde=LinearPDE(domain,numEquations=3)
275 mask=whereZero(x[0])
276 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
277 A=Tensor4(0,Function(domain))
278 A[0,:,0,:]=kronecker(3)
279 A[1,:,1,:]=kronecker(3)
280 A[2,:,2,:]=kronecker(3)
281 Y=Vector(0.,Function(domain))
282 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
283 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
284 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
285 pde.setValue(A=A,
286 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
287 Y=Y,
288 y=matrixmult(g_ex,domain.getNormal()))
289 # -------- get the solution ---------------------------
290 pde.setTolerance(SOLVER_TOL)
291 pde.setSolverMethod(pde.PCG,pde.JACOBI)
292 pde.setSolverPackage(pde.PASO)
293 u=pde.getSolution(verbose=SOLVER_VERBOSE)
294 # -------- test the solution ---------------------------
295 error=Lsup(u-u_ex)/Lsup(u_ex)
296 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
297 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
298 def test_solve(self):
299 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
300 x=Solution(domain).getX()
301 # --- set exact solution ----
302 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
303 # --- set exact gradient -----------
304 g_ex=Data(0.,(3,),Solution(domain))
305 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
306 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
307 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
308 # -------- test gradient --------------------------------
309 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
310 # -------- set-up PDE -----------------------------------
311 pde=LinearPDE(domain,numEquations=1)
312 mask=whereZero(x[0])
313 pde.setValue(r=u_ex,q=mask)
314 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
315 # -------- get the solution ---------------------------
316 pde.setTolerance(SOLVER_TOL)
317 pde.setSolverMethod(pde.PCG,pde.JACOBI)
318 pde.setSolverPackage(pde.PASO)
319 u=pde.getSolution(verbose=SOLVER_VERBOSE)
320 # -------- test the solution ---------------------------
321 error=Lsup(u-u_ex)/Lsup(u_ex)
322 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
323 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
324 def test_solve(self):
325 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
326 x=Solution(domain).getX()
327 # --- set exact solution ----
328 u_ex=Vector(0,Solution(domain))
329 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
330 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
331 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
332 # --- set exact gradient -----------
333 g_ex=Data(0.,(3,3),Solution(domain))
334 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
335 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
336 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
337 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
338 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
339 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
340 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
341 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
342 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
343 # -------- test gradient --------------------------------
344 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
345 # -------- set-up PDE -----------------------------------
346 pde=LinearPDE(domain,numEquations=3)
347 mask=whereZero(x[0])
348 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
349 Y=Vector(0.,Function(domain))
350 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
351 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
352 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
353 A=Tensor4(0,Function(domain))
354 A[0,:,0,:]=kronecker(3)
355 A[1,:,1,:]=kronecker(3)
356 A[2,:,2,:]=kronecker(3)
357 pde.setValue(A=A,
358 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
359 Y=Y-numarray.array([60.,20.,22.]),
360 y=matrixmult(g_ex,domain.getNormal()))
361 # -------- get the solution ---------------------------
362 pde.setTolerance(SOLVER_TOL)
363 pde.setSolverMethod(pde.PCG,pde.JACOBI)
364 pde.setSolverPackage(pde.PASO)
365 u=pde.getSolution(verbose=SOLVER_VERBOSE)
366 # -------- test the solution ---------------------------
367 error=Lsup(u-u_ex)/Lsup(u_ex)
368 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
369
370 if __name__ == '__main__':
371 suite = unittest.TestSuite()
372 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
373 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
374 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
375 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
376 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
377 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
378 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
379 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
380 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
381 s=unittest.TextTestRunner(verbosity=2).run(suite)
382 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26