/[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 4019 - (show annotations)
Thu Oct 11 08:12:55 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 38081 byte(s)
More tabbing errors,
range/xrange
...
1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2012 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development since 2012 by School of Earth Sciences
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 """
25 Test suite for AMG
26
27 :remark:
28
29 :var __author__: name of author
30 :var __licence__: licence agreement
31 :var __url__: url entry point on documentation
32 :var __version__: version
33 :var __date__: date of the version
34 """
35
36 __author__="Lutz Gross, l.gross@uq.edu.au"
37
38 import unittest, sys
39
40 from esys.escript import *
41 from esys.finley import Rectangle,Brick
42 from esys.escript.linearPDEs import LinearPDE, SolverOptions
43 import numpy
44 OPTIMIZE=True # and False
45 SOLVER_VERBOSE=True or False
46
47 MIN_MATRIX_SIZE=1
48 MIN_SPARSITY=1.
49 MIN_MATRIX_SIZE=None
50 MIN_SPARSITY=None
51 MAX_LEVEL=None
52 USE_AMG=True or False
53
54 try:
55 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
56 except KeyError:
57 FINLEY_TEST_DATA='.'
58 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
59
60 # number of elements in the spatial directions
61 NE_TOTAL=4096
62 #NE_TOTAL=4
63
64 class Test_AMG(unittest.TestCase):
65
66 def test_Poisson(self):
67 global USE_AMG
68 x=Solution(self.domain).getX()
69 # --- set exact solution ----
70 u_ex=Scalar(1,Solution(self.domain))
71 g_ex=Vector(0.,Solution(self.domain))
72 for i in range(self.domain.getDim()):
73 u_ex+=(i+1)*x[i]
74 g_ex[i]=(i+1)
75
76 # create PDE:
77 pde=LinearPDE(self.domain,numEquations=1)
78 pde.setSymmetryOn()
79 mask=whereZero(x[0])
80 pde.setValue(r=u_ex,q=mask)
81 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
82 # -------- get the solution ---------------------------
83 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
84 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
85 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
86 print("AMG is disabled for MPI builds")
87 USE_AMG=0
88 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
89 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
90 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
91 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
92 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
93
94 u=pde.getSolution()
95 # -------- test the solution ---------------------------
96 error=Lsup(u-u_ex)/Lsup(u_ex)
97 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
98
99 def test_PoissonWithDirectInterpolation(self):
100 global USE_AMG
101 x=Solution(self.domain).getX()
102 # --- set exact solution ----
103 u_ex=Scalar(1,Solution(self.domain))
104 g_ex=Vector(0.,Solution(self.domain))
105 for i in range(self.domain.getDim()):
106 u_ex+=(i+1)*x[i]
107 g_ex[i]=(i+1)
108
109 # create PDE:
110 pde=LinearPDE(self.domain,numEquations=1)
111 pde.setSymmetryOn()
112 mask=whereZero(x[0])
113 pde.setValue(r=u_ex,q=mask)
114 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
115 # -------- get the solution ---------------------------
116 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
117 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
118 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
119 print("AMG is disabled for MPI builds")
120 USE_AMG=0
121 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
122 pde.getSolverOptions().setNumPreSweeps(3)
123 pde.getSolverOptions().setNumPostSweeps(3)
124 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
125 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
126 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
127 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
128 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().DIRECT_INTERPOLATION)
129
130 u=pde.getSolution()
131 # -------- test the solution ---------------------------
132 error=Lsup(u-u_ex)/Lsup(u_ex)
133 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
134
135 def test_PoissonClassic(self):
136 global USE_AMG
137 x=Solution(self.domain).getX()
138 # --- set exact solution ----
139 u_ex=Scalar(1,Solution(self.domain))
140 g_ex=Vector(0.,Solution(self.domain))
141 for i in range(self.domain.getDim()):
142 u_ex+=(i+1)*x[i]
143 g_ex[i]=(i+1)
144
145 # create PDE:
146 pde=LinearPDE(self.domain,numEquations=1)
147 pde.setSymmetryOn()
148 mask=whereZero(x[0])
149 pde.setValue(r=u_ex,q=mask)
150 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
151 # -------- get the solution ---------------------------
152 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
153 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
154 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
155 print("AMG is disabled for MPI builds")
156 USE_AMG=0
157 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
158 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
159 pde.getSolverOptions().setNumPreSweeps(3)
160 pde.getSolverOptions().setNumPostSweeps(3)
161 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
162 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
163 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
164 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
165
166 u=pde.getSolution()
167 # -------- test the solution ---------------------------
168 error=Lsup(u-u_ex)/Lsup(u_ex)
169 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
170
171 def test_PoissonClassicWithFFCoupling(self):
172 global USE_AMG
173 x=Solution(self.domain).getX()
174 # --- set exact solution ----
175 u_ex=Scalar(1,Solution(self.domain))
176 g_ex=Vector(0.,Solution(self.domain))
177 for i in range(self.domain.getDim()):
178 u_ex+=(i+1)*x[i]
179 g_ex[i]=(i+1)
180
181 # create PDE:
182 pde=LinearPDE(self.domain,numEquations=1)
183 pde.setSymmetryOn()
184 mask=whereZero(x[0])
185 pde.setValue(r=u_ex,q=mask)
186 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
187 # -------- get the solution ---------------------------
188 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
189 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
190 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
191 print("AMG is disabled for MPI builds")
192 USE_AMG=0
193 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
194 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
195 pde.getSolverOptions().setNumPreSweeps(3)
196 pde.getSolverOptions().setNumPostSweeps(3)
197 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
198 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
199 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
200 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION_WITH_FF_COUPLING)
201
202 u=pde.getSolution()
203 # -------- test the solution ---------------------------
204 error=Lsup(u-u_ex)/Lsup(u_ex)
205 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
206
207 def test_PoissonSqueezedX(self):
208 global USE_AMG
209 x=self.domain.getX().copy()
210 x[0]*=0.5
211 self.domain.setX(x)
212 x=Solution(self.domain).getX()
213 # --- set exact solution ----
214 u_ex=Scalar(1,Solution(self.domain))
215 g_ex=Vector(0.,Solution(self.domain))
216 for i in range(self.domain.getDim()):
217 u_ex+=(i+1)*x[i]
218 g_ex[i]=(i+1)
219
220 # create PDE:
221 pde=LinearPDE(self.domain,numEquations=1)
222 pde.setSymmetryOn()
223 mask=whereZero(x[0])
224 pde.setValue(r=u_ex,q=mask)
225 pde.setValue(A=kronecker(self.domain),y=inner(g_ex,self.domain.getNormal()))
226 # -------- get the solution ---------------------------
227 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
228 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
229 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
230 print("AMG is disabled for MPI builds")
231 USE_AMG=0
232 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
233 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
234 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
235 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
236 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
237
238 u=pde.getSolution()
239 # -------- test the solution ---------------------------
240 error=Lsup(u-u_ex)/Lsup(u_ex)
241 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
242
243
244 def test_Poisson2(self):
245 global USE_AMG
246 x=Solution(self.domain).getX()
247 # --- set exact solution ----
248 u_ex=Data(1.,(2,),Solution(self.domain))
249 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
250 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
251 for i in range(self.domain.getDim()):
252 u_ex[0]+= 1*(i+1) *x[i]
253 g_ex[0,i]=1*(i+1)
254 u_ex[1]+= 2*(i+1)*x[i]
255 g_ex[1,i]=2*(i+1)
256 A[0,i,0,i]=1
257 A[1,i,1,i]=1
258
259 # create PDE:
260 pde=LinearPDE(self.domain,numEquations=2)
261 pde.setSymmetryOn()
262 mask=whereZero(x[0])*[1,1]
263 pde.setValue(r=u_ex,q=mask)
264 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
265 # -------- get the solution ---------------------------
266 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
267 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
268 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
269 print("AMG is disabled for MPI builds")
270 USE_AMG=0
271 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
272 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
273 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
274 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
275 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
276
277 u=pde.getSolution()
278 # -------- test the solution ---------------------------
279 error=Lsup(u-u_ex)/Lsup(u_ex)
280 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
281
282 def test_Poisson3(self):
283 global USE_AMG
284 x=Solution(self.domain).getX()
285 # --- set exact solution ----
286 u_ex=Data(1.,(3,),Solution(self.domain))
287 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
288 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
289 for i in range(self.domain.getDim()):
290 u_ex[0]+= 1*(i+1) *x[i]
291 g_ex[0,i]=1*(i+1)
292 u_ex[1]+= 2*(i+1)*x[i]
293 g_ex[1,i]=2*(i+1)
294 u_ex[2]+= 3*(i+1)*x[i]
295 g_ex[2,i]=3*(i+1)
296 A[0,i,0,i]=1
297 A[1,i,1,i]=1
298 A[2,i,2,i]=1
299
300 # create PDE:
301 pde=LinearPDE(self.domain,numEquations=3)
302 pde.setSymmetryOn()
303 mask=whereZero(x[0])*[1,1,1]
304 pde.setValue(r=u_ex,q=mask)
305 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
306 # -------- get the solution ---------------------------
307 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
308 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
309 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
310 print("AMG is disabled for MPI builds")
311 USE_AMG=0
312 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
313 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
314 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
315 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
316 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
317
318 u=pde.getSolution()
319 # -------- test the solution ---------------------------
320 error=Lsup(u-u_ex)/Lsup(u_ex)
321 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
322
323 def test_Poisson4(self):
324 global USE_AMG
325 x=Solution(self.domain).getX()
326 # --- set exact solution ----
327 u_ex=Data(1.,(4,),Solution(self.domain))
328 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
329 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
330 for i in range(self.domain.getDim()):
331 u_ex[0]+= 1*(i+1) *x[i]
332 g_ex[0,i]=1*(i+1)
333 u_ex[1]+= 2*(i+1)*x[i]
334 g_ex[1,i]=2*(i+1)
335 u_ex[2]+= 3*(i+1)*x[i]
336 g_ex[2,i]=3*(i+1)
337 u_ex[3]+= 4*(i+1)*x[i]
338 g_ex[3,i]=4*(i+1)
339 A[0,i,0,i]=1
340 A[1,i,1,i]=1
341 A[2,i,2,i]=1
342 A[3,i,3,i]=1
343
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,y=matrixmult(g_ex,self.domain.getNormal()))
350 # -------- get the solution ---------------------------
351 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
352 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
353 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
354 print("AMG is disabled for MPI builds")
355 USE_AMG=0
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
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 def test_Poisson2Classic(self):
367 global USE_AMG
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 range(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 # create PDE:
382 pde=LinearPDE(self.domain,numEquations=2)
383 pde.setSymmetryOn()
384 mask=whereZero(x[0])*[1,1]
385 pde.setValue(r=u_ex,q=mask)
386 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
387 # -------- get the solution ---------------------------
388 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
389 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
390 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
391 print("AMG is disabled for MPI builds")
392 USE_AMG=0
393 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
394 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
395 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
396 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
397 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
398 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
399
400 u=pde.getSolution()
401 # -------- test the solution ---------------------------
402 error=Lsup(u-u_ex)/Lsup(u_ex)
403 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
404
405 def test_Poisson3Classic(self):
406 global USE_AMG
407 x=Solution(self.domain).getX()
408 # --- set exact solution ----
409 u_ex=Data(1.,(3,),Solution(self.domain))
410 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
411 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
412 for i in range(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 A[0,i,0,i]=1
420 A[1,i,1,i]=1
421 A[2,i,2,i]=1
422
423 # create PDE:
424 pde=LinearPDE(self.domain,numEquations=3)
425 pde.setSymmetryOn()
426 mask=whereZero(x[0])*[1,1,1]
427 pde.setValue(r=u_ex,q=mask)
428 pde.setValue(A=A,y=matrixmult(g_ex,self.domain.getNormal()))
429 # -------- get the solution ---------------------------
430 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
431 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
432 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
433 print("AMG is disabled for MPI builds")
434 USE_AMG=0
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_Poisson4Classic(self):
448 global USE_AMG
449 x=Solution(self.domain).getX()
450 # --- set exact solution ----
451 u_ex=Data(1.,(4,),Solution(self.domain))
452 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
453 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
454 for i in range(self.domain.getDim()):
455 u_ex[0]+= 1*(i+1) *x[i]
456 g_ex[0,i]=1*(i+1)
457 u_ex[1]+= 2*(i+1)*x[i]
458 g_ex[1,i]=2*(i+1)
459 u_ex[2]+= 3*(i+1)*x[i]
460 g_ex[2,i]=3*(i+1)
461 u_ex[3]+= 4*(i+1)*x[i]
462 g_ex[3,i]=4*(i+1)
463 A[0,i,0,i]=1
464 A[1,i,1,i]=1
465 A[2,i,2,i]=1
466 A[3,i,3,i]=1
467
468 # create PDE:
469 pde=LinearPDE(self.domain,numEquations=4)
470 pde.setSymmetryOn()
471 mask=whereZero(x[0])*[1,1,1,1]
472 pde.setValue(r=u_ex,q=mask)
473 pde.setValue(A=A,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 and getEscriptParamInt('DISABLE_AMG',0):
478 print("AMG is disabled for MPI builds")
479 USE_AMG=0
480 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
481 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
482 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
483 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
484 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
485 pde.getSolverOptions().setAMGInterpolation(pde.getSolverOptions().CLASSIC_INTERPOLATION)
486
487 u=pde.getSolution()
488 # -------- test the solution ---------------------------
489 error=Lsup(u-u_ex)/Lsup(u_ex)
490 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
491
492 def test_WeakCoupled2(self):
493 global USE_AMG
494 x=Solution(self.domain).getX()
495 # --- set exact solution ----
496 u_ex=Data(1.,(2,),Solution(self.domain))
497 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
498 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
499 for i in range(self.domain.getDim()):
500 u_ex[0]+= 1*(i+1)*x[i]
501 g_ex[0,i]=1*(i+1)
502 u_ex[1]+= 2*(i+1)*x[i]
503 g_ex[1,i]=2*(i+1)
504 A[0,i,0,i]=1
505 A[1,i,1,i]=1
506
507 Y=Data(0.,(2,),Function(self.domain))
508 a=-1./2.*0.01
509 Y[0]= u_ex[0] + a * u_ex[1]
510 Y[1]=a * u_ex[0] + u_ex[1]
511 # create PDE:
512 pde=LinearPDE(self.domain,numEquations=2)
513 pde.setSymmetryOn()
514 mask=whereZero(x[0])*[1,1]
515 pde.setValue(r=u_ex,q=mask)
516 pde.setValue(A=A,
517 D=kronecker(2)+a*(numpy.ones((2,2))-kronecker(2)),
518 Y=Y,
519 y=matrixmult(g_ex,self.domain.getNormal()))
520 # -------- get the solution ---------------------------
521 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
522 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
523 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
524 print("AMG is disabled for MPI builds")
525 USE_AMG=0
526 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
527 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
528 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
529 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
530 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
531
532 u=pde.getSolution()
533 # -------- test the solution ---------------------------
534 error=Lsup(u-u_ex)/Lsup(u_ex)
535 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
536
537 def test_WeakCoupled3(self):
538 global USE_AMG
539 x=Solution(self.domain).getX()
540 # --- set exact solution ----
541 u_ex=Data(1.,(3,),Solution(self.domain))
542 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
543 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
544 for i in range(self.domain.getDim()):
545 u_ex[0]+= 1*(i+1)*x[i]
546 g_ex[0,i]=1*(i+1)
547 u_ex[1]+= 2*(i+1)*x[i]
548 g_ex[1,i]=2*(i+1)
549 u_ex[2]+= 3*(i+1)*x[i]
550 g_ex[2,i]=3*(i+1)
551 A[0,i,0,i]=1
552 A[1,i,1,i]=1
553 A[2,i,2,i]=1
554
555 a=-1./3.*0.01
556 Y=Data(0.,(3,),Function(self.domain))
557 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]
558 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]
559 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]
560 # create PDE:
561 pde=LinearPDE(self.domain,numEquations=3)
562 pde.setSymmetryOn()
563 mask=whereZero(x[0])*[1,1,1]
564 pde.setValue(r=u_ex,q=mask)
565 pde.setValue(A=A,
566 D=kronecker(3)+a*(numpy.ones((3,3))-kronecker(3)),
567 Y=Y,
568 y=matrixmult(g_ex,self.domain.getNormal()))
569 # -------- get the solution ---------------------------
570 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
571 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
572 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
573 print("AMG is disabled for MPI builds")
574 USE_AMG=0
575 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
576 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
577 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
578 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
579 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
580 u=pde.getSolution()
581 # -------- test the solution ---------------------------
582 error=Lsup(u-u_ex)/Lsup(u_ex)
583 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
584
585 def test_WeakCoupled4(self):
586 global USE_AMG
587 x=Solution(self.domain).getX()
588 # --- set exact solution ----
589 u_ex=Data(1.,(4,),Solution(self.domain))
590 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
591 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
592 for i in range(self.domain.getDim()):
593 u_ex[0]+= 1*(i+1)*x[i]
594 g_ex[0,i]=1*(i+1)
595 u_ex[1]+= 2*(i+1)*x[i]
596 g_ex[1,i]=2*(i+1)
597 u_ex[2]+= 3*(i+1)*x[i]
598 g_ex[2,i]=3*(i+1)
599 u_ex[3]+= 4*(i+1)*x[i]
600 g_ex[3,i]=4*(i+1)
601 A[0,i,0,i]=1
602 A[1,i,1,i]=1
603 A[2,i,2,i]=1
604 A[3,i,3,i]=1
605
606 a=-1./4.*0.01
607
608 Y=Data(0.,(4,),Function(self.domain))
609 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
610 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
611 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]+ a * u_ex[3]
612 Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ u_ex[3]
613 # create PDE:
614 pde=LinearPDE(self.domain,numEquations=4)
615 pde.setSymmetryOn()
616 mask=whereZero(x[0])*[1,1,1,1]
617 pde.setValue(r=u_ex,q=mask)
618 pde.setValue(A=A,
619 D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
620 Y=Y,
621 y=matrixmult(g_ex,self.domain.getNormal()))
622 # -------- get the solution ---------------------------
623 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
624 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
625 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
626 print("AMG is disabled for MPI builds")
627 USE_AMG=0
628 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
629 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
630 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
631 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
632 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
633 u=pde.getSolution()
634 # -------- test the solution ---------------------------
635 error=Lsup(u-u_ex)/Lsup(u_ex)
636 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
637
638 def test_StrongCoupled2(self):
639 global USE_AMG
640 x=Solution(self.domain).getX()
641 # --- set exact solution ----
642 u_ex=Data(1.,(2,),Solution(self.domain))
643 g_ex=Data(0.,(2,self.domain.getDim()), Solution(self.domain))
644 A=Data(0.,(2,self.domain.getDim(),2,self.domain.getDim()), Function(self.domain))
645 for i in range(self.domain.getDim()):
646 u_ex[0]+= 1*(i+1)*x[i]
647 g_ex[0,i]=1*(i+1)
648 u_ex[1]+= 2*(i+1)*x[i]
649 g_ex[1,i]=2*(i+1)
650 A[0,i,0,i]=1
651 A[1,i,1,i]=1
652
653 Y=Data(0.,(2,),Function(self.domain))
654 a=-1./2.*0.99
655 Y[0]= u_ex[0] + a * u_ex[1]
656 Y[1]=a * u_ex[0] + u_ex[1]
657 # create PDE:
658 pde=LinearPDE(self.domain,numEquations=2)
659 pde.setSymmetryOn()
660 mask=whereZero(x[0])*[1,1]
661 pde.setValue(r=u_ex,q=mask)
662 pde.setValue(A=A,
663 D=kronecker(2)+a*(numpy.ones((2,2))-kronecker(2)),
664 Y=Y,
665 y=matrixmult(g_ex,self.domain.getNormal()))
666 # -------- get the solution ---------------------------
667 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
668 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
669 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
670 print("AMG is disabled for MPI builds")
671 USE_AMG=0
672 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
673 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
674 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
675 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
676 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
677
678 u=pde.getSolution()
679 # -------- test the solution ---------------------------
680 error=Lsup(u-u_ex)/Lsup(u_ex)
681 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
682
683 def test_StrongCoupled3(self):
684 global USE_AMG
685 x=Solution(self.domain).getX()
686 # --- set exact solution ----
687 u_ex=Data(1.,(3,),Solution(self.domain))
688 g_ex=Data(0.,(3,self.domain.getDim()), Solution(self.domain))
689 A=Data(0.,(3,self.domain.getDim(),3,self.domain.getDim()), Function(self.domain))
690 for i in range(self.domain.getDim()):
691 u_ex[0]+= 1*(i+1)*x[i]
692 g_ex[0,i]=1*(i+1)
693 u_ex[1]+= 2*(i+1)*x[i]
694 g_ex[1,i]=2*(i+1)
695 u_ex[2]+= 3*(i+1)*x[i]
696 g_ex[2,i]=3*(i+1)
697 A[0,i,0,i]=1
698 A[1,i,1,i]=1
699 A[2,i,2,i]=1
700
701 a=-1./3.*0.99
702 Y=Data(0.,(3,),Function(self.domain))
703 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]
704 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]
705 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]
706 # create PDE:
707 pde=LinearPDE(self.domain,numEquations=3)
708 pde.setSymmetryOn()
709 mask=whereZero(x[0])*[1,1,1]
710 pde.setValue(r=u_ex,q=mask)
711 pde.setValue(A=A,
712 D=kronecker(3)+a*(numpy.ones((3,3))-kronecker(3)),
713 Y=Y,
714 y=matrixmult(g_ex,self.domain.getNormal()))
715 # -------- get the solution ---------------------------
716 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
717 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
718 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
719 print("AMG is disabled for MPI builds")
720 USE_AMG=0
721 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
722 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
723 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
724 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
725 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
726 u=pde.getSolution()
727 # -------- test the solution ---------------------------
728 error=Lsup(u-u_ex)/Lsup(u_ex)
729 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
730
731 def test_StrongCoupled4(self):
732 global USE_AMG
733 x=Solution(self.domain).getX()
734 # --- set exact solution ----
735 u_ex=Data(1.,(4,),Solution(self.domain))
736 g_ex=Data(0.,(4,self.domain.getDim()), Solution(self.domain))
737 A=Data(0.,(4,self.domain.getDim(),4,self.domain.getDim()), Function(self.domain))
738 for i in range(self.domain.getDim()):
739 u_ex[0]+= 1*(i+1)*x[i]
740 g_ex[0,i]=1*(i+1)
741 u_ex[1]+= 2*(i+1)*x[i]
742 g_ex[1,i]=2*(i+1)
743 u_ex[2]+= 3*(i+1)*x[i]
744 g_ex[2,i]=3*(i+1)
745 u_ex[3]+= 4*(i+1)*x[i]
746 g_ex[3,i]=4*(i+1)
747 A[0,i,0,i]=1
748 A[1,i,1,i]=1
749 A[2,i,2,i]=1
750 A[3,i,3,i]=1
751
752 a=-1./4.*0.99
753
754 Y=Data(0.,(4,),Function(self.domain))
755 Y[0]= u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
756 Y[1]= a * u_ex[0]+ u_ex[1]+ a * u_ex[2]+ a * u_ex[3]
757 Y[2]= a * u_ex[0]+ a * u_ex[1]+ u_ex[2]+ a * u_ex[3]
758 Y[3]= a * u_ex[0]+ a * u_ex[1]+ a * u_ex[2]+ u_ex[3]
759 # create PDE:
760 pde=LinearPDE(self.domain,numEquations=4)
761 pde.setSymmetryOn()
762 mask=whereZero(x[0])*[1,1,1,1]
763 pde.setValue(r=u_ex,q=mask)
764 pde.setValue(A=A,
765 D=kronecker(4)+a*(numpy.ones((4,4))-kronecker(4)),
766 Y=Y,
767 y=matrixmult(g_ex,self.domain.getNormal()))
768 # -------- get the solution ---------------------------
769 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
770 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
771 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
772 print("AMG is disabled for MPI builds")
773 USE_AMG=0
774 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
775 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
776 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
777 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
778 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
779 u=pde.getSolution()
780 # -------- test the solution ---------------------------
781 error=Lsup(u-u_ex)/Lsup(u_ex)
782 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
783
784
785 def test_Square(self):
786 global USE_AMG
787 # PDE constants
788 h1 = 0.5
789 h2 = 0.5
790 beta = 1.
791 alpha1 = 2.
792 alpha2 = 1.
793
794 # domain masks and domain-specific constants
795 x = Solution(self.domain).getX(); x0 = x[0]; x1 = x[1]
796 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
797 omega1 = 1-omega2
798 ratio = alpha1/alpha2
799 alpha = alpha1*omega1 + alpha2*omega2
800
801 # --- set exact solution ----
802 a1 = 1.
803 d1 = 1.
804 b1 = -d1*h2
805 c1 = -d1*h1
806
807 a2 = a1 - (1.-ratio)*d1*h1*h2
808 b2 = ratio*b1
809 c2 = ratio*c1
810 d2 = ratio*d1
811
812 u_ex = omega1*(a1 + b1*x0 + c1*x1 + d1*x0*x1) + \
813 omega2*(a2 + b2*x0 + c2*x1 + d2*x0*x1)
814
815 # create PDE:
816 pde = LinearPDE(self.domain,numEquations=1)
817
818 # set the value to that of the solution on the boundary
819 q = whereZero(x0) + whereZero(x1) + \
820 whereZero(sup(x0)-x0) + whereZero(sup(x1)-x1)
821 pde.setValue(q=q,r=u_ex)
822
823 # create X points in the centre of the grid elements
824 xe = Function(self.domain).getX()
825 x0 = xe[0]
826 x1 = xe[1]
827
828 # redefine omega so that apha is more precise on the diagonal (?)
829 omega2 = wherePositive(x0-h1)*wherePositive(x1-h2)
830 omega1 = 1-omega2
831 ratio = alpha1/alpha2
832 alpha = alpha1*omega1 + alpha2*omega2
833
834 # set up PDE coefficients
835 pde.setValue(A=alpha*kronecker(self.domain), D=beta, Y=beta*u_ex)
836 pde.setSymmetryOn()
837
838 # -------- get the solution ---------------------------
839 pde.getSolverOptions().setTolerance(self.SOLVER_TOL)
840 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
841 if USE_AMG and getEscriptParamInt('DISABLE_AMG',0):
842 print("AMG is disabled for MPI builds")
843 USE_AMG=0
844 if (USE_AMG): pde.getSolverOptions().setPreconditioner(SolverOptions.AMG)
845 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
846 if MIN_MATRIX_SIZE!= None: pde.getSolverOptions().setMinCoarseMatrixSize(MIN_MATRIX_SIZE)
847 if MIN_SPARSITY!=None: pde.getSolverOptions().setMinCoarseMatrixSparsity(MIN_SPARSITY)
848 if MAX_LEVEL!=None: pde.getSolverOptions().setLevelMax(MAX_LEVEL)
849 u = pde.getSolution()
850
851 # -------- test the solution ---------------------------
852 error=Lsup(u-u_ex)/Lsup(u_ex)
853 self.assertTrue(error<self.RES_TOL, "solution error %s is too big."%error)
854
855
856 class Test_AMGOnFinleyHex2DOrder1(Test_AMG):
857 RES_TOL=5.e-7
858 SOLVER_TOL=1.e-8
859 def setUp(self):
860 NE=int(float(NE_TOTAL)**(1./2.)+0.5)
861 self.domain = Rectangle(NE,NE,1, optimize=OPTIMIZE)
862 def tearDown(self):
863 del self.domain
864
865 class Test_AMGOnFinleyHex3DOrder1(Test_AMG):
866 RES_TOL=5.e-7
867 SOLVER_TOL=1.e-8
868 def setUp(self):
869 NE=int(float(NE_TOTAL)**(1./3.)+0.5)
870 self.domain = Brick(NE,NE,NE,1, optimize=OPTIMIZE)
871 def tearDown(self):
872 del self.domain
873 if __name__ == '__main__':
874 suite = unittest.TestSuite()
875 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex2DOrder1))
876 suite.addTest(unittest.makeSuite(Test_AMGOnFinleyHex3DOrder1))
877 # suite.addTest(Test_AMGOnFinleyHex3DOrder1("test_Poisson4"))
878 # suite.addTest(Test_AMGOnFinleyHex2DOrder1("test_WeakCoupled4"))
879
880 s=unittest.TextTestRunner(verbosity=2).run(suite)
881 if not s.wasSuccessful(): sys.exit(12)

  ViewVC Help
Powered by ViewVC 1.1.26