/[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 6676 - (show annotations)
Mon May 21 00:34:56 2018 UTC (3 years, 2 months ago) by uqagarro
File MIME type: text/x-python
File size: 36745 byte(s)
New version with gradient test corrected

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2018 by The University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Apache License, version 2.0
9 # http://www.apache.org/licenses/LICENSE-2.0
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-2018 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(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 %e >= %e"%(order,
233 res, self.TOLERANCE)).format("" if expanded else "un-"))
234
235 class Test_Speckley(unittest.TestCase):
236 TOLERANCE = 5e-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: %e != 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: %e != 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: %e >= %e"%(fs, order, sup(d[1]+0.5)-3*ranks, self.TOLERANCE))
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: %e != 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: %e != 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: %e >= %e"%(fs, order, sup(d[1]+0.5)-3*ranks, self.TOLERANCE))
270 self.assertLess(sup(d[2]+0.5) - 3, self.TOLERANCE,
271 "reduced->%s failure with order %d: %e != 3"%(fs, order, sup(d[2]+0.5)))
272
273 def test_Rectangle_Function_gradient_AntiSimetric(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=ContinuousFunction(dom).getX()
279 m=-2*1j+4
280 b=3+8*1j
281 f=Vector(0., Solution(dom))
282 sample_data=[[0.000000000+0.000000000j, -4.000000000+2.000000000j],[4.000000000-2.000000000j, 0.000000000+0.000000000j]]
283 sample=Data(sample_data, Function(dom), True)
284 sample[0,0]=0.000000000+0.000000000j
285 sample[0,1]=-4.000000000+2.000000000j
286 sample[1,0]=4.000000000-2.000000000j
287 sample[1,1]=0.000000000+0.000000000j
288 for i in range(dom.getDim()):
289 for j in range(dom.getDim()):
290 f[i]+=m*(i-j)*x[j]+b-i*j
291 g=grad(f)
292 self.assertLess(Lsup(g-sample), 1e-9,"single component failure, order %d%s, %e >= 1e-10"%(order,(" expanded" if expanded else ""), Lsup(g-sample)))
293
294 def test_Rectangle_Function_gradient_Simetric(self): #expanded and non-expanded
295 ranks = getMPISizeWorld()
296 for expanded in [True, False]:
297 for order in range(2,11):
298 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
299 x=ContinuousFunction(dom).getX()
300 m=-2*1j+4
301 b=3+8*1j
302 f=Vector(0., Solution(dom))
303 sample_data=[[0.000000000+0.000000000j, 4.000000000-2.000000000j],[4.000000000-2.000000000j, 8.000000000-4.000000000j]]
304 sample=Data(sample_data, Function(dom), True)
305 sample[0,0]=0.000000000+0.000000000j
306 sample[0,1]=4.000000000-2.000000000j
307 sample[1,0]=4.000000000-2.000000000j
308 sample[1,1]=8.000000000-4.000000000j
309 for i in range(dom.getDim()):
310 for j in range(dom.getDim()):
311 f[i]+=m*(i+j)*x[j]+b-i*j
312 g=grad(f)
313 self.assertLess(Lsup(g-sample), 1e-9,"single component failure, order %d%s, %e >= 1e-10"%(order,(" expanded" if expanded else ""), Lsup(g-sample)))
314
315 def test_Brick_Function_gradient_AntiSimetric(self): #expanded and non-expanded
316 ranks = getMPISizeWorld()
317 for expanded in [True, False]:
318 for order in range(2,11):
319 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
320 x=ContinuousFunction(dom).getX()
321 m=-2*1j+4
322 b=3+8*1j
323 f=Vector(0., Solution(dom))
324 sample_data=[[0.000000000+0.000000000j, -4.000000000+2.000000000j, -8.000000000+4.000000000j],[4.000000000-2.000000000j, 0.000000000+0.000000000j, -4.000000000+2.000000000j],[8.000000000-4.000000000j, 4.000000000-2.000000000j, 0.000000000+0.000000000j]]
325 sample=Data(sample_data, Function(dom), True)
326 sample[0,0]=0.000000000+0.000000000j
327 sample[0,1]=-4.000000000+2.000000000j
328 sample[0,2]=-8.000000000+4.000000000j
329 sample[1,0]=4.000000000-2.000000000j
330 sample[1,1]=0.000000000+0.000000000j
331 sample[1,2]=-4.000000000+2.000000000j
332 sample[2,0]=8.000000000-4.000000000j
333 sample[2,1]=4.000000000-2.000000000j
334 sample[2,2]=0.000000000+0.000000000j
335 for i in range(dom.getDim()):
336 for j in range(dom.getDim()):
337 f[i]+=m*(i-j)*x[j]+b-i*j
338 g=grad(f)
339 self.assertLess(Lsup(g-sample), 1e-9,"single component failure, order %d%s, %e >= 1e-10"%(order,(" expanded" if expanded else ""), Lsup(g-sample)))
340
341 def test_Brick_Function_gradient_Simetric(self): #expanded and non-expanded
342 ranks = getMPISizeWorld()
343 for expanded in [True, False]:
344 for order in range(2,11):
345 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
346 x=ContinuousFunction(dom).getX()
347 m=-2*1j+4
348 b=3+8*1j
349 f=Vector(0., Solution(dom))
350 sample_data=[[0.000000000+0.000000000j, 4.000000000-2.000000000j, 8.000000000-4.000000000j],[4.000000000-2.000000000j, 8.000000000-4.000000000j, 12.000000000-6.000000000j],[8.000000000-4.000000000j, 12.000000000-6.000000000j, 16.000000000-8.000000000j]]
351 sample=Data(sample_data, Function(dom), True)
352 sample[0,0]=0.000000000+0.000000000j
353 sample[0,1]=4.000000000-2.000000000j
354 sample[0,2]=8.000000000-4.000000000j
355 sample[1,0]=4.000000000-2.000000000j
356 sample[1,1]=8.000000000-4.000000000j
357 sample[1,2]=12.000000000-6.000000000j
358 sample[2,0]=8.000000000-4.000000000j
359 sample[2,1]=12.000000000-6.000000000j
360 sample[2,2]=16.000000000-8.000000000j
361 for i in range(dom.getDim()):
362 for j in range(dom.getDim()):
363 f[i]+=m*(i+j)*x[j]+b-i*j
364 g=grad(f)
365 self.assertLess(Lsup(g-sample), 1e-9,"single component failure, order %d%s, %e >= 1e-10"%(order,(" expanded" if expanded else ""), Lsup(g-sample)))
366
367
368 def test_Rectangle_Function_gradient(self): #expanded and non-expanded
369 ranks = getMPISizeWorld()
370 for expanded in [True, False]:
371 for order in range(2,11):
372 dom = Rectangle(order, 3, 3*ranks, d1=ranks)
373 x = Data(5, Function(dom), True)
374 self.assertLess(Lsup(grad(x)), 1e-10,
375 "single component failure, order %d%s, %e >= 1e-10"%(order,
376 (" expanded" if expanded else ""), Lsup(grad(x))))
377 for data in [[5,1], [-5,-1], [5,1,1e-5]]:
378 x = Data(data, Function(dom), True)
379 g = grad(x)
380 for n,d in enumerate(data):
381 self.assertLess(Lsup(g[n]), 1e-10,
382 "%d-component failure, order %d %sexpanded, %e >= 1e-10"%(len(data),
383 order, ("" if expanded else "un-"), Lsup(g[n])))
384
385 def test_Brick_Function_gradient(self):
386 ranks = getMPISizeWorld()
387 for expanded in [True, False]:
388 for order in range(2,11):
389 dom = Brick(order, 3, 3*ranks, 3, d1=ranks)
390 x = Data(5, Function(dom), True)
391 self.assertLess(Lsup(grad(x)), 1e-10,
392 "single component failure, order %d%s, %e >= 1e-10"%(order,
393 (" expanded" if expanded else ""), Lsup(grad(x))))
394 for data in [[5,1], [-5,-1], [5,1,1e-5]]:
395 x = Data(data, Function(dom), True)
396 g = grad(x)
397 for n,d in enumerate(data):
398 self.assertLess(Lsup(g[n]), 1e-10,
399 "%d-component failure, order %d %sexpanded, %e >= 1e-10"%(len(data),
400 order, ("" if expanded else "un-"), Lsup(g[n])))
401
402
403 def test_Rectangle_ContinuousFunction_gradient(self):
404 ranks = getMPISizeWorld()
405 for order in range(2,11):
406 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=100, l1=100)
407 X = dom.getX()
408 u = X[0] + X[1] + 1
409 v = Lsup(grad(u) - 1)
410 self.assertLess(v, 1e-10, "order %d, %e >= 1e-10, %s"%(order, v, str(grad(u)-1)))
411 for power in range(1, order+1):
412 for power2 in range(1, order+1):
413 a = X[0]**power * X[1]**power2
414 da = grad(a)
415 first = Lsup(da[0] - power*X[0]**(power-1) * X[1]**power2) \
416 /Lsup(power*X[0]**(power-1) * X[1]**power2)
417 second = Lsup(da[1] - power2*X[1]**(power2-1) * X[0]**power) \
418 /Lsup(power2*X[1]**(power2-1) * X[0]**power)
419 self.assertLess(first, 1e-10,
420 "order %d and degree %d,%d, %e >= 1e-9"%(order,
421 power, power2, first))
422 self.assertLess(second, 1e-10,
423 "order %d and degree %d,%d, %e >= 1e-9"%(order,
424 power, power2, second))
425
426 def test_Brick_ContinuousFunction_gradient(self):
427 ranks = getMPISizeWorld()
428 for order in range(2,11):
429 dom = Brick(order, 3, 3*ranks, 3, l0=100, l1=100, l2=100, d1=ranks)
430 X = dom.getX()
431 u = X[0] + X[1] + X[2] + 1
432 v = Lsup(grad(u) - 1)
433 self.assertLess(v, 1e-10, "order %d, %e >= 1e-10, %s"%(order, v,
434 str(grad(u)-1)))
435 for power1 in range(1, order+1, order//2):
436 for power2 in range(1, order+1, order//2):
437 for power3 in range(1, order+1, order//2):
438 a = X[0]**power1 * X[1]**power2 * X[2]**power3
439 da = grad(a)
440
441 temp = power1*X[0]**(power1-1) * X[1]**power2*X[2]**power3
442 first = Lsup(da[0] - temp) / Lsup(temp)
443 temp = power2*X[1]**(power2-1) * X[0]**power1*X[2]**power3
444 second = Lsup(da[1] - temp) / Lsup(temp)
445 temp = power3*X[2]**(power3-1) * X[0]**power1*X[1]**power2
446 third = Lsup(da[2] - temp) / Lsup(temp)
447
448 self.assertLess(first, 1e-10,
449 "order %d and degree %d,%d,%d, %e >= 1e-9"%(order,
450 power1, power2, power3, first))
451 self.assertLess(second, 1e-10,
452 "order %d and degree %d,%d,%d, %e >= 1e-9"%(order,
453 power1, power2, power3, second))
454 self.assertLess(third, 1e-10,
455 "order %d and degree %d,%d,%d, %e >= 1e-9"%(order,
456 power1, power2, power3, third))
457
458 def test_Rectangle_interpolation_continuous_noncontinuous_and_back(self):
459 ranks = getMPISizeWorld()
460 for order in range(2,11):
461 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=6, l1=6)
462 original = Data(5, Function(dom), True)
463 cont = interpolate(original, ContinuousFunction(dom))
464 func = interpolate(cont, Function(dom))
465 self.assertEqual(Lsup(original-func), 0,
466 "interpolation of constant, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
467 x = dom.getX()
468 original = x[0] + x[1] + 1
469 cont = interpolate(original, ContinuousFunction(dom))
470 func = interpolate(cont, Function(dom))
471 self.assertEqual(Lsup(original-func), 0,
472 "interpolation of expanded, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
473 original = whereZero(x[0]-2) + whereZero(x[1]-2)
474 cont = interpolate(original, ContinuousFunction(dom))
475 func = interpolate(cont, Function(dom))
476 self.assertEqual(Lsup(original-func), 0,
477 "interpolation of point, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
478
479 def xtest_Brick_interpolation_continuous_noncontinuous_and_back(self):
480 ranks = getMPISizeWorld()
481 for order in range(2,11):
482 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=ranks, l2=6, d1=ranks)
483 original = Data(5, Function(dom), True)
484 cont = interpolate(original, ContinuousFunction(dom))
485 func = interpolate(cont, Function(dom))
486 self.assertEqual(Lsup(original-func), 0,
487 "interpolation of constant, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
488 x = dom.getX()
489 original = x[0] + x[1] + x[2] + 1
490 cont = interpolate(original, ContinuousFunction(dom))
491 func = interpolate(cont, Function(dom))
492 self.assertEqual(Lsup(original-func), 0,
493 "interpolation of expanded, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
494 original = whereZero(x[0]-2) + whereZero(x[1]-2) + whereZero(x[2] - 2)
495 cont = interpolate(original, ContinuousFunction(dom))
496 func = interpolate(cont, Function(dom))
497 self.assertEqual(Lsup(original-func), 0,
498 "interpolation of point, order %d: original and final not equal, %e != 0"%(order, Lsup(original-func)))
499
500 def xtest_Rectangle_integration(self):
501 ranks = getMPISizeWorld()
502 for order in range(2,11):
503 size = 6
504 dom = Rectangle(order, 3, 3*ranks, d1=ranks, l0=size, l1=size)
505 X = dom.getX()
506 for k in range(1, order * 2):
507 for l in range(1, order * 2):
508 powered = X[0]**k * X[1]**l
509 integral = integrate(powered)
510 actual = size**(k+1) * size**(l+1) / ((k+1.)*(l+1.))
511 self.assertLess(abs(integral - actual)/actual, 1e-11,
512 "too much variance in integral result (order %d, degrees %d %d)"%(order, k, l))
513
514 def xtest_Brick_integration(self):
515 ranks = getMPISizeWorld()
516 for order in range(2,11):
517 size = 6
518 dom = Brick(order, 3, 3*ranks, 3, l0=6, l1=6, l2=6, d1=ranks)
519 X = dom.getX()
520 for k in [1, order, order*2 - 1]:
521 for l in [1, order, order*2 - 1]:
522 for m in [1, order, order*2 - 1]:
523 powered = X[0]**k * X[1]**l * X[2]**m
524 integral = integrate(powered)
525 actual = size**(k+1) * size**(l+1) * size**(m+1) \
526 /((k+1.)*(l+1.)*(m+1.))
527 res = abs(integral - actual)/actual
528 self.assertLess(res, 1e-11,
529 "too much variance in integral result (order %d, degrees %d %d, %e >= 1e-11)"%(order,
530 k, l, res))
531
532 @unittest.skipIf(getMPISizeWorld() == 1,
533 "only works with more than one rank")
534 def xtest_Rectangle_MPI_construction_single_dimensional(self):
535 ranks = getMPISizeWorld()
536 for order in range(2,11):
537 dom = Rectangle(order, 2, 2*ranks, l1=ranks, d1=ranks)
538 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]), ranks+1,
539 "invalid getX() for y-splits order %d"%order)
540
541 filt = whereZero(dom.getX()[1] - 1)
542 for i in range(2,ranks):
543 filt += whereZero(dom.getX()[1] - i)
544 if getMPIRankWorld() % 2:
545 filt *= -1
546 d = Vector(0, Function(dom))
547 d[0] = 1 * filt
548 d[1] = 10 * filt
549 X = interpolate(d, Function(dom))
550 res = interpolate(X, ContinuousFunction(dom))
551 val = Lsup(res)
552 self.assertEqual(val, 0,
553 "summation stage failure for y-splits in order %d"%order)
554
555 dom = Rectangle(2, 2*ranks, 2, l0=ranks, d0=ranks)
556 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]), ranks+1,
557 "invalid getX() for x-splits order %d"%order)
558 filt = whereZero(dom.getX()[0] - 1)
559 for i in range(2,getMPISizeWorld()):
560 filt += whereZero(dom.getX()[0] - i)
561 if getMPIRankWorld() % 2:
562 filt *= -1
563 d = Vector(0, Function(dom))
564 d[0] = 1 * filt
565 d[1] = 10 * filt
566 X = interpolate(d, Function(dom))
567 res = interpolate(X, ContinuousFunction(dom))
568 val = Lsup(res)
569 self.assertEqual(val, 0,
570 "summation stage failure for x-splits in order %d"%order)
571
572 X = interpolate(d+2, Function(dom))
573 res = interpolate(X, ContinuousFunction(dom))
574 val = Lsup(res-2)
575 self.assertEqual(val, 0,
576 "averaging stage failure for x-splits in order %d"%order)
577
578 @unittest.skipIf(getMPISizeWorld() != 4, "requires 4 ranks")
579 def xtest_Rectangle_MPI_construction_multi_dimensional(self):
580 ranks = getMPISizeWorld()
581 half = 2 #change if ranks != 4 (sqrt(ranks))
582 for order in range(2, 11):
583 dom = Rectangle(order, ranks, ranks, l0=half, l1=half,
584 d0=half, d1=half)
585 self.assertEqual(Lsup(dom.getX()[0]), half,
586 "invalid getX() for multidimensional splits in order %d"%order)
587 self.assertEqual(Lsup(dom.getX()[1]), half,
588 "invalid getX() for multidimensional splits in order %d"%order)
589 xfilt = whereZero(dom.getX()[0] - 1) + whereZero(dom.getX()[1] - 1)
590 for i in range(2,half):
591 xfilt += whereZero(dom.getX()[0] - i)
592 xfilt += whereZero(dom.getX()[1] - i)
593 xfilt = whereNonZero(xfilt)
594 if getMPIRankWorld() in [1,2]: #change if ranks != 4
595 xfilt *= -1
596 X = interpolate(xfilt, Function(dom))
597 res = interpolate(X, ContinuousFunction(dom))
598 val = Lsup(res)
599 self.assertEqual(val, 0,
600 "summation failure for mixed-splits in order %d"%order)
601 X = interpolate(xfilt+2, Function(dom))
602 res = interpolate(X, ContinuousFunction(dom))
603 val = Lsup(res-2)
604 self.assertEqual(val, 0,
605 "averaging failure for mixed-splits in order %d"%order)
606
607 @unittest.skipIf(getMPISizeWorld() == 1,
608 "only works with more than one rank")
609 def xtest_Brick_MPI_construction(self):
610 for order in range(2,11):
611 dom = Brick(order, 2*getMPISizeWorld(), 2, 2, d0=getMPISizeWorld())
612 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
613 3, "invalid Lsup(getX()) for x-split order %d"%order)
614 X = interpolate(dom.getX()[0], Function(dom)) \
615 + interpolate(dom.getX()[1], Function(dom))
616 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
617 - (dom.getX()[0] + dom.getX()[1]))
618 self.assertLess(val, 1e-10,
619 "interpolation failure for x-split order %d"%order)
620
621 dom = Brick(order, 2, 2*getMPISizeWorld(), 2, d1=getMPISizeWorld())
622 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
623 3, "invalid Lsup(getX()) for y-split order %d"%order)
624 X = interpolate(dom.getX()[0], Function(dom)) \
625 + interpolate(dom.getX()[1], Function(dom))
626 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
627 - (dom.getX()[0] + dom.getX()[1]))
628 self.assertLess(val, 1e-10,
629 "interpolation failure for y-split order %d"%order)
630
631 dom = Brick(order, 2, 2, 2*getMPISizeWorld(), d2=getMPISizeWorld())
632 self.assertEqual(Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
633 3, "invalid Lsup(getX()) for z-split order %d"%order)
634 X = interpolate(dom.getX()[0], Function(dom)) \
635 + interpolate(dom.getX()[1], Function(dom))
636 val = Lsup(interpolate(X, ContinuousFunction(dom)) \
637 - (dom.getX()[0] + dom.getX()[1]))
638 self.assertLess(val, 1e-10,
639 "interpolation failure for z-split order %d"%order)
640
641 @unittest.skipIf(getMPISizeWorld() == 1, "requires multiple MPI processes")
642 def xtest_Brick_singledim_subdivision(self):
643 ranks = getMPISizeWorld()
644 for dim in range(0,3):
645 label = ["x","y","z"][dim]
646 size = [2,2,2]
647 size[dim] *= ranks
648 lengths = [1,1,1]
649 lengths[dim] *= ranks
650 splits = [1,1,1]
651 splits[dim] *= ranks
652
653 for order in range(2, 11):
654 dom = Brick(order, size[0], size[1], size[2],
655 l0=lengths[0], l1=lengths[1], l2=lengths[2],
656 d0=splits[0], d1=splits[1], d2=splits[2])
657 self.assertEqual(Lsup(dom.getX()[1]+dom.getX()[0]+dom.getX()[2]),
658 ranks+2, "invalid getX() for %s-splits order %d"%(\
659 label, order))
660
661 filt = whereZero(dom.getX()[dim] - 1)
662 for i in range(2,ranks):
663 filt += whereZero(dom.getX()[0] - i)
664 if getMPIRankWorld() % 2:
665 filt *= -1
666 d = Vector(0, Function(dom))
667 d[0] = 1 * filt
668 d[1] = 10 * filt
669 d[2] = 100 * filt
670 X = interpolate(d, Function(dom))
671 res = interpolate(X, ContinuousFunction(dom))
672 val = Lsup(res)
673 self.assertEqual(val, 0,
674 "summation stage failure for %s-splits in order %d,"%(\
675 label, order))
676 X = interpolate(d+2, Function(dom))
677 res = interpolate(X, ContinuousFunction(dom))
678 val = Lsup(res-2)
679 self.assertEqual(val, 0,
680 "averaging stage failure for %s-splits in order %d,"%(\
681 label, order))
682
683 @unittest.skipIf(getMPISizeWorld() != 4, "requires 4 ranks exactly")
684 def xtest_Brick_multidim_subdivision(self):
685 ranks = getMPISizeWorld()
686 half = 2 #change if ranks != 4 (sqrt(ranks))
687 for order in range(2, 11):
688 for dom,dim1,dim2 in [
689 (Brick(order, 2*half, 2*half, 2, l0=half, l1=half,
690 d0=half, d1=half), 0,1),
691 (Brick(order, 2*half, 2, 2*half, l0=half, l2=half,
692 d0=half, d2=half), 0,2),
693 (Brick(order, 2, 2*half, 2*half, l1=half, l2=half,
694 d1=half, d2=half), 1,2)]:
695 self.assertEqual(ranks + 1,
696 Lsup(dom.getX()[0] + dom.getX()[1] + dom.getX()[2]),
697 "invalid getX() for multidimensional split " + \
698 "(dims %d,%d) in order %d"%(dim1,dim2,order))
699 xfilt = whereZero(dom.getX()[dim1] - 1) \
700 + whereZero(dom.getX()[dim2] - 1)
701 for i in range(2,half):
702 xfilt += whereZero(dom.getX()[dim1] - i)
703 xfilt += whereZero(dom.getX()[dim2] - i)
704 xfilt = whereNonZero(xfilt)
705 if getMPIRankWorld() in [1,2]: #change if ranks != 4
706 xfilt *= -1
707 X = interpolate(xfilt, Function(dom))
708 res = interpolate(X, ContinuousFunction(dom))
709 val = Lsup(res)
710 self.assertEqual(val, 0,
711 "summation failure for mixed-splits " \
712 + "(dims %d,%d) in order %d"%(dim1,dim2,order))
713 X = interpolate(xfilt+2, Function(dom))
714 res = interpolate(X, ContinuousFunction(dom))
715 val = Lsup(res-2)
716 self.assertEqual(val, 0,
717 "averaging failure for mixed-splits "\
718 + "(dims %d,%d) in order %d"%(dim1,dim2,order))
719
720 if __name__ == '__main__':
721 run_tests(__name__, exit_on_failure=True)
722

  ViewVC Help
Powered by ViewVC 1.1.26