/[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 3803 - (show annotations)
Fri Feb 3 06:16:30 2012 UTC (7 years, 8 months ago) by caltinay
File MIME type: text/x-python
File size: 25232 byte(s)
Fixes to ripley unit test setup

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

  ViewVC Help
Powered by ViewVC 1.1.26