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

  ViewVC Help
Powered by ViewVC 1.1.26