/[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 3449 - (show annotations)
Fri Jan 21 05:48:15 2011 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 35189 byte(s)
compile problem in AMG fixed.
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 # and False
44 SOLVER_VERBOSE=True and False
45
46 MIN_MATRIX_SIZE=1
47 MIN_SPARSITY=1.
48 MIN_MATRIX_SIZE=None
49 MIN_SPARSITY=None
50 MAX_LEVEL=None
51 USE_AMG=True or False
52
53 try:
54 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
55 except KeyError:
56 FINLEY_TEST_DATA='.'
57 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
58
59 # number of elements in the spatial directions
60 NE_TOTAL=4096
61 #NE_TOTAL=4
62
63 class Test_AMG(unittest.TestCase):
64
65 def test_Poisson(self):
66 x=Solution(self.domain).getX()
67 # --- set exact solution ----
68 u_ex=Scalar(1,Solution(self.domain))
69 g_ex=Vector(0.,Solution(self.domain))
70 for i in xrange(self.domain.getDim()):
71 u_ex+=(i+1)*x[i]
72 g_ex[i]=(i+1)
73
74 # create PDE:
75 pde=LinearPDE(self.domain,numEquations=1)
76 pde.setSymmetryOn()
77 mask=whereZero(x[0])
78 pde.setValue(r=u_ex,q=mask)
79 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
80 # -------- get the solution ---------------------------
81 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
82 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
83 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
84 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
85 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
86 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
87 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
88
89 u=pde.getSolution()
90 # -------- test the solution ---------------------------
91 error=Lsup(u-u_ex)/Lsup(u_ex)
92 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
93
94 def test_PoissonWithDirectInterpolation(self):
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().setNumPreSweeps(3)
114 pde.getSolverOptions().setNumPostSweeps(3)
115 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
116 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
117 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
118 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
119 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().DIRECT_INTERPOLATION)
120
121 u=pde.getSolution()
122 # -------- test the solution ---------------------------
123 error=Lsup(u-u_ex)/Lsup(u_ex)
124 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
125
126 def test_PoissonClassic(self):
127 x=Solution(self.domain).getX()
128 # --- set exact solution ----
129 u_ex=Scalar(1,Solution(self.domain))
130 g_ex=Vector(0.,Solution(self.domain))
131 for i in xrange(self.domain.getDim()):
132 u_ex+=(i+1)*x[i]
133 g_ex[i]=(i+1)
134
135 # create PDE:
136 pde=LinearPDE(self.domain,numEquations=1)
137 pde.setSymmetryOn()
138 mask=whereZero(x[0])
139 pde.setValue(r=u_ex,q=mask)
140 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
141 # -------- get the solution ---------------------------
142 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
143 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
144 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
145 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
146 pde.getSolverOptions().setNumPreSweeps(3)
147 pde.getSolverOptions().setNumPostSweeps(3)
148 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
149 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
150 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
151 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
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_PoissonClassicWithFFCoupling(self):
159 x=Solution(self.domain).getX()
160 # --- set exact solution ----
161 u_ex=Scalar(1,Solution(self.domain))
162 g_ex=Vector(0.,Solution(self.domain))
163 for i in xrange(self.domain.getDim()):
164 u_ex+=(i+1)*x[i]
165 g_ex[i]=(i+1)
166
167 # create PDE:
168 pde=LinearPDE(self.domain,numEquations=1)
169 pde.setSymmetryOn()
170 mask=whereZero(x[0])
171 pde.setValue(r=u_ex,q=mask)
172 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
173 # -------- get the solution ---------------------------
174 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
175 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
176 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
177 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
178 pde.getSolverOptions().setNumPreSweeps(3)
179 pde.getSolverOptions().setNumPostSweeps(3)
180 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
181 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
182 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
183 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
184
185 u=pde.getSolution()
186 # -------- test the solution ---------------------------
187 error=Lsup(u-u_ex)/Lsup(u_ex)
188 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
189
190 def test_PoissonSqueezedX(self):
191 x=self.domain.getX().copy()
192 x[0]*=0.5
193 self.domain.setX(x)
194 x=Solution(self.domain).getX()
195 # --- set exact solution ----
196 u_ex=Scalar(1,Solution(self.domain))
197 g_ex=Vector(0.,Solution(self.domain))
198 for i in xrange(self.domain.getDim()):
199 u_ex+=(i+1)*x[i]
200 g_ex[i]=(i+1)
201
202 # create PDE:
203 pde=LinearPDE(self.domain,numEquations=1)
204 pde.setSymmetryOn()
205 mask=whereZero(x[0])
206 pde.setValue(r=u_ex,q=mask)
207 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
208 # -------- get the solution ---------------------------
209 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
210 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
211 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
212 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
213 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
214 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
215 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
216
217 u=pde.getSolution()
218 # -------- test the solution ---------------------------
219 error=Lsup(u-u_ex)/Lsup(u_ex)
220 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
221
222
223 def test_Poisson2(self):
224 x=Solution(self.domain).getX()
225 # --- set exact solution ----
226 u_ex=Data(1.,(2,),Solution(self.domain))
227 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
228 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
229 for i in xrange(self.domain.getDim()):
230 u_ex[0]+= 1*(i+1) *x[i]
231 g_ex[0,i]=1*(i+1)
232 u_ex[1]+= 2*(i+1)*x[i]
233 g_ex[1,i]=2*(i+1)
234 A[0,i,0,i]=1
235 A[1,i,1,i]=1
236
237 # create PDE:
238 pde=LinearPDE(self.domain,numEquations=2)
239 pde.setSymmetryOn()
240 mask=whereZero(x[0])*[1,1]
241 pde.setValue(r=u_ex,q=mask)
242 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
243 # -------- get the solution ---------------------------
244 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
245 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
246 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
247 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
248 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
249 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
250 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
251
252 u=pde.getSolution()
253 # -------- test the solution ---------------------------
254 error=Lsup(u-u_ex)/Lsup(u_ex)
255 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
256
257 def test_Poisson3(self):
258 x=Solution(self.domain).getX()
259 # --- set exact solution ----
260 u_ex=Data(1.,(3,),Solution(self.domain))
261 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
262 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
263 for i in xrange(self.domain.getDim()):
264 u_ex[0]+= 1*(i+1) *x[i]
265 g_ex[0,i]=1*(i+1)
266 u_ex[1]+= 2*(i+1)*x[i]
267 g_ex[1,i]=2*(i+1)
268 u_ex[2]+= 3*(i+1)*x[i]
269 g_ex[2,i]=3*(i+1)
270 A[0,i,0,i]=1
271 A[1,i,1,i]=1
272 A[2,i,2,i]=1
273
274 # create PDE:
275 pde=LinearPDE(self.domain,numEquations=3)
276 pde.setSymmetryOn()
277 mask=whereZero(x[0])*[1,1,1]
278 pde.setValue(r=u_ex,q=mask)
279 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
280 # -------- get the solution ---------------------------
281 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
282 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
283 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
284 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
285 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
286 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
287 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
288
289 u=pde.getSolution()
290 # -------- test the solution ---------------------------
291 error=Lsup(u-u_ex)/Lsup(u_ex)
292 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
293
294 def test_Poisson4(self):
295 x=Solution(self.domain).getX()
296 # --- set exact solution ----
297 u_ex=Data(1.,(4,),Solution(self.domain))
298 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
299 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
300 for i in xrange(self.domain.getDim()):
301 u_ex[0]+= 1*(i+1) *x[i]
302 g_ex[0,i]=1*(i+1)
303 u_ex[1]+= 2*(i+1)*x[i]
304 g_ex[1,i]=2*(i+1)
305 u_ex[2]+= 3*(i+1)*x[i]
306 g_ex[2,i]=3*(i+1)
307 u_ex[3]+= 4*(i+1)*x[i]
308 g_ex[3,i]=4*(i+1)
309 A[0,i,0,i]=1
310 A[1,i,1,i]=1
311 A[2,i,2,i]=1
312 A[3,i,3,i]=1
313
314 # create PDE:
315 pde=LinearPDE(self.domain,numEquations=4)
316 pde.setSymmetryOn()
317 mask=whereZero(x[0])*[1,1,1,1]
318 pde.setValue(r=u_ex,q=mask)
319 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
320 # -------- get the solution ---------------------------
321 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
322 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
323 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
324 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
325 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
326 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
327 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
328
329 u=pde.getSolution()
330 # -------- test the solution ---------------------------
331 error=Lsup(u-u_ex)/Lsup(u_ex)
332 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
333 def test_Poisson2Classic(self):
334 x=Solution(self.domain).getX()
335 # --- set exact solution ----
336 u_ex=Data(1.,(2,),Solution(self.domain))
337 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
338 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
339 for i in xrange(self.domain.getDim()):
340 u_ex[0]+= 1*(i+1) *x[i]
341 g_ex[0,i]=1*(i+1)
342 u_ex[1]+= 2*(i+1)*x[i]
343 g_ex[1,i]=2*(i+1)
344 A[0,i,0,i]=1
345 A[1,i,1,i]=1
346
347 # create PDE:
348 pde=LinearPDE(self.domain,numEquations=2)
349 pde.setSymmetryOn()
350 mask=whereZero(x[0])*[1,1]
351 pde.setValue(r=u_ex,q=mask)
352 pde.setValue(A=A,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().setLevelMax(MAX_LEVEL)
361 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
362
363 u=pde.getSolution()
364 # -------- test the solution ---------------------------
365 error=Lsup(u-u_ex)/Lsup(u_ex)
366 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
367
368 def test_Poisson3Classic(self):
369 x=Solution(self.domain).getX()
370 # --- set exact solution ----
371 u_ex=Data(1.,(3,),Solution(self.domain))
372 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
373 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
374 for i in xrange(self.domain.getDim()):
375 u_ex[0]+= 1*(i+1) *x[i]
376 g_ex[0,i]=1*(i+1)
377 u_ex[1]+= 2*(i+1)*x[i]
378 g_ex[1,i]=2*(i+1)
379 u_ex[2]+= 3*(i+1)*x[i]
380 g_ex[2,i]=3*(i+1)
381 A[0,i,0,i]=1
382 A[1,i,1,i]=1
383 A[2,i,2,i]=1
384
385 # create PDE:
386 pde=LinearPDE(self.domain,numEquations=3)
387 pde.setSymmetryOn()
388 mask=whereZero(x[0])*[1,1,1]
389 pde.setValue(r=u_ex,q=mask)
390 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
391 # -------- get the solution ---------------------------
392 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
393 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
394 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
395 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
396 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
397 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
398 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
399 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
400
401 u=pde.getSolution()
402 # -------- test the solution ---------------------------
403 error=Lsup(u-u_ex)/Lsup(u_ex)
404 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
405
406 def test_Poisson4Classic(self):
407 x=Solution(self.domain).getX()
408 # --- set exact solution ----
409 u_ex=Data(1.,(4,),Solution(self.domain))
410 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
411 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
412 for i in xrange(self.domain.getDim()):
413 u_ex[0]+= 1*(i+1) *x[i]
414 g_ex[0,i]=1*(i+1)
415 u_ex[1]+= 2*(i+1)*x[i]
416 g_ex[1,i]=2*(i+1)
417 u_ex[2]+= 3*(i+1)*x[i]
418 g_ex[2,i]=3*(i+1)
419 u_ex[3]+= 4*(i+1)*x[i]
420 g_ex[3,i]=4*(i+1)
421 A[0,i,0,i]=1
422 A[1,i,1,i]=1
423 A[2,i,2,i]=1
424 A[3,i,3,i]=1
425
426 # create PDE:
427 pde=LinearPDE(self.domain,numEquations=4)
428 pde.setSymmetryOn()
429 mask=whereZero(x[0])*[1,1,1,1]
430 pde.setValue(r=u_ex,q=mask)
431 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
432 # -------- get the solution ---------------------------
433 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
434 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
435 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
436 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
437 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
438 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
439 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
440 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
441
442 u=pde.getSolution()
443 # -------- test the solution ---------------------------
444 error=Lsup(u-u_ex)/Lsup(u_ex)
445 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
446
447 def test_WeakCoupled2(self):
448 x=Solution(self.domain).getX()
449 # --- set exact solution ----
450 u_ex=Data(1.,(2,),Solution(self.domain))
451 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
452 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
453 for i in xrange(self.domain.getDim()):
454 u_ex[0]+= 1*(i+1)*x[i]
455 g_ex[0,i]=1*(i+1)
456 u_ex[1]+= 2*(i+1)*x[i]
457 g_ex[1,i]=2*(i+1)
458 A[0,i,0,i]=1
459 A[1,i,1,i]=1
460
461 Y=Data(0.,(2,),Function(self.domain))
462 a=-1./2.*0.01
463 Y[0]= u_ex[0] + a * u_ex[1]
464 Y[1]=a * u_ex[0] + u_ex[1]
465 # create PDE:
466 pde=LinearPDE(self.domain,numEquations=2)
467 pde.setSymmetryOn()
468 mask=whereZero(x[0])*[1,1]
469 pde.setValue(r=u_ex,q=mask)
470 pde.setValue(A=A,
471 D=kronecker(2)+a*(numpy.ones((2,2))-kronecker(2)),
472 Y=Y,
473 y=matrixmult(g_ex,self.domain.getNormal()))
474 # -------- get the solution ---------------------------
475 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
476 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
477 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
478 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
479 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
480 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
481 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
482
483 u=pde.getSolution()
484 # -------- test the solution ---------------------------
485 error=Lsup(u-u_ex)/Lsup(u_ex)
486 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
487
488 def test_WeakCoupled3(self):
489 x=Solution(self.domain).getX()
490 # --- set exact solution ----
491 u_ex=Data(1.,(3,),Solution(self.domain))
492 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
493 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
494 for i in xrange(self.domain.getDim()):
495 u_ex[0]+= 1*(i+1)*x[i]
496 g_ex[0,i]=1*(i+1)
497 u_ex[1]+= 2*(i+1)*x[i]
498 g_ex[1,i]=2*(i+1)
499 u_ex[2]+= 3*(i+1)*x[i]
500 g_ex[2,i]=3*(i+1)
501 A[0,i,0,i]=1
502 A[1,i,1,i]=1
503 A[2,i,2,i]=1
504
505 a=-1./3.*0.01
506 Y=Data(0.,(3,),Function(self.domain))
507 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]
508 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]
509 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]
510 # create PDE:
511 pde=LinearPDE(self.domain,numEquations=3)
512 pde.setSymmetryOn()
513 mask=whereZero(x[0])*[1,1,1]
514 pde.setValue(r=u_ex,q=mask)
515 pde.setValue(A=A,
516 D=kronecker(3)+a*(numpy.ones((3,3))-kronecker(3)),
517 Y=Y,
518 y=matrixmult(g_ex,self.domain.getNormal()))
519 # -------- get the solution ---------------------------
520 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
521 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
522 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
523 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
524 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
525 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
526 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
527 u=pde.getSolution()
528 # -------- test the solution ---------------------------
529 error=Lsup(u-u_ex)/Lsup(u_ex)
530 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
531
532 def test_WeakCoupled4(self):
533 x=Solution(self.domain).getX()
534 # --- set exact solution ----
535 u_ex=Data(1.,(4,),Solution(self.domain))
536 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
537 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
538 for i in xrange(self.domain.getDim()):
539 u_ex[0]+= 1*(i+1)*x[i]
540 g_ex[0,i]=1*(i+1)
541 u_ex[1]+= 2*(i+1)*x[i]
542 g_ex[1,i]=2*(i+1)
543 u_ex[2]+= 3*(i+1)*x[i]
544 g_ex[2,i]=3*(i+1)
545 u_ex[3]+= 4*(i+1)*x[i]
546 g_ex[3,i]=4*(i+1)
547 A[0,i,0,i]=1
548 A[1,i,1,i]=1
549 A[2,i,2,i]=1
550 A[3,i,3,i]=1
551
552 a=-1./4.*0.01
553
554 Y=Data(0.,(4,),Function(self.domain))
555 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
556 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
557 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]+ a * u_ex[3]
558 Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ u_ex[3]
559 # create PDE:
560 pde=LinearPDE(self.domain,numEquations=4)
561 pde.setSymmetryOn()
562 mask=whereZero(x[0])*[1,1,1,1]
563 pde.setValue(r=u_ex,q=mask)
564 pde.setValue(A=A,
565 D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
566 Y=Y,
567 y=matrixmult(g_ex,self.domain.getNormal()))
568 # -------- get the solution ---------------------------
569 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
570 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
571 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
572 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
573 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
574 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
575 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
576 u=pde.getSolution()
577 # -------- test the solution ---------------------------
578 error=Lsup(u-u_ex)/Lsup(u_ex)
579 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
580
581 def test_StrongCoupled2(self):
582 x=Solution(self.domain).getX()
583 # --- set exact solution ----
584 u_ex=Data(1.,(2,),Solution(self.domain))
585 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
586 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
587 for i in xrange(self.domain.getDim()):
588 u_ex[0]+= 1*(i+1)*x[i]
589 g_ex[0,i]=1*(i+1)
590 u_ex[1]+= 2*(i+1)*x[i]
591 g_ex[1,i]=2*(i+1)
592 A[0,i,0,i]=1
593 A[1,i,1,i]=1
594
595 Y=Data(0.,(2,),Function(self.domain))
596 a=-1./2.*0.99
597 Y[0]= u_ex[0] + a * u_ex[1]
598 Y[1]=a * u_ex[0] + u_ex[1]
599 # create PDE:
600 pde=LinearPDE(self.domain,numEquations=2)
601 pde.setSymmetryOn()
602 mask=whereZero(x[0])*[1,1]
603 pde.setValue(r=u_ex,q=mask)
604 pde.setValue(A=A,
605 D=kronecker(2)+a*(numpy.ones((2,2))-kronecker(2)),
606 Y=Y,
607 y=matrixmult(g_ex,self.domain.getNormal()))
608 # -------- get the solution ---------------------------
609 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
610 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
611 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
612 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
613 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
614 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
615 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
616
617 u=pde.getSolution()
618 # -------- test the solution ---------------------------
619 error=Lsup(u-u_ex)/Lsup(u_ex)
620 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
621
622 def test_StrongCoupled3(self):
623 x=Solution(self.domain).getX()
624 # --- set exact solution ----
625 u_ex=Data(1.,(3,),Solution(self.domain))
626 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
627 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
628 for i in xrange(self.domain.getDim()):
629 u_ex[0]+= 1*(i+1)*x[i]
630 g_ex[0,i]=1*(i+1)
631 u_ex[1]+= 2*(i+1)*x[i]
632 g_ex[1,i]=2*(i+1)
633 u_ex[2]+= 3*(i+1)*x[i]
634 g_ex[2,i]=3*(i+1)
635 A[0,i,0,i]=1
636 A[1,i,1,i]=1
637 A[2,i,2,i]=1
638
639 a=-1./3.*0.99
640 Y=Data(0.,(3,),Function(self.domain))
641 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]
642 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]
643 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]
644 # create PDE:
645 pde=LinearPDE(self.domain,numEquations=3)
646 pde.setSymmetryOn()
647 mask=whereZero(x[0])*[1,1,1]
648 pde.setValue(r=u_ex,q=mask)
649 pde.setValue(A=A,
650 D=kronecker(3)+a*(numpy.ones((3,3))-kronecker(3)),
651 Y=Y,
652 y=matrixmult(g_ex,self.domain.getNormal()))
653 # -------- get the solution ---------------------------
654 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
655 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
656 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
657 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
658 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
659 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
660 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
661 u=pde.getSolution()
662 # -------- test the solution ---------------------------
663 error=Lsup(u-u_ex)/Lsup(u_ex)
664 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
665
666 def test_StrongCoupled4(self):
667 x=Solution(self.domain).getX()
668 # --- set exact solution ----
669 u_ex=Data(1.,(4,),Solution(self.domain))
670 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
671 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
672 for i in xrange(self.domain.getDim()):
673 u_ex[0]+= 1*(i+1)*x[i]
674 g_ex[0,i]=1*(i+1)
675 u_ex[1]+= 2*(i+1)*x[i]
676 g_ex[1,i]=2*(i+1)
677 u_ex[2]+= 3*(i+1)*x[i]
678 g_ex[2,i]=3*(i+1)
679 u_ex[3]+= 4*(i+1)*x[i]
680 g_ex[3,i]=4*(i+1)
681 A[0,i,0,i]=1
682 A[1,i,1,i]=1
683 A[2,i,2,i]=1
684 A[3,i,3,i]=1
685
686 a=-1./4.*0.99
687
688 Y=Data(0.,(4,),Function(self.domain))
689 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
690 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
691 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]+ a * u_ex[3]
692 Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ u_ex[3]
693 # create PDE:
694 pde=LinearPDE(self.domain,numEquations=4)
695 pde.setSymmetryOn()
696 mask=whereZero(x[0])*[1,1,1,1]
697 pde.setValue(r=u_ex,q=mask)
698 pde.setValue(A=A,
699 D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
700 Y=Y,
701 y=matrixmult(g_ex,self.domain.getNormal()))
702 # -------- get the solution ---------------------------
703 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
704 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
705 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
706 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
707 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
708 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
709 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
710 u=pde.getSolution()
711 # -------- test the solution ---------------------------
712 error=Lsup(u-u_ex)/Lsup(u_ex)
713 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
714
715
716 def test_Square(self):
717 # PDE constants
718 h1 = 0.5
719 h2 = 0.5
720 beta = 1.
721 alpha1 = 2.
722 alpha2 = 1.
723
724 # domain masks and domain-specific constants
725 x = Solution(self.domain).getX(); x0 = x[0]; x1 = x[1]
726 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
727 omega1 = 1-omega2
728 ratio = alpha1/alpha2
729 alpha = alpha1*omega1 + alpha2*omega2
730
731 # --- set exact solution ----
732 a1 = 1.
733 d1 = 1.
734 b1 = -d1*h2
735 c1 = -d1*h1
736
737 a2 = a1 - (1.-ratio)*d1*h1*h2
738 b2 = ratio*b1
739 c2 = ratio*c1
740 d2 = ratio*d1
741
742 u_ex = omega1*(a1 + b1*x0 + c1*x1 + d1*x0*x1) + \
743 omega2*(a2 + b2*x0 + c2*x1 + d2*x0*x1)
744
745 # create PDE:
746 pde = LinearPDE(self.domain,numEquations=1)
747
748 # set the value to that of the solution on the boundary
749 q = whereZero(x0) + whereZero(x1) + \
750 whereZero(sup(x0)-x0) + whereZero(sup(x1)-x1)
751 pde.setValue(q=q,r=u_ex)
752
753 # create X points in the centre of the grid elements
754 xe = Function(self.domain).getX()
755 x0 = xe[0]
756 x1 = xe[1]
757
758 # redefine omega so that apha is more precise on the diagonal (?)
759 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
760 omega1 = 1-omega2
761 ratio = alpha1/alpha2
762 alpha = alpha1*omega1 + alpha2*omega2
763
764 # set up PDE coefficients
765 pde.setValue(A=alpha*kronecker(self.domain), D=beta, Y=beta*u_ex)
766 pde.setSymmetryOn()
767
768 # -------- get the solution ---------------------------
769 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
770 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
771 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
772 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
773 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
774 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
775 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
776 u = pde.getSolution()
777
778 # -------- test the solution ---------------------------
779 error=Lsup(u-u_ex)/Lsup(u_ex)
780 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
781
782
783 class Test_AMGOnFinleyHex2DOrder1(Test_AMG):
784 RES_TOL=5.e-7
785 SOLVER_TOL=1.e-8
786 def setUp(self):
787 NE=int(float(NE_TOTAL)**(1./2.)+0.5)
788 self.domain = Rectangle(NE,NE,1, optimize=OPTIMIZE)
789 def tearDown(self):
790 del self.domain
791
792 class Test_AMGOnFinleyHex3DOrder1(Test_AMG):
793 RES_TOL=5.e-7
794 SOLVER_TOL=1.e-8
795 def setUp(self):
796 NE=int(float(NE_TOTAL)**(1./3.)+0.5)
797 self.domain = Brick(NE,NE,NE,1, optimize=OPTIMIZE)
798 def tearDown(self):
799 del self.domain
800 if __name__ == '__main__':
801 suite = unittest.TestSuite()
802 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex2DOrder1))
803 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex3DOrder1))
804 # suite.addTest(Test_AMGOnFinleyHex3DOrder1("test_Poisson4"))
805 # suite.addTest(Test_AMGOnFinleyHex2DOrder1("test_WeakCoupled4"))
806
807 s=unittest.TextTestRunner(verbosity=2).run(suite)
808 if not s.wasSuccessful(): sys.exit(12)

  ViewVC Help
Powered by ViewVC 1.1.26