| 1 |
# |
| 2 |
# $Id$ |
| 3 |
# |
| 4 |
####################################################### |
| 5 |
# |
| 6 |
# Copyright 2003-2007 by ACceSS MNRF |
| 7 |
# Copyright 2007 by University of Queensland |
| 8 |
# |
| 9 |
# http://esscc.uq.edu.au |
| 10 |
# Primary Business: Queensland, Australia |
| 11 |
# Licensed under the Open Software License version 3.0 |
| 12 |
# http://www.opensource.org/licenses/osl-3.0.php |
| 13 |
# |
| 14 |
####################################################### |
| 15 |
# |
| 16 |
|
| 17 |
""" |
| 18 |
Test suite for the pdetools module |
| 19 |
|
| 20 |
The tests must be linked with a Domain class object in the setUp method: |
| 21 |
|
| 22 |
from esys.finley import Rectangle |
| 23 |
class Test_LinearPDEOnFinley(Test_LinearPDE): |
| 24 |
RES_TOL=1.e-8 |
| 25 |
def setUp(self): |
| 26 |
self.domain = Rectangle(10,10,2) |
| 27 |
def tearDown(self): |
| 28 |
del self.domain |
| 29 |
suite = unittest.TestSuite() |
| 30 |
suite.addTest(unittest.makeSuite(Test_LinearPDEOnFinley)) |
| 31 |
unittest.TextTestRunner(verbosity=2).run(suite) |
| 32 |
|
| 33 |
@var __author__: name of author |
| 34 |
@var __copyright__: copyrights |
| 35 |
@var __license__: licence agreement |
| 36 |
@var __url__: url entry point on documentation |
| 37 |
@var __version__: version |
| 38 |
@var __date__: date of the version |
| 39 |
""" |
| 40 |
|
| 41 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
| 42 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
| 43 |
http://www.access.edu.au |
| 44 |
Primary Business: Queensland, Australia""" |
| 45 |
__license__="""Licensed under the Open Software License version 3.0 |
| 46 |
http://www.opensource.org/licenses/osl-3.0.php""" |
| 47 |
__url__="http://www.iservo.edu.au/esys/escript" |
| 48 |
__version__="$Revision$" |
| 49 |
__date__="$Date$" |
| 50 |
|
| 51 |
import unittest |
| 52 |
from esys.escript import * |
| 53 |
from esys.escript.pdetools import Locator,Projector,TimeIntegrationManager,NoPDE,PCG, IterationHistory, ArithmeticTuple, GMRES |
| 54 |
|
| 55 |
class Test_pdetools_noLumping(unittest.TestCase): |
| 56 |
DEBUG=False |
| 57 |
VERBOSE=False |
| 58 |
def test_TimeIntegrationManager_scalar(self): |
| 59 |
t=0. |
| 60 |
dt=0.1 |
| 61 |
tm=TimeIntegrationManager(0.,p=1) |
| 62 |
while t<1.: |
| 63 |
t+=dt |
| 64 |
tm.checkin(dt,t) |
| 65 |
v_guess=tm.extrapolate(dt) |
| 66 |
self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong") |
| 67 |
|
| 68 |
def test_TimeIntegrationManager_vector(self): |
| 69 |
t=0. |
| 70 |
dt=0.3 |
| 71 |
tm=TimeIntegrationManager(0.,0.,p=1) |
| 72 |
while t<1.: |
| 73 |
t+=dt |
| 74 |
tm.checkin(dt,t,3*t) |
| 75 |
v_guess=tm.extrapolate(dt) |
| 76 |
e=max(abs(v_guess[0]-(tm.getTime()+dt)),abs(v_guess[1]-(tm.getTime()+dt)*3.)) |
| 77 |
self.failUnless(e<self.RES_TOL,"extrapolation is wrong") |
| 78 |
|
| 79 |
def test_Locator(self): |
| 80 |
x=self.domain.getX() |
| 81 |
l=Locator(self.domain,numarray.ones((self.domain.getDim(),))) |
| 82 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 83 |
|
| 84 |
l=Locator(ContinuousFunction(self.domain),numarray.ones((self.domain.getDim(),))) |
| 85 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 86 |
|
| 87 |
xx=l.getX() |
| 88 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 89 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong") |
| 90 |
xx=l(x) |
| 91 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 92 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector") |
| 93 |
xx=l(x[0]+x[1]) |
| 94 |
self.failUnless(isinstance(xx,float),"wrong scalar type") |
| 95 |
self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar") |
| 96 |
|
| 97 |
def test_Locator_withList(self): |
| 98 |
x=self.domain.getX() |
| 99 |
arg=[numarray.ones((self.domain.getDim(),)), numarray.zeros((self.domain.getDim(),))] |
| 100 |
l=Locator(self.domain,arg) |
| 101 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 102 |
|
| 103 |
l=Locator(ContinuousFunction(self.domain),arg) |
| 104 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 105 |
|
| 106 |
xx=l.getX() |
| 107 |
self.failUnless(isinstance(xx,list),"list expected") |
| 108 |
for i in range(len(xx)): |
| 109 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item"%i) |
| 110 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong"%i) |
| 111 |
xx=l(x) |
| 112 |
self.failUnless(isinstance(xx,list),"list expected (2)") |
| 113 |
for i in range(len(xx)): |
| 114 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item (2)"%i) |
| 115 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong (2)"%i) |
| 116 |
xx=l(x[0]+x[1]) |
| 117 |
self.failUnless(isinstance(xx,list),"list expected (3)") |
| 118 |
for i in range(len(xx)): |
| 119 |
self.failUnless(isinstance(xx[i],float),"wrong scalar type") |
| 120 |
self.failUnless(abs(xx[i]-(arg[i][0]+arg[i][1]))<self.RES_TOL,"value wrong scalar") |
| 121 |
|
| 122 |
def testProjector_rank0(self): |
| 123 |
x=ContinuousFunction(self.domain).getX() |
| 124 |
p=Projector(self.domain,reduce=False,fast=False) |
| 125 |
td_ref=x[0] |
| 126 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 127 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 128 |
|
| 129 |
def testProjector_rank1(self): |
| 130 |
x=ContinuousFunction(self.domain).getX() |
| 131 |
p=Projector(self.domain,reduce=False,fast=False) |
| 132 |
td_ref=x |
| 133 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 134 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 135 |
|
| 136 |
def testProjector_rank2(self): |
| 137 |
x=ContinuousFunction(self.domain).getX() |
| 138 |
p=Projector(self.domain,reduce=False,fast=False) |
| 139 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 140 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 141 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 142 |
|
| 143 |
def testProjector_rank3(self): |
| 144 |
x=ContinuousFunction(self.domain).getX() |
| 145 |
p=Projector(self.domain,reduce=False,fast=False) |
| 146 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 147 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 148 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 149 |
|
| 150 |
def testProjector_rank4(self): |
| 151 |
x=ContinuousFunction(self.domain).getX() |
| 152 |
p=Projector(self.domain,reduce=False,fast=False) |
| 153 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 154 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 155 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 156 |
|
| 157 |
|
| 158 |
def testProjector_rank0_reduced(self): |
| 159 |
x=ContinuousFunction(self.domain).getX() |
| 160 |
p=Projector(self.domain,reduce=True,fast=False) |
| 161 |
td_ref=x[0] |
| 162 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 163 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 164 |
|
| 165 |
def testProjector_rank1_reduced(self): |
| 166 |
x=ContinuousFunction(self.domain).getX() |
| 167 |
p=Projector(self.domain,reduce=True,fast=False) |
| 168 |
td_ref=x |
| 169 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 170 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 171 |
|
| 172 |
def testProjector_rank2_reduced(self): |
| 173 |
x=ContinuousFunction(self.domain).getX() |
| 174 |
p=Projector(self.domain,reduce=True,fast=False) |
| 175 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 176 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 177 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 178 |
|
| 179 |
def testProjector_rank3_reduced(self): |
| 180 |
x=ContinuousFunction(self.domain).getX() |
| 181 |
p=Projector(self.domain,reduce=True,fast=False) |
| 182 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 183 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 184 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 185 |
|
| 186 |
def testProjector_rank4_reduced(self): |
| 187 |
x=ContinuousFunction(self.domain).getX() |
| 188 |
p=Projector(self.domain,reduce=True,fast=False) |
| 189 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 190 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 191 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 192 |
|
| 193 |
def testProjector_rank0_with_reduced_input(self): |
| 194 |
x=ContinuousFunction(self.domain).getX() |
| 195 |
p=Projector(self.domain,reduce=False,fast=False) |
| 196 |
td_ref=x[0] |
| 197 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 198 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 199 |
|
| 200 |
def testProjector_rank1_with_reduced_input(self): |
| 201 |
x=ContinuousFunction(self.domain).getX() |
| 202 |
p=Projector(self.domain,reduce=False,fast=False) |
| 203 |
td_ref=x |
| 204 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 205 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 206 |
|
| 207 |
def testProjector_rank2_with_reduced_input(self): |
| 208 |
x=ContinuousFunction(self.domain).getX() |
| 209 |
p=Projector(self.domain,reduce=False,fast=False) |
| 210 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 211 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 212 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 213 |
|
| 214 |
def testProjector_rank3_with_reduced_input(self): |
| 215 |
x=ContinuousFunction(self.domain).getX() |
| 216 |
p=Projector(self.domain,reduce=False,fast=False) |
| 217 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 218 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 219 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 220 |
|
| 221 |
def testProjector_rank4_with_reduced_input(self): |
| 222 |
x=ContinuousFunction(self.domain).getX() |
| 223 |
p=Projector(self.domain,reduce=False,fast=False) |
| 224 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 225 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 226 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 227 |
|
| 228 |
|
| 229 |
def testProjector_rank0_reduced_with_reduced_input(self): |
| 230 |
x=ContinuousFunction(self.domain).getX() |
| 231 |
p=Projector(self.domain,reduce=True,fast=False) |
| 232 |
td_ref=1. |
| 233 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 234 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 235 |
|
| 236 |
def testProjector_rank1_reduced_with_reduced_input(self): |
| 237 |
x=ContinuousFunction(self.domain).getX() |
| 238 |
p=Projector(self.domain,reduce=True,fast=False) |
| 239 |
td_ref=numarray.array([1.,2.,3.]) |
| 240 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 241 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 242 |
|
| 243 |
def testProjector_rank2_reduced_with_reduced_input(self): |
| 244 |
x=ContinuousFunction(self.domain).getX() |
| 245 |
p=Projector(self.domain,reduce=True,fast=False) |
| 246 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 247 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 248 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 249 |
|
| 250 |
def testProjector_rank3_reduced_with_reduced_input(self): |
| 251 |
x=ContinuousFunction(self.domain).getX() |
| 252 |
p=Projector(self.domain,reduce=True,fast=False) |
| 253 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 254 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 255 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 256 |
|
| 257 |
def testProjector_rank4_reduced_with_reduced_input(self): |
| 258 |
x=ContinuousFunction(self.domain).getX() |
| 259 |
p=Projector(self.domain,reduce=True,fast=False) |
| 260 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 261 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 262 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 263 |
|
| 264 |
|
| 265 |
def test_NoPDE_scalar_missing_r(self): |
| 266 |
p=NoPDE(self.domain) |
| 267 |
x=self.domain.getX() |
| 268 |
msk=whereZero(x[0]) |
| 269 |
p.setValue(D=1.,Y=1.,q=msk) |
| 270 |
u=p.getSolution() |
| 271 |
u_ex=(1.-msk) |
| 272 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 273 |
|
| 274 |
def test_NoPDE_scalar_missing_Y(self): |
| 275 |
p=NoPDE(self.domain) |
| 276 |
x=self.domain.getX() |
| 277 |
msk=whereZero(x[0]) |
| 278 |
p.setValue(D=1.,q=msk,r=2.) |
| 279 |
u=p.getSolution() |
| 280 |
u_ex=msk*2. |
| 281 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 282 |
|
| 283 |
def test_NoPDE_scalar_constant(self): |
| 284 |
p=NoPDE(self.domain) |
| 285 |
x=self.domain.getX() |
| 286 |
msk=whereZero(x[0]) |
| 287 |
p.setValue(D=1.,Y=1.,q=msk,r=2.) |
| 288 |
u=p.getSolution() |
| 289 |
u_ex=(1.-msk)+msk*2. |
| 290 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 291 |
|
| 292 |
def test_NoPDE_scalar_variable(self): |
| 293 |
p=NoPDE(self.domain) |
| 294 |
x=self.domain.getX() |
| 295 |
msk=whereZero(x[0]) |
| 296 |
p.setValue(D=x[0],Y=2*x[0],q=msk,r=2.) |
| 297 |
u=p.getSolution() |
| 298 |
u_ex=2. |
| 299 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 300 |
|
| 301 |
def test_NoPDE_vector_missing_Y(self): |
| 302 |
p=NoPDE(self.domain) |
| 303 |
x=self.domain.getX() |
| 304 |
msk=whereZero(x[0])*[1.,0.] |
| 305 |
p.setValue(D=numarray.ones([2]),q=msk,r=2.) |
| 306 |
u=p.getSolution() |
| 307 |
u_ex=msk*2. |
| 308 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 309 |
|
| 310 |
def test_NoPDE_vector_missing_r(self): |
| 311 |
p=NoPDE(self.domain) |
| 312 |
x=self.domain.getX() |
| 313 |
msk=whereZero(x[0])*[1.,0.] |
| 314 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk) |
| 315 |
u=p.getSolution() |
| 316 |
u_ex=(1.-msk) |
| 317 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 318 |
|
| 319 |
def test_NoPDE_vector_constant(self): |
| 320 |
p=NoPDE(self.domain) |
| 321 |
x=self.domain.getX() |
| 322 |
msk=whereZero(x[0])*[1.,0.] |
| 323 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk,r=2.) |
| 324 |
u=p.getSolution() |
| 325 |
u_ex=(1.-msk)+msk*2. |
| 326 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 327 |
|
| 328 |
def test_NoPDE_vector_variable(self): |
| 329 |
p=NoPDE(self.domain) |
| 330 |
x=self.domain.getX() |
| 331 |
msk=whereZero(x[0])*[1.,0.] |
| 332 |
p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.) |
| 333 |
u=p.getSolution() |
| 334 |
u_ex=2. |
| 335 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 336 |
|
| 337 |
def testPCG(self): |
| 338 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 339 |
from math import sqrt |
| 340 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 341 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 342 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 343 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 344 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 345 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 346 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 347 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 348 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 349 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 350 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 351 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 352 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 353 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 354 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 355 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 356 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 357 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 358 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 359 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 360 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 361 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 362 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 363 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 364 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 365 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 366 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 367 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 368 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 369 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 370 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 371 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 372 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 373 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 374 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 375 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 376 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 377 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 378 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 379 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 380 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 381 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 382 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 383 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 384 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 385 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 386 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 387 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 388 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 389 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 390 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 391 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 392 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 393 |
0.314596738052894]) |
| 394 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 395 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 396 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 397 |
-8.7934289814322 ]) |
| 398 |
|
| 399 |
def Ap(x): |
| 400 |
return matrixmultiply(A,x) |
| 401 |
def Ms(b): |
| 402 |
out=zeros((size(b),),Float64) |
| 403 |
for i in xrange(size(b)): |
| 404 |
out[i]=b[i]/A[i,i] |
| 405 |
return out |
| 406 |
|
| 407 |
tol=1.e-4 |
| 408 |
x,r=PCG(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium,x=x_ref*1.5, iter_max=12) |
| 409 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 410 |
self.failUnless(Lsup(r-(b-matrixmultiply(A,x)))<=Lsup(b)*EPSILON*100.,"wrong solution") |
| 411 |
|
| 412 |
def testGMRES(self): |
| 413 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 414 |
from math import sqrt |
| 415 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 416 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 417 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 418 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 419 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 420 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 421 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 422 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 423 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 424 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 425 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 426 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 427 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 428 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 429 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 430 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 431 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 432 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 433 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 434 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 435 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 436 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 437 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 438 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 439 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 440 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 441 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 442 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 443 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 444 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 445 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 446 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 447 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 448 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 449 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 450 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 451 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 452 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 453 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 454 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 455 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 456 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 457 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 458 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 459 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 460 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 461 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 462 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 463 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 464 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 465 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 466 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 467 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 468 |
0.314596738052894]) |
| 469 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 470 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 471 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 472 |
-8.7934289814322 ]) |
| 473 |
|
| 474 |
def Ap(x): |
| 475 |
return matrixmultiply(A,x) |
| 476 |
def Ms(b): |
| 477 |
out=zeros((size(b),),Float64) |
| 478 |
for i in xrange(size(b)): |
| 479 |
out[i]=b[i]/A[i,i] |
| 480 |
return out |
| 481 |
|
| 482 |
tol=1.e-4 |
| 483 |
x=GMRES(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium2,x=x_ref*1.5, iter_max=12) |
| 484 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 485 |
|
| 486 |
def testArithmeticTuple(self): |
| 487 |
a=ArithmeticTuple(1.,2.) |
| 488 |
self.failUnless(len(a)==2,"wrong length") |
| 489 |
self.failUnless(a[0]==1.,"wrong first item") |
| 490 |
self.failUnless(a[1]==2.,"wrong second item") |
| 491 |
c=a*6. |
| 492 |
self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple") |
| 493 |
self.failUnless(len(c)==2,"c has wrong length") |
| 494 |
self.failUnless(c[0]==6.,"c has wrong first item") |
| 495 |
self.failUnless(c[1]==12.,"c has wrong second item") |
| 496 |
b=5.*a |
| 497 |
self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple") |
| 498 |
self.failUnless(len(b)==2,"b has wrong length") |
| 499 |
self.failUnless(b[0]==5.,"b has wrong first item") |
| 500 |
self.failUnless(b[1]==10.,"b has wrong second item") |
| 501 |
a+=ArithmeticTuple(3.,4.) |
| 502 |
self.failUnless(a[0]==4.,"wrong first item of inplace update") |
| 503 |
self.failUnless(a[1]==6.,"wrong second item of inplace update") |
| 504 |
|
| 505 |
|
| 506 |
|
| 507 |
class Test_pdetools(Test_pdetools_noLumping): |
| 508 |
def testProjector_rank0_fast_reduced(self): |
| 509 |
x=ContinuousFunction(self.domain).getX() |
| 510 |
h=Lsup(self.domain.getSize()) |
| 511 |
p=Projector(self.domain,reduce=True,fast=True) |
| 512 |
td_ref=x[0] |
| 513 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 514 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 515 |
|
| 516 |
def testProjector_rank1_fast_reduced(self): |
| 517 |
x=ContinuousFunction(self.domain).getX() |
| 518 |
h=Lsup(self.domain.getSize()) |
| 519 |
p=Projector(self.domain,reduce=True,fast=True) |
| 520 |
td_ref=x |
| 521 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 522 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 523 |
|
| 524 |
def testProjector_rank2_fast_reduced(self): |
| 525 |
x=ContinuousFunction(self.domain).getX() |
| 526 |
h=Lsup(self.domain.getSize()) |
| 527 |
p=Projector(self.domain,reduce=True,fast=True) |
| 528 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 529 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 530 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 531 |
|
| 532 |
def testProjector_rank3_fast_reduced(self): |
| 533 |
x=ContinuousFunction(self.domain).getX() |
| 534 |
h=Lsup(self.domain.getSize()) |
| 535 |
p=Projector(self.domain,reduce=True,fast=True) |
| 536 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 537 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 538 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 539 |
|
| 540 |
def testProjector_rank4_fast_reduced(self): |
| 541 |
x=ContinuousFunction(self.domain).getX() |
| 542 |
h=Lsup(self.domain.getSize()) |
| 543 |
p=Projector(self.domain,reduce=True,fast=True) |
| 544 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 545 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 546 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 547 |
|
| 548 |
def testProjector_rank0_fast_reduced_with_reduced_input(self): |
| 549 |
x=ContinuousFunction(self.domain).getX() |
| 550 |
h=Lsup(self.domain.getSize()) |
| 551 |
p=Projector(self.domain,reduce=True,fast=True) |
| 552 |
td_ref=1. |
| 553 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 554 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 555 |
|
| 556 |
def testProjector_rank1_fast_reduced_with_reduced_input(self): |
| 557 |
x=ContinuousFunction(self.domain).getX() |
| 558 |
h=Lsup(self.domain.getSize()) |
| 559 |
p=Projector(self.domain,reduce=True,fast=True) |
| 560 |
td_ref=numarray.array([1.,2.,3.]) |
| 561 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 562 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 563 |
|
| 564 |
def testProjector_rank2_fast_reduced_with_reduced_input(self): |
| 565 |
x=ContinuousFunction(self.domain).getX() |
| 566 |
h=Lsup(self.domain.getSize()) |
| 567 |
p=Projector(self.domain,reduce=True,fast=True) |
| 568 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 569 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 570 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 571 |
|
| 572 |
def testProjector_rank3_fast_reduced_with_reduced_input(self): |
| 573 |
x=ContinuousFunction(self.domain).getX() |
| 574 |
h=Lsup(self.domain.getSize()) |
| 575 |
p=Projector(self.domain,reduce=True,fast=True) |
| 576 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 577 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 578 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 579 |
|
| 580 |
def testProjector_rank4_fast_reduced_with_reduced_input(self): |
| 581 |
x=ContinuousFunction(self.domain).getX() |
| 582 |
h=Lsup(self.domain.getSize()) |
| 583 |
p=Projector(self.domain,reduce=True,fast=True) |
| 584 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 585 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 586 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 587 |
|
Parent Directory
| 
Revision Log