/[escript]/trunk/finley/test/python/run_amg.py
ViewVC logotype

Contents of /trunk/finley/test/python/run_amg.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3436 - (show annotations)
Mon Jan 10 02:06:07 2011 UTC (8 years, 7 months ago) by plaub
File MIME type: text/x-python
File size: 24686 byte(s)
Added/modified AMG tests:

  - "self.failUnless(...)" became "self.assertTrue(...)"
    because fail* is soon to be deprecated,
  - fixed logical error in self.assertTrue argument,
  - added another test, 'testSquare'.

Added some print statements to AMG.c in verbose mode to
identify which of the stopping conditions were fulfilled.

Changed default AMG parameters for 'setNumPreSweeps' and 
'setNumPostSweeps' to 1 (huge speed improvment).

P.S Beware of bugs in the above code; I have only proved 
it correct, not tried it... Jokes!


1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Test suite for AMG
25
26 :remark:
27
28 :var __author__: name of author
29 :var __licence__: licence agreement
30 :var __url__: url entry point on documentation
31 :var __version__: version
32 :var __date__: date of the version
33 """
34
35 __author__="Lutz Gross, l.gross@uq.edu.au"
36
37 import unittest, sys
38
39 from esys.escript import *
40 from esys.finley import Rectangle,Brick
41 from esys.escript.linearPDEs import LinearPDE, SolverOptions
42 import numpy
43 OPTIMIZE=True
44 SOLVER_VERBOSE=True
45
46 MIN_MATRIX_SIZE=None
47 MIN_SPARSITY=None
48 MAX_LEVEL=None
49 USE_AMG=True or False
50
51 try:
52 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
53 except KeyError:
54 FINLEY_TEST_DATA='.'
55 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
56
57 # number of elements in the spatial directions
58 NE_TOTAL=4096
59
60 class Test_AMG(unittest.TestCase):
61
62 def test_Poisson(self):
63 x=Solution(self.domain).getX()
64 # --- set exact solution ----
65 u_ex=Scalar(1,Solution(self.domain))
66 g_ex=Vector(0.,Solution(self.domain))
67 for i in xrange(self.domain.getDim()):
68 u_ex+=(i+1)*x[i]
69 g_ex[i]=(i+1)
70
71 # create PDE:
72 pde=LinearPDE(self.domain,numEquations=1)
73 pde.setSymmetryOn()
74 mask=whereZero(x[0])
75 pde.setValue(r=u_ex,q=mask)
76 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
77 # -------- get the solution ---------------------------
78 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
79 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
80 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
81 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
82 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
83 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
84 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
85
86 u=pde.getSolution()
87 # -------- test the solution ---------------------------
88 error=Lsup(u-u_ex)/Lsup(u_ex)
89 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
90
91 def test_PoissonSqueezedX(self):
92 x=self.domain.getX().copy()
93 x[0]*=0.5
94 self.domain.setX(x)
95 x=Solution(self.domain).getX()
96 # --- set exact solution ----
97 u_ex=Scalar(1,Solution(self.domain))
98 g_ex=Vector(0.,Solution(self.domain))
99 for i in xrange(self.domain.getDim()):
100 u_ex+=(i+1)*x[i]
101 g_ex[i]=(i+1)
102
103 # create PDE:
104 pde=LinearPDE(self.domain,numEquations=1)
105 pde.setSymmetryOn()
106 mask=whereZero(x[0])
107 pde.setValue(r=u_ex,q=mask)
108 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
109 # -------- get the solution ---------------------------
110 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
111 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
112 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
113 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
114 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
115 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
116 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
117
118 u=pde.getSolution()
119 # -------- test the solution ---------------------------
120 error=Lsup(u-u_ex)/Lsup(u_ex)
121 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
122
123
124 def test_Poisson2(self):
125 x=Solution(self.domain).getX()
126 # --- set exact solution ----
127 u_ex=Data(1.,(2,),Solution(self.domain))
128 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
129 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
130 for i in xrange(self.domain.getDim()):
131 u_ex[0]+= 1*(i+1) *x[i]
132 g_ex[0,i]=1*(i+1)
133 u_ex[1]+= 2*(i+1)*x[i]
134 g_ex[1,i]=2*(i+1)
135 A[0,i,0,i]=1
136 A[1,i,1,i]=1
137
138 # create PDE:
139 pde=LinearPDE(self.domain,numEquations=2)
140 pde.setSymmetryOn()
141 mask=whereZero(x[0])*[1,1]
142 pde.setValue(r=u_ex,q=mask)
143 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
144 # -------- get the solution ---------------------------
145 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
146 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
147 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
148 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
149 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
150 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
151 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
152
153 u=pde.getSolution()
154 # -------- test the solution ---------------------------
155 error=Lsup(u-u_ex)/Lsup(u_ex)
156 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
157
158 def test_Poisson3(self):
159 x=Solution(self.domain).getX()
160 # --- set exact solution ----
161 u_ex=Data(1.,(3,),Solution(self.domain))
162 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
163 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
164 for i in xrange(self.domain.getDim()):
165 u_ex[0]+= 1*(i+1) *x[i]
166 g_ex[0,i]=1*(i+1)
167 u_ex[1]+= 2*(i+1)*x[i]
168 g_ex[1,i]=2*(i+1)
169 u_ex[2]+= 3*(i+1)*x[i]
170 g_ex[2,i]=3*(i+1)
171 A[0,i,0,i]=1
172 A[1,i,1,i]=1
173 A[2,i,2,i]=1
174
175 # create PDE:
176 pde=LinearPDE(self.domain,numEquations=3)
177 pde.setSymmetryOn()
178 mask=whereZero(x[0])*[1,1,1]
179 pde.setValue(r=u_ex,q=mask)
180 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
181 # -------- get the solution ---------------------------
182 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
183 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
184 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
185 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
186 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
187 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
188 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
189
190 u=pde.getSolution()
191 # -------- test the solution ---------------------------
192 error=Lsup(u-u_ex)/Lsup(u_ex)
193 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
194
195 def test_Poisson4(self):
196 x=Solution(self.domain).getX()
197 # --- set exact solution ----
198 u_ex=Data(1.,(4,),Solution(self.domain))
199 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
200 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
201 for i in xrange(self.domain.getDim()):
202 u_ex[0]+= 1*(i+1) *x[i]
203 g_ex[0,i]=1*(i+1)
204 u_ex[1]+= 2*(i+1)*x[i]
205 g_ex[1,i]=2*(i+1)
206 u_ex[2]+= 3*(i+1)*x[i]
207 g_ex[2,i]=3*(i+1)
208 u_ex[3]+= 4*(i+1)*x[i]
209 g_ex[3,i]=4*(i+1)
210 A[0,i,0,i]=1
211 A[1,i,1,i]=1
212 A[2,i,2,i]=1
213 A[3,i,3,i]=1
214
215 # create PDE:
216 pde=LinearPDE(self.domain,numEquations=4)
217 pde.setSymmetryOn()
218 mask=whereZero(x[0])*[1,1,1,1]
219 pde.setValue(r=u_ex,q=mask)
220 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
221 # -------- get the solution ---------------------------
222 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
223 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
224 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
225 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
226 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
227 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
228 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
229
230 u=pde.getSolution()
231 # -------- test the solution ---------------------------
232 error=Lsup(u-u_ex)/Lsup(u_ex)
233 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
234
235 def test_WeakCoupled2(self):
236 x=Solution(self.domain).getX()
237 # --- set exact solution ----
238 u_ex=Data(1.,(2,),Solution(self.domain))
239 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
240 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
241 for i in xrange(self.domain.getDim()):
242 u_ex[0]+= 1*(i+1)*x[i]
243 g_ex[0,i]=1*(i+1)
244 u_ex[1]+= 2*(i+1)*x[i]
245 g_ex[1,i]=2*(i+1)
246 A[0,i,0,i]=1
247 A[1,i,1,i]=1
248
249 Y=Data(0.,(2,),Function(self.domain))
250 Y[0]= u_ex[1]
251 Y[1]=u_ex[0]
252 # create PDE:
253 pde=LinearPDE(self.domain,numEquations=2)
254 pde.setSymmetryOn()
255 mask=whereZero(x[0])*[1,1]
256 pde.setValue(r=u_ex,q=mask)
257 pde.setValue(A=A,
258 D=-0.01*(numpy.ones((2,2))-kronecker(2)),
259 Y=-0.01*Y,
260 y=matrixmult(g_ex,self.domain.getNormal()))
261 # -------- get the solution ---------------------------
262 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
263 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
264 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
265 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
266 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
267 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
268 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
269
270 u=pde.getSolution()
271 # -------- test the solution ---------------------------
272 error=Lsup(u-u_ex)/Lsup(u_ex)
273 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
274
275 def test_WeakCoupled3(self):
276 x=Solution(self.domain).getX()
277 # --- set exact solution ----
278 u_ex=Data(1.,(3,),Solution(self.domain))
279 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
280 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
281 for i in xrange(self.domain.getDim()):
282 u_ex[0]+= 1*(i+1)*x[i]
283 g_ex[0,i]=1*(i+1)
284 u_ex[1]+= 2*(i+1)*x[i]
285 g_ex[1,i]=2*(i+1)
286 u_ex[2]+= 3*(i+1)*x[i]
287 g_ex[2,i]=3*(i+1)
288 A[0,i,0,i]=1
289 A[1,i,1,i]=1
290 A[2,i,2,i]=1
291
292 Y=Data(0.,(3,),Function(self.domain))
293 Y[0]= u_ex[1]+u_ex[2]
294 Y[1]=u_ex[0]+ u_ex[2]
295 Y[2]=u_ex[0]+u_ex[1]
296 # create PDE:
297 pde=LinearPDE(self.domain,numEquations=3)
298 pde.setSymmetryOn()
299 mask=whereZero(x[0])*[1,1,1]
300 pde.setValue(r=u_ex,q=mask)
301 pde.setValue(A=A,
302 D=-0.01*(numpy.ones((3,3))-kronecker(3)),
303 Y=-0.01*Y,
304 y=matrixmult(g_ex,self.domain.getNormal()))
305 # -------- get the solution ---------------------------
306 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
307 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
308 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
309 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
310 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
311 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
312 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
313
314 u=pde.getSolution()
315 # -------- test the solution ---------------------------
316 error=Lsup(u-u_ex)/Lsup(u_ex)
317 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
318
319 def test_WeakCoupled4(self):
320 x=Solution(self.domain).getX()
321 # --- set exact solution ----
322 u_ex=Data(1.,(4,),Solution(self.domain))
323 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
324 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
325 for i in xrange(self.domain.getDim()):
326 u_ex[0]+= 1*(i+1)*x[i]
327 g_ex[0,i]=1*(i+1)
328 u_ex[1]+= 2*(i+1)*x[i]
329 g_ex[1,i]=2*(i+1)
330 u_ex[2]+= 3*(i+1)*x[i]
331 g_ex[2,i]=3*(i+1)
332 u_ex[3]+= 4*(i+1)*x[i]
333 g_ex[3,i]=4*(i+1)
334 A[0,i,0,i]=1
335 A[1,i,1,i]=1
336 A[2,i,2,i]=1
337 A[3,i,3,i]=1
338
339 Y=Data(0.,(4,),Function(self.domain))
340 Y[0]= u_ex[1]+u_ex[2]+u_ex[3]
341 Y[1]=u_ex[0]+ u_ex[2]+u_ex[3]
342 Y[2]=u_ex[0]+u_ex[1]+ u_ex[3]
343 Y[3]=u_ex[0]+u_ex[1]+u_ex[2]
344 # create PDE:
345 pde=LinearPDE(self.domain,numEquations=4)
346 pde.setSymmetryOn()
347 mask=whereZero(x[0])*[1,1,1,1]
348 pde.setValue(r=u_ex,q=mask)
349 pde.setValue(A=A,
350 D=-0.01*(numpy.ones((4,4))-kronecker(4)),
351 Y=-0.01*Y,
352 y=matrixmult(g_ex,self.domain.getNormal()))
353 # -------- get the solution ---------------------------
354 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
355 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
356 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
357 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
358 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
359 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
360 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
361
362 u=pde.getSolution()
363 # -------- test the solution ---------------------------
364 error=Lsup(u-u_ex)/Lsup(u_ex)
365 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
366
367 def test_StrongCoupled2(self):
368 x=Solution(self.domain).getX()
369 # --- set exact solution ----
370 u_ex=Data(1.,(2,),Solution(self.domain))
371 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
372 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
373 for i in xrange(self.domain.getDim()):
374 u_ex[0]+= 1*(i+1)*x[i]
375 g_ex[0,i]=1*(i+1)
376 u_ex[1]+= 2*(i+1)*x[i]
377 g_ex[1,i]=2*(i+1)
378 A[0,i,0,i]=1
379 A[1,i,1,i]=1
380
381 Y=Data(0.,(2,),Function(self.domain))
382 Y[0]= u_ex[1]
383 Y[1]=u_ex[0]
384 # create PDE:
385 pde=LinearPDE(self.domain,numEquations=2)
386 pde.setSymmetryOn()
387 mask=whereZero(x[0])*[1,1]
388 pde.setValue(r=u_ex,q=mask)
389 pde.setValue(A=A,
390 D=-20.*(numpy.ones((2,2))-kronecker(2)),
391 Y=-20.*Y,
392 y=matrixmult(g_ex,self.domain.getNormal()))
393 # -------- get the solution ---------------------------
394 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
395 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
396 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
397 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
398 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
399 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
400 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
401
402 u=pde.getSolution()
403 # -------- test the solution ---------------------------
404 error=Lsup(u-u_ex)/Lsup(u_ex)
405 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
406
407 def test_StrongCoupled3(self):
408 x=Solution(self.domain).getX()
409 # --- set exact solution ----
410 u_ex=Data(1.,(3,),Solution(self.domain))
411 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
412 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
413 for i in xrange(self.domain.getDim()):
414 u_ex[0]+= 1*(i+1)*x[i]
415 g_ex[0,i]=1*(i+1)
416 u_ex[1]+= 2*(i+1)*x[i]
417 g_ex[1,i]=2*(i+1)
418 u_ex[2]+= 3*(i+1)*x[i]
419 g_ex[2,i]=3*(i+1)
420 A[0,i,0,i]=1
421 A[1,i,1,i]=1
422 A[2,i,2,i]=1
423
424 Y=Data(0.,(3,),Function(self.domain))
425 Y[0]= u_ex[1]+u_ex[2]
426 Y[1]=u_ex[0]+ u_ex[2]
427 Y[2]=u_ex[0]+u_ex[1]
428 # create PDE:
429 pde=LinearPDE(self.domain,numEquations=3)
430 pde.setSymmetryOn()
431 mask=whereZero(x[0])*[1,1,1]
432 pde.setValue(r=u_ex,q=mask)
433 pde.setValue(A=A,
434 D=-20.*(numpy.ones((3,3))-kronecker(3)),
435 Y=-20.*Y,
436 y=matrixmult(g_ex,self.domain.getNormal()))
437 # -------- get the solution ---------------------------
438 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
439 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
440 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
441 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
442 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
443 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
444 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
445
446 u=pde.getSolution()
447 # -------- test the solution ---------------------------
448 error=Lsup(u-u_ex)/Lsup(u_ex)
449 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
450 def test_StrongCoupled4(self):
451 x=Solution(self.domain).getX()
452 # --- set exact solution ----
453 u_ex=Data(1.,(4,),Solution(self.domain))
454 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
455 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
456 for i in xrange(self.domain.getDim()):
457 u_ex[0]+= 1*(i+1)*x[i]
458 g_ex[0,i]=1*(i+1)
459 u_ex[1]+= 2*(i+1)*x[i]
460 g_ex[1,i]=2*(i+1)
461 u_ex[2]+= 3*(i+1)*x[i]
462 g_ex[2,i]=3*(i+1)
463 u_ex[3]+= 4*(i+1)*x[i]
464 g_ex[3,i]=4*(i+1)
465 A[0,i,0,i]=1
466 A[1,i,1,i]=1
467 A[2,i,2,i]=1
468 A[3,i,3,i]=1
469
470 Y=Data(0.,(4,),Function(self.domain))
471 Y[0]= u_ex[1]+u_ex[2]+u_ex[3]
472 Y[1]=u_ex[0]+ u_ex[2]+u_ex[3]
473 Y[2]=u_ex[0]+u_ex[1]+ u_ex[3]
474 Y[3]=u_ex[0]+u_ex[1]+u_ex[2]
475 # create PDE:
476 pde=LinearPDE(self.domain,numEquations=4)
477 pde.setSymmetryOn()
478 mask=whereZero(x[0])*[1,1,1,1]
479 pde.setValue(r=u_ex,q=mask)
480 pde.setValue(A=A,
481 D=-20.*(numpy.ones((4,4))-kronecker(4)),
482 Y=-20.*Y,
483 y=matrixmult(g_ex,self.domain.getNormal()))
484 # -------- get the solution ---------------------------
485 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
486 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
487 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
488 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
489 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
490 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
491 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
492
493 u=pde.getSolution()
494 # -------- test the solution ---------------------------
495 error=Lsup(u-u_ex)/Lsup(u_ex)
496 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
497
498 def test_Square(self):
499 # PDE constants
500 h1 = 0.5
501 h2 = 0.5
502 beta = 1.
503 alpha1 = 2.
504 alpha2 = 1.
505
506 # domain masks and domain-specific constants
507 x = Solution(self.domain).getX(); x0 = x[0]; x1 = x[1]
508 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
509 omega1 = 1-omega2
510 ratio = alpha1/alpha2
511 alpha = alpha1*omega1 + alpha2*omega2
512
513 # --- set exact solution ----
514 a1 = 1.
515 d1 = 1.
516 b1 = -d1*h2
517 c1 = -d1*h1
518
519 a2 = a1 - (1.-ratio)*d1*h1*h2
520 b2 = ratio*b1
521 c2 = ratio*c1
522 d2 = ratio*d1
523
524 u_ex = omega1*(a1 + b1*x0 + c1*x1 + d1*x0*x1) + \
525 omega2*(a2 + b2*x0 + c2*x1 + d2*x0*x1)
526
527 # create PDE:
528 pde = LinearPDE(self.domain,numEquations=1)
529
530 # set the value to that of the solution on the boundary
531 q = whereZero(x0) + whereZero(x1) + \
532 whereZero(sup(x0)-x0) + whereZero(sup(x1)-x1)
533 pde.setValue(q=q,r=u_ex)
534
535 # create X points in the centre of the grid elements
536 xe = Function(self.domain).getX()
537 x0 = xe[0]
538 x1 = xe[1]
539
540 # redefine omega so that apha is more precise on the diagonal (?)
541 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
542 omega1 = 1-omega2
543 ratio = alpha1/alpha2
544 alpha = alpha1*omega1 + alpha2*omega2
545
546 # set up PDE coefficients
547 pde.setValue(A=alpha*kronecker(self.domain), D=beta, Y=beta*u_ex)
548 pde.setSymmetryOn()
549
550 # -------- get the solution ---------------------------
551 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
552 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
553 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
554 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
555 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
556 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
557 if MAX_LEVEL!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
558
559 u = pde.getSolution()
560
561 # -------- test the solution ---------------------------
562 error=Lsup(u-u_ex)/Lsup(u_ex)
563 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
564
565
566 class Test_AMGOnFinleyHex2DOrder1(Test_AMG):
567 RES_TOL=5.e-7
568 SOLVER_TOL=1.e-8
569 def setUp(self):
570 NE=int(float(NE_TOTAL)**(1./2.)+0.5)
571 self.domain = Rectangle(NE,NE,1, optimize=OPTIMIZE)
572 def tearDown(self):
573 del self.domain
574
575 class Test_AMGOnFinleyHex3DOrder1(Test_AMG):
576 RES_TOL=5.e-7
577 SOLVER_TOL=1.e-8
578 def setUp(self):
579 NE=int(float(NE_TOTAL)**(1./3.)+0.5)
580 self.domain = Brick(NE,NE,NE,1, optimize=OPTIMIZE)
581 def tearDown(self):
582 del self.domain
583 if __name__ == '__main__':
584 suite = unittest.TestSuite()
585 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex2DOrder1))
586 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex3DOrder1))
587 #suite.addTest(Test_AMGOnFinleyHex2DOrder1("test_PoissonSqueezedX"))
588
589 s=unittest.TextTestRunner(verbosity=2).run(suite)
590 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26