/[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 2344 - (show annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 45175 byte(s)
Change __url__ to launchpad site

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

  ViewVC Help
Powered by ViewVC 1.1.26