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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4938 - (show annotations)
Wed May 14 01:13:23 2014 UTC (5 years, 3 months ago) by jfenwick
File MIME type: text/x-python
File size: 67560 byte(s)
Modify unit tests to read their classes from
esys.escriptcore.utestselect

Change the line in that file to switch between unittest and unittest2


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

  ViewVC Help
Powered by ViewVC 1.1.26