/[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 3892 - (show annotations)
Tue Apr 10 08:57:23 2012 UTC (6 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 67375 byte(s)
Merged changes across from the attempt2 branch.
This version builds and passes python2 tests.
It also passes most python3 tests.



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

  ViewVC Help
Powered by ViewVC 1.1.26