/[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 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 61295 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



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,Ms,dot, atol=0, rtol=tol,x=x_ref*1.5, 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,Ms,dot, atol=0, rtol=tol,x=x_ref*1.5, 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 return matrixmultiply(A,x)
549 def Ms(b):
550 out=zeros((size(b),),Float64)
551 for i in xrange(size(b)):
552 out[i]=b[i]/A[i,i]
553 return out
554
555 tol=1.e-5
556 x=TFQMR(b*1.,Ap,Ms,dot, atol=0, rtol=tol,x=x_ref*1.5, iter_max=12)
557 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
558
559 def testGMRES(self):
560 from numarray import array,matrixmultiply, zeros, dot, size, Float64
561 from math import sqrt
562 A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
563 5.834798554135237e-01, -3.704394311709722e+00,
564 5.765369186984777e+00, -1.309786358737351e+01,
565 2.522087134507148e+01, -3.393956279045637e+01,
566 1.046856914770830e+02, -2.447764190849540e+02],
567 [ -2.391895572674098e-01, 1.256797283910693e+02,
568 -9.188270412920813e-01, 1.300169538880688e+00,
569 -5.353714719231424e-01, 2.674709444667012e+00,
570 -1.116097841269580e+01, 2.801193427514478e+01,
571 -3.877806125898224e+01, 3.063505753648256e+01],
572 [ 5.834798554135237e-01, -9.188270412920813e-01,
573 6.240841811806843e+01, -8.176289504109282e-01,
574 1.447935098417076e-01, -9.721424148655324e-01,
575 6.713551574117577e-01, -3.656297654168375e+00,
576 7.015141656913973e+00, -4.195525932156250e+01],
577 [ -3.704394311709722e+00, 1.300169538880688e+00,
578 -8.176289504109282e-01, 3.604980536782198e+01,
579 -6.241238423759328e-01, 1.142345320047869e+00,
580 -3.438816797096519e+00, 5.854857481367470e+00,
581 -4.524311288596452e+00, 1.136590280389803e+01],
582 [ 5.765369186984777e+00, -5.353714719231424e-01,
583 1.447935098417076e-01, -6.241238423759328e-01,
584 2.953997190215862e+01, -9.474729233464712e-01,
585 1.883516378345809e+00, -1.906274765704230e+00,
586 4.401859671778645e+00, -1.064573816075257e+01],
587 [ -1.309786358737351e+01, 2.674709444667012e+00,
588 -9.721424148655324e-01, 1.142345320047869e+00,
589 -9.474729233464712e-01, 2.876998216302979e+01,
590 -4.853065259692995e-01, 7.088596468102618e-01,
591 -8.972224295152829e-01, 5.228606946522749e+00],
592 [ 2.522087134507148e+01, -1.116097841269580e+01,
593 6.713551574117577e-01, -3.438816797096519e+00,
594 1.883516378345809e+00, -4.853065259692995e-01,
595 5.121175860935919e+01, -3.523133115905478e-01,
596 1.782136702229135e+00, -1.560849559916187e+00],
597 [ -3.393956279045637e+01, 2.801193427514478e+01,
598 -3.656297654168375e+00, 5.854857481367470e+00,
599 -1.906274765704230e+00, 7.088596468102618e-01,
600 -3.523133115905478e-01, 8.411681423853814e+01,
601 -5.238590858177903e-01, 1.515872114883926e+00],
602 [ 1.046856914770830e+02, -3.877806125898224e+01,
603 7.015141656913973e+00, -4.524311288596452e+00,
604 4.401859671778645e+00, -8.972224295152829e-01,
605 1.782136702229135e+00, -5.238590858177903e-01,
606 1.797889693808014e+02, -8.362340479938084e-01],
607 [ -2.447764190849540e+02, 3.063505753648256e+01,
608 -4.195525932156250e+01, 1.136590280389803e+01,
609 -1.064573816075257e+01, 5.228606946522749e+00,
610 -1.560849559916187e+00, 1.515872114883926e+00,
611 -8.362340479938084e-01, 3.833719335346630e+02]])
612 x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
613 0.807186823427233, 0.48950999450145, 0.995486532098031,
614 0.351243009576568, 0.704352576819321, 0.850648989740204,
615 0.314596738052894])
616 b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
617 30.344553414038817, 15.247917439094513, 24.060664905403492,
618 27.210293789825833, 47.122067744075842, 199.267136417856847,
619 -8.7934289814322 ])
620
621 def Ap(x):
622 b=matrixmultiply(A,x)
623 for i in xrange(size(b)):
624 b[i]/=A[i,i]
625 return b
626
627 tol=1.e-4
628 for i in xrange(size(b)): b[i]/=A[i,i]
629 x=GMRES(b*1.,Ap,dot,atol=0, rtol=tol,x=x_ref*1.5, iter_max=12)
630 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
631
632 def testNewtonGMRES(self):
633 from numarray import array,matrixmultiply, zeros, dot, size, Float64
634 from math import sqrt
635 class LL(Defect):
636 def __init__(self,*kwargs):
637 super(LL, self).__init__(*kwargs)
638 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
639 5.834798554135237e-01, -3.704394311709722e+00,
640 5.765369186984777e+00, -1.309786358737351e+01,
641 2.522087134507148e+01, -3.393956279045637e+01,
642 1.046856914770830e+02, -2.447764190849540e+02],
643 [ -2.391895572674098e-01, 1.256797283910693e+02,
644 -9.188270412920813e-01, 1.300169538880688e+00,
645 -5.353714719231424e-01, 2.674709444667012e+00,
646 -1.116097841269580e+01, 2.801193427514478e+01,
647 -3.877806125898224e+01, 3.063505753648256e+01],
648 [ 5.834798554135237e-01, -9.188270412920813e-01,
649 6.240841811806843e+01, -8.176289504109282e-01,
650 1.447935098417076e-01, -9.721424148655324e-01,
651 6.713551574117577e-01, -3.656297654168375e+00,
652 7.015141656913973e+00, -4.195525932156250e+01],
653 [ -3.704394311709722e+00, 1.300169538880688e+00,
654 -8.176289504109282e-01, 3.604980536782198e+01,
655 -6.241238423759328e-01, 1.142345320047869e+00,
656 -3.438816797096519e+00, 5.854857481367470e+00,
657 -4.524311288596452e+00, 1.136590280389803e+01],
658 [ 5.765369186984777e+00, -5.353714719231424e-01,
659 1.447935098417076e-01, -6.241238423759328e-01,
660 2.953997190215862e+01, -9.474729233464712e-01,
661 1.883516378345809e+00, -1.906274765704230e+00,
662 4.401859671778645e+00, -1.064573816075257e+01],
663 [ -1.309786358737351e+01, 2.674709444667012e+00,
664 -9.721424148655324e-01, 1.142345320047869e+00,
665 -9.474729233464712e-01, 2.876998216302979e+01,
666 -4.853065259692995e-01, 7.088596468102618e-01,
667 -8.972224295152829e-01, 5.228606946522749e+00],
668 [ 2.522087134507148e+01, -1.116097841269580e+01,
669 6.713551574117577e-01, -3.438816797096519e+00,
670 1.883516378345809e+00, -4.853065259692995e-01,
671 5.121175860935919e+01, -3.523133115905478e-01,
672 1.782136702229135e+00, -1.560849559916187e+00],
673 [ -3.393956279045637e+01, 2.801193427514478e+01,
674 -3.656297654168375e+00, 5.854857481367470e+00,
675 -1.906274765704230e+00, 7.088596468102618e-01,
676 -3.523133115905478e-01, 8.411681423853814e+01,
677 -5.238590858177903e-01, 1.515872114883926e+00],
678 [ 1.046856914770830e+02, -3.877806125898224e+01,
679 7.015141656913973e+00, -4.524311288596452e+00,
680 4.401859671778645e+00, -8.972224295152829e-01,
681 1.782136702229135e+00, -5.238590858177903e-01,
682 1.797889693808014e+02, -8.362340479938084e-01],
683 [ -2.447764190849540e+02, 3.063505753648256e+01,
684 -4.195525932156250e+01, 1.136590280389803e+01,
685 -1.064573816075257e+01, 5.228606946522749e+00,
686 -1.560849559916187e+00, 1.515872114883926e+00,
687 -8.362340479938084e-01, 3.833719335346630e+02]])
688 self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
689 0.807186823427233, 0.48950999450145, 0.995486532098031,
690 0.351243009576568, 0.704352576819321, 0.850648989740204,
691 0.314596738052894])
692 self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
693 30.344553414038817, 15.247917439094513, 24.060664905403492,
694 27.210293789825833, 47.122067744075842, 199.267136417856847,
695 -8.7934289814322 ])
696 def eval(self,x):
697 out=matrixmultiply(self.A,x)-self.b
698 for i in xrange(size(self.b)):
699 out[i]/=self.A[i,i]
700 return out
701 def bilinearform(self,x0,x1):
702 return dot(x0,x1)
703
704 tol=1.e-8
705 ll=LL()
706 x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
707 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
708
709 def testNewtonGMRES(self):
710 from numarray import array,matrixmultiply, zeros, dot, size, Float64
711 from math import sqrt
712 class LL(Defect):
713 def __init__(self,*kwargs):
714 super(LL, self).__init__(*kwargs)
715 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
716 5.834798554135237e-01, -3.704394311709722e+00,
717 5.765369186984777e+00, -1.309786358737351e+01,
718 2.522087134507148e+01, -3.393956279045637e+01,
719 1.046856914770830e+02, -2.447764190849540e+02],
720 [ -2.391895572674098e-01, 1.256797283910693e+02,
721 -9.188270412920813e-01, 1.300169538880688e+00,
722 -5.353714719231424e-01, 2.674709444667012e+00,
723 -1.116097841269580e+01, 2.801193427514478e+01,
724 -3.877806125898224e+01, 3.063505753648256e+01],
725 [ 5.834798554135237e-01, -9.188270412920813e-01,
726 6.240841811806843e+01, -8.176289504109282e-01,
727 1.447935098417076e-01, -9.721424148655324e-01,
728 6.713551574117577e-01, -3.656297654168375e+00,
729 7.015141656913973e+00, -4.195525932156250e+01],
730 [ -3.704394311709722e+00, 1.300169538880688e+00,
731 -8.176289504109282e-01, 3.604980536782198e+01,
732 -6.241238423759328e-01, 1.142345320047869e+00,
733 -3.438816797096519e+00, 5.854857481367470e+00,
734 -4.524311288596452e+00, 1.136590280389803e+01],
735 [ 5.765369186984777e+00, -5.353714719231424e-01,
736 1.447935098417076e-01, -6.241238423759328e-01,
737 2.953997190215862e+01, -9.474729233464712e-01,
738 1.883516378345809e+00, -1.906274765704230e+00,
739 4.401859671778645e+00, -1.064573816075257e+01],
740 [ -1.309786358737351e+01, 2.674709444667012e+00,
741 -9.721424148655324e-01, 1.142345320047869e+00,
742 -9.474729233464712e-01, 2.876998216302979e+01,
743 -4.853065259692995e-01, 7.088596468102618e-01,
744 -8.972224295152829e-01, 5.228606946522749e+00],
745 [ 2.522087134507148e+01, -1.116097841269580e+01,
746 6.713551574117577e-01, -3.438816797096519e+00,
747 1.883516378345809e+00, -4.853065259692995e-01,
748 5.121175860935919e+01, -3.523133115905478e-01,
749 1.782136702229135e+00, -1.560849559916187e+00],
750 [ -3.393956279045637e+01, 2.801193427514478e+01,
751 -3.656297654168375e+00, 5.854857481367470e+00,
752 -1.906274765704230e+00, 7.088596468102618e-01,
753 -3.523133115905478e-01, 8.411681423853814e+01,
754 -5.238590858177903e-01, 1.515872114883926e+00],
755 [ 1.046856914770830e+02, -3.877806125898224e+01,
756 7.015141656913973e+00, -4.524311288596452e+00,
757 4.401859671778645e+00, -8.972224295152829e-01,
758 1.782136702229135e+00, -5.238590858177903e-01,
759 1.797889693808014e+02, -8.362340479938084e-01],
760 [ -2.447764190849540e+02, 3.063505753648256e+01,
761 -4.195525932156250e+01, 1.136590280389803e+01,
762 -1.064573816075257e+01, 5.228606946522749e+00,
763 -1.560849559916187e+00, 1.515872114883926e+00,
764 -8.362340479938084e-01, 3.833719335346630e+02]])
765 self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
766 0.807186823427233, 0.48950999450145, 0.995486532098031,
767 0.351243009576568, 0.704352576819321, 0.850648989740204,
768 0.314596738052894])
769 self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
770 30.344553414038817, 15.247917439094513, 24.060664905403492,
771 27.210293789825833, 47.122067744075842, 199.267136417856847,
772 -8.7934289814322 ])
773 def eval(self,x):
774 out=matrixmultiply(self.A,x)-self.b
775 for i in xrange(size(self.b)):
776 out[i]/=self.A[i,i]
777 return out
778 def bilinearform(self,x0,x1):
779 return dot(x0,x1)
780
781 tol=1.e-8
782 ll=LL()
783 x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
784 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
785
786 def testHomogeneousSaddlePointProblem_PCG(self):
787 from numarray import array,matrixmultiply, zeros, dot, size, Float64
788 from math import sqrt
789 class LL(HomogeneousSaddlePointProblem):
790 def initialize(self):
791 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
792 5.834798554135237e-01, -3.704394311709722e+00,
793 5.765369186984777e+00, -1.309786358737351e+01,
794 2.522087134507148e+01, -3.393956279045637e+01,
795 1.046856914770830e+02, -2.447764190849540e+02],
796 [ -2.391895572674098e-01, 1.256797283910693e+02,
797 -9.188270412920813e-01, 1.300169538880688e+00,
798 -5.353714719231424e-01, 2.674709444667012e+00,
799 -1.116097841269580e+01, 2.801193427514478e+01,
800 -3.877806125898224e+01, 3.063505753648256e+01],
801 [ 5.834798554135237e-01, -9.188270412920813e-01,
802 6.240841811806843e+01, -8.176289504109282e-01,
803 1.447935098417076e-01, -9.721424148655324e-01,
804 6.713551574117577e-01, -3.656297654168375e+00,
805 7.015141656913973e+00, -4.195525932156250e+01],
806 [ -3.704394311709722e+00, 1.300169538880688e+00,
807 -8.176289504109282e-01, 3.604980536782198e+01,
808 -6.241238423759328e-01, 1.142345320047869e+00,
809 -3.438816797096519e+00, 5.854857481367470e+00,
810 -4.524311288596452e+00, 1.136590280389803e+01],
811 [ 5.765369186984777e+00, -5.353714719231424e-01,
812 1.447935098417076e-01, -6.241238423759328e-01,
813 2.953997190215862e+01, -9.474729233464712e-01,
814 1.883516378345809e+00, -1.906274765704230e+00,
815 4.401859671778645e+00, -1.064573816075257e+01],
816 [ -1.309786358737351e+01, 2.674709444667012e+00,
817 -9.721424148655324e-01, 1.142345320047869e+00,
818 -9.474729233464712e-01, 2.876998216302979e+01,
819 -4.853065259692995e-01, 7.088596468102618e-01,
820 -8.972224295152829e-01, 5.228606946522749e+00],
821 [ 2.522087134507148e+01, -1.116097841269580e+01,
822 6.713551574117577e-01, -3.438816797096519e+00,
823 1.883516378345809e+00, -4.853065259692995e-01,
824 5.121175860935919e+01, -3.523133115905478e-01,
825 1.782136702229135e+00, -1.560849559916187e+00],
826 [ -3.393956279045637e+01, 2.801193427514478e+01,
827 -3.656297654168375e+00, 5.854857481367470e+00,
828 -1.906274765704230e+00, 7.088596468102618e-01,
829 -3.523133115905478e-01, 8.411681423853814e+01,
830 -5.238590858177903e-01, 1.515872114883926e+00],
831 [ 1.046856914770830e+02, -3.877806125898224e+01,
832 7.015141656913973e+00, -4.524311288596452e+00,
833 4.401859671778645e+00, -8.972224295152829e-01,
834 1.782136702229135e+00, -5.238590858177903e-01,
835 1.797889693808014e+02, -8.362340479938084e-01],
836 [ -2.447764190849540e+02, 3.063505753648256e+01,
837 -4.195525932156250e+01, 1.136590280389803e+01,
838 -1.064573816075257e+01, 5.228606946522749e+00,
839 -1.560849559916187e+00, 1.515872114883926e+00,
840 -8.362340479938084e-01, 3.833719335346630e+02]])
841 self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
842 0.253012436539738, -0.346223308561905, 0.2425508275422,
843 -0.194695862196008, 0.09451439391473, 0.302961126826511,
844 -0.236043777597633] )
845
846 self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
847 [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
848 [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
849 [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
850 [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
851 [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
852 [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
853 [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
854 [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
855 [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
856 self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
857
858 self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
859 6.512480714694167 , -5.727371407390975 , 4.802494840775022,
860 -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
861 -118.602464657076439])
862
863 self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
864 [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
865 [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
866 def inner_pBv(self,p,Bv):
867 return dot(p,Bv)
868 def inner_p(self,p0,p1):
869 return dot(p0,p1)
870 def inner_v(self,v0,v1):
871 return dot(v0,v1)
872 def B(self,v):
873 return matrixmultiply(transpose(self.Bt),v)
874 def solve_A(self,u,p):
875 out=self.b-matrixmultiply(self.A,u)-matrixmultiply(self.Bt,p)
876 # for i in xrange(size(self.b)): out[i]/=self.A[i,i]
877 return solve_linear_equations(self.A,out)+u*self.getSubProblemTolerance()
878 def solve_prec(self,p):
879 out=zeros((size(p),),Float64)
880 for i in xrange(size(p)): out[i]=p[i]*self.Sinv[i,i]
881 return out
882
883
884 tol=1.e-8
885 ll=LL()
886 ll.initialize()
887 ll.setTolerance(tol)
888 # ll.setSubToleranceReductionFactor(0.1)
889 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)
890 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
891 self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
892
893 def testHomogeneousSaddlePointProblem_GMRES(self):
894 from numarray import array,matrixmultiply, zeros, dot, size, Float64
895 from math import sqrt
896 class LL(HomogeneousSaddlePointProblem):
897 def initialize(self):
898 self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
899 5.834798554135237e-01, -3.704394311709722e+00,
900 5.765369186984777e+00, -1.309786358737351e+01,
901 2.522087134507148e+01, -3.393956279045637e+01,
902 1.046856914770830e+02, -2.447764190849540e+02],
903 [ -2.391895572674098e-01, 1.256797283910693e+02,
904 -9.188270412920813e-01, 1.300169538880688e+00,
905 -5.353714719231424e-01, 2.674709444667012e+00,
906 -1.116097841269580e+01, 2.801193427514478e+01,
907 -3.877806125898224e+01, 3.063505753648256e+01],
908 [ 5.834798554135237e-01, -9.188270412920813e-01,
909 6.240841811806843e+01, -8.176289504109282e-01,
910 1.447935098417076e-01, -9.721424148655324e-01,
911 6.713551574117577e-01, -3.656297654168375e+00,
912 7.015141656913973e+00, -4.195525932156250e+01],
913 [ -3.704394311709722e+00, 1.300169538880688e+00,
914 -8.176289504109282e-01, 3.604980536782198e+01,
915 -6.241238423759328e-01, 1.142345320047869e+00,
916 -3.438816797096519e+00, 5.854857481367470e+00,
917 -4.524311288596452e+00, 1.136590280389803e+01],
918 [ 5.765369186984777e+00, -5.353714719231424e-01,
919 1.447935098417076e-01, -6.241238423759328e-01,
920 2.953997190215862e+01, -9.474729233464712e-01,
921 1.883516378345809e+00, -1.906274765704230e+00,
922 4.401859671778645e+00, -1.064573816075257e+01],
923 [ -1.309786358737351e+01, 2.674709444667012e+00,
924 -9.721424148655324e-01, 1.142345320047869e+00,
925 -9.474729233464712e-01, 2.876998216302979e+01,
926 -4.853065259692995e-01, 7.088596468102618e-01,
927 -8.972224295152829e-01, 5.228606946522749e+00],
928 [ 2.522087134507148e+01, -1.116097841269580e+01,
929 6.713551574117577e-01, -3.438816797096519e+00,
930 1.883516378345809e+00, -4.853065259692995e-01,
931 5.121175860935919e+01, -3.523133115905478e-01,
932 1.782136702229135e+00, -1.560849559916187e+00],
933 [ -3.393956279045637e+01, 2.801193427514478e+01,
934 -3.656297654168375e+00, 5.854857481367470e+00,
935 -1.906274765704230e+00, 7.088596468102618e-01,
936 -3.523133115905478e-01, 8.411681423853814e+01,
937 -5.238590858177903e-01, 1.515872114883926e+00],
938 [ 1.046856914770830e+02, -3.877806125898224e+01,
939 7.015141656913973e+00, -4.524311288596452e+00,
940 4.401859671778645e+00, -8.972224295152829e-01,
941 1.782136702229135e+00, -5.238590858177903e-01,
942 1.797889693808014e+02, -8.362340479938084e-01],
943 [ -2.447764190849540e+02, 3.063505753648256e+01,
944 -4.195525932156250e+01, 1.136590280389803e+01,
945 -1.064573816075257e+01, 5.228606946522749e+00,
946 -1.560849559916187e+00, 1.515872114883926e+00,
947 -8.362340479938084e-01, 3.833719335346630e+02]])
948 self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
949 0.253012436539738, -0.346223308561905, 0.2425508275422,
950 -0.194695862196008, 0.09451439391473, 0.302961126826511,
951 -0.236043777597633] )
952
953 self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
954 [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
955 [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
956 [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
957 [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
958 [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
959 [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
960 [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
961 [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
962 [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
963 self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
964
965 self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
966 6.512480714694167 , -5.727371407390975 , 4.802494840775022,
967 -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
968 -118.602464657076439])
969
970 self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
971 [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
972 [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
973 def inner_pBv(self,p,Bv):
974 return dot(p,Bv)
975 def inner_p(self,p0,p1):
976 return dot(p0,p1)
977 def inner_v(self,v0,v1):
978 return dot(v0,v1)
979 def B(self,v):
980 return matrixmultiply(transpose(self.Bt),v)
981 def solve_A(self,u,p):
982 out=self.b-matrixmultiply(self.A,u)-matrixmultiply(self.Bt,p)
983 return solve_linear_equations(self.A,out)+u*self.getSubProblemTolerance()
984 def solve_prec(self,p):
985 out=zeros((size(p),),Float64)
986 for i in xrange(size(p)): out[i]=p[i]*self.Sinv[i,i]
987 # out=matrixmultiply(self.Sinv,p)
988 return out
989
990
991 tol=1.e-8
992 ll=LL()
993 ll.initialize()
994 ll.setTolerance(tol)
995 # ll.setSubToleranceReductionFactor(0.1)
996 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)
997 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
998 self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
999
1000 def testArithmeticTuple(self):
1001 a=ArithmeticTuple(1.,2.)
1002 self.failUnless(len(a)==2,"wrong length")
1003 self.failUnless(a[0]==1.,"wrong first item")
1004 self.failUnless(a[1]==2.,"wrong second item")
1005 c=a*6.
1006 self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple")
1007 self.failUnless(len(c)==2,"c has wrong length")
1008 self.failUnless(c[0]==6.,"c has wrong first item")
1009 self.failUnless(c[1]==12.,"c has wrong second item")
1010 b=5.*a
1011 self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple")
1012 self.failUnless(len(b)==2,"b has wrong length")
1013 self.failUnless(b[0]==5.,"b has wrong first item")
1014 self.failUnless(b[1]==10.,"b has wrong second item")
1015 a+=ArithmeticTuple(3.,4.)
1016 self.failUnless(a[0]==4.,"wrong first item of inplace update")
1017 self.failUnless(a[1]==6.,"wrong second item of inplace update")
1018
1019
1020
1021 class Test_pdetools(Test_pdetools_noLumping):
1022 def testProjector_rank0_fast_reduced(self):
1023 x=ContinuousFunction(self.domain).getX()
1024 h=Lsup(self.domain.getSize())
1025 p=Projector(self.domain,reduce=True,fast=True)
1026 td_ref=x[0]
1027 td=p(td_ref.interpolate(Function(self.domain)))
1028 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1029
1030 def testProjector_rank1_fast_reduced(self):
1031 x=ContinuousFunction(self.domain).getX()
1032 h=Lsup(self.domain.getSize())
1033 p=Projector(self.domain,reduce=True,fast=True)
1034 td_ref=x
1035 td=p(td_ref.interpolate(Function(self.domain)))
1036 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1037
1038 def testProjector_rank2_fast_reduced(self):
1039 x=ContinuousFunction(self.domain).getX()
1040 h=Lsup(self.domain.getSize())
1041 p=Projector(self.domain,reduce=True,fast=True)
1042 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
1043 td=p(td_ref.interpolate(Function(self.domain)))
1044 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1045
1046 def testProjector_rank3_fast_reduced(self):
1047 x=ContinuousFunction(self.domain).getX()
1048 h=Lsup(self.domain.getSize())
1049 p=Projector(self.domain,reduce=True,fast=True)
1050 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
1051 td=p(td_ref.interpolate(Function(self.domain)))
1052 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1053
1054 def testProjector_rank4_fast_reduced(self):
1055 x=ContinuousFunction(self.domain).getX()
1056 h=Lsup(self.domain.getSize())
1057 p=Projector(self.domain,reduce=True,fast=True)
1058 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
1059 td=p(td_ref.interpolate(Function(self.domain)))
1060 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1061
1062 def testProjector_rank0_fast_reduced_with_reduced_input(self):
1063 x=ContinuousFunction(self.domain).getX()
1064 h=Lsup(self.domain.getSize())
1065 p=Projector(self.domain,reduce=True,fast=True)
1066 td_ref=1.
1067 td=p(Data(td_ref,ReducedFunction(self.domain)))
1068 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1069
1070 def testProjector_rank1_fast_reduced_with_reduced_input(self):
1071 x=ContinuousFunction(self.domain).getX()
1072 h=Lsup(self.domain.getSize())
1073 p=Projector(self.domain,reduce=True,fast=True)
1074 td_ref=numarray.array([1.,2.,3.])
1075 td=p(Data(td_ref,ReducedFunction(self.domain)))
1076 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1077
1078 def testProjector_rank2_fast_reduced_with_reduced_input(self):
1079 x=ContinuousFunction(self.domain).getX()
1080 h=Lsup(self.domain.getSize())
1081 p=Projector(self.domain,reduce=True,fast=True)
1082 td_ref=numarray.array([[11.,12.],[21,22.]])
1083 td=p(Data(td_ref,ReducedFunction(self.domain)))
1084 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1085
1086 def testProjector_rank3_fast_reduced_with_reduced_input(self):
1087 x=ContinuousFunction(self.domain).getX()
1088 h=Lsup(self.domain.getSize())
1089 p=Projector(self.domain,reduce=True,fast=True)
1090 td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
1091 td=p(Data(td_ref,ReducedFunction(self.domain)))
1092 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1093
1094 def testProjector_rank4_fast_reduced_with_reduced_input(self):
1095 x=ContinuousFunction(self.domain).getX()
1096 h=Lsup(self.domain.getSize())
1097 p=Projector(self.domain,reduce=True,fast=True)
1098 td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
1099 td=p(Data(td_ref,ReducedFunction(self.domain)))
1100 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1101

  ViewVC Help
Powered by ViewVC 1.1.26