/[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 1787 - (show annotations)
Mon Sep 15 01:36:34 2008 UTC (12 years, 9 months ago) by artak
File MIME type: text/x-python
File size: 45294 byte(s)
MINRES solver is added to escript. Additional 16 tests are added to run_simplesolve for MINRES and TFQMR solvers
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 # Tell about how many MPI CPUs and OpenMP threads
74 printParallelThreadCounts()
75 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
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.failUnless(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.setTolerance(SOLVER_TOL)
94 pde.setSolverMethod(pde.BICGSTAB,pde.JACOBI)
95 pde.setSolverPackage(pde.PASO)
96 u=pde.getSolution(verbose=SOLVER_VERBOSE)
97 # -------- test the solution ---------------------------
98 error=Lsup(u-u_ex)/Lsup(u_ex)
99 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
100 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
101 def test_solve(self):
102 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
103 x=Solution(domain).getX()
104 # --- set exact solution ----
105 u_ex=Scalar(0,Solution(domain))
106 u_ex=1.+2.*x[0]+3.*x[1]
107 # --- set exact gradient -----------
108 g_ex=Data(0.,(2,),Solution(domain))
109 g_ex[0]=2.
110 g_ex[1]=3.
111 # -------- test gradient --------------------------------
112 g=grad(u_ex)
113 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
114 # -------- set-up PDE -----------------------------------
115 pde=LinearPDE(domain,numEquations=1)
116 mask=whereZero(x[0])
117 pde.setValue(r=u_ex,q=mask)
118 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
119 # -------- get the solution ---------------------------
120 pde.setTolerance(SOLVER_TOL)
121 pde.setSolverMethod(pde.PCG,pde.JACOBI)
122 pde.setSolverPackage(pde.PASO)
123 u=pde.getSolution(verbose=SOLVER_VERBOSE)
124 # -------- test the solution ---------------------------
125 error=Lsup(u-u_ex)/Lsup(u_ex)
126 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
127 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
128 def test_solve(self):
129 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
130 x=Solution(domain).getX()
131 # --- set exact solution ----
132 u_ex=Vector(0,Solution(domain))
133 u_ex[0]=1.+2.*x[0]+3.*x[1]
134 u_ex[1]=-1.+3.*x[0]+2.*x[1]
135 # --- set exact gradient -----------
136 g_ex=Data(0.,(2,2),Solution(domain))
137 g_ex[0,0]=2.
138 g_ex[0,1]=3.
139 g_ex[1,0]=3.
140 g_ex[1,1]=2.
141 # -------- test gradient --------------------------------
142 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
143 # -------- set-up PDE -----------------------------------
144 pde=LinearPDE(domain,numEquations=2)
145 mask=whereZero(x[0])
146 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
147 A=Tensor4(0,Function(domain))
148 A[0,:,0,:]=kronecker(2)
149 A[1,:,1,:]=kronecker(2)
150 Y=Vector(0.,Function(domain))
151 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
152 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
153 pde.setValue(A=A,
154 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
155 Y=Y,
156 y=matrixmult(g_ex,domain.getNormal()))
157 # -------- get the solution ---------------------------
158 pde.setTolerance(SOLVER_TOL)
159 pde.setSolverMethod(pde.PCG,pde.JACOBI)
160 pde.setSolverPackage(pde.PASO)
161 u=pde.getSolution(verbose=SOLVER_VERBOSE)
162 # -------- test the solution ---------------------------
163 error=Lsup(u-u_ex)/Lsup(u_ex)
164 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
165 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
166 def test_solve(self):
167 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
168 x=Solution(domain).getX()
169 # --- set exact solution ----
170 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
171 # --- set exact gradient -----------
172 g_ex=Data(0.,(2,),Solution(domain))
173 g_ex[0]=2.+8.*x[0]+5.*x[1]
174 g_ex[1]=3.+5.*x[0]+12.*x[1]
175 # -------- test gradient --------------------------------
176 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
177 # -------- set-up PDE -----------------------------------
178 pde=LinearPDE(domain,numEquations=1)
179 mask=whereZero(x[0])
180 pde.setValue(r=u_ex,q=mask)
181 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
182 # -------- get the solution ---------------------------
183 pde.setTolerance(SOLVER_TOL)
184 pde.setSolverMethod(pde.PCG,pde.JACOBI)
185 pde.setSolverPackage(pde.PASO)
186 u=pde.getSolution(verbose=SOLVER_VERBOSE)
187 # -------- test the solution ---------------------------
188 error=Lsup(u-u_ex)/Lsup(u_ex)
189 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
190 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
191 def test_solve(self):
192 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
193 x=Solution(domain).getX()
194 # --- set exact solution ----
195 u_ex=Vector(0,Solution(domain))
196 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
197 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
198 # --- set exact gradient -----------
199 g_ex=Data(0.,(2,2),Solution(domain))
200 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
201 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
202 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
203 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
204 # -------- test gradient --------------------------------
205 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
206 # -------- set-up PDE -----------------------------------
207 pde=LinearPDE(domain,numEquations=2)
208 mask=whereZero(x[0])
209 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
210 A=Tensor4(0,Function(domain))
211 A[0,:,0,:]=kronecker(2)
212 A[1,:,1,:]=kronecker(2)
213 Y=Vector(0.,Function(domain))
214 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
215 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
216 pde.setValue(A=A,
217 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
218 Y=Y-[20.,10.],
219 y=matrixmult(g_ex,domain.getNormal()))
220 # -------- get the solution ---------------------------
221 pde.setTolerance(SOLVER_TOL)
222 pde.setSolverMethod(pde.PCG,pde.JACOBI)
223 pde.setSolverPackage(pde.PASO)
224 u=pde.getSolution(verbose=SOLVER_VERBOSE)
225 # -------- test the solution ---------------------------
226 error=Lsup(u-u_ex)/Lsup(u_ex)
227 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
228 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
229 def test_solve(self):
230 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
231 x=Solution(domain).getX()
232 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
233 # --- set exact gradient -----------
234 g_ex=Data(0.,(3,),Solution(domain))
235 g_ex[0]=2.
236 g_ex[1]=3.
237 g_ex[2]=4.
238 # -------- test gradient --------------------------------
239 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
240 # -------- set-up PDE -----------------------------------
241 pde=LinearPDE(domain,numEquations=1)
242 mask=whereZero(x[0])
243 pde.setValue(r=u_ex,q=mask)
244 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
245 # -------- get the solution ---------------------------
246 pde.setTolerance(SOLVER_TOL)
247 pde.setSolverMethod(pde.PCG,pde.JACOBI)
248 pde.setSolverPackage(pde.PASO)
249 u=pde.getSolution(verbose=SOLVER_VERBOSE)
250 # -------- test the solution ---------------------------
251 error=Lsup(u-u_ex)/Lsup(u_ex)
252 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
253 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
254 def test_solve(self):
255 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
256 x=Solution(domain).getX()
257 # --- set exact solution ----
258 u_ex=Vector(0,Solution(domain))
259 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
260 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
261 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
262 # --- set exact gradient -----------
263 g_ex=Data(0.,(3,3),Solution(domain))
264 g_ex[0,0]=2.
265 g_ex[0,1]=3.
266 g_ex[0,2]=4.
267 g_ex[1,0]=4.
268 g_ex[1,1]=1.
269 g_ex[1,2]=-2.
270 g_ex[2,0]=8.
271 g_ex[2,1]=4.
272 g_ex[2,2]=5.
273 # -------- test gradient --------------------------------
274 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
275 # -------- set-up PDE -----------------------------------
276 pde=LinearPDE(domain,numEquations=3)
277 mask=whereZero(x[0])
278 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
279 A=Tensor4(0,Function(domain))
280 A[0,:,0,:]=kronecker(3)
281 A[1,:,1,:]=kronecker(3)
282 A[2,:,2,:]=kronecker(3)
283 Y=Vector(0.,Function(domain))
284 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
285 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
286 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
287 pde.setValue(A=A,
288 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
289 Y=Y,
290 y=matrixmult(g_ex,domain.getNormal()))
291 # -------- get the solution ---------------------------
292 pde.setTolerance(SOLVER_TOL)
293 pde.setSolverMethod(pde.PCG,pde.JACOBI)
294 pde.setSolverPackage(pde.PASO)
295 u=pde.getSolution(verbose=SOLVER_VERBOSE)
296 # -------- test the solution ---------------------------
297 error=Lsup(u-u_ex)/Lsup(u_ex)
298 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
299 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
300 def test_solve(self):
301 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
302 x=Solution(domain).getX()
303 # --- set exact solution ----
304 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
305 # --- set exact gradient -----------
306 g_ex=Data(0.,(3,),Solution(domain))
307 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
308 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
309 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
310 # -------- test gradient --------------------------------
311 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
312 # -------- set-up PDE -----------------------------------
313 pde=LinearPDE(domain,numEquations=1)
314 mask=whereZero(x[0])
315 pde.setValue(r=u_ex,q=mask)
316 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
317 # -------- get the solution ---------------------------
318 pde.setTolerance(SOLVER_TOL)
319 pde.setSolverMethod(pde.PCG,pde.JACOBI)
320 pde.setSolverPackage(pde.PASO)
321 u=pde.getSolution(verbose=SOLVER_VERBOSE)
322 # -------- test the solution ---------------------------
323 error=Lsup(u-u_ex)/Lsup(u_ex)
324 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
325 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
326 def test_solve(self):
327 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
328 x=Solution(domain).getX()
329 # --- set exact solution ----
330 u_ex=Vector(0,Solution(domain))
331 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
332 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
333 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
334 # --- set exact gradient -----------
335 g_ex=Data(0.,(3,3),Solution(domain))
336 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
337 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
338 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
339 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
340 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
341 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
342 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
343 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
344 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
345 # -------- test gradient --------------------------------
346 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
347 # -------- set-up PDE -----------------------------------
348 pde=LinearPDE(domain,numEquations=3)
349 mask=whereZero(x[0])
350 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
351 Y=Vector(0.,Function(domain))
352 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
353 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
354 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
355 A=Tensor4(0,Function(domain))
356 A[0,:,0,:]=kronecker(3)
357 A[1,:,1,:]=kronecker(3)
358 A[2,:,2,:]=kronecker(3)
359 pde.setValue(A=A,
360 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
361 Y=Y-numarray.array([60.,20.,22.]),
362 y=matrixmult(g_ex,domain.getNormal()))
363 # -------- get the solution ---------------------------
364 pde.setTolerance(SOLVER_TOL)
365 pde.setSolverMethod(pde.PCG,pde.JACOBI)
366 pde.setSolverPackage(pde.PASO)
367 u=pde.getSolution(verbose=SOLVER_VERBOSE)
368 # -------- test the solution ---------------------------
369 error=Lsup(u-u_ex)/Lsup(u_ex)
370 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
371
372 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
373 def test_solve(self):
374 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
375 x=Solution(domain).getX()
376 # --- set exact solution ----
377 u_ex=Scalar(0,Solution(domain))
378 u_ex=1.+2.*x[0]+3.*x[1]
379 # --- set exact gradient -----------
380 g_ex=Data(0.,(2,),Solution(domain))
381 g_ex[0]=2.
382 g_ex[1]=3.
383 # -------- test gradient --------------------------------
384 g=grad(u_ex)
385 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
386 # -------- set-up PDE -----------------------------------
387 pde=LinearPDE(domain,numEquations=1)
388 mask=whereZero(x[0])
389 pde.setValue(r=u_ex,q=mask)
390 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
391 # -------- get the solution ---------------------------
392 pde.setTolerance(SOLVER_TOL)
393 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
394 pde.setSolverPackage(pde.PASO)
395 u=pde.getSolution(verbose=SOLVER_VERBOSE)
396 # -------- test the solution ---------------------------
397 error=Lsup(u-u_ex)/Lsup(u_ex)
398 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
399
400 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
401 def test_solve(self):
402 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
403 x=Solution(domain).getX()
404 # --- set exact solution ----
405 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
406 # --- set exact gradient -----------
407 g_ex=Data(0.,(2,),Solution(domain))
408 g_ex[0]=2.+8.*x[0]+5.*x[1]
409 g_ex[1]=3.+5.*x[0]+12.*x[1]
410 # -------- test gradient --------------------------------
411 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
412 # -------- set-up PDE -----------------------------------
413 pde=LinearPDE(domain,numEquations=1)
414 mask=whereZero(x[0])
415 pde.setValue(r=u_ex,q=mask)
416 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
417 # -------- get the solution ---------------------------
418 pde.setTolerance(SOLVER_TOL)
419 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
420 pde.setSolverPackage(pde.PASO)
421 u=pde.getSolution(verbose=SOLVER_VERBOSE)
422 # -------- test the solution ---------------------------
423 error=Lsup(u-u_ex)/Lsup(u_ex)
424 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
425
426 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
427 def test_solve(self):
428 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
429 x=Solution(domain).getX()
430 # --- set exact solution ----
431 u_ex=Vector(0,Solution(domain))
432 u_ex[0]=1.+2.*x[0]+3.*x[1]
433 u_ex[1]=-1.+3.*x[0]+2.*x[1]
434 # --- set exact gradient -----------
435 g_ex=Data(0.,(2,2),Solution(domain))
436 g_ex[0,0]=2.
437 g_ex[0,1]=3.
438 g_ex[1,0]=3.
439 g_ex[1,1]=2.
440 # -------- test gradient --------------------------------
441 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
442 # -------- set-up PDE -----------------------------------
443 pde=LinearPDE(domain,numEquations=2)
444 mask=whereZero(x[0])
445 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
446 A=Tensor4(0,Function(domain))
447 A[0,:,0,:]=kronecker(2)
448 A[1,:,1,:]=kronecker(2)
449 Y=Vector(0.,Function(domain))
450 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
451 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
452 pde.setValue(A=A,
453 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
454 Y=Y,
455 y=matrixmult(g_ex,domain.getNormal()))
456 # -------- get the solution ---------------------------
457 pde.setTolerance(SOLVER_TOL)
458 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
459 pde.setSolverPackage(pde.PASO)
460 u=pde.getSolution(verbose=SOLVER_VERBOSE)
461 # -------- test the solution ---------------------------
462 error=Lsup(u-u_ex)/Lsup(u_ex)
463 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
464
465 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
466 def test_solve(self):
467 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
468 x=Solution(domain).getX()
469 # --- set exact solution ----
470 u_ex=Vector(0,Solution(domain))
471 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
472 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
473 # --- set exact gradient -----------
474 g_ex=Data(0.,(2,2),Solution(domain))
475 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
476 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
477 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
478 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
479 # -------- test gradient --------------------------------
480 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
481 # -------- set-up PDE -----------------------------------
482 pde=LinearPDE(domain,numEquations=2)
483 mask=whereZero(x[0])
484 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
485 A=Tensor4(0,Function(domain))
486 A[0,:,0,:]=kronecker(2)
487 A[1,:,1,:]=kronecker(2)
488 Y=Vector(0.,Function(domain))
489 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
490 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
491 pde.setValue(A=A,
492 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
493 Y=Y-[20.,10.],
494 y=matrixmult(g_ex,domain.getNormal()))
495 # -------- get the solution ---------------------------
496 pde.setTolerance(SOLVER_TOL)
497 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
498 pde.setSolverPackage(pde.PASO)
499 u=pde.getSolution(verbose=SOLVER_VERBOSE)
500 # -------- test the solution ---------------------------
501 error=Lsup(u-u_ex)/Lsup(u_ex)
502 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
503
504 class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
505 def test_solve(self):
506 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
507 x=Solution(domain).getX()
508 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
509 # --- set exact gradient -----------
510 g_ex=Data(0.,(3,),Solution(domain))
511 g_ex[0]=2.
512 g_ex[1]=3.
513 g_ex[2]=4.
514 # -------- test gradient --------------------------------
515 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
516 # -------- set-up PDE -----------------------------------
517 pde=LinearPDE(domain,numEquations=1)
518 mask=whereZero(x[0])
519 pde.setValue(r=u_ex,q=mask)
520 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
521 # -------- get the solution ---------------------------
522 pde.setTolerance(SOLVER_TOL)
523 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
524 pde.setSolverPackage(pde.PASO)
525 u=pde.getSolution(verbose=SOLVER_VERBOSE)
526 # -------- test the solution ---------------------------
527 error=Lsup(u-u_ex)/Lsup(u_ex)
528 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
529
530 class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
531 def test_solve(self):
532 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
533 x=Solution(domain).getX()
534 # --- set exact solution ----
535 u_ex=Vector(0,Solution(domain))
536 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
537 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
538 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
539 # --- set exact gradient -----------
540 g_ex=Data(0.,(3,3),Solution(domain))
541 g_ex[0,0]=2.
542 g_ex[0,1]=3.
543 g_ex[0,2]=4.
544 g_ex[1,0]=4.
545 g_ex[1,1]=1.
546 g_ex[1,2]=-2.
547 g_ex[2,0]=8.
548 g_ex[2,1]=4.
549 g_ex[2,2]=5.
550 # -------- test gradient --------------------------------
551 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
552 # -------- set-up PDE -----------------------------------
553 pde=LinearPDE(domain,numEquations=3)
554 mask=whereZero(x[0])
555 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
556 A=Tensor4(0,Function(domain))
557 A[0,:,0,:]=kronecker(3)
558 A[1,:,1,:]=kronecker(3)
559 A[2,:,2,:]=kronecker(3)
560 Y=Vector(0.,Function(domain))
561 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
562 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
563 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
564 pde.setValue(A=A,
565 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
566 Y=Y,
567 y=matrixmult(g_ex,domain.getNormal()))
568 # -------- get the solution ---------------------------
569 pde.setTolerance(SOLVER_TOL)
570 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
571 pde.setSolverPackage(pde.PASO)
572 u=pde.getSolution(verbose=SOLVER_VERBOSE)
573 # -------- test the solution ---------------------------
574 error=Lsup(u-u_ex)/Lsup(u_ex)
575 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
576
577 class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
578 def test_solve(self):
579 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
580 x=Solution(domain).getX()
581 # --- set exact solution ----
582 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
583 # --- set exact gradient -----------
584 g_ex=Data(0.,(3,),Solution(domain))
585 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
586 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
587 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
588 # -------- test gradient --------------------------------
589 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
590 # -------- set-up PDE -----------------------------------
591 pde=LinearPDE(domain,numEquations=1)
592 mask=whereZero(x[0])
593 pde.setValue(r=u_ex,q=mask)
594 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
595 # -------- get the solution ---------------------------
596 pde.setTolerance(SOLVER_TOL)
597 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
598 pde.setSolverPackage(pde.PASO)
599 u=pde.getSolution(verbose=SOLVER_VERBOSE)
600 # -------- test the solution ---------------------------
601 error=Lsup(u-u_ex)/Lsup(u_ex)
602 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
603
604 class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
605 def test_solve(self):
606 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
607 x=Solution(domain).getX()
608 # --- set exact solution ----
609 u_ex=Vector(0,Solution(domain))
610 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
611 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
612 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
613 # --- set exact gradient -----------
614 g_ex=Data(0.,(3,3),Solution(domain))
615 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
616 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
617 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
618 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
619 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
620 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
621 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
622 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
623 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
624 # -------- test gradient --------------------------------
625 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
626 # -------- set-up PDE -----------------------------------
627 pde=LinearPDE(domain,numEquations=3)
628 mask=whereZero(x[0])
629 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
630 Y=Vector(0.,Function(domain))
631 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
632 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
633 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
634 A=Tensor4(0,Function(domain))
635 A[0,:,0,:]=kronecker(3)
636 A[1,:,1,:]=kronecker(3)
637 A[2,:,2,:]=kronecker(3)
638 pde.setValue(A=A,
639 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
640 Y=Y-numarray.array([60.,20.,22.]),
641 y=matrixmult(g_ex,domain.getNormal()))
642 # -------- get the solution ---------------------------
643 pde.setTolerance(SOLVER_TOL)
644 pde.setSolverMethod(pde.TFQMR,pde.JACOBI)
645 pde.setSolverPackage(pde.PASO)
646 u=pde.getSolution(verbose=SOLVER_VERBOSE)
647 # -------- test the solution ---------------------------
648 error=Lsup(u-u_ex)/Lsup(u_ex)
649 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
650
651
652 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
653 def test_solve(self):
654 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
655 x=Solution(domain).getX()
656 # --- set exact solution ----
657 u_ex=Scalar(0,Solution(domain))
658 u_ex=1.+2.*x[0]+3.*x[1]
659 # --- set exact gradient -----------
660 g_ex=Data(0.,(2,),Solution(domain))
661 g_ex[0]=2.
662 g_ex[1]=3.
663 # -------- test gradient --------------------------------
664 g=grad(u_ex)
665 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
666 # -------- set-up PDE -----------------------------------
667 pde=LinearPDE(domain,numEquations=1)
668 mask=whereZero(x[0])
669 pde.setValue(r=u_ex,q=mask)
670 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
671 # -------- get the solution ---------------------------
672 pde.setTolerance(SOLVER_TOL)
673 pde.setSolverMethod(pde.MINRES,pde.JACOBI)
674 pde.setSolverPackage(pde.PASO)
675
676 u=pde.getSolution(verbose=SOLVER_VERBOSE)
677 # -------- test the solution ---------------------------
678 error=Lsup(u-u_ex)/Lsup(u_ex)
679 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
680
681 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
682 def test_solve(self):
683 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
684 x=Solution(domain).getX()
685 # --- set exact solution ----
686 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
687 # --- set exact gradient -----------
688 g_ex=Data(0.,(2,),Solution(domain))
689 g_ex[0]=2.+8.*x[0]+5.*x[1]
690 g_ex[1]=3.+5.*x[0]+12.*x[1]
691 # -------- test gradient --------------------------------
692 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
693 # -------- set-up PDE -----------------------------------
694 pde=LinearPDE(domain,numEquations=1)
695 mask=whereZero(x[0])
696 pde.setValue(r=u_ex,q=mask)
697 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
698 # -------- get the solution ---------------------------
699 pde.setTolerance(SOLVER_TOL)
700 pde.setSolverMethod(pde.MINRES,pde.JACOBI)
701 pde.setSolverPackage(pde.PASO)
702 u=pde.getSolution(verbose=SOLVER_VERBOSE)
703 # -------- test the solution ---------------------------
704 error=Lsup(u-u_ex)/Lsup(u_ex)
705 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
706
707 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
708 def test_solve(self):
709 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
710 x=Solution(domain).getX()
711 # --- set exact solution ----
712 u_ex=Vector(0,Solution(domain))
713 u_ex[0]=1.+2.*x[0]+3.*x[1]
714 u_ex[1]=-1.+3.*x[0]+2.*x[1]
715 # --- set exact gradient -----------
716 g_ex=Data(0.,(2,2),Solution(domain))
717 g_ex[0,0]=2.
718 g_ex[0,1]=3.
719 g_ex[1,0]=3.
720 g_ex[1,1]=2.
721 # -------- test gradient --------------------------------
722 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
723 # -------- set-up PDE -----------------------------------
724 pde=LinearPDE(domain,numEquations=2)
725 mask=whereZero(x[0])
726 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
727 A=Tensor4(0,Function(domain))
728 A[0,:,0,:]=kronecker(2)
729 A[1,:,1,:]=kronecker(2)
730 Y=Vector(0.,Function(domain))
731 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
732 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
733 pde.setValue(A=A,
734 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
735 Y=Y,
736 y=matrixmult(g_ex,domain.getNormal()))
737 # -------- get the solution ---------------------------
738 pde.setTolerance(SOLVER_TOL)
739 pde.setSolverMethod(pde.MINRES,pde.JACOBI)
740 pde.setSolverPackage(pde.PASO)
741 u=pde.getSolution(verbose=SOLVER_VERBOSE)
742 # -------- test the solution ---------------------------
743 error=Lsup(u-u_ex)/Lsup(u_ex)
744 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
745
746 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_RILU(unittest.TestCase):
747 def test_solve(self):
748 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
749 x=Solution(domain).getX()
750 # --- set exact solution ----
751 u_ex=Vector(0,Solution(domain))
752 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
753 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
754 # --- set exact gradient -----------
755 g_ex=Data(0.,(2,2),Solution(domain))
756 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
757 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
758 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
759 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
760 # -------- test gradient --------------------------------
761 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
762 # -------- set-up PDE -----------------------------------
763 pde=LinearPDE(domain,numEquations=2)
764 mask=whereZero(x[0])
765 pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
766 A=Tensor4(0,Function(domain))
767 A[0,:,0,:]=kronecker(2)
768 A[1,:,1,:]=kronecker(2)
769 Y=Vector(0.,Function(domain))
770 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
771 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
772 pde.setValue(A=A,
773 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
774 Y=Y-[20.,10.],
775 y=matrixmult(g_ex,domain.getNormal()))
776 # -------- get the solution ---------------------------
777 pde.setTolerance(SOLVER_TOL)
778 pde.setSolverMethod(pde.MINRES,pde.RILU)
779 pde.setSolverPackage(pde.PASO)
780 u=pde.getSolution(verbose=SOLVER_VERBOSE)
781 # -------- test the solution ---------------------------
782 error=Lsup(u-u_ex)/Lsup(u_ex)
783 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
784
785 class SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
786 def test_solve(self):
787 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
788 x=Solution(domain).getX()
789 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
790 # --- set exact gradient -----------
791 g_ex=Data(0.,(3,),Solution(domain))
792 g_ex[0]=2.
793 g_ex[1]=3.
794 g_ex[2]=4.
795 # -------- test gradient --------------------------------
796 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
797 # -------- set-up PDE -----------------------------------
798 pde=LinearPDE(domain,numEquations=1)
799 mask=whereZero(x[0])
800 pde.setValue(r=u_ex,q=mask)
801 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
802 # -------- get the solution ---------------------------
803 pde.setTolerance(SOLVER_TOL)
804 pde.setSolverMethod(pde.MINRES,pde.JACOBI)
805 pde.setSolverPackage(pde.PASO)
806 u=pde.getSolution(verbose=SOLVER_VERBOSE)
807 # -------- test the solution ---------------------------
808 error=Lsup(u-u_ex)/Lsup(u_ex)
809 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
810
811 class SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
812 def test_solve(self):
813 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
814 x=Solution(domain).getX()
815 # --- set exact solution ----
816 u_ex=Vector(0,Solution(domain))
817 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
818 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
819 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
820 # --- set exact gradient -----------
821 g_ex=Data(0.,(3,3),Solution(domain))
822 g_ex[0,0]=2.
823 g_ex[0,1]=3.
824 g_ex[0,2]=4.
825 g_ex[1,0]=4.
826 g_ex[1,1]=1.
827 g_ex[1,2]=-2.
828 g_ex[2,0]=8.
829 g_ex[2,1]=4.
830 g_ex[2,2]=5.
831 # -------- test gradient --------------------------------
832 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
833 # -------- set-up PDE -----------------------------------
834 pde=LinearPDE(domain,numEquations=3)
835 mask=whereZero(x[0])
836 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
837 A=Tensor4(0,Function(domain))
838 A[0,:,0,:]=kronecker(3)
839 A[1,:,1,:]=kronecker(3)
840 A[2,:,2,:]=kronecker(3)
841 Y=Vector(0.,Function(domain))
842 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
843 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
844 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
845 pde.setValue(A=A,
846 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
847 Y=Y,
848 y=matrixmult(g_ex,domain.getNormal()))
849 # -------- get the solution ---------------------------
850 pde.setTolerance(SOLVER_TOL)
851 pde.setSolverMethod(pde.MINRES,pde.JACOBI)
852 pde.setSolverPackage(pde.PASO)
853 u=pde.getSolution(verbose=SOLVER_VERBOSE)
854 # -------- test the solution ---------------------------
855 error=Lsup(u-u_ex)/Lsup(u_ex)
856 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
857
858 class SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_RILU(unittest.TestCase):
859 def test_solve(self):
860 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
861 x=Solution(domain).getX()
862 # --- set exact solution ----
863 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
864 # --- set exact gradient -----------
865 g_ex=Data(0.,(3,),Solution(domain))
866 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
867 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
868 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
869 # -------- test gradient --------------------------------
870 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
871 # -------- set-up PDE -----------------------------------
872 pde=LinearPDE(domain,numEquations=1)
873 mask=whereZero(x[0])
874 pde.setValue(r=u_ex,q=mask)
875 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
876 # -------- get the solution ---------------------------
877 pde.setTolerance(SOLVER_TOL)
878 pde.setSolverMethod(pde.MINRES,pde.RILU)
879 pde.setSolverPackage(pde.PASO)
880 u=pde.getSolution(verbose=SOLVER_VERBOSE)
881 # -------- test the solution ---------------------------
882 error=Lsup(u-u_ex)/Lsup(u_ex)
883 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
884
885 class SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_RILU(unittest.TestCase):
886 def test_solve(self):
887 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
888 x=Solution(domain).getX()
889 # --- set exact solution ----
890 u_ex=Vector(0,Solution(domain))
891 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
892 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
893 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
894 # --- set exact gradient -----------
895 g_ex=Data(0.,(3,3),Solution(domain))
896 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
897 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
898 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
899 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
900 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
901 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
902 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
903 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
904 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
905 # -------- test gradient --------------------------------
906 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
907 # -------- set-up PDE -----------------------------------
908 pde=LinearPDE(domain,numEquations=3)
909 mask=whereZero(x[0])
910 pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
911 Y=Vector(0.,Function(domain))
912 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
913 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
914 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
915 A=Tensor4(0,Function(domain))
916 A[0,:,0,:]=kronecker(3)
917 A[1,:,1,:]=kronecker(3)
918 A[2,:,2,:]=kronecker(3)
919 pde.setValue(A=A,
920 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
921 Y=Y-numarray.array([60.,20.,22.]),
922 y=matrixmult(g_ex,domain.getNormal()))
923 # -------- get the solution ---------------------------
924 pde.setTolerance(SOLVER_TOL)
925 pde.setSolverMethod(pde.MINRES,pde.RILU)
926 pde.setSolverPackage(pde.PASO)
927 u=pde.getSolution(verbose=SOLVER_VERBOSE)
928 # -------- test the solution ---------------------------
929 error=Lsup(u-u_ex)/Lsup(u_ex)
930 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
931
932 if __name__ == '__main__':
933 suite = unittest.TestSuite()
934 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
935 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
936 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
937 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
938 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
939 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
940 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
941 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
942 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
943 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_Jacobi))
944 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_Jacobi))
945 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_Jacobi))
946 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_Jacobi))
947 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_Jacobi))
948 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_Jacobi))
949 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_Jacobi))
950 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_Jacobi))
951 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi))
952 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_Jacobi))
953 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi))
954 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_RILU))
955 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_Jacobi))
956 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi))
957 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_RILU))
958 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_RILU))
959
960 s=unittest.TextTestRunner(verbosity=2).run(suite)
961 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26