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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4135 - (show annotations)
Mon Jan 14 04:54:14 2013 UTC (6 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 22998 byte(s)
Fixed test file.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Test suite for the linearPDE and pdetools test on ripley
25
26 :remark:
27
28 :var __author__: name of author
29 :var __licence__: licence agreement
30 :var __url__: url entry point on documentation
31 :var __version__: version
32 :var __date__: date of the version
33 """
34
35 __author__="Lutz Gross, l.gross@uq.edu.au"
36
37 import unittest, sys
38 from esys.escript import *
39 from esys.ripley import Rectangle,Brick
40 from esys.escript.linearPDEs import LinearPDE, SolverOptions
41 import numpy
42 SOLVER_VERBOSE=False
43
44 try:
45 RIPLEY_TEST_DATA=os.environ['RIPLEY_TEST_DATA']
46 except KeyError:
47 RIPLEY_TEST_DATA='.'
48
49 SOLVER_TOL=1.e-8
50 REL_TOL=1.e-6
51
52 FAC_DIAG=1.
53 FAC_OFFDIAG=-0.4
54
55 # number of elements in the spatial directions
56 NE0=12
57 NE1=12
58 NE2=8
59 mpiSize=getMPISizeWorld()
60 for x in [int(sqrt(mpiSize)),2,3,5,7,1]:
61 NX=x
62 NY=mpiSize//x
63 if NX*NY == mpiSize:
64 break
65
66 for x in [(int(mpiSize**(1/3.)),int(mpiSize**(1/3.))),(2,3),(2,2),(1,2),(1,1)]:
67 NXb=x[0]
68 NYb=x[1]
69 NZb=mpiSize//(x[0]*x[1])
70 if NXb*NYb*NZb == mpiSize:
71 break
72
73 class SimpleSolve_Rectangle_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
74 def test_solve(self):
75 # Tell about how many MPI CPUs and OpenMP threads
76 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
77 x=Solution(domain).getX()
78 # --- set exact solution ----
79 u_ex=Scalar(0,Solution(domain))
80 u_ex=1.+2.*x[0]+3.*x[1]
81 # --- set exact gradient -----------
82 g_ex=Data(0.,(2,),Solution(domain))
83 g_ex[0]=2.
84 g_ex[1]=3.
85 # -------- test gradient --------------------------------
86 g=grad(u_ex)
87 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
88 # -------- set-up PDE -----------------------------------
89 pde=LinearPDE(domain,numEquations=1)
90 mask=whereZero(x[0])
91 pde.setValue(r=u_ex,q=mask)
92 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
93 # -------- get the solution ---------------------------
94 pde.getSolverOptions().setTolerance(SOLVER_TOL)
95 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
96 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
97 pde.getSolverOptions().setPackage(SolverOptions.PASO)
98 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
99 u=pde.getSolution()
100 # -------- test the solution ---------------------------
101 error=Lsup(u-u_ex)
102 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
103 class SimpleSolve_Rectangle_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
104 def test_solve(self):
105 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
106 x=Solution(domain).getX()
107 # --- set exact solution ----
108 u_ex=Scalar(0,Solution(domain))
109 u_ex=1.+2.*x[0]+3.*x[1]
110 # --- set exact gradient -----------
111 g_ex=Data(0.,(2,),Solution(domain))
112 g_ex[0]=2.
113 g_ex[1]=3.
114 # -------- test gradient --------------------------------
115 g=grad(u_ex)
116 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
117 # -------- set-up PDE -----------------------------------
118 pde=LinearPDE(domain,numEquations=1)
119 mask=whereZero(x[0])
120 pde.setValue(r=u_ex,q=mask)
121 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
122 # -------- get the solution ---------------------------
123 pde.getSolverOptions().setTolerance(SOLVER_TOL)
124 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
125 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
126 pde.getSolverOptions().setPackage(SolverOptions.PASO)
127 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
128 u=pde.getSolution()
129 # -------- test the solution ---------------------------
130 error=Lsup(u-u_ex)
131 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
132 class SimpleSolve_Rectangle_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
133 def test_solve(self):
134 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
135 x=Solution(domain).getX()
136 # --- set exact solution ----
137 u_ex=Vector(0,Solution(domain))
138 u_ex[0]=1.+2.*x[0]+3.*x[1]
139 u_ex[1]=-1.+3.*x[0]+2.*x[1]
140 # --- set exact gradient -----------
141 g_ex=Data(0.,(2,2),Solution(domain))
142 g_ex[0,0]=2.
143 g_ex[0,1]=3.
144 g_ex[1,0]=3.
145 g_ex[1,1]=2.
146 # -------- test gradient --------------------------------
147 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
148 # -------- set-up PDE -----------------------------------
149 pde=LinearPDE(domain,numEquations=2)
150 mask=whereZero(x[0])
151 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
152 A=Tensor4(0,Function(domain))
153 A[0,:,0,:]=kronecker(2)
154 A[1,:,1,:]=kronecker(2)
155 Y=Vector(0.,Function(domain))
156 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
157 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
158 pde.setValue(A=A,
159 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
160 Y=Y,
161 y=matrixmult(g_ex,domain.getNormal()))
162 # -------- get the solution ---------------------------
163 pde.getSolverOptions().setTolerance(SOLVER_TOL)
164 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
165 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
166 pde.getSolverOptions().setPackage(SolverOptions.PASO)
167 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
168 u=pde.getSolution()
169 # -------- test the solution ---------------------------
170 error=Lsup(u-u_ex)
171 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
172 class SimpleSolve_Brick_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
173 def test_solve(self):
174 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
175 x=Solution(domain).getX()
176 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
177 # --- set exact gradient -----------
178 g_ex=Data(0.,(3,),Solution(domain))
179 g_ex[0]=2.
180 g_ex[1]=3.
181 g_ex[2]=4.
182 # -------- test gradient --------------------------------
183 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
184 # -------- set-up PDE -----------------------------------
185 pde=LinearPDE(domain,numEquations=1)
186 mask=whereZero(x[0])
187 pde.setValue(r=u_ex,q=mask)
188 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
189 # -------- get the solution ---------------------------
190 pde.getSolverOptions().setTolerance(SOLVER_TOL)
191 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
192 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
193 pde.getSolverOptions().setPackage(SolverOptions.PASO)
194 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
195 u=pde.getSolution()
196 # -------- test the solution ---------------------------
197 error=Lsup(u-u_ex)
198 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
199 class SimpleSolve_Brick_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
200 def test_solve(self):
201 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
202 x=Solution(domain).getX()
203 # --- set exact solution ----
204 u_ex=Vector(0,Solution(domain))
205 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
206 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
207 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
208 # --- set exact gradient -----------
209 g_ex=Data(0.,(3,3),Solution(domain))
210 g_ex[0,0]=2.
211 g_ex[0,1]=3.
212 g_ex[0,2]=4.
213 g_ex[1,0]=4.
214 g_ex[1,1]=1.
215 g_ex[1,2]=-2.
216 g_ex[2,0]=8.
217 g_ex[2,1]=4.
218 g_ex[2,2]=5.
219 # -------- test gradient --------------------------------
220 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
221 # -------- set-up PDE -----------------------------------
222 pde=LinearPDE(domain,numEquations=3)
223 mask=whereZero(x[0])
224 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
225 A=Tensor4(0,Function(domain))
226 A[0,:,0,:]=kronecker(3)
227 A[1,:,1,:]=kronecker(3)
228 A[2,:,2,:]=kronecker(3)
229 Y=Vector(0.,Function(domain))
230 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
231 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
232 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
233 pde.setValue(A=A,
234 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
235 Y=Y,
236 y=matrixmult(g_ex,domain.getNormal()))
237 # -------- get the solution ---------------------------
238 pde.getSolverOptions().setTolerance(SOLVER_TOL)
239 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
240 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
241 pde.getSolverOptions().setPackage(SolverOptions.PASO)
242 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
243 u=pde.getSolution()
244 # -------- test the solution ---------------------------
245 error=Lsup(u-u_ex)
246 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
247
248 class SimpleSolve_Rectangle_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
249 def test_solve(self):
250 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
251 x=Solution(domain).getX()
252 # --- set exact solution ----
253 u_ex=Scalar(0,Solution(domain))
254 u_ex=1.+2.*x[0]+3.*x[1]
255 # --- set exact gradient -----------
256 g_ex=Data(0.,(2,),Solution(domain))
257 g_ex[0]=2.
258 g_ex[1]=3.
259 # -------- test gradient --------------------------------
260 g=grad(u_ex)
261 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
262 # -------- set-up PDE -----------------------------------
263 pde=LinearPDE(domain,numEquations=1)
264 mask=whereZero(x[0])
265 pde.setValue(r=u_ex,q=mask)
266 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
267 # -------- get the solution ---------------------------
268 pde.getSolverOptions().setTolerance(SOLVER_TOL)
269 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
270 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
271 pde.getSolverOptions().setPackage(SolverOptions.PASO)
272 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
273 u=pde.getSolution()
274 # -------- test the solution ---------------------------
275 error=Lsup(u-u_ex)
276 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
277
278 class SimpleSolve_Rectangle_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
279 def test_solve(self):
280 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
281 x=Solution(domain).getX()
282 # --- set exact solution ----
283 u_ex=Vector(0,Solution(domain))
284 u_ex[0]=1.+2.*x[0]+3.*x[1]
285 u_ex[1]=-1.+3.*x[0]+2.*x[1]
286 # --- set exact gradient -----------
287 g_ex=Data(0.,(2,2),Solution(domain))
288 g_ex[0,0]=2.
289 g_ex[0,1]=3.
290 g_ex[1,0]=3.
291 g_ex[1,1]=2.
292 # -------- test gradient --------------------------------
293 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
294 # -------- set-up PDE -----------------------------------
295 pde=LinearPDE(domain,numEquations=2)
296 mask=whereZero(x[0])
297 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
298 A=Tensor4(0,Function(domain))
299 A[0,:,0,:]=kronecker(2)
300 A[1,:,1,:]=kronecker(2)
301 Y=Vector(0.,Function(domain))
302 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
303 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
304 pde.setValue(A=A,
305 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
306 Y=Y,
307 y=matrixmult(g_ex,domain.getNormal()))
308 # -------- get the solution ---------------------------
309 pde.getSolverOptions().setTolerance(SOLVER_TOL)
310 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
311 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
312 pde.getSolverOptions().setPackage(SolverOptions.PASO)
313 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
314 u=pde.getSolution()
315 # -------- test the solution ---------------------------
316 error=Lsup(u-u_ex)
317 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
318
319 class SimpleSolve_Brick_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
320 def test_solve(self):
321 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
322 x=Solution(domain).getX()
323 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
324 # --- set exact gradient -----------
325 g_ex=Data(0.,(3,),Solution(domain))
326 g_ex[0]=2.
327 g_ex[1]=3.
328 g_ex[2]=4.
329 # -------- test gradient --------------------------------
330 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
331 # -------- set-up PDE -----------------------------------
332 pde=LinearPDE(domain,numEquations=1)
333 mask=whereZero(x[0])
334 pde.setValue(r=u_ex,q=mask)
335 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
336 # -------- get the solution ---------------------------
337 pde.getSolverOptions().setTolerance(SOLVER_TOL)
338 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
339 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
340 pde.getSolverOptions().setPackage(SolverOptions.PASO)
341 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
342 u=pde.getSolution()
343 # -------- test the solution ---------------------------
344 error=Lsup(u-u_ex)
345 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
346
347 class SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
348 def test_solve(self):
349 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
350 x=Solution(domain).getX()
351 # --- set exact solution ----
352 u_ex=Scalar(0,Solution(domain))
353 u_ex=1.+2.*x[0]+3.*x[1]
354 # --- set exact gradient -----------
355 g_ex=Data(0.,(2,),Solution(domain))
356 g_ex[0]=2.
357 g_ex[1]=3.
358 # -------- test gradient --------------------------------
359 g=grad(u_ex)
360 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
361 # -------- set-up PDE -----------------------------------
362 pde=LinearPDE(domain,numEquations=1)
363 mask=whereZero(x[0])
364 pde.setValue(r=u_ex,q=mask)
365 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
366 # -------- get the solution ---------------------------
367 pde.getSolverOptions().setTolerance(SOLVER_TOL)
368 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
369 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
370 pde.getSolverOptions().setPackage(SolverOptions.PASO)
371
372 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
373 u=pde.getSolution()
374 # -------- test the solution ---------------------------
375 error=Lsup(u-u_ex)
376 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
377
378 class SimpleSolve_Rectangle_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
379 def test_solve(self):
380 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
381 x=Solution(domain).getX()
382 # --- set exact solution ----
383 u_ex=Vector(0,Solution(domain))
384 u_ex[0]=1.+2.*x[0]+3.*x[1]
385 u_ex[1]=-1.+3.*x[0]+2.*x[1]
386 # --- set exact gradient -----------
387 g_ex=Data(0.,(2,2),Solution(domain))
388 g_ex[0,0]=2.
389 g_ex[0,1]=3.
390 g_ex[1,0]=3.
391 g_ex[1,1]=2.
392 # -------- test gradient --------------------------------
393 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
394 # -------- set-up PDE -----------------------------------
395 pde=LinearPDE(domain,numEquations=2)
396 mask=whereZero(x[0])
397 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
398 A=Tensor4(0,Function(domain))
399 A[0,:,0,:]=kronecker(2)
400 A[1,:,1,:]=kronecker(2)
401 Y=Vector(0.,Function(domain))
402 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
403 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
404 pde.setValue(A=A,
405 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
406 Y=Y,
407 y=matrixmult(g_ex,domain.getNormal()))
408 # -------- get the solution ---------------------------
409 pde.getSolverOptions().setTolerance(SOLVER_TOL)
410 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
411 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
412 pde.getSolverOptions().setPackage(SolverOptions.PASO)
413 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
414 u=pde.getSolution()
415 # -------- test the solution ---------------------------
416 error=Lsup(u-u_ex)
417 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
418 class SimpleSolve_Brick_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
419 def test_solve(self):
420 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
421 x=Solution(domain).getX()
422 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
423 # --- set exact gradient -----------
424 g_ex=Data(0.,(3,),Solution(domain))
425 g_ex[0]=2.
426 g_ex[1]=3.
427 g_ex[2]=4.
428 # -------- test gradient --------------------------------
429 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
430 # -------- set-up PDE -----------------------------------
431 pde=LinearPDE(domain,numEquations=1)
432 mask=whereZero(x[0])
433 pde.setValue(r=u_ex,q=mask)
434 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
435 # -------- get the solution ---------------------------
436 pde.getSolverOptions().setTolerance(SOLVER_TOL)
437 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
438 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
439 pde.getSolverOptions().setPackage(SolverOptions.PASO)
440 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
441 u=pde.getSolution()
442 # -------- test the solution ---------------------------
443 error=Lsup(u-u_ex)
444 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
445
446 class SimpleSolve_Brick_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
447 def test_solve(self):
448 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
449 x=Solution(domain).getX()
450 # --- set exact solution ----
451 u_ex=Vector(0,Solution(domain))
452 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
453 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
454 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
455 # --- set exact gradient -----------
456 g_ex=Data(0.,(3,3),Solution(domain))
457 g_ex[0,0]=2.
458 g_ex[0,1]=3.
459 g_ex[0,2]=4.
460 g_ex[1,0]=4.
461 g_ex[1,1]=1.
462 g_ex[1,2]=-2.
463 g_ex[2,0]=8.
464 g_ex[2,1]=4.
465 g_ex[2,2]=5.
466 # -------- test gradient --------------------------------
467 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
468 # -------- set-up PDE -----------------------------------
469 pde=LinearPDE(domain,numEquations=3)
470 mask=whereZero(x[0])
471 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
472 A=Tensor4(0,Function(domain))
473 A[0,:,0,:]=kronecker(3)
474 A[1,:,1,:]=kronecker(3)
475 A[2,:,2,:]=kronecker(3)
476 Y=Vector(0.,Function(domain))
477 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
478 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
479 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
480 pde.setValue(A=A,
481 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
482 Y=Y,
483 y=matrixmult(g_ex,domain.getNormal()))
484 # -------- get the solution ---------------------------
485 pde.getSolverOptions().setTolerance(SOLVER_TOL)
486 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
487 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
488 pde.getSolverOptions().setPackage(SolverOptions.PASO)
489 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
490 u=pde.getSolution()
491 # -------- test the solution ---------------------------
492 error=Lsup(u-u_ex)
493 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
494
495
496 if __name__ == '__main__':
497 suite = unittest.TestSuite()
498 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_BICGSTAB_Jacobi))
499 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_PCG_Jacobi))
500 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SystemPDE_Paso_PCG_Jacobi))
501 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_PCG_Jacobi))
502 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SystemPDE_Paso_PCG_Jacobi))
503 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_TFQMR_Jacobi))
504 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SystemPDE_Paso_TFQMR_Jacobi))
505 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_TFQMR_Jacobi))
506 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi))
507 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi))
508 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_MINRES_Jacobi))
509 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SystemPDE_Paso_MINRES_Jacobi))
510
511 s=unittest.TextTestRunner(verbosity=2).run(suite)
512 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26