/[escript]/trunk/speckley/test/python/run_specialOnSpeckley.py
ViewVC logotype

Contents of /trunk/speckley/test/python/run_specialOnSpeckley.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 31235 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2015 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 from __future__ import division, print_function
18
19 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 import os
27 import sys
28 import esys.escriptcore.utestselect as unittest
29 from esys.escriptcore.testing import *
30 from esys.escript import *
31 from esys.speckley import Rectangle, Brick, speckleycpp
32
33 class Test_Speckley_Assemblers(unittest.TestCase):
34 TOLERANCE = 1e-10
35
36 def test_Brick_XY_single(self):
37 ranks = getMPISizeWorld()
38 for expanded in (True, False):
39 for order in range(2,11):
40 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
41 Y = Data(3, Function(dom), expanded)
42 X = Data(0, (3,), Function(dom), expanded)
43 X[0] = dom.getX()[0]
44 X[1] = dom.getX()[1]
45 X[2] = dom.getX()[2]
46 f = Data(0, Solution(dom), expanded)
47
48 dom.addToRHS(f, [("Y",Y)],
49 dom.createAssembler("DefaultAssembler", []))
50 dom.addToRHS(f, [("X", X)],
51 dom.createAssembler("DefaultAssembler", []))
52 #nuke the boundary back to zero since it's not dealt with here
53 for dim in range(3):
54 f *= whereNegative(dom.getX()[dim]-6)
55 res = Lsup(f)
56 self.assertLess(res, self.TOLERANCE,
57 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
58 res, self.TOLERANCE)).format("" if expanded else "un-"))
59
60 def test_Rectangle_XY_single(self):
61 ranks = getMPISizeWorld()
62 for expanded in (True, False):
63 for order in range(2,11):
64 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=6, l1=6)
65 Y = Data(2, Function(dom), expanded)
66 X = Data(0, (2,), Function(dom), expanded)
67 X[0] = dom.getX()[0]
68 X[1] = dom.getX()[1]
69 f = Data(0, Solution(dom), expanded)
70
71 dom.addToRHS(f, [("Y",Y)],
72 dom.createAssembler("DefaultAssembler", []))
73 dom.addToRHS(f, [("X", X)],
74 dom.createAssembler("DefaultAssembler", []))
75 #nuke the boundary back to zero since it's not dealt with here
76 for dim in range(2):
77 f *= whereNegative(dom.getX()[dim]-6)
78 res = Lsup(f)
79 self.assertLess(res, self.TOLERANCE,
80 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
81 res, self.TOLERANCE)).format("" if expanded else "un-"))
82
83 def test_Brick_XY_system(self):
84 ranks = getMPISizeWorld()
85 for expanded in (True, False):
86 for order in range(2,11):
87 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
88 Y = Data(1, (3,), Function(dom), expanded)
89 X = Data(0, (3,3), Function(dom), expanded)
90 X[0,0] = dom.getX()[0]
91 X[1,1] = dom.getX()[1]
92 X[2,2] = dom.getX()[2]
93
94 f = Data(0, (3,), Solution(dom), expanded)
95
96 dom.addToRHS(f, [("Y",Y)],
97 dom.createAssembler("DefaultAssembler", []))
98 dom.addToRHS(f, [("X", X)],
99 dom.createAssembler("DefaultAssembler", []))
100 #nuke the boundary back to zero since it's not dealt with here
101 for dim in range(3):
102 f *= whereNegative(dom.getX()[dim]-6)
103 res = Lsup(f)
104 self.assertLess(res, self.TOLERANCE,
105 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
106 res, self.TOLERANCE)).format("" if expanded else "un-"))
107
108 def test_Rectangle_XY_system(self):
109 ranks = getMPISizeWorld()
110 for expanded in (True, False):
111 for order in range(2,11):
112 dom = Rectangle(order, 3, 3*ranks, l0=6, l1=6, d1=ranks)
113 Y = Data(1, (2,), Function(dom), expanded)
114 X = Data(0, (2,2), Function(dom), expanded)
115 X[0,0] = dom.getX()[0]
116 X[1,1] = dom.getX()[1]
117
118 f = Data(0, (2,), Solution(dom), expanded)
119
120 dom.addToRHS(f, [("Y",Y)],
121 dom.createAssembler("DefaultAssembler", []))
122 dom.addToRHS(f, [("X", X)],
123 dom.createAssembler("DefaultAssembler", []))
124 #nuke the boundary back to zero since it's not dealt with here
125 f *= whereNegative(dom.getX()[0]-6)*whereNegative(dom.getX()[1]-6)
126 res = Lsup(f)
127 self.assertLess(res, self.TOLERANCE,
128 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
129 res, self.TOLERANCE)).format("" if expanded else "un-"))
130
131 def test_Brick_Du_Y_single(self):
132 #solves for u in Du = Y, where D = 1, Y = 2
133 ranks = getMPISizeWorld()
134 for expanded in (True, False):
135 for order in range(2,11):
136 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
137 D = Data(1, ContinuousFunction(dom), expanded)
138 Y = Data(2, ContinuousFunction(dom), expanded)
139
140 Y = interpolate(Y, Function(dom))
141 D = interpolate(D, Function(dom))
142
143 rhs = Data(0, Solution(dom), expanded)
144 lhs = Data(0, Solution(dom), expanded)
145
146 dom.addToRHS(rhs, [("Y",Y)],
147 dom.createAssembler("DefaultAssembler", []))
148 dom.addToRHS(lhs, [("D",D)],
149 dom.createAssembler("DefaultAssembler", []))
150
151 res = Lsup((rhs/lhs)-2)
152 self.assertLess(res, self.TOLERANCE,
153 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
154 res, self.TOLERANCE)).format("" if expanded else "un-"))
155
156 def test_Rectangle_Du_Y_single(self):
157 #solves for u in Du = Y, where D = 1, Y = 2
158 ranks = getMPISizeWorld()
159 for expanded in (True, False):
160 for order in range(2,11):
161 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
162 D = Data(1, ContinuousFunction(dom), expanded)
163 Y = Data(2, ContinuousFunction(dom), expanded)
164
165 Y = interpolate(Y, Function(dom))
166 D = interpolate(D, Function(dom))
167
168 rhs = Data(0, Solution(dom), expanded)
169 lhs = Data(0, Solution(dom), expanded)
170
171 dom.addToRHS(rhs, [("Y",Y)],
172 dom.createAssembler("DefaultAssembler", []))
173 dom.addToRHS(lhs, [("D",D)],
174 dom.createAssembler("DefaultAssembler", []))
175
176 res = Lsup((rhs/lhs)-2)
177 self.assertLess(res, self.TOLERANCE,
178 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
179 res, self.TOLERANCE)).format("" if expanded else "un-"))
180
181 def test_Brick_Du_Y_system(self):
182 #solves for u in Du = Y, where D = [1,2], Y = [2,4]
183 ranks = getMPISizeWorld()
184 for expanded in (True, False):
185 for order in range(2,11):
186 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
187 D = Data(1, (2,), ContinuousFunction(dom), expanded)
188 D[1] = 2
189 Y = Data(2, (2,), ContinuousFunction(dom), expanded)
190 Y[1] = 4
191
192 Y = interpolate(Y, Function(dom))
193 D = interpolate(D, Function(dom))
194
195 rhs = Data(0, (2,), Solution(dom), expanded)
196 lhs = Data(0, (2,), Solution(dom), expanded)
197
198 dom.addToRHS(rhs, [("Y",Y)],
199 dom.createAssembler("DefaultAssembler", []))
200 dom.addToRHS(lhs, [("D",D)],
201 dom.createAssembler("DefaultAssembler", []))
202
203 res = Lsup((rhs/lhs)-2)
204 self.assertLess(res, self.TOLERANCE,
205 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
206 res, self.TOLERANCE)).format("" if expanded else "un-"))
207
208 def test_Rectangle_Du_Y_system(self):
209 #solves for u in Du = Y, where D = [1,2], Y = [2,4]
210 ranks = getMPISizeWorld()
211 for expanded in (True, False):
212 for order in range(2,11):
213 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
214 D = Data(1, (2,), ContinuousFunction(dom), expanded)
215 D[1] = 2
216 Y = Data(2, (2,), ContinuousFunction(dom), expanded)
217 Y[1] = 4
218
219 Y = interpolate(Y, Function(dom))
220 D = interpolate(D, Function(dom))
221
222 rhs = Data(0, (2,), Solution(dom), expanded)
223 lhs = Data(0, (2,), Solution(dom), expanded)
224
225 dom.addToRHS(rhs, [("Y",Y)],
226 dom.createAssembler("DefaultAssembler", []))
227 dom.addToRHS(lhs, [("D",D)],
228 dom.createAssembler("DefaultAssembler", []))
229
230 res = Lsup((rhs/lhs)-2)
231 self.assertLess(res, self.TOLERANCE,
232 ("assembly for {0}expanded order %d failed with %g >= %g"%(order,
233 res, self.TOLERANCE)).format("" if expanded else "un-"))
234
235 class Test_Speckley(unittest.TestCase):
236 TOLERANCE = 1e-10
237 def test_Rectangle_ReducedFunction(self):
238 ranks = getMPISizeWorld()
239 for order in range(2, 11):
240 dom = Rectangle(order, 3, 3*ranks, l0=3, l1=3*ranks, d1=ranks)
241 X = dom.getX()
242 redData = interpolate(X, ReducedFunction(dom))
243 data = [(interpolate(redData, ReducedFunction(dom)), "ReducedFunction"),
244 (interpolate(redData, Function(dom)), "Function"),
245 (interpolate(redData, ContinuousFunction(dom)), "ContinuousFunction")]
246 for d, fs in data:
247 self.assertLess(inf(d-[0.5]*2), self.TOLERANCE,
248 "reduced->%s failure with order %d: %g != 0"%(fs, order, inf(d-[0.5]*2)))
249 self.assertLess(sup(d[0]+0.5) - 3, self.TOLERANCE,
250 "reduced->%s failure with order %d: %g != 3"%(fs, order, sup(d[0]+0.5)))
251 self.assertLess(sup(d[1]+0.5) - 3*ranks, self.TOLERANCE,
252 "reduced->%s failure with order %d: %g != %g"%(fs, order, sup(d[1]+0.5), 3*ranks))
253
254 def test_Brick_ReducedFunction(self):
255 ranks = getMPISizeWorld()
256 for order in range(2, 11):
257 dom = Brick(order, 3, 3*ranks, 3, l0=3, l1=3*ranks, l2=3, d1=ranks)
258 X = dom.getX()
259 redData = interpolate(X, ReducedFunction(dom))
260 data = [(interpolate(redData, ReducedFunction(dom)), "ReducedFunction"),
261 (interpolate(redData, Function(dom)), "Function"),
262 (interpolate(redData, ContinuousFunction(dom)), "ContinuousFunction")]
263 for d, fs in data:
264 self.assertLess(inf(d-[0.5]*3), self.TOLERANCE,
265 "reduced->%s failure with order %d: %g != 0"%(fs, order, inf(d-[0.5]*3)))
266 self.assertLess(sup(d[0]+0.5) - 3, self.TOLERANCE,
267 "reduced->%s failure with order %d: %g != 3"%(fs, order, sup(d[0]+0.5)))
268 self.assertLess(sup(d[1]+0.5) - 3*ranks, self.TOLERANCE,
269 "reduced->%s failure with order %d: %g != %g"%(fs, order, sup(d[1]+0.5), 3*ranks))
270 self.assertLess(sup(d[2]+0.5) - 3, self.TOLERANCE,
271 "reduced->%s failure with order %d: %g != 3"%(fs, order, sup(d[2]+0.5)))
272
273 def test_Rectangle_Function_gradient(self): #expanded and non-expanded
274 ranks = getMPISizeWorld()
275 for expanded in [True, False]:
276 for order in range(2,11):
277 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
278 x = Data(5, Function(dom), True)
279 self.assertLess(Lsup(grad(x)), 1e-10,
280 "single component failure, order %d%s, %g >= 1e-10"%(order,
281 (" expanded" if expanded else ""), Lsup(grad(x))))
282 for data in [[5,1], [-5,-1], [5,1,1e-5]]:
283 x = Data(data, Function(dom), True)
284 g = grad(x)
285 for n,d in enumerate(data):
286 self.assertLess(Lsup(g[n]), 1e-10,
287 "%d-component failure, order %d %sexpanded, %g >= 1e-10"%(len(data),
288 order, ("" if expanded else "un-"), Lsup(g[n])))
289
290 def test_Brick_Function_gradient(self):
291 ranks = getMPISizeWorld()
292 for expanded in [True, False]:
293 for order in range(2,11):
294 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
295 x = Data(5, Function(dom), True)
296 self.assertLess(Lsup(grad(x)), 1e-10,
297 "single component failure, order %d%s, %g >= 1e-10"%(order,
298 (" expanded" if expanded else ""), Lsup(grad(x))))
299 for data in [[5,1], [-5,-1], [5,1,1e-5]]:
300 x = Data(data, Function(dom), True)
301 g = grad(x)
302 for n,d in enumerate(data):
303 self.assertLess(Lsup(g[n]), 1e-10,
304 "%d-component failure, order %d %sexpanded, %g >= 1e-10"%(len(data),
305 order, ("" if expanded else "un-"), Lsup(g[n])))
306
307
308 def test_Rectangle_ContinuousFunction_gradient(self):
309 ranks = getMPISizeWorld()
310 for order in range(2,11):
311 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=100, l1=100)
312 X = dom.getX()
313 u = X[0] + X[1] + 1
314 v = Lsup(grad(u) - 1)
315 self.assertLess(v, 1e-10, "order %d, %g >= 1e-10, %s"%(order, v, str(grad(u)-1)))
316 for power in range(1, order+1):
317 for power2 in range(1, order+1):
318 a = X[0]**power * X[1]**power2
319 da = grad(a)
320 first = Lsup(da[0] - power*X[0]**(power-1) * X[1]**power2) \
321 /Lsup(power*X[0]**(power-1) * X[1]**power2)
322 second = Lsup(da[1] - power2*X[1]**(power2-1) * X[0]**power) \
323 /Lsup(power2*X[1]**(power2-1) * X[0]**power)
324 self.assertLess(first, 1e-10,
325 "order %d and degree %d,%d, %g >= 1e-9"%(order,
326 power, power2, first))
327 self.assertLess(second, 1e-10,
328 "order %d and degree %d,%d, %g >= 1e-9"%(order,
329 power, power2, second))
330
331 def test_Brick_ContinuousFunction_gradient(self):
332 ranks = getMPISizeWorld()
333 for order in range(2,11):
334 dom = Brick(order, 3, 3*ranks, 3, l0=100, l1=100, l2=100, d1=ranks)
335 X = dom.getX()
336 u = X[0] + X[1] + X[2] + 1
337 v = Lsup(grad(u) - 1)
338 self.assertLess(v, 1e-10, "order %d, %g >= 1e-10, %s"%(order, v,
339 str(grad(u)-1)))
340 for power1 in range(1, order+1, order//2):
341 for power2 in range(1, order+1, order//2):
342 for power3 in range(1, order+1, order//2):
343 a = X[0]**power1 * X[1]**power2 * X[2]**power3
344 da = grad(a)
345
346 temp = power1*X[0]**(power1-1) * X[1]**power2*X[2]**power3
347 first = Lsup(da[0] - temp) / Lsup(temp)
348 temp = power2*X[1]**(power2-1) * X[0]**power1*X[2]**power3
349 second = Lsup(da[1] - temp) / Lsup(temp)
350 temp = power3*X[2]**(power3-1) * X[0]**power1*X[1]**power2
351 third = Lsup(da[2] - temp) / Lsup(temp)
352
353 self.assertLess(first, 1e-10,
354 "order %d and degree %d,%d,%d, %g >= 1e-9"%(order,
355 power1, power2, power3, first))
356 self.assertLess(second, 1e-10,
357 "order %d and degree %d,%d,%d, %g >= 1e-9"%(order,
358 power1, power2, power3, second))
359 self.assertLess(third, 1e-10,
360 "order %d and degree %d,%d,%d, %g >= 1e-9"%(order,
361 power1, power2, power3, third))
362
363 def test_Rectangle_interpolation_continuous_noncontinuous_and_back(self):
364 ranks = getMPISizeWorld()
365 for order in range(2,11):
366 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=6, l1=6)
367 original = Data(5, Function(dom), True)
368 cont = interpolate(original, ContinuousFunction(dom))
369 func = interpolate(cont, Function(dom))
370 self.assertEqual(Lsup(original-func), 0,
371 "interpolation of constant, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
372 x = dom.getX()
373 original = x[0] + x[1] + 1
374 cont = interpolate(original, ContinuousFunction(dom))
375 func = interpolate(cont, Function(dom))
376 self.assertEqual(Lsup(original-func), 0,
377 "interpolation of expanded, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
378 original = whereZero(x[0]-2) + whereZero(x[1]-2)
379 cont = interpolate(original, ContinuousFunction(dom))
380 func = interpolate(cont, Function(dom))
381 self.assertEqual(Lsup(original-func), 0,
382 "interpolation of point, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
383
384 def test_Brick_interpolation_continuous_noncontinuous_and_back(self):
385 ranks = getMPISizeWorld()
386 for order in range(2,11):
387 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=ranks, l2=6, d1=ranks)
388 original = Data(5, Function(dom), True)
389 cont = interpolate(original, ContinuousFunction(dom))
390 func = interpolate(cont, Function(dom))
391 self.assertEqual(Lsup(original-func), 0,
392 "interpolation of constant, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
393 x = dom.getX()
394 original = x[0] + x[1] + x[2] + 1
395 cont = interpolate(original, ContinuousFunction(dom))
396 func = interpolate(cont, Function(dom))
397 self.assertEqual(Lsup(original-func), 0,
398 "interpolation of expanded, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
399 original = whereZero(x[0]-2) + whereZero(x[1]-2) + whereZero(x[2] - 2)
400 cont = interpolate(original, ContinuousFunction(dom))
401 func = interpolate(cont, Function(dom))
402 self.assertEqual(Lsup(original-func), 0,
403 "interpolation of point, order %d: original and final not equal, %g != 0"%(order, Lsup(original-func)))
404
405 def test_Rectangle_integration(self):
406 ranks = getMPISizeWorld()
407 for order in range(2,11):
408 size = 6
409 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=size, l1=size)
410 X = dom.getX()
411 for k in range(1, order * 2):
412 for l in range(1, order * 2):
413 powered = X[0]**k * X[1]**l
414 integral = integrate(powered)
415 actual = size**(k+1) * size**(l+1) / ((k+1.)*(l+1.))
416 self.assertLess(abs(integral - actual)/actual, 1e-11,
417 "too much variance in integral result (order %d, degrees %d %d)"%(order, k, l))
418
419 def test_Brick_integration(self):
420 ranks = getMPISizeWorld()
421 for order in range(2,11):
422 size = 6
423 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
424 X = dom.getX()
425 for k in [1, order, order*2 - 1]:
426 for l in [1, order, order*2 - 1]:
427 for m in [1, order, order*2 - 1]:
428 powered = X[0]**k * X[1]**l * X[2]**m
429 integral = integrate(powered)
430 actual = size**(k+1) * size**(l+1) * size**(m+1) \
431 /((k+1.)*(l+1.)*(m+1.))
432 res = abs(integral - actual)/actual
433 self.assertLess(res, 1e-11,
434 "too much variance in integral result (order %d, degrees %d %d, %g >= 1e-11)"%(order,
435 k, l, res))
436
437 @unittest.skipIf(getMPISizeWorld() == 1,
438 "only works with more than one rank")
439 def test_Rectangle_MPI_construction_single_dimensional(self):
440 ranks = getMPISizeWorld()
441 for order in range(2,11):
442 dom = Rectangle(order, 2, 2*ranks, l1=ranks, d1=ranks)
443 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]), ranks+1,
444 "invalid getX() for y-splits order %d"%order)
445
446 filt = whereZero(dom.getX()[1] - 1)
447 for i in range(2,ranks):
448 filt += whereZero(dom.getX()[1] - i)
449 if getMPIRankWorld() % 2:
450 filt *= -1
451 d = Vector(0, Function(dom))
452 d[0] = 1 * filt
453 d[1] = 10 * filt
454 X = interpolate(d, Function(dom))
455 res = interpolate(X, ContinuousFunction(dom))
456 val = Lsup(res)
457 self.assertEqual(val, 0,
458 "summation stage failure for y-splits in order %d"%order)
459
460 dom = Rectangle(2, 2*ranks, 2, l0=ranks, d0=ranks)
461 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]), ranks+1,
462 "invalid getX() for x-splits order %d"%order)
463 filt = whereZero(dom.getX()[0] - 1)
464 for i in range(2,getMPISizeWorld()):
465 filt += whereZero(dom.getX()[0] - i)
466 if getMPIRankWorld() % 2:
467 filt *= -1
468 d = Vector(0, Function(dom))
469 d[0] = 1 * filt
470 d[1] = 10 * filt
471 X = interpolate(d, Function(dom))
472 res = interpolate(X, ContinuousFunction(dom))
473 val = Lsup(res)
474 self.assertEqual(val, 0,
475 "summation stage failure for x-splits in order %d"%order)
476
477 X = interpolate(d+2, Function(dom))
478 res = interpolate(X, ContinuousFunction(dom))
479 val = Lsup(res-2)
480 self.assertEqual(val, 0,
481 "averaging stage failure for x-splits in order %d"%order)
482
483 @unittest.skipIf(getMPISizeWorld() != 4, "requires 4 ranks")
484 def test_Rectangle_MPI_construction_multi_dimensional(self):
485 ranks = getMPISizeWorld()
486 half = 2 #change if ranks != 4 (sqrt(ranks))
487 for order in range(2, 11):
488 dom = Rectangle(order, ranks, ranks, l0=half, l1=half,
489 d0=half, d1=half)
490 self.assertEqual(Lsup(dom.getX()[0]), half,
491 "invalid getX() for multidimensional splits in order %d"%order)
492 self.assertEqual(Lsup(dom.getX()[1]), half,
493 "invalid getX() for multidimensional splits in order %d"%order)
494 xfilt = whereZero(dom.getX()[0] - 1) + whereZero(dom.getX()[1] - 1)
495 for i in range(2,half):
496 xfilt += whereZero(dom.getX()[0] - i)
497 xfilt += whereZero(dom.getX()[1] - i)
498 xfilt = whereNonZero(xfilt)
499 if getMPIRankWorld() in [1,2]: #change if ranks != 4
500 xfilt *= -1
501 X = interpolate(xfilt, Function(dom))
502 res = interpolate(X, ContinuousFunction(dom))
503 val = Lsup(res)
504 self.assertEqual(val, 0,
505 "summation failure for mixed-splits in order %d"%order)
506 X = interpolate(xfilt+2, Function(dom))
507 res = interpolate(X, ContinuousFunction(dom))
508 val = Lsup(res-2)
509 self.assertEqual(val, 0,
510 "averaging failure for mixed-splits in order %d"%order)
511
512 @unittest.skipIf(getMPISizeWorld() == 1,
513 "only works with more than one rank")
514 def test_Brick_MPI_construction(self):
515 for order in range(2,11):
516 dom = Brick(order, 2*getMPISizeWorld(), 2, 2, d0=getMPISizeWorld())
517 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
518 3, "invalid Lsup(getX()) for x-split order %d"%order)
519 X = interpolate(dom.getX()[0], Function(dom)) \
520 + interpolate(dom.getX()[1], Function(dom))
521 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
522 - (dom.getX()[0] + dom.getX()[1]))
523 self.assertLess(val, 1e-10,
524 "interpolation failure for x-split order %d"%order)
525
526 dom = Brick(order, 2, 2*getMPISizeWorld(), 2, d1=getMPISizeWorld())
527 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
528 3, "invalid Lsup(getX()) for y-split order %d"%order)
529 X = interpolate(dom.getX()[0], Function(dom)) \
530 + interpolate(dom.getX()[1], Function(dom))
531 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
532 - (dom.getX()[0] + dom.getX()[1]))
533 self.assertLess(val, 1e-10,
534 "interpolation failure for y-split order %d"%order)
535
536 dom = Brick(order, 2, 2, 2*getMPISizeWorld(), d2=getMPISizeWorld())
537 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
538 3, "invalid Lsup(getX()) for z-split order %d"%order)
539 X = interpolate(dom.getX()[0], Function(dom)) \
540 + interpolate(dom.getX()[1], Function(dom))
541 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
542 - (dom.getX()[0] + dom.getX()[1]))
543 self.assertLess(val, 1e-10,
544 "interpolation failure for z-split order %d"%order)
545
546 @unittest.skipIf(getMPISizeWorld() == 1, "requires multiple MPI processes")
547 def test_Brick_singledim_subdivision(self):
548 ranks = getMPISizeWorld()
549 for dim in range(0,3):
550 label = ["x","y","z"][dim]
551 size = [2,2,2]
552 size[dim] *= ranks
553 lengths = [1,1,1]
554 lengths[dim] *= ranks
555 splits = [1,1,1]
556 splits[dim] *= ranks
557
558 for order in range(2, 11):
559 dom = Brick(order, size[0], size[1], size[2],
560 l0=lengths[0], l1=lengths[1], l2=lengths[2],
561 d0=splits[0], d1=splits[1], d2=splits[2])
562 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]+dom.getX()[2]),
563 ranks+2, "invalid getX() for %s-splits order %d"%(\
564 label, order))
565
566 filt = whereZero(dom.getX()[dim] - 1)
567 for i in range(2,ranks):
568 filt += whereZero(dom.getX()[0] - i)
569 if getMPIRankWorld() % 2:
570 filt *= -1
571 d = Vector(0, Function(dom))
572 d[0] = 1 * filt
573 d[1] = 10 * filt
574 d[2] = 100 * filt
575 X = interpolate(d, Function(dom))
576 res = interpolate(X, ContinuousFunction(dom))
577 val = Lsup(res)
578 self.assertEqual(val, 0,
579 "summation stage failure for %s-splits in order %d,"%(\
580 label, order))
581 X = interpolate(d+2, Function(dom))
582 res = interpolate(X, ContinuousFunction(dom))
583 val = Lsup(res-2)
584 self.assertEqual(val, 0,
585 "averaging stage failure for %s-splits in order %d,"%(\
586 label, order))
587
588 @unittest.skipIf(getMPISizeWorld() != 4, "requires 4 ranks exactly")
589 def test_Brick_multidim_subdivision(self):
590 ranks = getMPISizeWorld()
591 half = 2 #change if ranks != 4 (sqrt(ranks))
592 for order in range(2, 11):
593 for dom,dim1,dim2 in [
594 (Brick(order, 2*half, 2*half, 2, l0=half, l1=half,
595 d0=half, d1=half), 0,1),
596 (Brick(order, 2*half, 2, 2*half, l0=half, l2=half,
597 d0=half, d2=half), 0,2),
598 (Brick(order, 2, 2*half, 2*half, l1=half, l2=half,
599 d1=half, d2=half), 1,2)]:
600 self.assertEqual(ranks + 1,
601 Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
602 "invalid getX() for multidimensional split " + \
603 "(dims %d,%d) in order %d"%(dim1,dim2,order))
604 xfilt = whereZero(dom.getX()[dim1] - 1) \
605 + whereZero(dom.getX()[dim2] - 1)
606 for i in range(2,half):
607 xfilt += whereZero(dom.getX()[dim1] - i)
608 xfilt += whereZero(dom.getX()[dim2] - i)
609 xfilt = whereNonZero(xfilt)
610 if getMPIRankWorld() in [1,2]: #change if ranks != 4
611 xfilt *= -1
612 X = interpolate(xfilt, Function(dom))
613 res = interpolate(X, ContinuousFunction(dom))
614 val = Lsup(res)
615 self.assertEqual(val, 0,
616 "summation failure for mixed-splits " \
617 + "(dims %d,%d) in order %d"%(dim1,dim2,order))
618 X = interpolate(xfilt+2, Function(dom))
619 res = interpolate(X, ContinuousFunction(dom))
620 val = Lsup(res-2)
621 self.assertEqual(val, 0,
622 "averaging failure for mixed-splits "\
623 + "(dims %d,%d) in order %d"%(dim1,dim2,order))
624
625 if __name__ == '__main__':
626 run_tests(__name__, exit_on_failure=True)
627

  ViewVC Help
Powered by ViewVC 1.1.26