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

Annotation of /trunk/escript/test/python/test_pdetools.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2264 - (hide annotations)
Wed Feb 11 06:48:28 2009 UTC (11 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 65819 byte(s)
a new darcy flux solver.
1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # 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 ksteube 1312 #
12 ksteube 1809 ########################################################
13 gross 525
14 ksteube 1809 __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 gross 525 """
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 artak 1695 RES_TOL=1.e-8
30 gross 525 def setUp(self):
31     self.domain = Rectangle(10,10,2)
32 gross 798 def tearDown(self):
33     del self.domain
34 gross 525 suite = unittest.TestSuite()
35     suite.addTest(unittest.makeSuite(Test_LinearPDEOnFinley))
36     unittest.TextTestRunner(verbosity=2).run(suite)
37    
38 gross 637 @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 gross 525 """
45    
46     __author__="Lutz Gross, l.gross@uq.edu.au"
47    
48     import unittest
49     from esys.escript import *
50 gross 2100 from esys.escript.pdetools import Locator,Projector,TimeIntegrationManager,NoPDE,PCG, ArithmeticTuple, GMRES, MINRES, TFQMR, HomogeneousSaddlePointProblem
51 gross 1878 from esys.escript.pdetools import Defect, NewtonGMRES
52 gross 2100 from numarray.linear_algebra import solve_linear_equations
53 gross 525
54 gross 855 class Test_pdetools_noLumping(unittest.TestCase):
55 gross 525 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 gross 798 self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong")
66 gross 525
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 gross 798 self.failUnless(e<self.RES_TOL,"extrapolation is wrong")
77 gross 525
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 gross 798 self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong")
89 gross 525 xx=l(x)
90     self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type")
91 gross 798 self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector")
92 gross 525 xx=l(x[0]+x[1])
93     self.failUnless(isinstance(xx,float),"wrong scalar type")
94 gross 798 self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar")
95 gross 880
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 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
127 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
134 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
141 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
148 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
155 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
163 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
170 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
177 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
184 gross 525
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 gross 798 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
191 gross 525
192 gross 1072 def testProjector_rank0_with_reduced_input(self):
193 gross 525 x=ContinuousFunction(self.domain).getX()
194 gross 1072 p=Projector(self.domain,reduce=False,fast=False)
195 gross 525 td_ref=x[0]
196     td=p(td_ref.interpolate(Function(self.domain)))
197 gross 1072 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
198 gross 525
199 gross 1072 def testProjector_rank1_with_reduced_input(self):
200 gross 525 x=ContinuousFunction(self.domain).getX()
201 gross 1072 p=Projector(self.domain,reduce=False,fast=False)
202 gross 525 td_ref=x
203     td=p(td_ref.interpolate(Function(self.domain)))
204 gross 1072 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
205 gross 525
206 gross 1072 def testProjector_rank2_with_reduced_input(self):
207 gross 525 x=ContinuousFunction(self.domain).getX()
208 gross 1072 p=Projector(self.domain,reduce=False,fast=False)
209 gross 525 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
210     td=p(td_ref.interpolate(Function(self.domain)))
211 gross 1072 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
212 gross 525
213 gross 1072 def testProjector_rank3_with_reduced_input(self):
214 gross 525 x=ContinuousFunction(self.domain).getX()
215 gross 1072 p=Projector(self.domain,reduce=False,fast=False)
216 gross 525 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
217     td=p(td_ref.interpolate(Function(self.domain)))
218 gross 1072 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
219 gross 525
220 gross 1072 def testProjector_rank4_with_reduced_input(self):
221 gross 525 x=ContinuousFunction(self.domain).getX()
222 gross 1072 p=Projector(self.domain,reduce=False,fast=False)
223 gross 525 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 gross 1072 self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong")
226 gross 525
227 gross 1072
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 gross 525 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 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
272 gross 525
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 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
281 gross 525
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 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
290 gross 525
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 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
299 gross 525
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 gross 854 p.setValue(D=numarray.ones([2]),q=msk,r=2.)
305 gross 525 u=p.getSolution()
306     u_ex=msk*2.
307 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
308 gross 525
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 gross 854 p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk)
314 gross 525 u=p.getSolution()
315     u_ex=(1.-msk)
316 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
317 gross 525
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 gross 854 p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk,r=2.)
323 gross 525 u=p.getSolution()
324     u_ex=(1.-msk)+msk*2.
325 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
326 gross 525
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 artak 1695 p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.)
332 gross 525 u=p.getSolution()
333     u_ex=2.
334 gross 798 self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong")
335 gross 2156 #=====
336 ksteube 1312 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 gross 855
398 ksteube 1312 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 gross 2264 x,r,a_norm=PCG(b*1.,Ap,x_ref*0.,Ms,dot, atol=0, rtol=tol, iter_max=12)
408 ksteube 1312 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 gross 2100 def testMINRES(self):
412 gross 1468 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 gross 2156 x=MINRES(b*1.,Ap,x_ref*0,Ms,dot, atol=0, rtol=tol, iter_max=12)
483 gross 1468 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
484    
485 gross 2100 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 gross 2156 out=matrixmultiply(A,x)
549     for i in xrange(size(x)):
550     out[i]/=A[i,i]
551 gross 2100 return out
552    
553     tol=1.e-5
554 gross 2156 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 gross 2100 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 gross 2156 x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12)
629 gross 2100 self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
630    
631 gross 2264 def testGMRES_P_R(self):
632     from numarray import array,matrixmultiply, zeros, dot, size, Float64
633     from math import sqrt
634     A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
635     5.834798554135237e-01, -3.704394311709722e+00,
636     5.765369186984777e+00, -1.309786358737351e+01,
637     2.522087134507148e+01, -3.393956279045637e+01,
638     1.046856914770830e+02, -2.447764190849540e+02],
639     [ -2.391895572674098e-01, 1.256797283910693e+02,
640     -9.188270412920813e-01, 1.300169538880688e+00,
641     -5.353714719231424e-01, 2.674709444667012e+00,
642     -1.116097841269580e+01, 2.801193427514478e+01,
643     -3.877806125898224e+01, 3.063505753648256e+01],
644     [ 5.834798554135237e-01, -9.188270412920813e-01,
645     6.240841811806843e+01, -8.176289504109282e-01,
646     1.447935098417076e-01, -9.721424148655324e-01,
647     6.713551574117577e-01, -3.656297654168375e+00,
648     7.015141656913973e+00, -4.195525932156250e+01],
649     [ -3.704394311709722e+00, 1.300169538880688e+00,
650     -8.176289504109282e-01, 3.604980536782198e+01,
651     -6.241238423759328e-01, 1.142345320047869e+00,
652     -3.438816797096519e+00, 5.854857481367470e+00,
653     -4.524311288596452e+00, 1.136590280389803e+01],
654     [ 5.765369186984777e+00, -5.353714719231424e-01,
655     1.447935098417076e-01, -6.241238423759328e-01,
656     2.953997190215862e+01, -9.474729233464712e-01,
657     1.883516378345809e+00, -1.906274765704230e+00,
658     4.401859671778645e+00, -1.064573816075257e+01],
659     [ -1.309786358737351e+01, 2.674709444667012e+00,
660     -9.721424148655324e-01, 1.142345320047869e+00,
661     -9.474729233464712e-01, 2.876998216302979e+01,
662     -4.853065259692995e-01, 7.088596468102618e-01,
663     -8.972224295152829e-01, 5.228606946522749e+00],
664     [ 2.522087134507148e+01, -1.116097841269580e+01,
665     6.713551574117577e-01, -3.438816797096519e+00,
666     1.883516378345809e+00, -4.853065259692995e-01,
667     5.121175860935919e+01, -3.523133115905478e-01,
668     1.782136702229135e+00, -1.560849559916187e+00],
669     [ -3.393956279045637e+01, 2.801193427514478e+01,
670     -3.656297654168375e+00, 5.854857481367470e+00,
671     -1.906274765704230e+00, 7.088596468102618e-01,
672     -3.523133115905478e-01, 8.411681423853814e+01,
673     -5.238590858177903e-01, 1.515872114883926e+00],
674     [ 1.046856914770830e+02, -3.877806125898224e+01,
675     7.015141656913973e+00, -4.524311288596452e+00,
676     4.401859671778645e+00, -8.972224295152829e-01,
677     1.782136702229135e+00, -5.238590858177903e-01,
678     1.797889693808014e+02, -8.362340479938084e-01],
679     [ -2.447764190849540e+02, 3.063505753648256e+01,
680     -4.195525932156250e+01, 1.136590280389803e+01,
681     -1.064573816075257e+01, 5.228606946522749e+00,
682     -1.560849559916187e+00, 1.515872114883926e+00,
683     -8.362340479938084e-01, 3.833719335346630e+02]])
684     x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
685     0.807186823427233, 0.48950999450145, 0.995486532098031,
686     0.351243009576568, 0.704352576819321, 0.850648989740204,
687     0.314596738052894])
688     b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
689     30.344553414038817, 15.247917439094513, 24.060664905403492,
690     27.210293789825833, 47.122067744075842, 199.267136417856847,
691     -8.7934289814322 ])
692    
693     def Ap(x):
694     return matrixmultiply(A,x)
695     def P_Rp(x):
696     out=zeros(size(x))
697     for i in xrange(size(b)):
698     out[i]=x[i]/A[i,i]
699     return out
700    
701     tol=1.e-4
702     x=GMRES(b,Ap,x_ref*0,dot,atol=0, rtol=tol, iter_max=12,P_R=P_Rp)
703     self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution")
704    
705 gross 1878 def testNewtonGMRES(self):
706     from numarray import array,matrixmultiply, zeros, dot, size, Float64
707     from math import sqrt
708     class LL(Defect):
709     def __init__(self,*kwargs):
710     super(LL, self).__init__(*kwargs)
711     self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
712     5.834798554135237e-01, -3.704394311709722e+00,
713     5.765369186984777e+00, -1.309786358737351e+01,
714     2.522087134507148e+01, -3.393956279045637e+01,
715     1.046856914770830e+02, -2.447764190849540e+02],
716     [ -2.391895572674098e-01, 1.256797283910693e+02,
717     -9.188270412920813e-01, 1.300169538880688e+00,
718     -5.353714719231424e-01, 2.674709444667012e+00,
719     -1.116097841269580e+01, 2.801193427514478e+01,
720     -3.877806125898224e+01, 3.063505753648256e+01],
721     [ 5.834798554135237e-01, -9.188270412920813e-01,
722     6.240841811806843e+01, -8.176289504109282e-01,
723     1.447935098417076e-01, -9.721424148655324e-01,
724     6.713551574117577e-01, -3.656297654168375e+00,
725     7.015141656913973e+00, -4.195525932156250e+01],
726     [ -3.704394311709722e+00, 1.300169538880688e+00,
727     -8.176289504109282e-01, 3.604980536782198e+01,
728     -6.241238423759328e-01, 1.142345320047869e+00,
729     -3.438816797096519e+00, 5.854857481367470e+00,
730     -4.524311288596452e+00, 1.136590280389803e+01],
731     [ 5.765369186984777e+00, -5.353714719231424e-01,
732     1.447935098417076e-01, -6.241238423759328e-01,
733     2.953997190215862e+01, -9.474729233464712e-01,
734     1.883516378345809e+00, -1.906274765704230e+00,
735     4.401859671778645e+00, -1.064573816075257e+01],
736     [ -1.309786358737351e+01, 2.674709444667012e+00,
737     -9.721424148655324e-01, 1.142345320047869e+00,
738     -9.474729233464712e-01, 2.876998216302979e+01,
739     -4.853065259692995e-01, 7.088596468102618e-01,
740     -8.972224295152829e-01, 5.228606946522749e+00],
741     [ 2.522087134507148e+01, -1.116097841269580e+01,
742     6.713551574117577e-01, -3.438816797096519e+00,
743     1.883516378345809e+00, -4.853065259692995e-01,
744     5.121175860935919e+01, -3.523133115905478e-01,
745     1.782136702229135e+00, -1.560849559916187e+00],
746     [ -3.393956279045637e+01, 2.801193427514478e+01,
747     -3.656297654168375e+00, 5.854857481367470e+00,
748     -1.906274765704230e+00, 7.088596468102618e-01,
749     -3.523133115905478e-01, 8.411681423853814e+01,
750     -5.238590858177903e-01, 1.515872114883926e+00],
751     [ 1.046856914770830e+02, -3.877806125898224e+01,
752     7.015141656913973e+00, -4.524311288596452e+00,
753     4.401859671778645e+00, -8.972224295152829e-01,
754     1.782136702229135e+00, -5.238590858177903e-01,
755     1.797889693808014e+02, -8.362340479938084e-01],
756     [ -2.447764190849540e+02, 3.063505753648256e+01,
757     -4.195525932156250e+01, 1.136590280389803e+01,
758     -1.064573816075257e+01, 5.228606946522749e+00,
759     -1.560849559916187e+00, 1.515872114883926e+00,
760     -8.362340479938084e-01, 3.833719335346630e+02]])
761     self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
762     0.807186823427233, 0.48950999450145, 0.995486532098031,
763     0.351243009576568, 0.704352576819321, 0.850648989740204,
764     0.314596738052894])
765     self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
766     30.344553414038817, 15.247917439094513, 24.060664905403492,
767     27.210293789825833, 47.122067744075842, 199.267136417856847,
768     -8.7934289814322 ])
769     def eval(self,x):
770     out=matrixmultiply(self.A,x)-self.b
771     for i in xrange(size(self.b)):
772     out[i]/=self.A[i,i]
773     return out
774     def bilinearform(self,x0,x1):
775     return dot(x0,x1)
776    
777     tol=1.e-8
778     ll=LL()
779     x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
780     self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
781    
782 gross 2100 def testNewtonGMRES(self):
783     from numarray import array,matrixmultiply, zeros, dot, size, Float64
784     from math import sqrt
785     class LL(Defect):
786     def __init__(self,*kwargs):
787     super(LL, self).__init__(*kwargs)
788     self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
789     5.834798554135237e-01, -3.704394311709722e+00,
790     5.765369186984777e+00, -1.309786358737351e+01,
791     2.522087134507148e+01, -3.393956279045637e+01,
792     1.046856914770830e+02, -2.447764190849540e+02],
793     [ -2.391895572674098e-01, 1.256797283910693e+02,
794     -9.188270412920813e-01, 1.300169538880688e+00,
795     -5.353714719231424e-01, 2.674709444667012e+00,
796     -1.116097841269580e+01, 2.801193427514478e+01,
797     -3.877806125898224e+01, 3.063505753648256e+01],
798     [ 5.834798554135237e-01, -9.188270412920813e-01,
799     6.240841811806843e+01, -8.176289504109282e-01,
800     1.447935098417076e-01, -9.721424148655324e-01,
801     6.713551574117577e-01, -3.656297654168375e+00,
802     7.015141656913973e+00, -4.195525932156250e+01],
803     [ -3.704394311709722e+00, 1.300169538880688e+00,
804     -8.176289504109282e-01, 3.604980536782198e+01,
805     -6.241238423759328e-01, 1.142345320047869e+00,
806     -3.438816797096519e+00, 5.854857481367470e+00,
807     -4.524311288596452e+00, 1.136590280389803e+01],
808     [ 5.765369186984777e+00, -5.353714719231424e-01,
809     1.447935098417076e-01, -6.241238423759328e-01,
810     2.953997190215862e+01, -9.474729233464712e-01,
811     1.883516378345809e+00, -1.906274765704230e+00,
812     4.401859671778645e+00, -1.064573816075257e+01],
813     [ -1.309786358737351e+01, 2.674709444667012e+00,
814     -9.721424148655324e-01, 1.142345320047869e+00,
815     -9.474729233464712e-01, 2.876998216302979e+01,
816     -4.853065259692995e-01, 7.088596468102618e-01,
817     -8.972224295152829e-01, 5.228606946522749e+00],
818     [ 2.522087134507148e+01, -1.116097841269580e+01,
819     6.713551574117577e-01, -3.438816797096519e+00,
820     1.883516378345809e+00, -4.853065259692995e-01,
821     5.121175860935919e+01, -3.523133115905478e-01,
822     1.782136702229135e+00, -1.560849559916187e+00],
823     [ -3.393956279045637e+01, 2.801193427514478e+01,
824     -3.656297654168375e+00, 5.854857481367470e+00,
825     -1.906274765704230e+00, 7.088596468102618e-01,
826     -3.523133115905478e-01, 8.411681423853814e+01,
827     -5.238590858177903e-01, 1.515872114883926e+00],
828     [ 1.046856914770830e+02, -3.877806125898224e+01,
829     7.015141656913973e+00, -4.524311288596452e+00,
830     4.401859671778645e+00, -8.972224295152829e-01,
831     1.782136702229135e+00, -5.238590858177903e-01,
832     1.797889693808014e+02, -8.362340479938084e-01],
833     [ -2.447764190849540e+02, 3.063505753648256e+01,
834     -4.195525932156250e+01, 1.136590280389803e+01,
835     -1.064573816075257e+01, 5.228606946522749e+00,
836     -1.560849559916187e+00, 1.515872114883926e+00,
837     -8.362340479938084e-01, 3.833719335346630e+02]])
838     self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401,
839     0.807186823427233, 0.48950999450145, 0.995486532098031,
840     0.351243009576568, 0.704352576819321, 0.850648989740204,
841     0.314596738052894])
842     self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201,
843     30.344553414038817, 15.247917439094513, 24.060664905403492,
844     27.210293789825833, 47.122067744075842, 199.267136417856847,
845     -8.7934289814322 ])
846     def eval(self,x):
847     out=matrixmultiply(self.A,x)-self.b
848     for i in xrange(size(self.b)):
849     out[i]/=self.A[i,i]
850     return out
851     def bilinearform(self,x0,x1):
852     return dot(x0,x1)
853    
854     tol=1.e-8
855     ll=LL()
856     x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE)
857     self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution")
858    
859     def testHomogeneousSaddlePointProblem_PCG(self):
860     from numarray import array,matrixmultiply, zeros, dot, size, Float64
861     from math import sqrt
862     class LL(HomogeneousSaddlePointProblem):
863     def initialize(self):
864     self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
865     5.834798554135237e-01, -3.704394311709722e+00,
866     5.765369186984777e+00, -1.309786358737351e+01,
867     2.522087134507148e+01, -3.393956279045637e+01,
868     1.046856914770830e+02, -2.447764190849540e+02],
869     [ -2.391895572674098e-01, 1.256797283910693e+02,
870     -9.188270412920813e-01, 1.300169538880688e+00,
871     -5.353714719231424e-01, 2.674709444667012e+00,
872     -1.116097841269580e+01, 2.801193427514478e+01,
873     -3.877806125898224e+01, 3.063505753648256e+01],
874     [ 5.834798554135237e-01, -9.188270412920813e-01,
875     6.240841811806843e+01, -8.176289504109282e-01,
876     1.447935098417076e-01, -9.721424148655324e-01,
877     6.713551574117577e-01, -3.656297654168375e+00,
878     7.015141656913973e+00, -4.195525932156250e+01],
879     [ -3.704394311709722e+00, 1.300169538880688e+00,
880     -8.176289504109282e-01, 3.604980536782198e+01,
881     -6.241238423759328e-01, 1.142345320047869e+00,
882     -3.438816797096519e+00, 5.854857481367470e+00,
883     -4.524311288596452e+00, 1.136590280389803e+01],
884     [ 5.765369186984777e+00, -5.353714719231424e-01,
885     1.447935098417076e-01, -6.241238423759328e-01,
886     2.953997190215862e+01, -9.474729233464712e-01,
887     1.883516378345809e+00, -1.906274765704230e+00,
888     4.401859671778645e+00, -1.064573816075257e+01],
889     [ -1.309786358737351e+01, 2.674709444667012e+00,
890     -9.721424148655324e-01, 1.142345320047869e+00,
891     -9.474729233464712e-01, 2.876998216302979e+01,
892     -4.853065259692995e-01, 7.088596468102618e-01,
893     -8.972224295152829e-01, 5.228606946522749e+00],
894     [ 2.522087134507148e+01, -1.116097841269580e+01,
895     6.713551574117577e-01, -3.438816797096519e+00,
896     1.883516378345809e+00, -4.853065259692995e-01,
897     5.121175860935919e+01, -3.523133115905478e-01,
898     1.782136702229135e+00, -1.560849559916187e+00],
899     [ -3.393956279045637e+01, 2.801193427514478e+01,
900     -3.656297654168375e+00, 5.854857481367470e+00,
901     -1.906274765704230e+00, 7.088596468102618e-01,
902     -3.523133115905478e-01, 8.411681423853814e+01,
903     -5.238590858177903e-01, 1.515872114883926e+00],
904     [ 1.046856914770830e+02, -3.877806125898224e+01,
905     7.015141656913973e+00, -4.524311288596452e+00,
906     4.401859671778645e+00, -8.972224295152829e-01,
907     1.782136702229135e+00, -5.238590858177903e-01,
908     1.797889693808014e+02, -8.362340479938084e-01],
909     [ -2.447764190849540e+02, 3.063505753648256e+01,
910     -4.195525932156250e+01, 1.136590280389803e+01,
911     -1.064573816075257e+01, 5.228606946522749e+00,
912     -1.560849559916187e+00, 1.515872114883926e+00,
913     -8.362340479938084e-01, 3.833719335346630e+02]])
914     self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
915     0.253012436539738, -0.346223308561905, 0.2425508275422,
916     -0.194695862196008, 0.09451439391473, 0.302961126826511,
917     -0.236043777597633] )
918    
919     self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
920     [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
921     [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
922     [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
923     [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
924     [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
925     [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
926     [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
927     [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
928     [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
929     self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
930    
931     self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
932     6.512480714694167 , -5.727371407390975 , 4.802494840775022,
933     -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
934     -118.602464657076439])
935    
936     self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
937     [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
938     [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
939 gross 2251 def inner_pBv(self,p,v):
940     return dot(p,matrixmultiply(transpose(self.Bt),v))
941 gross 2100 def inner_p(self,p0,p1):
942     return dot(p0,p1)
943 gross 2251 def norm_v(self,v):
944     return sqrt(dot(v,v))
945     def getV(self,p,v):
946     out=self.b-matrixmultiply(self.Bt,p)
947     return solve_linear_equations(self.A,out)+v*self.getSubProblemTolerance()
948     def norm_Bv(self,v):
949     t=matrixmultiply(transpose(self.Bt),v)
950     return sqrt(dot(t,t))
951     def solve_AinvBt(self,p):
952     out=matrixmultiply(self.Bt,p)
953     return solve_linear_equations(self.A,out)
954     def solve_precB(self,v):
955     out=matrixmultiply(transpose(self.Bt),v)
956     for i in xrange(size(out)): out[i]*=self.Sinv[i,i]
957     return out
958 gross 2100
959     tol=1.e-8
960     ll=LL()
961     ll.initialize()
962     ll.setTolerance(tol)
963     # ll.setSubToleranceReductionFactor(0.1)
964 gross 2251 x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, show_details=False, usePCG=True, iter_restart=20,max_correction_steps=3)
965 gross 2100 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
966     self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
967    
968     def testHomogeneousSaddlePointProblem_GMRES(self):
969     from numarray import array,matrixmultiply, zeros, dot, size, Float64
970     from math import sqrt
971     class LL(HomogeneousSaddlePointProblem):
972     def initialize(self):
973     self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01,
974     5.834798554135237e-01, -3.704394311709722e+00,
975     5.765369186984777e+00, -1.309786358737351e+01,
976     2.522087134507148e+01, -3.393956279045637e+01,
977     1.046856914770830e+02, -2.447764190849540e+02],
978     [ -2.391895572674098e-01, 1.256797283910693e+02,
979     -9.188270412920813e-01, 1.300169538880688e+00,
980     -5.353714719231424e-01, 2.674709444667012e+00,
981     -1.116097841269580e+01, 2.801193427514478e+01,
982     -3.877806125898224e+01, 3.063505753648256e+01],
983     [ 5.834798554135237e-01, -9.188270412920813e-01,
984     6.240841811806843e+01, -8.176289504109282e-01,
985     1.447935098417076e-01, -9.721424148655324e-01,
986     6.713551574117577e-01, -3.656297654168375e+00,
987     7.015141656913973e+00, -4.195525932156250e+01],
988     [ -3.704394311709722e+00, 1.300169538880688e+00,
989     -8.176289504109282e-01, 3.604980536782198e+01,
990     -6.241238423759328e-01, 1.142345320047869e+00,
991     -3.438816797096519e+00, 5.854857481367470e+00,
992     -4.524311288596452e+00, 1.136590280389803e+01],
993     [ 5.765369186984777e+00, -5.353714719231424e-01,
994     1.447935098417076e-01, -6.241238423759328e-01,
995     2.953997190215862e+01, -9.474729233464712e-01,
996     1.883516378345809e+00, -1.906274765704230e+00,
997     4.401859671778645e+00, -1.064573816075257e+01],
998     [ -1.309786358737351e+01, 2.674709444667012e+00,
999     -9.721424148655324e-01, 1.142345320047869e+00,
1000     -9.474729233464712e-01, 2.876998216302979e+01,
1001     -4.853065259692995e-01, 7.088596468102618e-01,
1002     -8.972224295152829e-01, 5.228606946522749e+00],
1003     [ 2.522087134507148e+01, -1.116097841269580e+01,
1004     6.713551574117577e-01, -3.438816797096519e+00,
1005     1.883516378345809e+00, -4.853065259692995e-01,
1006     5.121175860935919e+01, -3.523133115905478e-01,
1007     1.782136702229135e+00, -1.560849559916187e+00],
1008     [ -3.393956279045637e+01, 2.801193427514478e+01,
1009     -3.656297654168375e+00, 5.854857481367470e+00,
1010     -1.906274765704230e+00, 7.088596468102618e-01,
1011     -3.523133115905478e-01, 8.411681423853814e+01,
1012     -5.238590858177903e-01, 1.515872114883926e+00],
1013     [ 1.046856914770830e+02, -3.877806125898224e+01,
1014     7.015141656913973e+00, -4.524311288596452e+00,
1015     4.401859671778645e+00, -8.972224295152829e-01,
1016     1.782136702229135e+00, -5.238590858177903e-01,
1017     1.797889693808014e+02, -8.362340479938084e-01],
1018     [ -2.447764190849540e+02, 3.063505753648256e+01,
1019     -4.195525932156250e+01, 1.136590280389803e+01,
1020     -1.064573816075257e+01, 5.228606946522749e+00,
1021     -1.560849559916187e+00, 1.515872114883926e+00,
1022     -8.362340479938084e-01, 3.833719335346630e+02]])
1023     self.x_ref=array([ 0.100225501676291, -0.308862704993209, 0.064097238997721,
1024     0.253012436539738, -0.346223308561905, 0.2425508275422,
1025     -0.194695862196008, 0.09451439391473, 0.302961126826511,
1026     -0.236043777597633] )
1027    
1028     self.Bt=array([[ 0.01627853113636 ,0.06688235764255 , 0.004870689484614],
1029     [ 0.062879587145773 ,0.038798770300146, 0.022155850155616],
1030     [ 0.09312121957248 ,0.110244632756116, 0.14053347386784 ],
1031     [ 0.059000597728388 ,0.090986953740106, 0.035316011834982],
1032     [ 0.091209362659698 ,0.13205572801294 , 0.069462874306956],
1033     [ 0.077790176986096 ,0.133626423045765, 0.011149969846981],
1034     [ 0.01407283482513 ,0.094910926488907, 0.133498532648644],
1035     [ 0.025728916673085 ,0.102542818811672, 0.13657268163218 ],
1036     [ 0.071254288170748 ,0.071738715618163, 0.078005951991733],
1037     [ 0.049463014576779 ,0.103559223780991, 0.003356415647637]])
1038     self.p_ref = array([ 2.580984952252628 ,4.054090902056985, 0.935138168128546])
1039    
1040     self.b=array([ 123.322775367582238, -51.556206655564573 , 16.220697868056913,
1041     6.512480714694167 , -5.727371407390975 , 4.802494840775022,
1042     -4.171606044721161 , -1.862366353566293 ,74.850226163257105,
1043     -118.602464657076439])
1044    
1045     self.Sinv=array([[ 9313.705360982807179,-5755.536981691270739, 806.289245589733696],
1046     [-5755.536981691271649, 4606.321002756208145,-1630.50619635660928 ],
1047     [ 806.289245589733468,-1630.506196356609053, 2145.65035816388945 ]])
1048 gross 2251 def inner_pBv(self,p,v):
1049     return dot(p,matrixmultiply(transpose(self.Bt),v))
1050 gross 2100 def inner_p(self,p0,p1):
1051     return dot(p0,p1)
1052 gross 2251 def norm_v(self,v):
1053     return sqrt(dot(v,v))
1054     def getV(self,p,v):
1055     out=self.b-matrixmultiply(self.Bt,p)
1056     return solve_linear_equations(self.A,out)+v*self.getSubProblemTolerance()
1057     def norm_Bv(self,v):
1058     t=matrixmultiply(transpose(self.Bt),v)
1059     return sqrt(dot(t,t))
1060     def solve_AinvBt(self,p):
1061     out=matrixmultiply(self.Bt,p)
1062     return solve_linear_equations(self.A,out)
1063     def solve_precB(self,v):
1064     out=matrixmultiply(transpose(self.Bt),v)
1065     for i in xrange(size(out)): out[i]*=self.Sinv[i,i]
1066     return out
1067 gross 2100
1068     tol=1.e-8
1069     ll=LL()
1070     ll.initialize()
1071     ll.setTolerance(tol)
1072     # ll.setSubToleranceReductionFactor(0.1)
1073 gross 2251 x,p=ll.solve(ll.x_ref*1.20,ll.p_ref*(-2),max_iter=20, verbose=False, show_details=False, usePCG=False, iter_restart=20,max_correction_steps=3)
1074 gross 2100 self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong x solution")
1075     self.failUnless(Lsup(p-ll.p_ref)<=Lsup(ll.p_ref)*tol*10.,"wrong p solution")
1076    
1077 gross 1331 def testArithmeticTuple(self):
1078     a=ArithmeticTuple(1.,2.)
1079     self.failUnless(len(a)==2,"wrong length")
1080     self.failUnless(a[0]==1.,"wrong first item")
1081     self.failUnless(a[1]==2.,"wrong second item")
1082     c=a*6.
1083     self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple")
1084     self.failUnless(len(c)==2,"c has wrong length")
1085     self.failUnless(c[0]==6.,"c has wrong first item")
1086     self.failUnless(c[1]==12.,"c has wrong second item")
1087     b=5.*a
1088     self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple")
1089     self.failUnless(len(b)==2,"b has wrong length")
1090     self.failUnless(b[0]==5.,"b has wrong first item")
1091     self.failUnless(b[1]==10.,"b has wrong second item")
1092     a+=ArithmeticTuple(3.,4.)
1093     self.failUnless(a[0]==4.,"wrong first item of inplace update")
1094     self.failUnless(a[1]==6.,"wrong second item of inplace update")
1095 ksteube 1312
1096 gross 1331
1097    
1098 gross 855 class Test_pdetools(Test_pdetools_noLumping):
1099 gross 1072 def testProjector_rank0_fast_reduced(self):
1100 gross 855 x=ContinuousFunction(self.domain).getX()
1101     h=Lsup(self.domain.getSize())
1102 gross 1072 p=Projector(self.domain,reduce=True,fast=True)
1103 gross 855 td_ref=x[0]
1104     td=p(td_ref.interpolate(Function(self.domain)))
1105     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1106    
1107 gross 1072 def testProjector_rank1_fast_reduced(self):
1108 gross 855 x=ContinuousFunction(self.domain).getX()
1109     h=Lsup(self.domain.getSize())
1110 gross 1072 p=Projector(self.domain,reduce=True,fast=True)
1111 gross 855 td_ref=x
1112     td=p(td_ref.interpolate(Function(self.domain)))
1113     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1114    
1115 gross 1072 def testProjector_rank2_fast_reduced(self):
1116 gross 855 x=ContinuousFunction(self.domain).getX()
1117     h=Lsup(self.domain.getSize())
1118 gross 1072 p=Projector(self.domain,reduce=True,fast=True)
1119 gross 855 td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1])
1120     td=p(td_ref.interpolate(Function(self.domain)))
1121     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1122    
1123 gross 1072 def testProjector_rank3_fast_reduced(self):
1124 gross 855 x=ContinuousFunction(self.domain).getX()
1125     h=Lsup(self.domain.getSize())
1126 gross 1072 p=Projector(self.domain,reduce=True,fast=True)
1127 gross 855 td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1])
1128     td=p(td_ref.interpolate(Function(self.domain)))
1129     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1130    
1131 gross 1072 def testProjector_rank4_fast_reduced(self):
1132 gross 855 x=ContinuousFunction(self.domain).getX()
1133     h=Lsup(self.domain.getSize())
1134 gross 1072 p=Projector(self.domain,reduce=True,fast=True)
1135 gross 855 td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1])
1136     td=p(td_ref.interpolate(Function(self.domain)))
1137     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1138 gross 1072
1139     def testProjector_rank0_fast_reduced_with_reduced_input(self):
1140     x=ContinuousFunction(self.domain).getX()
1141     h=Lsup(self.domain.getSize())
1142     p=Projector(self.domain,reduce=True,fast=True)
1143     td_ref=1.
1144     td=p(Data(td_ref,ReducedFunction(self.domain)))
1145     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1146    
1147     def testProjector_rank1_fast_reduced_with_reduced_input(self):
1148     x=ContinuousFunction(self.domain).getX()
1149     h=Lsup(self.domain.getSize())
1150     p=Projector(self.domain,reduce=True,fast=True)
1151     td_ref=numarray.array([1.,2.,3.])
1152     td=p(Data(td_ref,ReducedFunction(self.domain)))
1153     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1154    
1155     def testProjector_rank2_fast_reduced_with_reduced_input(self):
1156     x=ContinuousFunction(self.domain).getX()
1157     h=Lsup(self.domain.getSize())
1158     p=Projector(self.domain,reduce=True,fast=True)
1159     td_ref=numarray.array([[11.,12.],[21,22.]])
1160     td=p(Data(td_ref,ReducedFunction(self.domain)))
1161     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1162    
1163     def testProjector_rank3_fast_reduced_with_reduced_input(self):
1164     x=ContinuousFunction(self.domain).getX()
1165     h=Lsup(self.domain.getSize())
1166     p=Projector(self.domain,reduce=True,fast=True)
1167     td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]])
1168     td=p(Data(td_ref,ReducedFunction(self.domain)))
1169     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1170    
1171     def testProjector_rank4_fast_reduced_with_reduced_input(self):
1172     x=ContinuousFunction(self.domain).getX()
1173     h=Lsup(self.domain.getSize())
1174     p=Projector(self.domain,reduce=True,fast=True)
1175     td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]])
1176     td=p(Data(td_ref,ReducedFunction(self.domain)))
1177     self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong")
1178    

  ViewVC Help
Powered by ViewVC 1.1.26