/[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 2768 - (show annotations)
Mon Nov 23 23:56:14 2009 UTC (9 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 66682 byte(s)
Dodge problem in NoPDE class

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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__="https://launchpad.net/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, ArithmeticTuple, GMRES, MINRES, TFQMR, HomogeneousSaddlePointProblem
51 from esys.escript.pdetools import Defect, NewtonGMRES
52 from numpy.linalg import solve as solve_linear_equations
53
54 class Test_pdetools_noLumping(unittest.TestCase):
55 DEBUG=False
56 VERBOSE=False
57 def test_TimeIntegrationManager_scalar(self):
58 t=0.
59 dt=0.1
60 tm=TimeIntegrationManager(0.,p=1)
61 while t<1.:
62 t+=dt
63 tm.checkin(dt,t)
64 v_guess=tm.extrapolate(dt)
65 self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong")
66
67 def test_TimeIntegrationManager_vector(self):
68 t=0.
69 dt=0.3
70 tm=TimeIntegrationManager(0.,0.,p=1)
71 while t<1.:
72 t+=dt
73 tm.checkin(dt,t,3*t)
74 v_guess=tm.extrapolate(dt)
75 e=max(abs(v_guess[0]-(tm.getTime()+dt)),abs(v_guess[1]-(tm.getTime()+dt)*3.))
76 self.failUnless(e<self.RES_TOL,"extrapolation is wrong")
77
78 def test_Locator(self):
79 x=self.domain.getX()
80 l=Locator(self.domain,numpy.ones((self.domain.getDim(),)))
81 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain")
82
83 l=Locator(ContinuousFunction(self.domain),numpy.ones((self.domain.getDim(),)))
84 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space")
85
86 xx=l.getX()
87 self.failUnless(isinstance(xx,numpy.ndarray),"wrong vector type")
88 self.failUnless(Lsup(xx-numpy.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong")
89 xx=l(x)
90 self.failUnless(isinstance(xx,numpy.ndarray),"wrong vector type")
91 self.failUnless(Lsup(xx-numpy.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector")
92 xx=l(x[0]+x[1])
93 self.failUnless(isinstance(xx,float),"wrong scalar type")
94 self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar")
95
96 # now with interpolation:
97 l=Locator(Function(self.domain),numpy.ones((self.domain.getDim(),)))
98 x2=Function(self.domain).getX()
99 xx=l(x)
100 self.failUnless(isinstance(xx,numpy.ndarray),"wrong vector type")
101 self.failUnless(Lsup(xx-l(x2))<self.RES_TOL,"location wrong")
102 xx=l(x[0]+x[1])
103 self.failUnless(isinstance(xx,float),"wrong scalar type")
104 self.failUnless(abs(xx-l(x2[0])-l(x2[1]))<self.RES_TOL,"value wrong scalar")
105
106
107 def test_Locator_withList(self):
108 x=self.domain.getX()
109 arg=[numpy.ones((self.domain.getDim(),)), numpy.zeros((self.domain.getDim(),))]
110 l=Locator(self.domain,arg)
111 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain")
112
113 l=Locator(ContinuousFunction(self.domain),arg)
114 self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space")
115
116 xx=l.getX()
117 self.failUnless(isinstance(xx,list),"list expected")
118 for i in range(len(xx)):
119 self.failUnless(isinstance(xx[i],numpy.ndarray),"vector expected for %s item"%i)
120 self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong"%i)
121 xx=l(x)
122 self.failUnless(isinstance(xx,list),"list expected (2)")
123 for i in range(len(xx)):
124 self.failUnless(isinstance(xx[i],numpy.ndarray),"vector expected for %s item (2)"%i)
125 self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong (2)"%i)
126 xx=l(x[0]+x[1])
127 self.failUnless(isinstance(xx,list),"list expected (3)")
128 for i in range(len(xx)):
129 self.failUnless(isinstance(xx[i],float),"wrong scalar type")
130 self.failUnless(abs(xx[i]-(arg[i][0]+arg[i][1]))<self.RES_TOL,"value wrong scalar")
131
132 # now with interpolation:
133 l=Locator(Function(self.domain),arg)
134 self.failUnless(Function(self.domain)==l.getFunctionSpace(),"wrong function space")
135 xx=l(x)
136 x2=Function(self.domain).getX()
137 self.failUnless(isinstance(xx,list),"list expected (2)")
138 for i in range(len(xx)):
139 self.failUnless(isinstance(xx[i],numpy.ndarray),"vector expected for %s item (2)"%i)
140 self.failUnless(Lsup(xx[i]-l(x2)[i])<self.RES_TOL,"%s-th location is wrong (2)"%i)
141 xx=l(x[0]+x[1])
142 self.failUnless(isinstance(xx,list),"list expected (3)")
143 for i in range(len(xx)):
144 self.failUnless(isinstance(xx[i],float),"wrong scalar type")
145 self.failUnless(abs(xx[i]-(l(x2[0])[i]+l(x2[1])[i]))<self.RES_TOL,"value wrong scalar")
146
147
148 def testProjector_rank0(self):
149 x=ContinuousFunction(self.domain).getX()
150 p=Projector(self.domain,reduce=False,fast=False)
151 td_ref=x[0]
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 def testProjector_rank1(self):
156 x=ContinuousFunction(self.domain).getX()
157 p=Projector(self.domain,reduce=False,fast=False)
158 td_ref=x
159 td=p(td_ref.interpolate(Function(self.domain)))
160 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
161
162 def testProjector_rank2(self):
163 x=ContinuousFunction(self.domain).getX()
164 p=Projector(self.domain,reduce=False,fast=False)
165 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
166 td=p(td_ref.interpolate(Function(self.domain)))
167 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
168
169 def testProjector_rank3(self):
170 x=ContinuousFunction(self.domain).getX()
171 p=Projector(self.domain,reduce=False,fast=False)
172 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
173 td=p(td_ref.interpolate(Function(self.domain)))
174 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
175
176 def testProjector_rank4(self):
177 x=ContinuousFunction(self.domain).getX()
178 p=Projector(self.domain,reduce=False,fast=False)
179 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
180 td=p(td_ref.interpolate(Function(self.domain)))
181 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
182
183
184 def testProjector_rank0_reduced(self):
185 x=ContinuousFunction(self.domain).getX()
186 p=Projector(self.domain,reduce=True,fast=False)
187 td_ref=x[0]
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_rank1_reduced(self):
192 x=ContinuousFunction(self.domain).getX()
193 p=Projector(self.domain,reduce=True,fast=False)
194 td_ref=x
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_rank2_reduced(self):
199 x=ContinuousFunction(self.domain).getX()
200 p=Projector(self.domain,reduce=True,fast=False)
201 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
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_rank3_reduced(self):
206 x=ContinuousFunction(self.domain).getX()
207 p=Projector(self.domain,reduce=True,fast=False)
208 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(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_rank4_reduced(self):
213 x=ContinuousFunction(self.domain).getX()
214 p=Projector(self.domain,reduce=True,fast=False)
215 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(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_rank0_with_reduced_input(self):
220 x=ContinuousFunction(self.domain).getX()
221 p=Projector(self.domain,reduce=False,fast=False)
222 td_ref=x[0]
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 def testProjector_rank1_with_reduced_input(self):
227 x=ContinuousFunction(self.domain).getX()
228 p=Projector(self.domain,reduce=False,fast=False)
229 td_ref=x
230 td=p(td_ref.interpolate(Function(self.domain)))
231 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
232
233 def testProjector_rank2_with_reduced_input(self):
234 x=ContinuousFunction(self.domain).getX()
235 p=Projector(self.domain,reduce=False,fast=False)
236 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
237 td=p(td_ref.interpolate(Function(self.domain)))
238 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
239
240 def testProjector_rank3_with_reduced_input(self):
241 x=ContinuousFunction(self.domain).getX()
242 p=Projector(self.domain,reduce=False,fast=False)
243 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
244 td=p(td_ref.interpolate(Function(self.domain)))
245 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
246
247 def testProjector_rank4_with_reduced_input(self):
248 x=ContinuousFunction(self.domain).getX()
249 p=Projector(self.domain,reduce=False,fast=False)
250 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
251 td=p(td_ref.interpolate(Function(self.domain)))
252 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
253
254
255 def testProjector_rank0_reduced_with_reduced_input(self):
256 x=ContinuousFunction(self.domain).getX()
257 p=Projector(self.domain,reduce=True,fast=False)
258 td_ref=1.
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 def testProjector_rank1_reduced_with_reduced_input(self):
263 x=ContinuousFunction(self.domain).getX()
264 p=Projector(self.domain,reduce=True,fast=False)
265 td_ref=numpy.array([1.,2.,3.])
266 td=p(Data(td_ref,ReducedFunction(self.domain)))
267 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
268
269 def testProjector_rank2_reduced_with_reduced_input(self):
270 x=ContinuousFunction(self.domain).getX()
271 p=Projector(self.domain,reduce=True,fast=False)
272 td_ref=numpy.array([[11.,12.],[21,22.]])
273 td=p(Data(td_ref,ReducedFunction(self.domain)))
274 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
275
276 def testProjector_rank3_reduced_with_reduced_input(self):
277 x=ContinuousFunction(self.domain).getX()
278 p=Projector(self.domain,reduce=True,fast=False)
279 td_ref=numpy.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
280 td=p(Data(td_ref,ReducedFunction(self.domain)))
281 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
282
283 def testProjector_rank4_reduced_with_reduced_input(self):
284 x=ContinuousFunction(self.domain).getX()
285 p=Projector(self.domain,reduce=True,fast=False)
286 td_ref=numpy.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
287 td=p(Data(td_ref,ReducedFunction(self.domain)))
288 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
289
290
291 def test_NoPDE_scalar_missing_r(self):
292 p=NoPDE(self.domain)
293 x=self.domain.getX()
294 msk=whereZero(x[0])
295 p.setValue(D=1.,Y=1.,q=msk)
296 u=p.getSolution()
297 u_ex=(1.-msk)
298 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
299
300 def test_NoPDE_scalar_missing_Y(self):
301 p=NoPDE(self.domain)
302 x=self.domain.getX()
303 msk=whereZero(x[0])
304 p.setValue(D=1.,q=msk,r=2.)
305 u=p.getSolution()
306 u_ex=msk*2.
307 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
308
309 def test_NoPDE_scalar_constant(self):
310 p=NoPDE(self.domain)
311 x=self.domain.getX()
312 msk=whereZero(x[0])
313 p.setValue(D=1.,Y=1.,q=msk,r=2.)
314 u=p.getSolution()
315 u_ex=(1.-msk)+msk*2.
316 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
317
318 def test_NoPDE_scalar_variable(self):
319 p=NoPDE(self.domain)
320 x=self.domain.getX()
321 msk=whereZero(x[0])
322 p.setValue(D=10,Y=2*x[0],q=msk,r=2.)
323 u=p.getSolution()
324 u_ex=2.
325 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
326
327 def test_NoPDE_vector_missing_Y(self):
328 p=NoPDE(self.domain)
329 x=self.domain.getX()
330 msk=whereZero(x[0])*[1.,0.]
331 p.setValue(D=numpy.ones([2]),q=msk,r=2.)
332 u=p.getSolution()
333 u_ex=msk*2.
334 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
335
336 def test_NoPDE_vector_missing_r(self):
337 p=NoPDE(self.domain)
338 x=self.domain.getX()
339 msk=whereZero(x[0])*[1.,0.]
340 p.setValue(D=numpy.ones([2]),Y=numpy.ones([2]),q=msk)
341 u=p.getSolution()
342 u_ex=(1.-msk)
343 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
344
345 def test_NoPDE_vector_constant(self):
346 p=NoPDE(self.domain)
347 x=self.domain.getX()
348 msk=whereZero(x[0])*[1.,0.]
349 p.setValue(D=numpy.ones([2]),Y=numpy.ones([2]),q=msk,r=2.)
350 u=p.getSolution()
351 u_ex=(1.-msk)+msk*2.
352 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
353
354 def test_NoPDE_vector_variable(self):
355 p=NoPDE(self.domain)
356 x=self.domain.getX()
357 msk=whereZero(x[0])*[1.,0.]
358 p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.)
359 u=p.getSolution()
360 u_ex=2.
361 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
362 #=====
363 def testPCG(self):
364 from numpy import array, dot, zeros, size, float64
365 from math import sqrt
366 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
367 5.834798554135237e-01, -3.704394311709722e+00,
368 5.765369186984777e+00, -1.309786358737351e+01,
369 2.522087134507148e+01, -3.393956279045637e+01,
370 1.046856914770830e+02, -2.447764190849540e+02],
371 [ -2.391895572674098e-01, 1.256797283910693e+02,
372 -9.188270412920813e-01, 1.300169538880688e+00,
373 -5.353714719231424e-01, 2.674709444667012e+00,
374 -1.116097841269580e+01, 2.801193427514478e+01,
375 -3.877806125898224e+01, 3.063505753648256e+01],
376 [ 5.834798554135237e-01, -9.188270412920813e-01,
377 6.240841811806843e+01, -8.176289504109282e-01,
378 1.447935098417076e-01, -9.721424148655324e-01,
379 6.713551574117577e-01, -3.656297654168375e+00,
380 7.015141656913973e+00, -4.195525932156250e+01],
381 [ -3.704394311709722e+00, 1.300169538880688e+00,
382 -8.176289504109282e-01, 3.604980536782198e+01,
383 -6.241238423759328e-01, 1.142345320047869e+00,
384 -3.438816797096519e+00, 5.854857481367470e+00,
385 -4.524311288596452e+00, 1.136590280389803e+01],
386 [ 5.765369186984777e+00, -5.353714719231424e-01,
387 1.447935098417076e-01, -6.241238423759328e-01,
388 2.953997190215862e+01, -9.474729233464712e-01,
389 1.883516378345809e+00, -1.906274765704230e+00,
390 4.401859671778645e+00, -1.064573816075257e+01],
391 [ -1.309786358737351e+01, 2.674709444667012e+00,
392 -9.721424148655324e-01, 1.142345320047869e+00,
393 -9.474729233464712e-01, 2.876998216302979e+01,
394 -4.853065259692995e-01, 7.088596468102618e-01,
395 -8.972224295152829e-01, 5.228606946522749e+00],
396 [ 2.522087134507148e+01, -1.116097841269580e+01,
397 6.713551574117577e-01, -3.438816797096519e+00,
398 1.883516378345809e+00, -4.853065259692995e-01,
399 5.121175860935919e+01, -3.523133115905478e-01,
400 1.782136702229135e+00, -1.560849559916187e+00],
401 [ -3.393956279045637e+01, 2.801193427514478e+01,
402 -3.656297654168375e+00, 5.854857481367470e+00,
403 -1.906274765704230e+00, 7.088596468102618e-01,
404 -3.523133115905478e-01, 8.411681423853814e+01,
405 -5.238590858177903e-01, 1.515872114883926e+00],
406 [ 1.046856914770830e+02, -3.877806125898224e+01,
407 7.015141656913973e+00, -4.524311288596452e+00,
408 4.401859671778645e+00, -8.972224295152829e-01,
409 1.782136702229135e+00, -5.238590858177903e-01,
410 1.797889693808014e+02, -8.362340479938084e-01],
411 [ -2.447764190849540e+02, 3.063505753648256e+01,
412 -4.195525932156250e+01, 1.136590280389803e+01,
413 -1.064573816075257e+01, 5.228606946522749e+00,
414 -1.560849559916187e+00, 1.515872114883926e+00,
415 -8.362340479938084e-01, 3.833719335346630e+02]])
416 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
417 0.807186823427233, 0.48950999450145, 0.995486532098031,
418 0.351243009576568, 0.704352576819321, 0.850648989740204,
419 0.314596738052894])
420 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
421 30.344553414038817, 15.247917439094513, 24.060664905403492,
422 27.210293789825833, 47.122067744075842, 199.267136417856847,
423 -8.7934289814322 ])
424
425 def Ap(x):
426 return dot(A,x)
427 def Ms(b):
428 out=zeros((b.size,),float64)
429 for i in xrange(size(b)):
430 out[i]=b[i]/A[i,i]
431 return out
432
433 tol=1.e-4
434 x,r,a_norm=PCG(b*1.,Ap,x_ref*0.,Ms,dot, atol=0, rtol=tol, iter_max=12)
435 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
436 self.failUnless(Lsup(r-(b-dot(A,x)))<=Lsup(b)*EPSILON*100.,"wrong solution")
437
438 def testMINRES(self):
439 from numpy import array, dot, zeros, size, float64
440 from math import sqrt
441 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
442 5.834798554135237e-01, -3.704394311709722e+00,
443 5.765369186984777e+00, -1.309786358737351e+01,
444 2.522087134507148e+01, -3.393956279045637e+01,
445 1.046856914770830e+02, -2.447764190849540e+02],
446 [ -2.391895572674098e-01, 1.256797283910693e+02,
447 -9.188270412920813e-01, 1.300169538880688e+00,
448 -5.353714719231424e-01, 2.674709444667012e+00,
449 -1.116097841269580e+01, 2.801193427514478e+01,
450 -3.877806125898224e+01, 3.063505753648256e+01],
451 [ 5.834798554135237e-01, -9.188270412920813e-01,
452 6.240841811806843e+01, -8.176289504109282e-01,
453 1.447935098417076e-01, -9.721424148655324e-01,
454 6.713551574117577e-01, -3.656297654168375e+00,
455 7.015141656913973e+00, -4.195525932156250e+01],
456 [ -3.704394311709722e+00, 1.300169538880688e+00,
457 -8.176289504109282e-01, 3.604980536782198e+01,
458 -6.241238423759328e-01, 1.142345320047869e+00,
459 -3.438816797096519e+00, 5.854857481367470e+00,
460 -4.524311288596452e+00, 1.136590280389803e+01],
461 [ 5.765369186984777e+00, -5.353714719231424e-01,
462 1.447935098417076e-01, -6.241238423759328e-01,
463 2.953997190215862e+01, -9.474729233464712e-01,
464 1.883516378345809e+00, -1.906274765704230e+00,
465 4.401859671778645e+00, -1.064573816075257e+01],
466 [ -1.309786358737351e+01, 2.674709444667012e+00,
467 -9.721424148655324e-01, 1.142345320047869e+00,
468 -9.474729233464712e-01, 2.876998216302979e+01,
469 -4.853065259692995e-01, 7.088596468102618e-01,
470 -8.972224295152829e-01, 5.228606946522749e+00],
471 [ 2.522087134507148e+01, -1.116097841269580e+01,
472 6.713551574117577e-01, -3.438816797096519e+00,
473 1.883516378345809e+00, -4.853065259692995e-01,
474 5.121175860935919e+01, -3.523133115905478e-01,
475 1.782136702229135e+00, -1.560849559916187e+00],
476 [ -3.393956279045637e+01, 2.801193427514478e+01,
477 -3.656297654168375e+00, 5.854857481367470e+00,
478 -1.906274765704230e+00, 7.088596468102618e-01,
479 -3.523133115905478e-01, 8.411681423853814e+01,
480 -5.238590858177903e-01, 1.515872114883926e+00],
481 [ 1.046856914770830e+02, -3.877806125898224e+01,
482 7.015141656913973e+00, -4.524311288596452e+00,
483 4.401859671778645e+00, -8.972224295152829e-01,
484 1.782136702229135e+00, -5.238590858177903e-01,
485 1.797889693808014e+02, -8.362340479938084e-01],
486 [ -2.447764190849540e+02, 3.063505753648256e+01,
487 -4.195525932156250e+01, 1.136590280389803e+01,
488 -1.064573816075257e+01, 5.228606946522749e+00,
489 -1.560849559916187e+00, 1.515872114883926e+00,
490 -8.362340479938084e-01, 3.833719335346630e+02]])
491 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
492 0.807186823427233, 0.48950999450145, 0.995486532098031,
493 0.351243009576568, 0.704352576819321, 0.850648989740204,
494 0.314596738052894])
495 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
496 30.344553414038817, 15.247917439094513, 24.060664905403492,
497 27.210293789825833, 47.122067744075842, 199.267136417856847,
498 -8.7934289814322 ])
499
500 def Ap(x):
501 return dot(A,x)
502 def Ms(b):
503 out=zeros((size(b),),float64)
504 for i in xrange(size(b)):
505 out[i]=b[i]/A[i,i]
506 return out
507
508 tol=1.e-4
509 x=MINRES(b*1.,Ap,x_ref*0,Ms,dot, atol=0, rtol=tol, iter_max=12)
510 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
511
512 def testTFQMR(self):
513 from numpy import array, dot, zeros, size, float64
514 from math import sqrt
515 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
516 5.834798554135237e-01, -3.704394311709722e+00,
517 5.765369186984777e+00, -1.309786358737351e+01,
518 2.522087134507148e+01, -3.393956279045637e+01,
519 1.046856914770830e+02, -2.447764190849540e+02],
520 [ -2.391895572674098e-01, 1.256797283910693e+02,
521 -9.188270412920813e-01, 1.300169538880688e+00,
522 -5.353714719231424e-01, 2.674709444667012e+00,
523 -1.116097841269580e+01, 2.801193427514478e+01,
524 -3.877806125898224e+01, 3.063505753648256e+01],
525 [ 5.834798554135237e-01, -9.188270412920813e-01,
526 6.240841811806843e+01, -8.176289504109282e-01,
527 1.447935098417076e-01, -9.721424148655324e-01,
528 6.713551574117577e-01, -3.656297654168375e+00,
529 7.015141656913973e+00, -4.195525932156250e+01],
530 [ -3.704394311709722e+00, 1.300169538880688e+00,
531 -8.176289504109282e-01, 3.604980536782198e+01,
532 -6.241238423759328e-01, 1.142345320047869e+00,
533 -3.438816797096519e+00, 5.854857481367470e+00,
534 -4.524311288596452e+00, 1.136590280389803e+01],
535 [ 5.765369186984777e+00, -5.353714719231424e-01,
536 1.447935098417076e-01, -6.241238423759328e-01,
537 2.953997190215862e+01, -9.474729233464712e-01,
538 1.883516378345809e+00, -1.906274765704230e+00,
539 4.401859671778645e+00, -1.064573816075257e+01],
540 [ -1.309786358737351e+01, 2.674709444667012e+00,
541 -9.721424148655324e-01, 1.142345320047869e+00,
542 -9.474729233464712e-01, 2.876998216302979e+01,
543 -4.853065259692995e-01, 7.088596468102618e-01,
544 -8.972224295152829e-01, 5.228606946522749e+00],
545 [ 2.522087134507148e+01, -1.116097841269580e+01,
546 6.713551574117577e-01, -3.438816797096519e+00,
547 1.883516378345809e+00, -4.853065259692995e-01,
548 5.121175860935919e+01, -3.523133115905478e-01,
549 1.782136702229135e+00, -1.560849559916187e+00],
550 [ -3.393956279045637e+01, 2.801193427514478e+01,
551 -3.656297654168375e+00, 5.854857481367470e+00,
552 -1.906274765704230e+00, 7.088596468102618e-01,
553 -3.523133115905478e-01, 8.411681423853814e+01,
554 -5.238590858177903e-01, 1.515872114883926e+00],
555 [ 1.046856914770830e+02, -3.877806125898224e+01,
556 7.015141656913973e+00, -4.524311288596452e+00,
557 4.401859671778645e+00, -8.972224295152829e-01,
558 1.782136702229135e+00, -5.238590858177903e-01,
559 1.797889693808014e+02, -8.362340479938084e-01],
560 [ -2.447764190849540e+02, 3.063505753648256e+01,
561 -4.195525932156250e+01, 1.136590280389803e+01,
562 -1.064573816075257e+01, 5.228606946522749e+00,
563 -1.560849559916187e+00, 1.515872114883926e+00,
564 -8.362340479938084e-01, 3.833719335346630e+02]])
565 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
566 0.807186823427233, 0.48950999450145, 0.995486532098031,
567 0.351243009576568, 0.704352576819321, 0.850648989740204,
568 0.314596738052894])
569 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
570 30.344553414038817, 15.247917439094513, 24.060664905403492,
571 27.210293789825833, 47.122067744075842, 199.267136417856847,
572 -8.7934289814322 ])
573
574 def Ap(x):
575 out=dot(A,x)
576 for i in xrange(size(x)):
577 out[i]/=A[i,i]
578 return out
579
580 tol=1.e-5
581 for i in xrange(size(b)): b[i]/=A[i,i]
582 x=TFQMR(b,Ap,x_ref*0,dot, atol=0, rtol=tol, iter_max=12)
583 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
584
585 def testGMRES(self):
586 from numpy import array, dot, zeros, size, float64
587 from math import sqrt
588 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
589 5.834798554135237e-01, -3.704394311709722e+00,
590 5.765369186984777e+00, -1.309786358737351e+01,
591 2.522087134507148e+01, -3.393956279045637e+01,
592 1.046856914770830e+02, -2.447764190849540e+02],
593 [ -2.391895572674098e-01, 1.256797283910693e+02,
594 -9.188270412920813e-01, 1.300169538880688e+00,
595 -5.353714719231424e-01, 2.674709444667012e+00,
596 -1.116097841269580e+01, 2.801193427514478e+01,
597 -3.877806125898224e+01, 3.063505753648256e+01],
598 [ 5.834798554135237e-01, -9.188270412920813e-01,
599 6.240841811806843e+01, -8.176289504109282e-01,
600 1.447935098417076e-01, -9.721424148655324e-01,
601 6.713551574117577e-01, -3.656297654168375e+00,
602 7.015141656913973e+00, -4.195525932156250e+01],
603 [ -3.704394311709722e+00, 1.300169538880688e+00,
604 -8.176289504109282e-01, 3.604980536782198e+01,
605 -6.241238423759328e-01, 1.142345320047869e+00,
606 -3.438816797096519e+00, 5.854857481367470e+00,
607 -4.524311288596452e+00, 1.136590280389803e+01],
608 [ 5.765369186984777e+00, -5.353714719231424e-01,
609 1.447935098417076e-01, -6.241238423759328e-01,
610 2.953997190215862e+01, -9.474729233464712e-01,
611 1.883516378345809e+00, -1.906274765704230e+00,
612 4.401859671778645e+00, -1.064573816075257e+01],
613 [ -1.309786358737351e+01, 2.674709444667012e+00,
614 -9.721424148655324e-01, 1.142345320047869e+00,
615 -9.474729233464712e-01, 2.876998216302979e+01,
616 -4.853065259692995e-01, 7.088596468102618e-01,
617 -8.972224295152829e-01, 5.228606946522749e+00],
618 [ 2.522087134507148e+01, -1.116097841269580e+01,
619 6.713551574117577e-01, -3.438816797096519e+00,
620 1.883516378345809e+00, -4.853065259692995e-01,
621 5.121175860935919e+01, -3.523133115905478e-01,
622 1.782136702229135e+00, -1.560849559916187e+00],
623 [ -3.393956279045637e+01, 2.801193427514478e+01,
624 -3.656297654168375e+00, 5.854857481367470e+00,
625 -1.906274765704230e+00, 7.088596468102618e-01,
626 -3.523133115905478e-01, 8.411681423853814e+01,
627 -5.238590858177903e-01, 1.515872114883926e+00],
628 [ 1.046856914770830e+02, -3.877806125898224e+01,
629 7.015141656913973e+00, -4.524311288596452e+00,
630 4.401859671778645e+00, -8.972224295152829e-01,
631 1.782136702229135e+00, -5.238590858177903e-01,
632 1.797889693808014e+02, -8.362340479938084e-01],
633 [ -2.447764190849540e+02, 3.063505753648256e+01,
634 -4.195525932156250e+01, 1.136590280389803e+01,
635 -1.064573816075257e+01, 5.228606946522749e+00,
636 -1.560849559916187e+00, 1.515872114883926e+00,
637 -8.362340479938084e-01, 3.833719335346630e+02]])
638 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
639 0.807186823427233, 0.48950999450145, 0.995486532098031,
640 0.351243009576568, 0.704352576819321, 0.850648989740204,
641 0.314596738052894])
642 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
643 30.344553414038817, 15.247917439094513, 24.060664905403492,
644 27.210293789825833, 47.122067744075842, 199.267136417856847,
645 -8.7934289814322 ])
646
647 def Ap(x):
648 b=dot(A,x)
649 for i in xrange(size(b)):
650 b[i]/=A[i,i]
651 return b
652
653 tol=1.e-4
654 for i in xrange(size(b)): b[i]/=A[i,i]
655 x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12)
656 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
657
658 def testGMRES_P_R(self):
659 from numpy import array, dot, zeros, size, float64
660 from math import sqrt
661 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
662 5.834798554135237e-01, -3.704394311709722e+00,
663 5.765369186984777e+00, -1.309786358737351e+01,
664 2.522087134507148e+01, -3.393956279045637e+01,
665 1.046856914770830e+02, -2.447764190849540e+02],
666 [ -2.391895572674098e-01, 1.256797283910693e+02,
667 -9.188270412920813e-01, 1.300169538880688e+00,
668 -5.353714719231424e-01, 2.674709444667012e+00,
669 -1.116097841269580e+01, 2.801193427514478e+01,
670 -3.877806125898224e+01, 3.063505753648256e+01],
671 [ 5.834798554135237e-01, -9.188270412920813e-01,
672 6.240841811806843e+01, -8.176289504109282e-01,
673 1.447935098417076e-01, -9.721424148655324e-01,
674 6.713551574117577e-01, -3.656297654168375e+00,
675 7.015141656913973e+00, -4.195525932156250e+01],
676 [ -3.704394311709722e+00, 1.300169538880688e+00,
677 -8.176289504109282e-01, 3.604980536782198e+01,
678 -6.241238423759328e-01, 1.142345320047869e+00,
679 -3.438816797096519e+00, 5.854857481367470e+00,
680 -4.524311288596452e+00, 1.136590280389803e+01],
681 [ 5.765369186984777e+00, -5.353714719231424e-01,
682 1.447935098417076e-01, -6.241238423759328e-01,
683 2.953997190215862e+01, -9.474729233464712e-01,
684 1.883516378345809e+00, -1.906274765704230e+00,
685 4.401859671778645e+00, -1.064573816075257e+01],
686 [ -1.309786358737351e+01, 2.674709444667012e+00,
687 -9.721424148655324e-01, 1.142345320047869e+00,
688 -9.474729233464712e-01, 2.876998216302979e+01,
689 -4.853065259692995e-01, 7.088596468102618e-01,
690 -8.972224295152829e-01, 5.228606946522749e+00],
691 [ 2.522087134507148e+01, -1.116097841269580e+01,
692 6.713551574117577e-01, -3.438816797096519e+00,
693 1.883516378345809e+00, -4.853065259692995e-01,
694 5.121175860935919e+01, -3.523133115905478e-01,
695 1.782136702229135e+00, -1.560849559916187e+00],
696 [ -3.393956279045637e+01, 2.801193427514478e+01,
697 -3.656297654168375e+00, 5.854857481367470e+00,
698 -1.906274765704230e+00, 7.088596468102618e-01,
699 -3.523133115905478e-01, 8.411681423853814e+01,
700 -5.238590858177903e-01, 1.515872114883926e+00],
701 [ 1.046856914770830e+02, -3.877806125898224e+01,
702 7.015141656913973e+00, -4.524311288596452e+00,
703 4.401859671778645e+00, -8.972224295152829e-01,
704 1.782136702229135e+00, -5.238590858177903e-01,
705 1.797889693808014e+02, -8.362340479938084e-01],
706 [ -2.447764190849540e+02, 3.063505753648256e+01,
707 -4.195525932156250e+01, 1.136590280389803e+01,
708 -1.064573816075257e+01, 5.228606946522749e+00,
709 -1.560849559916187e+00, 1.515872114883926e+00,
710 -8.362340479938084e-01, 3.833719335346630e+02]])
711 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
712 0.807186823427233, 0.48950999450145, 0.995486532098031,
713 0.351243009576568, 0.704352576819321, 0.850648989740204,
714 0.314596738052894])
715 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
716 30.344553414038817, 15.247917439094513, 24.060664905403492,
717 27.210293789825833, 47.122067744075842, 199.267136417856847,
718 -8.7934289814322 ])
719
720 def Ap(x):
721 return dot(A,x)
722 def P_Rp(x):
723 out=zeros(size(x), float64)
724 for i in xrange(size(x)):
725 out[i]=x[i]/A[i,i]
726 return out
727
728 tol=1.e-4
729 x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12,P_R=P_Rp)
730 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
731
732 def testNewtonGMRES(self):
733 from numpy import array, dot, zeros, size, float64
734 from math import sqrt
735 class LL(Defect):
736 def __init__(self,*kwargs):
737 super(LL, self).__init__(*kwargs)
738 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
739 5.834798554135237e-01, -3.704394311709722e+00,
740 5.765369186984777e+00, -1.309786358737351e+01,
741 2.522087134507148e+01, -3.393956279045637e+01,
742 1.046856914770830e+02, -2.447764190849540e+02],
743 [ -2.391895572674098e-01, 1.256797283910693e+02,
744 -9.188270412920813e-01, 1.300169538880688e+00,
745 -5.353714719231424e-01, 2.674709444667012e+00,
746 -1.116097841269580e+01, 2.801193427514478e+01,
747 -3.877806125898224e+01, 3.063505753648256e+01],
748 [ 5.834798554135237e-01, -9.188270412920813e-01,
749 6.240841811806843e+01, -8.176289504109282e-01,
750 1.447935098417076e-01, -9.721424148655324e-01,
751 6.713551574117577e-01, -3.656297654168375e+00,
752 7.015141656913973e+00, -4.195525932156250e+01],
753 [ -3.704394311709722e+00, 1.300169538880688e+00,
754 -8.176289504109282e-01, 3.604980536782198e+01,
755 -6.241238423759328e-01, 1.142345320047869e+00,
756 -3.438816797096519e+00, 5.854857481367470e+00,
757 -4.524311288596452e+00, 1.136590280389803e+01],
758 [ 5.765369186984777e+00, -5.353714719231424e-01,
759 1.447935098417076e-01, -6.241238423759328e-01,
760 2.953997190215862e+01, -9.474729233464712e-01,
761 1.883516378345809e+00, -1.906274765704230e+00,
762 4.401859671778645e+00, -1.064573816075257e+01],
763 [ -1.309786358737351e+01, 2.674709444667012e+00,
764 -9.721424148655324e-01, 1.142345320047869e+00,
765 -9.474729233464712e-01, 2.876998216302979e+01,
766 -4.853065259692995e-01, 7.088596468102618e-01,
767 -8.972224295152829e-01, 5.228606946522749e+00],
768 [ 2.522087134507148e+01, -1.116097841269580e+01,
769 6.713551574117577e-01, -3.438816797096519e+00,
770 1.883516378345809e+00, -4.853065259692995e-01,
771 5.121175860935919e+01, -3.523133115905478e-01,
772 1.782136702229135e+00, -1.560849559916187e+00],
773 [ -3.393956279045637e+01, 2.801193427514478e+01,
774 -3.656297654168375e+00, 5.854857481367470e+00,
775 -1.906274765704230e+00, 7.088596468102618e-01,
776 -3.523133115905478e-01, 8.411681423853814e+01,
777 -5.238590858177903e-01, 1.515872114883926e+00],
778 [ 1.046856914770830e+02, -3.877806125898224e+01,
779 7.015141656913973e+00, -4.524311288596452e+00,
780 4.401859671778645e+00, -8.972224295152829e-01,
781 1.782136702229135e+00, -5.238590858177903e-01,
782 1.797889693808014e+02, -8.362340479938084e-01],
783 [ -2.447764190849540e+02, 3.063505753648256e+01,
784 -4.195525932156250e+01, 1.136590280389803e+01,
785 -1.064573816075257e+01, 5.228606946522749e+00,
786 -1.560849559916187e+00, 1.515872114883926e+00,
787 -8.362340479938084e-01, 3.833719335346630e+02]])
788 self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
789 0.807186823427233, 0.48950999450145, 0.995486532098031,
790 0.351243009576568, 0.704352576819321, 0.850648989740204,
791 0.314596738052894])
792 self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
793 30.344553414038817, 15.247917439094513, 24.060664905403492,
794 27.210293789825833, 47.122067744075842, 199.267136417856847,
795 -8.7934289814322 ])
796 def eval(self,x):
797 out=dot(self.A,x)-self.b
798 for i in xrange(size(self.b)):
799 out[i]/=self.A[i,i]
800 return out
801 def bilinearform(self,x0,x1):
802 return dot(x0,x1)
803
804 tol=1.e-8
805 ll=LL()
806 x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
807 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
808
809 def testNewtonGMRES(self):
810 from numpy import array, dot, zeros, size, float64
811 from math import sqrt
812 class LL(Defect):
813 def __init__(self,*kwargs):
814 super(LL, self).__init__(*kwargs)
815 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
816 5.834798554135237e-01, -3.704394311709722e+00,
817 5.765369186984777e+00, -1.309786358737351e+01,
818 2.522087134507148e+01, -3.393956279045637e+01,
819 1.046856914770830e+02, -2.447764190849540e+02],
820 [ -2.391895572674098e-01, 1.256797283910693e+02,
821 -9.188270412920813e-01, 1.300169538880688e+00,
822 -5.353714719231424e-01, 2.674709444667012e+00,
823 -1.116097841269580e+01, 2.801193427514478e+01,
824 -3.877806125898224e+01, 3.063505753648256e+01],
825 [ 5.834798554135237e-01, -9.188270412920813e-01,
826 6.240841811806843e+01, -8.176289504109282e-01,
827 1.447935098417076e-01, -9.721424148655324e-01,
828 6.713551574117577e-01, -3.656297654168375e+00,
829 7.015141656913973e+00, -4.195525932156250e+01],
830 [ -3.704394311709722e+00, 1.300169538880688e+00,
831 -8.176289504109282e-01, 3.604980536782198e+01,
832 -6.241238423759328e-01, 1.142345320047869e+00,
833 -3.438816797096519e+00, 5.854857481367470e+00,
834 -4.524311288596452e+00, 1.136590280389803e+01],
835 [ 5.765369186984777e+00, -5.353714719231424e-01,
836 1.447935098417076e-01, -6.241238423759328e-01,
837 2.953997190215862e+01, -9.474729233464712e-01,
838 1.883516378345809e+00, -1.906274765704230e+00,
839 4.401859671778645e+00, -1.064573816075257e+01],
840 [ -1.309786358737351e+01, 2.674709444667012e+00,
841 -9.721424148655324e-01, 1.142345320047869e+00,
842 -9.474729233464712e-01, 2.876998216302979e+01,
843 -4.853065259692995e-01, 7.088596468102618e-01,
844 -8.972224295152829e-01, 5.228606946522749e+00],
845 [ 2.522087134507148e+01, -1.116097841269580e+01,
846 6.713551574117577e-01, -3.438816797096519e+00,
847 1.883516378345809e+00, -4.853065259692995e-01,
848 5.121175860935919e+01, -3.523133115905478e-01,
849 1.782136702229135e+00, -1.560849559916187e+00],
850 [ -3.393956279045637e+01, 2.801193427514478e+01,
851 -3.656297654168375e+00, 5.854857481367470e+00,
852 -1.906274765704230e+00, 7.088596468102618e-01,
853 -3.523133115905478e-01, 8.411681423853814e+01,
854 -5.238590858177903e-01, 1.515872114883926e+00],
855 [ 1.046856914770830e+02, -3.877806125898224e+01,
856 7.015141656913973e+00, -4.524311288596452e+00,
857 4.401859671778645e+00, -8.972224295152829e-01,
858 1.782136702229135e+00, -5.238590858177903e-01,
859 1.797889693808014e+02, -8.362340479938084e-01],
860 [ -2.447764190849540e+02, 3.063505753648256e+01,
861 -4.195525932156250e+01, 1.136590280389803e+01,
862 -1.064573816075257e+01, 5.228606946522749e+00,
863 -1.560849559916187e+00, 1.515872114883926e+00,
864 -8.362340479938084e-01, 3.833719335346630e+02]])
865 self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
866 0.807186823427233, 0.48950999450145, 0.995486532098031,
867 0.351243009576568, 0.704352576819321, 0.850648989740204,
868 0.314596738052894])
869 self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
870 30.344553414038817, 15.247917439094513, 24.060664905403492,
871 27.210293789825833, 47.122067744075842, 199.267136417856847,
872 -8.7934289814322 ])
873 def eval(self,x):
874 out=dot(self.A,x)-self.b
875 for i in xrange(size(self.b)):
876 out[i]/=self.A[i,i]
877 return out
878 def bilinearform(self,x0,x1):
879 return dot(x0,x1)
880
881 tol=1.e-8
882 ll=LL()
883 x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
884 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
885
886 def testHomogeneousSaddlePointProblem_PCG(self):
887 from numpy import array, dot, zeros, size, float64
888 from math import sqrt
889 class LL(HomogeneousSaddlePointProblem):
890 def initialize(self):
891 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
892 5.834798554135237e-01, -3.704394311709722e+00,
893 5.765369186984777e+00, -1.309786358737351e+01,
894 2.522087134507148e+01, -3.393956279045637e+01,
895 1.046856914770830e+02, -2.447764190849540e+02],
896 [ -2.391895572674098e-01, 1.256797283910693e+02,
897 -9.188270412920813e-01, 1.300169538880688e+00,
898 -5.353714719231424e-01, 2.674709444667012e+00,
899 -1.116097841269580e+01, 2.801193427514478e+01,
900 -3.877806125898224e+01, 3.063505753648256e+01],
901 [ 5.834798554135237e-01, -9.188270412920813e-01,
902 6.240841811806843e+01, -8.176289504109282e-01,
903 1.447935098417076e-01, -9.721424148655324e-01,
904 6.713551574117577e-01, -3.656297654168375e+00,
905 7.015141656913973e+00, -4.195525932156250e+01],
906 [ -3.704394311709722e+00, 1.300169538880688e+00,
907 -8.176289504109282e-01, 3.604980536782198e+01,
908 -6.241238423759328e-01, 1.142345320047869e+00,
909 -3.438816797096519e+00, 5.854857481367470e+00,
910 -4.524311288596452e+00, 1.136590280389803e+01],
911 [ 5.765369186984777e+00, -5.353714719231424e-01,
912 1.447935098417076e-01, -6.241238423759328e-01,
913 2.953997190215862e+01, -9.474729233464712e-01,
914 1.883516378345809e+00, -1.906274765704230e+00,
915 4.401859671778645e+00, -1.064573816075257e+01],
916 [ -1.309786358737351e+01, 2.674709444667012e+00,
917 -9.721424148655324e-01, 1.142345320047869e+00,
918 -9.474729233464712e-01, 2.876998216302979e+01,
919 -4.853065259692995e-01, 7.088596468102618e-01,
920 -8.972224295152829e-01, 5.228606946522749e+00],
921 [ 2.522087134507148e+01, -1.116097841269580e+01,
922 6.713551574117577e-01, -3.438816797096519e+00,
923 1.883516378345809e+00, -4.853065259692995e-01,
924 5.121175860935919e+01, -3.523133115905478e-01,
925 1.782136702229135e+00, -1.560849559916187e+00],
926 [ -3.393956279045637e+01, 2.801193427514478e+01,
927 -3.656297654168375e+00, 5.854857481367470e+00,
928 -1.906274765704230e+00, 7.088596468102618e-01,
929 -3.523133115905478e-01, 8.411681423853814e+01,
930 -5.238590858177903e-01, 1.515872114883926e+00],
931 [ 1.046856914770830e+02, -3.877806125898224e+01,
932 7.015141656913973e+00, -4.524311288596452e+00,
933 4.401859671778645e+00, -8.972224295152829e-01,
934 1.782136702229135e+00, -5.238590858177903e-01,
935 1.797889693808014e+02, -8.362340479938084e-01],
936 [ -2.447764190849540e+02, 3.063505753648256e+01,
937 -4.195525932156250e+01, 1.136590280389803e+01,
938 -1.064573816075257e+01, 5.228606946522749e+00,
939 -1.560849559916187e+00, 1.515872114883926e+00,
940 -8.362340479938084e-01, 3.833719335346630e+02]])
941 self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
942 0.253012436539738, -0.346223308561905, 0.2425508275422,
943 -0.194695862196008, 0.09451439391473, 0.302961126826511,
944 -0.236043777597633] )
945
946 self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
947 [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
948 [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
949 [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
950 [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
951 [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
952 [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
953 [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
954 [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
955 [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
956 self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
957
958 self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
959 6.512480714694167 , -5.727371407390975 , 4.802494840775022,
960 -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
961 -118.602464657076439])
962
963 self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
964 [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
965 [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
966 def inner_pBv(self,p,Bv):
967 return dot(p,Bv)
968 def Bv(self,v, tol):
969 return dot(transpose(self.Bt),v)
970 def inner_p(self,p0,p1):
971 return dot(p0,p1)
972 def norm_v(self,v):
973 return sqrt(dot(v,v))
974 def getDV(self,p,v, tol):
975 dv=solve_linear_equations(self.A, self.b-dot(self.Bt,p)-dot(self.A,v))
976 return dv*(1+tol)
977 def norm_Bv(self,Bv):
978 return sqrt(dot(Bv,Bv))
979 def solve_AinvBt(self,p, tol):
980 out=solve_linear_equations(self.A, dot(self.Bt,p))
981 return out*(1.+tol)
982 def solve_prec(self,Bv, tol):
983 out=Bv*1.
984 for i in xrange(size(out)): out[i]*=self.Sinv[i,i]
985 return out*(1-tol)
986
987 tol=1.e-8
988 ll=LL()
989 ll.initialize()
990 ll.setTolerance(tol)
991 # ll.setSubToleranceReductionFactor(0.1)
992 x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, usePCG=True, iter_restart=20,max_correction_steps=10)
993 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
994 self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
995
996 def testHomogeneousSaddlePointProblem_GMRES(self):
997 from numpy import array, prod, dot, zeros, size, float64
998 from math import sqrt
999 class LL(HomogeneousSaddlePointProblem):
1000 def initialize(self):
1001 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
1002 5.834798554135237e-01, -3.704394311709722e+00,
1003 5.765369186984777e+00, -1.309786358737351e+01,
1004 2.522087134507148e+01, -3.393956279045637e+01,
1005 1.046856914770830e+02, -2.447764190849540e+02],
1006 [ -2.391895572674098e-01, 1.256797283910693e+02,
1007 -9.188270412920813e-01, 1.300169538880688e+00,
1008 -5.353714719231424e-01, 2.674709444667012e+00,
1009 -1.116097841269580e+01, 2.801193427514478e+01,
1010 -3.877806125898224e+01, 3.063505753648256e+01],
1011 [ 5.834798554135237e-01, -9.188270412920813e-01,
1012 6.240841811806843e+01, -8.176289504109282e-01,
1013 1.447935098417076e-01, -9.721424148655324e-01,
1014 6.713551574117577e-01, -3.656297654168375e+00,
1015 7.015141656913973e+00, -4.195525932156250e+01],
1016 [ -3.704394311709722e+00, 1.300169538880688e+00,
1017 -8.176289504109282e-01, 3.604980536782198e+01,
1018 -6.241238423759328e-01, 1.142345320047869e+00,
1019 -3.438816797096519e+00, 5.854857481367470e+00,
1020 -4.524311288596452e+00, 1.136590280389803e+01],
1021 [ 5.765369186984777e+00, -5.353714719231424e-01,
1022 1.447935098417076e-01, -6.241238423759328e-01,
1023 2.953997190215862e+01, -9.474729233464712e-01,
1024 1.883516378345809e+00, -1.906274765704230e+00,
1025 4.401859671778645e+00, -1.064573816075257e+01],
1026 [ -1.309786358737351e+01, 2.674709444667012e+00,
1027 -9.721424148655324e-01, 1.142345320047869e+00,
1028 -9.474729233464712e-01, 2.876998216302979e+01,
1029 -4.853065259692995e-01, 7.088596468102618e-01,
1030 -8.972224295152829e-01, 5.228606946522749e+00],
1031 [ 2.522087134507148e+01, -1.116097841269580e+01,
1032 6.713551574117577e-01, -3.438816797096519e+00,
1033 1.883516378345809e+00, -4.853065259692995e-01,
1034 5.121175860935919e+01, -3.523133115905478e-01,
1035 1.782136702229135e+00, -1.560849559916187e+00],
1036 [ -3.393956279045637e+01, 2.801193427514478e+01,
1037 -3.656297654168375e+00, 5.854857481367470e+00,
1038 -1.906274765704230e+00, 7.088596468102618e-01,
1039 -3.523133115905478e-01, 8.411681423853814e+01,
1040 -5.238590858177903e-01, 1.515872114883926e+00],
1041 [ 1.046856914770830e+02, -3.877806125898224e+01,
1042 7.015141656913973e+00, -4.524311288596452e+00,
1043 4.401859671778645e+00, -8.972224295152829e-01,
1044 1.782136702229135e+00, -5.238590858177903e-01,
1045 1.797889693808014e+02, -8.362340479938084e-01],
1046 [ -2.447764190849540e+02, 3.063505753648256e+01,
1047 -4.195525932156250e+01, 1.136590280389803e+01,
1048 -1.064573816075257e+01, 5.228606946522749e+00,
1049 -1.560849559916187e+00, 1.515872114883926e+00,
1050 -8.362340479938084e-01, 3.833719335346630e+02]])
1051 self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
1052 0.253012436539738, -0.346223308561905, 0.2425508275422,
1053 -0.194695862196008, 0.09451439391473, 0.302961126826511,
1054 -0.236043777597633] )
1055
1056 self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
1057 [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
1058 [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
1059 [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
1060 [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
1061 [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
1062 [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
1063 [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
1064 [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
1065 [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
1066 self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
1067
1068 self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
1069 6.512480714694167 , -5.727371407390975 , 4.802494840775022,
1070 -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
1071 -118.602464657076439])
1072
1073 self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
1074 [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
1075 [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
1076 def inner_pBv(self,p,Bv):
1077 return dot(p,Bv)
1078 def Bv(self,v, tol):
1079 return dot(transpose(self.Bt),v)
1080 def inner_p(self,p0,p1):
1081 return dot(p0,p1)
1082 def norm_v(self,v):
1083 return sqrt(dot(v,v))
1084 def getDV(self,p,v, tol):
1085 dv=solve_linear_equations(self.A, self.b-dot(self.Bt,p)-dot(self.A,v))
1086 return dv*(1+tol)
1087 def norm_Bv(self,Bv):
1088 return sqrt(dot(Bv,Bv))
1089 def solve_AinvBt(self,p, tol):
1090 out=solve_linear_equations(self.A, dot(self.Bt,p))
1091 return out*(1.+tol)
1092 def solve_prec(self,Bv, tol):
1093 out=Bv*1.
1094 for i in xrange(size(out)): out[i]*=self.Sinv[i,i]
1095 return out*(1-tol)
1096
1097 tol=1.e-8
1098 ll=LL()
1099 ll.initialize()
1100 ll.setTolerance(tol)
1101 # ll.setSubToleranceReductionFactor(0.1)
1102 x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, usePCG=False, iter_restart=20,max_correction_steps=10)
1103 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
1104 self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
1105
1106 def testArithmeticTuple(self):
1107 a=ArithmeticTuple(1.,2.)
1108 self.failUnless(len(a)==2,"wrong length")
1109 self.failUnless(a[0]==1.,"wrong first item")
1110 self.failUnless(a[1]==2.,"wrong second item")
1111 c=a*6.
1112 self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple")
1113 self.failUnless(len(c)==2,"c has wrong length")
1114 self.failUnless(c[0]==6.,"c has wrong first item")
1115 self.failUnless(c[1]==12.,"c has wrong second item")
1116 b=5.*a
1117 self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple")
1118 self.failUnless(len(b)==2,"b has wrong length")
1119 self.failUnless(b[0]==5.,"b has wrong first item")
1120 self.failUnless(b[1]==10.,"b has wrong second item")
1121 a+=ArithmeticTuple(3.,4.)
1122 self.failUnless(a[0]==4.,"wrong first item of inplace update")
1123 self.failUnless(a[1]==6.,"wrong second item of inplace update")
1124
1125
1126
1127 class Test_pdetools(Test_pdetools_noLumping):
1128 def testProjector_rank0_fast_reduced(self):
1129 x=ContinuousFunction(self.domain).getX()
1130 h=Lsup(self.domain.getSize())
1131 p=Projector(self.domain,reduce=True,fast=True)
1132 td_ref=x[0]
1133 td=p(td_ref.interpolate(Function(self.domain)))
1134 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1135
1136 def testProjector_rank1_fast_reduced(self):
1137 x=ContinuousFunction(self.domain).getX()
1138 h=Lsup(self.domain.getSize())
1139 p=Projector(self.domain,reduce=True,fast=True)
1140 td_ref=x
1141 td=p(td_ref.interpolate(Function(self.domain)))
1142 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1143
1144 def testProjector_rank2_fast_reduced(self):
1145 x=ContinuousFunction(self.domain).getX()
1146 h=Lsup(self.domain.getSize())
1147 p=Projector(self.domain,reduce=True,fast=True)
1148 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
1149 td=p(td_ref.interpolate(Function(self.domain)))
1150 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1151
1152 def testProjector_rank3_fast_reduced(self):
1153 x=ContinuousFunction(self.domain).getX()
1154 h=Lsup(self.domain.getSize())
1155 p=Projector(self.domain,reduce=True,fast=True)
1156 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
1157 td=p(td_ref.interpolate(Function(self.domain)))
1158 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1159
1160 def testProjector_rank4_fast_reduced(self):
1161 x=ContinuousFunction(self.domain).getX()
1162 h=Lsup(self.domain.getSize())
1163 p=Projector(self.domain,reduce=True,fast=True)
1164 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
1165 td=p(td_ref.interpolate(Function(self.domain)))
1166 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1167
1168 def testProjector_rank0_fast_reduced_with_reduced_input(self):
1169 x=ContinuousFunction(self.domain).getX()
1170 h=Lsup(self.domain.getSize())
1171 p=Projector(self.domain,reduce=True,fast=True)
1172 td_ref=1.
1173 td=p(Data(td_ref,ReducedFunction(self.domain)))
1174 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1175
1176 def testProjector_rank1_fast_reduced_with_reduced_input(self):
1177 x=ContinuousFunction(self.domain).getX()
1178 h=Lsup(self.domain.getSize())
1179 p=Projector(self.domain,reduce=True,fast=True)
1180 td_ref=numpy.array([1.,2.,3.])
1181 td=p(Data(td_ref,ReducedFunction(self.domain)))
1182 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1183
1184 def testProjector_rank2_fast_reduced_with_reduced_input(self):
1185 x=ContinuousFunction(self.domain).getX()
1186 h=Lsup(self.domain.getSize())
1187 p=Projector(self.domain,reduce=True,fast=True)
1188 td_ref=numpy.array([[11.,12.],[21,22.]])
1189 td=p(Data(td_ref,ReducedFunction(self.domain)))
1190 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1191
1192 def testProjector_rank3_fast_reduced_with_reduced_input(self):
1193 x=ContinuousFunction(self.domain).getX()
1194 h=Lsup(self.domain.getSize())
1195 p=Projector(self.domain,reduce=True,fast=True)
1196 td_ref=numpy.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
1197 td=p(Data(td_ref,ReducedFunction(self.domain)))
1198 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1199
1200 def testProjector_rank4_fast_reduced_with_reduced_input(self):
1201 x=ContinuousFunction(self.domain).getX()
1202 h=Lsup(self.domain.getSize())
1203 p=Projector(self.domain,reduce=True,fast=True)
1204 td_ref=numpy.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
1205 td=p(Data(td_ref,ReducedFunction(self.domain)))
1206 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1207

  ViewVC Help
Powered by ViewVC 1.1.26