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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5478 - (hide annotations)
Wed Feb 18 01:57:49 2015 UTC (4 years, 8 months ago) by sshaw
File MIME type: text/x-python
File size: 31249 byte(s)
introducing shiny new reduced function spaces to speckley
1 sshaw 5123
2     ##############################################################################
3     #
4 jfenwick 5448 # Copyright (c) 2003-2015 by University of Queensland
5 sshaw 5123 # 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 jfenwick 5448 __copyright__="""Copyright (c) 2003-2015 by University of Queensland
21 sshaw 5123 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 sshaw 5139 ranks = getMPISizeWorld()
39 sshaw 5123 for expanded in (True, False):
40     for order in range(2,11):
41 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
42 sshaw 5123 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 sshaw 5132 ranks = getMPISizeWorld()
63 sshaw 5123 for expanded in (True, False):
64     for order in range(2,11):
65 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=6, l1=6)
66 sshaw 5123 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 sshaw 5139 ranks = getMPISizeWorld()
86 sshaw 5123 for expanded in (True, False):
87     for order in range(2,11):
88 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
89 sshaw 5123 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 sshaw 5132 ranks = getMPISizeWorld()
111 sshaw 5123 for expanded in (True, False):
112     for order in range(2,11):
113 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, l0=6, l1=6, d1=ranks)
114 sshaw 5123 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 sshaw 5139 ranks = getMPISizeWorld()
135 sshaw 5123 for expanded in (True, False):
136     for order in range(2,11):
137 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
138 sshaw 5123 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 sshaw 5132 dom.addToRHS(rhs, [("Y",Y)],
148     dom.createAssembler("DefaultAssembler", []))
149     dom.addToRHS(lhs, [("D",D)],
150     dom.createAssembler("DefaultAssembler", []))
151 sshaw 5123
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 sshaw 5132 ranks = getMPISizeWorld()
160 sshaw 5123 for expanded in (True, False):
161     for order in range(2,11):
162 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
163 sshaw 5123 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 sshaw 5132 dom.addToRHS(rhs, [("Y",Y)],
173     dom.createAssembler("DefaultAssembler", []))
174     dom.addToRHS(lhs, [("D",D)],
175     dom.createAssembler("DefaultAssembler", []))
176 sshaw 5123
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 sshaw 5139 ranks = getMPISizeWorld()
185 sshaw 5123 for expanded in (True, False):
186     for order in range(2,11):
187 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
188 sshaw 5123 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 sshaw 5132 dom.addToRHS(rhs, [("Y",Y)],
200     dom.createAssembler("DefaultAssembler", []))
201     dom.addToRHS(lhs, [("D",D)],
202     dom.createAssembler("DefaultAssembler", []))
203 sshaw 5123
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 sshaw 5132 ranks = getMPISizeWorld()
212 sshaw 5123 for expanded in (True, False):
213     for order in range(2,11):
214 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
215 sshaw 5123 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 sshaw 5132 dom.addToRHS(rhs, [("Y",Y)],
227     dom.createAssembler("DefaultAssembler", []))
228     dom.addToRHS(lhs, [("D",D)],
229     dom.createAssembler("DefaultAssembler", []))
230 sshaw 5123
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 sshaw 5478 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 sshaw 5123 def test_Rectangle_Function_gradient(self): #expanded and non-expanded
275 sshaw 5132 ranks = getMPISizeWorld()
276 sshaw 5123 for expanded in [True, False]:
277     for order in range(2,11):
278 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
279 sshaw 5123 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 sshaw 5139 ranks = getMPISizeWorld()
293 sshaw 5123 for expanded in [True, False]:
294     for order in range(2,11):
295 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
296 sshaw 5123 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 sshaw 5132 ranks = getMPISizeWorld()
311 sshaw 5123 for order in range(2,11):
312 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=100, l1=100)
313 sshaw 5123 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 sshaw 5132 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 sshaw 5123
332     def test_Brick_ContinuousFunction_gradient(self):
333 sshaw 5139 ranks = getMPISizeWorld()
334 sshaw 5123 for order in range(2,11):
335 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, l0=100, l1=100, l2=100, d1=ranks)
336 sshaw 5123 X = dom.getX()
337     u = X[0] + X[1] + X[2] + 1
338     v = Lsup(grad(u) - 1)
339 sshaw 5132 self.assertLess(v, 1e-10, "order %d, %g >= 1e-10, %s"%(order, v,
340     str(grad(u)-1)))
341 sshaw 5123 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 sshaw 5132 ranks = getMPISizeWorld()
366 sshaw 5123 for order in range(2,11):
367 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=6, l1=6)
368 sshaw 5123 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 sshaw 5139 ranks = getMPISizeWorld()
387 sshaw 5123 for order in range(2,11):
388 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=ranks, l2=6, d1=ranks)
389 sshaw 5123 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 sshaw 5132 ranks = getMPISizeWorld()
408 sshaw 5123 for order in range(2,11):
409     size = 6
410 sshaw 5132 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=size, l1=size)
411 sshaw 5123 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 sshaw 5132 self.assertLess(abs(integral - actual)/actual, 1e-11,
418     "too much variance in integral result (order %d, degrees %d %d)"%(order, k, l))
419 sshaw 5123
420     def test_Brick_integration(self):
421 sshaw 5139 ranks = getMPISizeWorld()
422 sshaw 5123 for order in range(2,11):
423     size = 6
424 sshaw 5139 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
425 sshaw 5123 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 sshaw 5132 actual = size**(k+1) * size**(l+1) * size**(m+1) \
432     /((k+1.)*(l+1.)*(m+1.))
433 sshaw 5123 res = abs(integral - actual)/actual
434 sshaw 5132 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 sshaw 5123
438 sshaw 5132 @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 sshaw 5123 for order in range(2,11):
443 sshaw 5132 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 sshaw 5123
447 sshaw 5132 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 sshaw 5123
461 sshaw 5132 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 sshaw 5123
478 sshaw 5132 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 sshaw 5139 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 sshaw 5132
513 sshaw 5139 @unittest.skipIf(getMPISizeWorld() == 1,
514 sshaw 5132 "only works with more than one rank")
515 sshaw 5123 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 sshaw 5132 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 sshaw 5123
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 sshaw 5132 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 sshaw 5123
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 sshaw 5132 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 sshaw 5139
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 sshaw 5123
567 sshaw 5139 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 sshaw 5170 (Brick(order, 2*half, 2*half, 2, l0=half, l1=half,
596 sshaw 5139 d0=half, d1=half), 0,1),
597 sshaw 5170 (Brick(order, 2*half, 2, 2*half, l0=half, l2=half,
598 sshaw 5139 d0=half, d2=half), 0,2),
599 sshaw 5170 (Brick(order, 2, 2*half, 2*half, l1=half, l2=half,
600 sshaw 5139 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 sshaw 5123 if __name__ == '__main__':
627     run_tests(__name__, exit_on_failure=True)
628    

  ViewVC Help
Powered by ViewVC 1.1.26