/[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 2156 - (show annotations)
Mon Dec 15 05:09:02 2008 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 61262 byte(s)
some modifications to the iterative solver to make them easier to use. 
There are also improved versions of the Darcy flux solver and the incompressible solver.


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

  ViewVC Help
Powered by ViewVC 1.1.26