/[escript]/branches/arrexp_trunk2098/escript/test/python/test_pdetools.py
ViewVC logotype

Contents of /branches/arrexp_trunk2098/escript/test/python/test_pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2145 - (show annotations)
Wed Dec 10 00:48:53 2008 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 32853 byte(s)
Branch commit.
This version of the python files passes unit tests compiled with 
usenumarray=no and usenumarray=yes
That is, it does not use any functions returning 
boost::python::numeric::array


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

  ViewVC Help
Powered by ViewVC 1.1.26