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

  ViewVC Help
Powered by ViewVC 1.1.26