| 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, IterationHistory, ArithmeticTuple, GMRES |
| 51 |
from esys.escript.pdetools import Defect, NewtonGMRES |
| 52 |
|
| 53 |
class Test_pdetools_noLumping(unittest.TestCase): |
| 54 |
DEBUG=False |
| 55 |
VERBOSE=False |
| 56 |
def test_TimeIntegrationManager_scalar(self): |
| 57 |
t=0. |
| 58 |
dt=0.1 |
| 59 |
tm=TimeIntegrationManager(0.,p=1) |
| 60 |
while t<1.: |
| 61 |
t+=dt |
| 62 |
tm.checkin(dt,t) |
| 63 |
v_guess=tm.extrapolate(dt) |
| 64 |
self.failUnless(abs(v_guess-(tm.getTime()+dt))<self.RES_TOL,"extrapolation is wrong") |
| 65 |
|
| 66 |
def test_TimeIntegrationManager_vector(self): |
| 67 |
t=0. |
| 68 |
dt=0.3 |
| 69 |
tm=TimeIntegrationManager(0.,0.,p=1) |
| 70 |
while t<1.: |
| 71 |
t+=dt |
| 72 |
tm.checkin(dt,t,3*t) |
| 73 |
v_guess=tm.extrapolate(dt) |
| 74 |
e=max(abs(v_guess[0]-(tm.getTime()+dt)),abs(v_guess[1]-(tm.getTime()+dt)*3.)) |
| 75 |
self.failUnless(e<self.RES_TOL,"extrapolation is wrong") |
| 76 |
|
| 77 |
def test_Locator(self): |
| 78 |
x=self.domain.getX() |
| 79 |
l=Locator(self.domain,numarray.ones((self.domain.getDim(),))) |
| 80 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 81 |
|
| 82 |
l=Locator(ContinuousFunction(self.domain),numarray.ones((self.domain.getDim(),))) |
| 83 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 84 |
|
| 85 |
xx=l.getX() |
| 86 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 87 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"location wrong") |
| 88 |
xx=l(x) |
| 89 |
self.failUnless(isinstance(xx,numarray.NumArray),"wrong vector type") |
| 90 |
self.failUnless(Lsup(xx-numarray.ones((self.domain.getDim(),)))<self.RES_TOL,"value wrong vector") |
| 91 |
xx=l(x[0]+x[1]) |
| 92 |
self.failUnless(isinstance(xx,float),"wrong scalar type") |
| 93 |
self.failUnless(abs(xx-2.)<self.RES_TOL,"value wrong scalar") |
| 94 |
|
| 95 |
def test_Locator_withList(self): |
| 96 |
x=self.domain.getX() |
| 97 |
arg=[numarray.ones((self.domain.getDim(),)), numarray.zeros((self.domain.getDim(),))] |
| 98 |
l=Locator(self.domain,arg) |
| 99 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space from domain") |
| 100 |
|
| 101 |
l=Locator(ContinuousFunction(self.domain),arg) |
| 102 |
self.failUnless(ContinuousFunction(self.domain)==l.getFunctionSpace(),"wrong function space") |
| 103 |
|
| 104 |
xx=l.getX() |
| 105 |
self.failUnless(isinstance(xx,list),"list expected") |
| 106 |
for i in range(len(xx)): |
| 107 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item"%i) |
| 108 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong"%i) |
| 109 |
xx=l(x) |
| 110 |
self.failUnless(isinstance(xx,list),"list expected (2)") |
| 111 |
for i in range(len(xx)): |
| 112 |
self.failUnless(isinstance(xx[i],numarray.NumArray),"vector expected for %s item (2)"%i) |
| 113 |
self.failUnless(Lsup(xx[i]-arg[i])<self.RES_TOL,"%s-th location is wrong (2)"%i) |
| 114 |
xx=l(x[0]+x[1]) |
| 115 |
self.failUnless(isinstance(xx,list),"list expected (3)") |
| 116 |
for i in range(len(xx)): |
| 117 |
self.failUnless(isinstance(xx[i],float),"wrong scalar type") |
| 118 |
self.failUnless(abs(xx[i]-(arg[i][0]+arg[i][1]))<self.RES_TOL,"value wrong scalar") |
| 119 |
|
| 120 |
def testProjector_rank0(self): |
| 121 |
x=ContinuousFunction(self.domain).getX() |
| 122 |
p=Projector(self.domain,reduce=False,fast=False) |
| 123 |
td_ref=x[0] |
| 124 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 125 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 126 |
|
| 127 |
def testProjector_rank1(self): |
| 128 |
x=ContinuousFunction(self.domain).getX() |
| 129 |
p=Projector(self.domain,reduce=False,fast=False) |
| 130 |
td_ref=x |
| 131 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 132 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 133 |
|
| 134 |
def testProjector_rank2(self): |
| 135 |
x=ContinuousFunction(self.domain).getX() |
| 136 |
p=Projector(self.domain,reduce=False,fast=False) |
| 137 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 138 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 139 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 140 |
|
| 141 |
def testProjector_rank3(self): |
| 142 |
x=ContinuousFunction(self.domain).getX() |
| 143 |
p=Projector(self.domain,reduce=False,fast=False) |
| 144 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 145 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 146 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 147 |
|
| 148 |
def testProjector_rank4(self): |
| 149 |
x=ContinuousFunction(self.domain).getX() |
| 150 |
p=Projector(self.domain,reduce=False,fast=False) |
| 151 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 152 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 153 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 154 |
|
| 155 |
|
| 156 |
def testProjector_rank0_reduced(self): |
| 157 |
x=ContinuousFunction(self.domain).getX() |
| 158 |
p=Projector(self.domain,reduce=True,fast=False) |
| 159 |
td_ref=x[0] |
| 160 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 161 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 162 |
|
| 163 |
def testProjector_rank1_reduced(self): |
| 164 |
x=ContinuousFunction(self.domain).getX() |
| 165 |
p=Projector(self.domain,reduce=True,fast=False) |
| 166 |
td_ref=x |
| 167 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 168 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 169 |
|
| 170 |
def testProjector_rank2_reduced(self): |
| 171 |
x=ContinuousFunction(self.domain).getX() |
| 172 |
p=Projector(self.domain,reduce=True,fast=False) |
| 173 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 174 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 175 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 176 |
|
| 177 |
def testProjector_rank3_reduced(self): |
| 178 |
x=ContinuousFunction(self.domain).getX() |
| 179 |
p=Projector(self.domain,reduce=True,fast=False) |
| 180 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 181 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 182 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 183 |
|
| 184 |
def testProjector_rank4_reduced(self): |
| 185 |
x=ContinuousFunction(self.domain).getX() |
| 186 |
p=Projector(self.domain,reduce=True,fast=False) |
| 187 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 188 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 189 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 190 |
|
| 191 |
def testProjector_rank0_with_reduced_input(self): |
| 192 |
x=ContinuousFunction(self.domain).getX() |
| 193 |
p=Projector(self.domain,reduce=False,fast=False) |
| 194 |
td_ref=x[0] |
| 195 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 196 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 197 |
|
| 198 |
def testProjector_rank1_with_reduced_input(self): |
| 199 |
x=ContinuousFunction(self.domain).getX() |
| 200 |
p=Projector(self.domain,reduce=False,fast=False) |
| 201 |
td_ref=x |
| 202 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 203 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 204 |
|
| 205 |
def testProjector_rank2_with_reduced_input(self): |
| 206 |
x=ContinuousFunction(self.domain).getX() |
| 207 |
p=Projector(self.domain,reduce=False,fast=False) |
| 208 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 209 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 210 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 211 |
|
| 212 |
def testProjector_rank3_with_reduced_input(self): |
| 213 |
x=ContinuousFunction(self.domain).getX() |
| 214 |
p=Projector(self.domain,reduce=False,fast=False) |
| 215 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 216 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 217 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 218 |
|
| 219 |
def testProjector_rank4_with_reduced_input(self): |
| 220 |
x=ContinuousFunction(self.domain).getX() |
| 221 |
p=Projector(self.domain,reduce=False,fast=False) |
| 222 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 223 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 224 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 225 |
|
| 226 |
|
| 227 |
def testProjector_rank0_reduced_with_reduced_input(self): |
| 228 |
x=ContinuousFunction(self.domain).getX() |
| 229 |
p=Projector(self.domain,reduce=True,fast=False) |
| 230 |
td_ref=1. |
| 231 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 232 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 233 |
|
| 234 |
def testProjector_rank1_reduced_with_reduced_input(self): |
| 235 |
x=ContinuousFunction(self.domain).getX() |
| 236 |
p=Projector(self.domain,reduce=True,fast=False) |
| 237 |
td_ref=numarray.array([1.,2.,3.]) |
| 238 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 239 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 240 |
|
| 241 |
def testProjector_rank2_reduced_with_reduced_input(self): |
| 242 |
x=ContinuousFunction(self.domain).getX() |
| 243 |
p=Projector(self.domain,reduce=True,fast=False) |
| 244 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 245 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 246 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 247 |
|
| 248 |
def testProjector_rank3_reduced_with_reduced_input(self): |
| 249 |
x=ContinuousFunction(self.domain).getX() |
| 250 |
p=Projector(self.domain,reduce=True,fast=False) |
| 251 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 252 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 253 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 254 |
|
| 255 |
def testProjector_rank4_reduced_with_reduced_input(self): |
| 256 |
x=ContinuousFunction(self.domain).getX() |
| 257 |
p=Projector(self.domain,reduce=True,fast=False) |
| 258 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 259 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 260 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*self.RES_TOL,"value wrong") |
| 261 |
|
| 262 |
|
| 263 |
def test_NoPDE_scalar_missing_r(self): |
| 264 |
p=NoPDE(self.domain) |
| 265 |
x=self.domain.getX() |
| 266 |
msk=whereZero(x[0]) |
| 267 |
p.setValue(D=1.,Y=1.,q=msk) |
| 268 |
u=p.getSolution() |
| 269 |
u_ex=(1.-msk) |
| 270 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 271 |
|
| 272 |
def test_NoPDE_scalar_missing_Y(self): |
| 273 |
p=NoPDE(self.domain) |
| 274 |
x=self.domain.getX() |
| 275 |
msk=whereZero(x[0]) |
| 276 |
p.setValue(D=1.,q=msk,r=2.) |
| 277 |
u=p.getSolution() |
| 278 |
u_ex=msk*2. |
| 279 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 280 |
|
| 281 |
def test_NoPDE_scalar_constant(self): |
| 282 |
p=NoPDE(self.domain) |
| 283 |
x=self.domain.getX() |
| 284 |
msk=whereZero(x[0]) |
| 285 |
p.setValue(D=1.,Y=1.,q=msk,r=2.) |
| 286 |
u=p.getSolution() |
| 287 |
u_ex=(1.-msk)+msk*2. |
| 288 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 289 |
|
| 290 |
def test_NoPDE_scalar_variable(self): |
| 291 |
p=NoPDE(self.domain) |
| 292 |
x=self.domain.getX() |
| 293 |
msk=whereZero(x[0]) |
| 294 |
p.setValue(D=x[0],Y=2*x[0],q=msk,r=2.) |
| 295 |
u=p.getSolution() |
| 296 |
u_ex=2. |
| 297 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 298 |
|
| 299 |
def test_NoPDE_vector_missing_Y(self): |
| 300 |
p=NoPDE(self.domain) |
| 301 |
x=self.domain.getX() |
| 302 |
msk=whereZero(x[0])*[1.,0.] |
| 303 |
p.setValue(D=numarray.ones([2]),q=msk,r=2.) |
| 304 |
u=p.getSolution() |
| 305 |
u_ex=msk*2. |
| 306 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 307 |
|
| 308 |
def test_NoPDE_vector_missing_r(self): |
| 309 |
p=NoPDE(self.domain) |
| 310 |
x=self.domain.getX() |
| 311 |
msk=whereZero(x[0])*[1.,0.] |
| 312 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk) |
| 313 |
u=p.getSolution() |
| 314 |
u_ex=(1.-msk) |
| 315 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 316 |
|
| 317 |
def test_NoPDE_vector_constant(self): |
| 318 |
p=NoPDE(self.domain) |
| 319 |
x=self.domain.getX() |
| 320 |
msk=whereZero(x[0])*[1.,0.] |
| 321 |
p.setValue(D=numarray.ones([2]),Y=numarray.ones([2]),q=msk,r=2.) |
| 322 |
u=p.getSolution() |
| 323 |
u_ex=(1.-msk)+msk*2. |
| 324 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 325 |
|
| 326 |
def test_NoPDE_vector_variable(self): |
| 327 |
p=NoPDE(self.domain) |
| 328 |
x=self.domain.getX() |
| 329 |
msk=whereZero(x[0])*[1.,0.] |
| 330 |
p.setValue(D=x[:2]+1,Y=2*(x[:2]+1),q=msk,r=2.) |
| 331 |
u=p.getSolution() |
| 332 |
u_ex=2. |
| 333 |
self.failUnless(Lsup(u_ex-u)<Lsup(u_ex)*self.RES_TOL,"value wrong") |
| 334 |
|
| 335 |
def testPCG(self): |
| 336 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 337 |
from math import sqrt |
| 338 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 339 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 340 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 341 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 342 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 343 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 344 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 345 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 346 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 347 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 348 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 349 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 350 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 351 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 352 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 353 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 354 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 355 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 356 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 357 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 358 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 359 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 360 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 361 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 362 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 363 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 364 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 365 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 366 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 367 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 368 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 369 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 370 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 371 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 372 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 373 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 374 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 375 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 376 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 377 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 378 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 379 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 380 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 381 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 382 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 383 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 384 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 385 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 386 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 387 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 388 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 389 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 390 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 391 |
0.314596738052894]) |
| 392 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 393 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 394 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 395 |
-8.7934289814322 ]) |
| 396 |
|
| 397 |
def Ap(x): |
| 398 |
return matrixmultiply(A,x) |
| 399 |
def Ms(b): |
| 400 |
out=zeros((size(b),),Float64) |
| 401 |
for i in xrange(size(b)): |
| 402 |
out[i]=b[i]/A[i,i] |
| 403 |
return out |
| 404 |
|
| 405 |
tol=1.e-4 |
| 406 |
x,r=PCG(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium,x=x_ref*1.5, iter_max=12) |
| 407 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 408 |
self.failUnless(Lsup(r-(b-matrixmultiply(A,x)))<=Lsup(b)*EPSILON*100.,"wrong solution") |
| 409 |
|
| 410 |
def testGMRES(self): |
| 411 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 412 |
from math import sqrt |
| 413 |
A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 414 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 415 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 416 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 417 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 418 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 419 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 420 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 421 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 422 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 423 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 424 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 425 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 426 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 427 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 428 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 429 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 430 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 431 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 432 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 433 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 434 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 435 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 436 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 437 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 438 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 439 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 440 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 441 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 442 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 443 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 444 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 445 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 446 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 447 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 448 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 449 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 450 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 451 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 452 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 453 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 454 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 455 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 456 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 457 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 458 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 459 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 460 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 461 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 462 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 463 |
x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 464 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 465 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 466 |
0.314596738052894]) |
| 467 |
b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 468 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 469 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 470 |
-8.7934289814322 ]) |
| 471 |
|
| 472 |
def Ap(x): |
| 473 |
return matrixmultiply(A,x) |
| 474 |
def Ms(b): |
| 475 |
out=zeros((size(b),),Float64) |
| 476 |
for i in xrange(size(b)): |
| 477 |
out[i]=b[i]/A[i,i] |
| 478 |
return out |
| 479 |
|
| 480 |
tol=1.e-4 |
| 481 |
x=GMRES(b*1.,Ap,Ms,dot, IterationHistory(tol).stoppingcriterium2,x=x_ref*1.5, iter_max=12) |
| 482 |
self.failUnless(Lsup(x-x_ref)<=Lsup(x_ref)*tol*10.,"wrong solution") |
| 483 |
|
| 484 |
def testNewtonGMRES(self): |
| 485 |
from numarray import array,matrixmultiply, zeros, dot, size, Float64 |
| 486 |
from math import sqrt |
| 487 |
class LL(Defect): |
| 488 |
def __init__(self,*kwargs): |
| 489 |
super(LL, self).__init__(*kwargs) |
| 490 |
self.A=array([[ 4.752141253159452e+02, -2.391895572674098e-01, |
| 491 |
5.834798554135237e-01, -3.704394311709722e+00, |
| 492 |
5.765369186984777e+00, -1.309786358737351e+01, |
| 493 |
2.522087134507148e+01, -3.393956279045637e+01, |
| 494 |
1.046856914770830e+02, -2.447764190849540e+02], |
| 495 |
[ -2.391895572674098e-01, 1.256797283910693e+02, |
| 496 |
-9.188270412920813e-01, 1.300169538880688e+00, |
| 497 |
-5.353714719231424e-01, 2.674709444667012e+00, |
| 498 |
-1.116097841269580e+01, 2.801193427514478e+01, |
| 499 |
-3.877806125898224e+01, 3.063505753648256e+01], |
| 500 |
[ 5.834798554135237e-01, -9.188270412920813e-01, |
| 501 |
6.240841811806843e+01, -8.176289504109282e-01, |
| 502 |
1.447935098417076e-01, -9.721424148655324e-01, |
| 503 |
6.713551574117577e-01, -3.656297654168375e+00, |
| 504 |
7.015141656913973e+00, -4.195525932156250e+01], |
| 505 |
[ -3.704394311709722e+00, 1.300169538880688e+00, |
| 506 |
-8.176289504109282e-01, 3.604980536782198e+01, |
| 507 |
-6.241238423759328e-01, 1.142345320047869e+00, |
| 508 |
-3.438816797096519e+00, 5.854857481367470e+00, |
| 509 |
-4.524311288596452e+00, 1.136590280389803e+01], |
| 510 |
[ 5.765369186984777e+00, -5.353714719231424e-01, |
| 511 |
1.447935098417076e-01, -6.241238423759328e-01, |
| 512 |
2.953997190215862e+01, -9.474729233464712e-01, |
| 513 |
1.883516378345809e+00, -1.906274765704230e+00, |
| 514 |
4.401859671778645e+00, -1.064573816075257e+01], |
| 515 |
[ -1.309786358737351e+01, 2.674709444667012e+00, |
| 516 |
-9.721424148655324e-01, 1.142345320047869e+00, |
| 517 |
-9.474729233464712e-01, 2.876998216302979e+01, |
| 518 |
-4.853065259692995e-01, 7.088596468102618e-01, |
| 519 |
-8.972224295152829e-01, 5.228606946522749e+00], |
| 520 |
[ 2.522087134507148e+01, -1.116097841269580e+01, |
| 521 |
6.713551574117577e-01, -3.438816797096519e+00, |
| 522 |
1.883516378345809e+00, -4.853065259692995e-01, |
| 523 |
5.121175860935919e+01, -3.523133115905478e-01, |
| 524 |
1.782136702229135e+00, -1.560849559916187e+00], |
| 525 |
[ -3.393956279045637e+01, 2.801193427514478e+01, |
| 526 |
-3.656297654168375e+00, 5.854857481367470e+00, |
| 527 |
-1.906274765704230e+00, 7.088596468102618e-01, |
| 528 |
-3.523133115905478e-01, 8.411681423853814e+01, |
| 529 |
-5.238590858177903e-01, 1.515872114883926e+00], |
| 530 |
[ 1.046856914770830e+02, -3.877806125898224e+01, |
| 531 |
7.015141656913973e+00, -4.524311288596452e+00, |
| 532 |
4.401859671778645e+00, -8.972224295152829e-01, |
| 533 |
1.782136702229135e+00, -5.238590858177903e-01, |
| 534 |
1.797889693808014e+02, -8.362340479938084e-01], |
| 535 |
[ -2.447764190849540e+02, 3.063505753648256e+01, |
| 536 |
-4.195525932156250e+01, 1.136590280389803e+01, |
| 537 |
-1.064573816075257e+01, 5.228606946522749e+00, |
| 538 |
-1.560849559916187e+00, 1.515872114883926e+00, |
| 539 |
-8.362340479938084e-01, 3.833719335346630e+02]]) |
| 540 |
self.x_ref=array([ 0.41794207085296, 0.031441086046563, 0.882801683420401, |
| 541 |
0.807186823427233, 0.48950999450145, 0.995486532098031, |
| 542 |
0.351243009576568, 0.704352576819321, 0.850648989740204, |
| 543 |
0.314596738052894]) |
| 544 |
self.b=array([ 182.911023960262952, -1.048322041992754, 44.181293875206201, |
| 545 |
30.344553414038817, 15.247917439094513, 24.060664905403492, |
| 546 |
27.210293789825833, 47.122067744075842, 199.267136417856847, |
| 547 |
-8.7934289814322 ]) |
| 548 |
def eval(self,x): |
| 549 |
out=matrixmultiply(self.A,x)-self.b |
| 550 |
for i in xrange(size(self.b)): |
| 551 |
out[i]/=self.A[i,i] |
| 552 |
return out |
| 553 |
def bilinearform(self,x0,x1): |
| 554 |
return dot(x0,x1) |
| 555 |
|
| 556 |
tol=1.e-8 |
| 557 |
ll=LL() |
| 558 |
x=NewtonGMRES(LL(),ll.x_ref*0., iter_max=100, sub_iter_max=20, atol=0,rtol=tol, verbose=self.VERBOSE) |
| 559 |
self.failUnless(Lsup(x-ll.x_ref)<=Lsup(ll.x_ref)*tol*10.,"wrong solution") |
| 560 |
|
| 561 |
def testArithmeticTuple(self): |
| 562 |
a=ArithmeticTuple(1.,2.) |
| 563 |
self.failUnless(len(a)==2,"wrong length") |
| 564 |
self.failUnless(a[0]==1.,"wrong first item") |
| 565 |
self.failUnless(a[1]==2.,"wrong second item") |
| 566 |
c=a*6. |
| 567 |
self.failUnless(isinstance(c,ArithmeticTuple),"c is not an instance of ArithmeticTuple") |
| 568 |
self.failUnless(len(c)==2,"c has wrong length") |
| 569 |
self.failUnless(c[0]==6.,"c has wrong first item") |
| 570 |
self.failUnless(c[1]==12.,"c has wrong second item") |
| 571 |
b=5.*a |
| 572 |
self.failUnless(isinstance(b,ArithmeticTuple),"b is not an instance of ArithmeticTuple") |
| 573 |
self.failUnless(len(b)==2,"b has wrong length") |
| 574 |
self.failUnless(b[0]==5.,"b has wrong first item") |
| 575 |
self.failUnless(b[1]==10.,"b has wrong second item") |
| 576 |
a+=ArithmeticTuple(3.,4.) |
| 577 |
self.failUnless(a[0]==4.,"wrong first item of inplace update") |
| 578 |
self.failUnless(a[1]==6.,"wrong second item of inplace update") |
| 579 |
|
| 580 |
|
| 581 |
|
| 582 |
class Test_pdetools(Test_pdetools_noLumping): |
| 583 |
def testProjector_rank0_fast_reduced(self): |
| 584 |
x=ContinuousFunction(self.domain).getX() |
| 585 |
h=Lsup(self.domain.getSize()) |
| 586 |
p=Projector(self.domain,reduce=True,fast=True) |
| 587 |
td_ref=x[0] |
| 588 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 589 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 590 |
|
| 591 |
def testProjector_rank1_fast_reduced(self): |
| 592 |
x=ContinuousFunction(self.domain).getX() |
| 593 |
h=Lsup(self.domain.getSize()) |
| 594 |
p=Projector(self.domain,reduce=True,fast=True) |
| 595 |
td_ref=x |
| 596 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 597 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 598 |
|
| 599 |
def testProjector_rank2_fast_reduced(self): |
| 600 |
x=ContinuousFunction(self.domain).getX() |
| 601 |
h=Lsup(self.domain.getSize()) |
| 602 |
p=Projector(self.domain,reduce=True,fast=True) |
| 603 |
td_ref=[[11.,12.],[21,22.]]*(x[0]+x[1]) |
| 604 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 605 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 606 |
|
| 607 |
def testProjector_rank3_fast_reduced(self): |
| 608 |
x=ContinuousFunction(self.domain).getX() |
| 609 |
h=Lsup(self.domain.getSize()) |
| 610 |
p=Projector(self.domain,reduce=True,fast=True) |
| 611 |
td_ref=[[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]*(x[0]+x[1]) |
| 612 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 613 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 614 |
|
| 615 |
def testProjector_rank4_fast_reduced(self): |
| 616 |
x=ContinuousFunction(self.domain).getX() |
| 617 |
h=Lsup(self.domain.getSize()) |
| 618 |
p=Projector(self.domain,reduce=True,fast=True) |
| 619 |
td_ref=[[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]*(x[0]+x[1]) |
| 620 |
td=p(td_ref.interpolate(Function(self.domain))) |
| 621 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 622 |
|
| 623 |
def testProjector_rank0_fast_reduced_with_reduced_input(self): |
| 624 |
x=ContinuousFunction(self.domain).getX() |
| 625 |
h=Lsup(self.domain.getSize()) |
| 626 |
p=Projector(self.domain,reduce=True,fast=True) |
| 627 |
td_ref=1. |
| 628 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 629 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 630 |
|
| 631 |
def testProjector_rank1_fast_reduced_with_reduced_input(self): |
| 632 |
x=ContinuousFunction(self.domain).getX() |
| 633 |
h=Lsup(self.domain.getSize()) |
| 634 |
p=Projector(self.domain,reduce=True,fast=True) |
| 635 |
td_ref=numarray.array([1.,2.,3.]) |
| 636 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 637 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 638 |
|
| 639 |
def testProjector_rank2_fast_reduced_with_reduced_input(self): |
| 640 |
x=ContinuousFunction(self.domain).getX() |
| 641 |
h=Lsup(self.domain.getSize()) |
| 642 |
p=Projector(self.domain,reduce=True,fast=True) |
| 643 |
td_ref=numarray.array([[11.,12.],[21,22.]]) |
| 644 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 645 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 646 |
|
| 647 |
def testProjector_rank3_fast_reduced_with_reduced_input(self): |
| 648 |
x=ContinuousFunction(self.domain).getX() |
| 649 |
h=Lsup(self.domain.getSize()) |
| 650 |
p=Projector(self.domain,reduce=True,fast=True) |
| 651 |
td_ref=numarray.array([[[111.,112.],[121,122.]],[[211.,212.],[221,222.]]]) |
| 652 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 653 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 654 |
|
| 655 |
def testProjector_rank4_fast_reduced_with_reduced_input(self): |
| 656 |
x=ContinuousFunction(self.domain).getX() |
| 657 |
h=Lsup(self.domain.getSize()) |
| 658 |
p=Projector(self.domain,reduce=True,fast=True) |
| 659 |
td_ref=numarray.array([[[[1111.,1112.],[1121,1122.]],[[1211.,1212.],[1221,1222.]]],[[[2111.,2112.],[2121,2122.]],[[2211.,2212.],[2221,2222.]]]]) |
| 660 |
td=p(Data(td_ref,ReducedFunction(self.domain))) |
| 661 |
self.failUnless(Lsup(td-td_ref)<Lsup(td_ref)*h,"value wrong") |
| 662 |
|
Parent Directory
| 
Revision Log