/[escript]/trunk/escript/test/python/test_pdetools.py
ViewVC logotype

Contents of /trunk/escript/test/python/test_pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1695 - (show annotations)
Fri Aug 8 04:28:28 2008 UTC (11 years, 2 months ago) by artak
File MIME type: text/x-python
File size: 27748 byte(s)
test_NoPDE_vector_variable is changed to avoid 0/0 division
1 #
2 # $Id$
3 #
4 #######################################################
5 #
6 # Copyright 2003-2007 by ACceSS MNRF
7 # Copyright 2007 by University of Queensland
8 #
9 # http://esscc.uq.edu.au
10 # Primary Business: Queensland, Australia
11 # Licensed under the Open Software License version 3.0
12 # http://www.opensource.org/licenses/osl-3.0.php
13 #
14 #######################################################
15 #
16
17 """
18 Test suite for the pdetools module
19
20 The tests must be linked with a Domain class object in the setUp method:
21
22 from esys.finley import Rectangle
23 class Test_LinearPDEOnFinley(Test_LinearPDE):
24 RES_TOL=1.e-8
25 def setUp(self):
26 self.domain = Rectangle(10,10,2)
27 def tearDown(self):
28 del self.domain
29 suite = unittest.TestSuite()
30 suite.addTest(unittest.makeSuite(Test_LinearPDEOnFinley))
31 unittest.TextTestRunner(verbosity=2).run(suite)
32
33 @var __author__: name of author
34 @var __copyright__: copyrights
35 @var __license__: licence agreement
36 @var __url__: url entry point on documentation
37 @var __version__: version
38 @var __date__: date of the version
39 """
40
41 __author__="Lutz Gross, l.gross@uq.edu.au"
42 __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
43 http://www.access.edu.au
44 Primary Business: Queensland, Australia"""
45 __license__="""Licensed under the Open Software License version 3.0
46 http://www.opensource.org/licenses/osl-3.0.php"""
47 __url__="http://www.iservo.edu.au/esys/escript"
48 __version__="$Revision$"
49 __date__="$Date$"
50
51 import unittest
52 from esys.escript import *
53 from esys.escript.pdetools import Locator,Projector,TimeIntegrationManager,NoPDE,PCG, IterationHistory, ArithmeticTuple, GMRES
54
55 class Test_pdetools_noLumping(unittest.TestCase):
56 DEBUG=False
57 VERBOSE=False
58 def test_TimeIntegrationManager_scalar(self):
59 t=0.
60 dt=0.1
61 tm=TimeIntegrationManager(0.,p=1)
62 while t<1.:
63 t+=dt
64 tm.checkin(dt,t)
65 v_guess=tm.extrapolate(dt)
66 self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong")
67
68 def test_TimeIntegrationManager_vector(self):
69 t=0.
70 dt=0.3
71 tm=TimeIntegrationManager(0.,0.,p=1)
72 while t<1.:
73 t+=dt
74 tm.checkin(dt,t,3*t)
75 v_guess=tm.extrapolate(dt)
76 e=max(abs(v_guess[0]-(tm.getTime()+dt)),abs(v_guess[1]-(tm.getTime()+dt)*3.))
77 self.failUnless(e<self.RES_TOL,"extrapolation is wrong")
78
79 def test_Locator(self):
80 x=self.domain.getX()
81 l=Locator(self.domain,numarray.ones((self.domain.getDim(),)))
82 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain")
83
84 l=Locator(ContinuousFunction(self.domain),numarray.ones((self.domain.getDim(),)))
85 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space")
86
87 xx=l.getX()
88 self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type")
89 self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong")
90 xx=l(x)
91 self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type")
92 self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector")
93 xx=l(x[0]+x[1])
94 self.failUnless(isinstance(xx,float),"wrong scalar type")
95 self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar")
96
97 def test_Locator_withList(self):
98 x=self.domain.getX()
99 arg=[numarray.ones((self.domain.getDim(),)), numarray.zeros((self.domain.getDim(),))]
100 l=Locator(self.domain,arg)
101 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain")
102
103 l=Locator(ContinuousFunction(self.domain),arg)
104 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space")
105
106 xx=l.getX()
107 self.failUnless(isinstance(xx,list),"list expected")
108 for i in range(len(xx)):
109 self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item"%i)
110 self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong"%i)
111 xx=l(x)
112 self.failUnless(isinstance(xx,list),"list expected (2)")
113 for i in range(len(xx)):
114 self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item (2)"%i)
115 self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong (2)"%i)
116 xx=l(x[0]+x[1])
117 self.failUnless(isinstance(xx,list),"list expected (3)")
118 for i in range(len(xx)):
119 self.failUnless(isinstance(xx[i],float),"wrong scalar type")
120 self.failUnless(abs(xx[i]-(arg[i][0]+arg[i][1]))<self.RES_TOL,"value wrong scalar")
121
122 def testProjector_rank0(self):
123 x=ContinuousFunction(self.domain).getX()
124 p=Projector(self.domain,reduce=False,fast=False)
125 td_ref=x[0]
126 td=p(td_ref.interpolate(Function(self.domain)))
127 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
128
129 def testProjector_rank1(self):
130 x=ContinuousFunction(self.domain).getX()
131 p=Projector(self.domain,reduce=False,fast=False)
132 td_ref=x
133 td=p(td_ref.interpolate(Function(self.domain)))
134 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
135
136 def testProjector_rank2(self):
137 x=ContinuousFunction(self.domain).getX()
138 p=Projector(self.domain,reduce=False,fast=False)
139 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
140 td=p(td_ref.interpolate(Function(self.domain)))
141 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
142
143 def testProjector_rank3(self):
144 x=ContinuousFunction(self.domain).getX()
145 p=Projector(self.domain,reduce=False,fast=False)
146 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
147 td=p(td_ref.interpolate(Function(self.domain)))
148 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
149
150 def testProjector_rank4(self):
151 x=ContinuousFunction(self.domain).getX()
152 p=Projector(self.domain,reduce=False,fast=False)
153 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
154 td=p(td_ref.interpolate(Function(self.domain)))
155 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
156
157
158 def testProjector_rank0_reduced(self):
159 x=ContinuousFunction(self.domain).getX()
160 p=Projector(self.domain,reduce=True,fast=False)
161 td_ref=x[0]
162 td=p(td_ref.interpolate(Function(self.domain)))
163 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
164
165 def testProjector_rank1_reduced(self):
166 x=ContinuousFunction(self.domain).getX()
167 p=Projector(self.domain,reduce=True,fast=False)
168 td_ref=x
169 td=p(td_ref.interpolate(Function(self.domain)))
170 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
171
172 def testProjector_rank2_reduced(self):
173 x=ContinuousFunction(self.domain).getX()
174 p=Projector(self.domain,reduce=True,fast=False)
175 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
176 td=p(td_ref.interpolate(Function(self.domain)))
177 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
178
179 def testProjector_rank3_reduced(self):
180 x=ContinuousFunction(self.domain).getX()
181 p=Projector(self.domain,reduce=True,fast=False)
182 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
183 td=p(td_ref.interpolate(Function(self.domain)))
184 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
185
186 def testProjector_rank4_reduced(self):
187 x=ContinuousFunction(self.domain).getX()
188 p=Projector(self.domain,reduce=True,fast=False)
189 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
190 td=p(td_ref.interpolate(Function(self.domain)))
191 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
192
193 def testProjector_rank0_with_reduced_input(self):
194 x=ContinuousFunction(self.domain).getX()
195 p=Projector(self.domain,reduce=False,fast=False)
196 td_ref=x[0]
197 td=p(td_ref.interpolate(Function(self.domain)))
198 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
199
200 def testProjector_rank1_with_reduced_input(self):
201 x=ContinuousFunction(self.domain).getX()
202 p=Projector(self.domain,reduce=False,fast=False)
203 td_ref=x
204 td=p(td_ref.interpolate(Function(self.domain)))
205 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
206
207 def testProjector_rank2_with_reduced_input(self):
208 x=ContinuousFunction(self.domain).getX()
209 p=Projector(self.domain,reduce=False,fast=False)
210 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
211 td=p(td_ref.interpolate(Function(self.domain)))
212 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
213
214 def testProjector_rank3_with_reduced_input(self):
215 x=ContinuousFunction(self.domain).getX()
216 p=Projector(self.domain,reduce=False,fast=False)
217 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
218 td=p(td_ref.interpolate(Function(self.domain)))
219 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
220
221 def testProjector_rank4_with_reduced_input(self):
222 x=ContinuousFunction(self.domain).getX()
223 p=Projector(self.domain,reduce=False,fast=False)
224 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
225 td=p(td_ref.interpolate(Function(self.domain)))
226 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
227
228
229 def testProjector_rank0_reduced_with_reduced_input(self):
230 x=ContinuousFunction(self.domain).getX()
231 p=Projector(self.domain,reduce=True,fast=False)
232 td_ref=1.
233 td=p(Data(td_ref,ReducedFunction(self.domain)))
234 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
235
236 def testProjector_rank1_reduced_with_reduced_input(self):
237 x=ContinuousFunction(self.domain).getX()
238 p=Projector(self.domain,reduce=True,fast=False)
239 td_ref=numarray.array([1.,2.,3.])
240 td=p(Data(td_ref,ReducedFunction(self.domain)))
241 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
242
243 def testProjector_rank2_reduced_with_reduced_input(self):
244 x=ContinuousFunction(self.domain).getX()
245 p=Projector(self.domain,reduce=True,fast=False)
246 td_ref=numarray.array([[11.,12.],[21,22.]])
247 td=p(Data(td_ref,ReducedFunction(self.domain)))
248 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
249
250 def testProjector_rank3_reduced_with_reduced_input(self):
251 x=ContinuousFunction(self.domain).getX()
252 p=Projector(self.domain,reduce=True,fast=False)
253 td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
254 td=p(Data(td_ref,ReducedFunction(self.domain)))
255 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
256
257 def testProjector_rank4_reduced_with_reduced_input(self):
258 x=ContinuousFunction(self.domain).getX()
259 p=Projector(self.domain,reduce=True,fast=False)
260 td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
261 td=p(Data(td_ref,ReducedFunction(self.domain)))
262 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
263
264
265 def test_NoPDE_scalar_missing_r(self):
266 p=NoPDE(self.domain)
267 x=self.domain.getX()
268 msk=whereZero(x[0])
269 p.setValue(D=1.,Y=1.,q=msk)
270 u=p.getSolution()
271 u_ex=(1.-msk)
272 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
273
274 def test_NoPDE_scalar_missing_Y(self):
275 p=NoPDE(self.domain)
276 x=self.domain.getX()
277 msk=whereZero(x[0])
278 p.setValue(D=1.,q=msk,r=2.)
279 u=p.getSolution()
280 u_ex=msk*2.
281 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
282
283 def test_NoPDE_scalar_constant(self):
284 p=NoPDE(self.domain)
285 x=self.domain.getX()
286 msk=whereZero(x[0])
287 p.setValue(D=1.,Y=1.,q=msk,r=2.)
288 u=p.getSolution()
289 u_ex=(1.-msk)+msk*2.
290 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
291
292 def test_NoPDE_scalar_variable(self):
293 p=NoPDE(self.domain)
294 x=self.domain.getX()
295 msk=whereZero(x[0])
296 p.setValue(D=x[0],Y=2*x[0],q=msk,r=2.)
297 u=p.getSolution()
298 u_ex=2.
299 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
300
301 def test_NoPDE_vector_missing_Y(self):
302 p=NoPDE(self.domain)
303 x=self.domain.getX()
304 msk=whereZero(x[0])*[1.,0.]
305 p.setValue(D=numarray.ones([2]),q=msk,r=2.)
306 u=p.getSolution()
307 u_ex=msk*2.
308 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
309
310 def test_NoPDE_vector_missing_r(self):
311 p=NoPDE(self.domain)
312 x=self.domain.getX()
313 msk=whereZero(x[0])*[1.,0.]
314 p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk)
315 u=p.getSolution()
316 u_ex=(1.-msk)
317 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
318
319 def test_NoPDE_vector_constant(self):
320 p=NoPDE(self.domain)
321 x=self.domain.getX()
322 msk=whereZero(x[0])*[1.,0.]
323 p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk,r=2.)
324 u=p.getSolution()
325 u_ex=(1.-msk)+msk*2.
326 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
327
328 def test_NoPDE_vector_variable(self):
329 p=NoPDE(self.domain)
330 x=self.domain.getX()
331 msk=whereZero(x[0])*[1.,0.]
332 p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.)
333 u=p.getSolution()
334 u_ex=2.
335 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
336
337 def testPCG(self):
338 from numarray import array,matrixmultiply, zeros, dot, size, Float64
339 from math import sqrt
340 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
341 5.834798554135237e-01, -3.704394311709722e+00,
342 5.765369186984777e+00, -1.309786358737351e+01,
343 2.522087134507148e+01, -3.393956279045637e+01,
344 1.046856914770830e+02, -2.447764190849540e+02],
345 [ -2.391895572674098e-01, 1.256797283910693e+02,
346 -9.188270412920813e-01, 1.300169538880688e+00,
347 -5.353714719231424e-01, 2.674709444667012e+00,
348 -1.116097841269580e+01, 2.801193427514478e+01,
349 -3.877806125898224e+01, 3.063505753648256e+01],
350 [ 5.834798554135237e-01, -9.188270412920813e-01,
351 6.240841811806843e+01, -8.176289504109282e-01,
352 1.447935098417076e-01, -9.721424148655324e-01,
353 6.713551574117577e-01, -3.656297654168375e+00,
354 7.015141656913973e+00, -4.195525932156250e+01],
355 [ -3.704394311709722e+00, 1.300169538880688e+00,
356 -8.176289504109282e-01, 3.604980536782198e+01,
357 -6.241238423759328e-01, 1.142345320047869e+00,
358 -3.438816797096519e+00, 5.854857481367470e+00,
359 -4.524311288596452e+00, 1.136590280389803e+01],
360 [ 5.765369186984777e+00, -5.353714719231424e-01,
361 1.447935098417076e-01, -6.241238423759328e-01,
362 2.953997190215862e+01, -9.474729233464712e-01,
363 1.883516378345809e+00, -1.906274765704230e+00,
364 4.401859671778645e+00, -1.064573816075257e+01],
365 [ -1.309786358737351e+01, 2.674709444667012e+00,
366 -9.721424148655324e-01, 1.142345320047869e+00,
367 -9.474729233464712e-01, 2.876998216302979e+01,
368 -4.853065259692995e-01, 7.088596468102618e-01,
369 -8.972224295152829e-01, 5.228606946522749e+00],
370 [ 2.522087134507148e+01, -1.116097841269580e+01,
371 6.713551574117577e-01, -3.438816797096519e+00,
372 1.883516378345809e+00, -4.853065259692995e-01,
373 5.121175860935919e+01, -3.523133115905478e-01,
374 1.782136702229135e+00, -1.560849559916187e+00],
375 [ -3.393956279045637e+01, 2.801193427514478e+01,
376 -3.656297654168375e+00, 5.854857481367470e+00,
377 -1.906274765704230e+00, 7.088596468102618e-01,
378 -3.523133115905478e-01, 8.411681423853814e+01,
379 -5.238590858177903e-01, 1.515872114883926e+00],
380 [ 1.046856914770830e+02, -3.877806125898224e+01,
381 7.015141656913973e+00, -4.524311288596452e+00,
382 4.401859671778645e+00, -8.972224295152829e-01,
383 1.782136702229135e+00, -5.238590858177903e-01,
384 1.797889693808014e+02, -8.362340479938084e-01],
385 [ -2.447764190849540e+02, 3.063505753648256e+01,
386 -4.195525932156250e+01, 1.136590280389803e+01,
387 -1.064573816075257e+01, 5.228606946522749e+00,
388 -1.560849559916187e+00, 1.515872114883926e+00,
389 -8.362340479938084e-01, 3.833719335346630e+02]])
390 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
391 0.807186823427233, 0.48950999450145, 0.995486532098031,
392 0.351243009576568, 0.704352576819321, 0.850648989740204,
393 0.314596738052894])
394 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
395 30.344553414038817, 15.247917439094513, 24.060664905403492,
396 27.210293789825833, 47.122067744075842, 199.267136417856847,
397 -8.7934289814322 ])
398
399 def Ap(x):
400 return matrixmultiply(A,x)
401 def Ms(b):
402 out=zeros((size(b),),Float64)
403 for i in xrange(size(b)):
404 out[i]=b[i]/A[i,i]
405 return out
406
407 tol=1.e-4
408 x,r=PCG(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium,x=x_ref*1.5, iter_max=12)
409 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
410 self.failUnless(Lsup(r-(b-matrixmultiply(A,x)))<=Lsup(b)*EPSILON*100.,"wrong solution")
411
412 def testGMRES(self):
413 from numarray import array,matrixmultiply, zeros, dot, size, Float64
414 from math import sqrt
415 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
416 5.834798554135237e-01, -3.704394311709722e+00,
417 5.765369186984777e+00, -1.309786358737351e+01,
418 2.522087134507148e+01, -3.393956279045637e+01,
419 1.046856914770830e+02, -2.447764190849540e+02],
420 [ -2.391895572674098e-01, 1.256797283910693e+02,
421 -9.188270412920813e-01, 1.300169538880688e+00,
422 -5.353714719231424e-01, 2.674709444667012e+00,
423 -1.116097841269580e+01, 2.801193427514478e+01,
424 -3.877806125898224e+01, 3.063505753648256e+01],
425 [ 5.834798554135237e-01, -9.188270412920813e-01,
426 6.240841811806843e+01, -8.176289504109282e-01,
427 1.447935098417076e-01, -9.721424148655324e-01,
428 6.713551574117577e-01, -3.656297654168375e+00,
429 7.015141656913973e+00, -4.195525932156250e+01],
430 [ -3.704394311709722e+00, 1.300169538880688e+00,
431 -8.176289504109282e-01, 3.604980536782198e+01,
432 -6.241238423759328e-01, 1.142345320047869e+00,
433 -3.438816797096519e+00, 5.854857481367470e+00,
434 -4.524311288596452e+00, 1.136590280389803e+01],
435 [ 5.765369186984777e+00, -5.353714719231424e-01,
436 1.447935098417076e-01, -6.241238423759328e-01,
437 2.953997190215862e+01, -9.474729233464712e-01,
438 1.883516378345809e+00, -1.906274765704230e+00,
439 4.401859671778645e+00, -1.064573816075257e+01],
440 [ -1.309786358737351e+01, 2.674709444667012e+00,
441 -9.721424148655324e-01, 1.142345320047869e+00,
442 -9.474729233464712e-01, 2.876998216302979e+01,
443 -4.853065259692995e-01, 7.088596468102618e-01,
444 -8.972224295152829e-01, 5.228606946522749e+00],
445 [ 2.522087134507148e+01, -1.116097841269580e+01,
446 6.713551574117577e-01, -3.438816797096519e+00,
447 1.883516378345809e+00, -4.853065259692995e-01,
448 5.121175860935919e+01, -3.523133115905478e-01,
449 1.782136702229135e+00, -1.560849559916187e+00],
450 [ -3.393956279045637e+01, 2.801193427514478e+01,
451 -3.656297654168375e+00, 5.854857481367470e+00,
452 -1.906274765704230e+00, 7.088596468102618e-01,
453 -3.523133115905478e-01, 8.411681423853814e+01,
454 -5.238590858177903e-01, 1.515872114883926e+00],
455 [ 1.046856914770830e+02, -3.877806125898224e+01,
456 7.015141656913973e+00, -4.524311288596452e+00,
457 4.401859671778645e+00, -8.972224295152829e-01,
458 1.782136702229135e+00, -5.238590858177903e-01,
459 1.797889693808014e+02, -8.362340479938084e-01],
460 [ -2.447764190849540e+02, 3.063505753648256e+01,
461 -4.195525932156250e+01, 1.136590280389803e+01,
462 -1.064573816075257e+01, 5.228606946522749e+00,
463 -1.560849559916187e+00, 1.515872114883926e+00,
464 -8.362340479938084e-01, 3.833719335346630e+02]])
465 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
466 0.807186823427233, 0.48950999450145, 0.995486532098031,
467 0.351243009576568, 0.704352576819321, 0.850648989740204,
468 0.314596738052894])
469 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
470 30.344553414038817, 15.247917439094513, 24.060664905403492,
471 27.210293789825833, 47.122067744075842, 199.267136417856847,
472 -8.7934289814322 ])
473
474 def Ap(x):
475 return matrixmultiply(A,x)
476 def Ms(b):
477 out=zeros((size(b),),Float64)
478 for i in xrange(size(b)):
479 out[i]=b[i]/A[i,i]
480 return out
481
482 tol=1.e-4
483 x=GMRES(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium2,x=x_ref*1.5, iter_max=12)
484 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
485
486 def testArithmeticTuple(self):
487 a=ArithmeticTuple(1.,2.)
488 self.failUnless(len(a)==2,"wrong length")
489 self.failUnless(a[0]==1.,"wrong first item")
490 self.failUnless(a[1]==2.,"wrong second item")
491 c=a*6.
492 self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple")
493 self.failUnless(len(c)==2,"c has wrong length")
494 self.failUnless(c[0]==6.,"c has wrong first item")
495 self.failUnless(c[1]==12.,"c has wrong second item")
496 b=5.*a
497 self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple")
498 self.failUnless(len(b)==2,"b has wrong length")
499 self.failUnless(b[0]==5.,"b has wrong first item")
500 self.failUnless(b[1]==10.,"b has wrong second item")
501 a+=ArithmeticTuple(3.,4.)
502 self.failUnless(a[0]==4.,"wrong first item of inplace update")
503 self.failUnless(a[1]==6.,"wrong second item of inplace update")
504
505
506
507 class Test_pdetools(Test_pdetools_noLumping):
508 def testProjector_rank0_fast_reduced(self):
509 x=ContinuousFunction(self.domain).getX()
510 h=Lsup(self.domain.getSize())
511 p=Projector(self.domain,reduce=True,fast=True)
512 td_ref=x[0]
513 td=p(td_ref.interpolate(Function(self.domain)))
514 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
515
516 def testProjector_rank1_fast_reduced(self):
517 x=ContinuousFunction(self.domain).getX()
518 h=Lsup(self.domain.getSize())
519 p=Projector(self.domain,reduce=True,fast=True)
520 td_ref=x
521 td=p(td_ref.interpolate(Function(self.domain)))
522 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
523
524 def testProjector_rank2_fast_reduced(self):
525 x=ContinuousFunction(self.domain).getX()
526 h=Lsup(self.domain.getSize())
527 p=Projector(self.domain,reduce=True,fast=True)
528 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
529 td=p(td_ref.interpolate(Function(self.domain)))
530 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
531
532 def testProjector_rank3_fast_reduced(self):
533 x=ContinuousFunction(self.domain).getX()
534 h=Lsup(self.domain.getSize())
535 p=Projector(self.domain,reduce=True,fast=True)
536 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
537 td=p(td_ref.interpolate(Function(self.domain)))
538 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
539
540 def testProjector_rank4_fast_reduced(self):
541 x=ContinuousFunction(self.domain).getX()
542 h=Lsup(self.domain.getSize())
543 p=Projector(self.domain,reduce=True,fast=True)
544 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
545 td=p(td_ref.interpolate(Function(self.domain)))
546 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
547
548 def testProjector_rank0_fast_reduced_with_reduced_input(self):
549 x=ContinuousFunction(self.domain).getX()
550 h=Lsup(self.domain.getSize())
551 p=Projector(self.domain,reduce=True,fast=True)
552 td_ref=1.
553 td=p(Data(td_ref,ReducedFunction(self.domain)))
554 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
555
556 def testProjector_rank1_fast_reduced_with_reduced_input(self):
557 x=ContinuousFunction(self.domain).getX()
558 h=Lsup(self.domain.getSize())
559 p=Projector(self.domain,reduce=True,fast=True)
560 td_ref=numarray.array([1.,2.,3.])
561 td=p(Data(td_ref,ReducedFunction(self.domain)))
562 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
563
564 def testProjector_rank2_fast_reduced_with_reduced_input(self):
565 x=ContinuousFunction(self.domain).getX()
566 h=Lsup(self.domain.getSize())
567 p=Projector(self.domain,reduce=True,fast=True)
568 td_ref=numarray.array([[11.,12.],[21,22.]])
569 td=p(Data(td_ref,ReducedFunction(self.domain)))
570 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
571
572 def testProjector_rank3_fast_reduced_with_reduced_input(self):
573 x=ContinuousFunction(self.domain).getX()
574 h=Lsup(self.domain.getSize())
575 p=Projector(self.domain,reduce=True,fast=True)
576 td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
577 td=p(Data(td_ref,ReducedFunction(self.domain)))
578 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
579
580 def testProjector_rank4_fast_reduced_with_reduced_input(self):
581 x=ContinuousFunction(self.domain).getX()
582 h=Lsup(self.domain.getSize())
583 p=Projector(self.domain,reduce=True,fast=True)
584 td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
585 td=p(Data(td_ref,ReducedFunction(self.domain)))
586 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
587

  ViewVC Help
Powered by ViewVC 1.1.26