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

  ViewVC Help
Powered by ViewVC 1.1.26