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 linearPDE and pdetools test on finley |
19 |
|
20 |
@remark: |
21 |
|
22 |
@var __author__: name of author |
23 |
@var __licence__: licence agreement |
24 |
@var __url__: url entry point on documentation |
25 |
@var __version__: version |
26 |
@var __date__: date of the version |
27 |
""" |
28 |
|
29 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
30 |
http://www.access.edu.au |
31 |
Primary Business: Queensland, Australia""" |
32 |
__license__="""Licensed under the Open Software License version 3.0 |
33 |
http://www.opensource.org/licenses/osl-3.0.php""" |
34 |
__author__="Lutz Gross, l.gross@uq.edu.au" |
35 |
__url__="http://www.iservo.edu.au/esys/escript" |
36 |
__version__="$Revision: 859 $" |
37 |
__date__="$Date: 2006-09-26 12:19:18 +1000 (Tue, 26 Sep 2006) $" |
38 |
|
39 |
|
40 |
import unittest |
41 |
|
42 |
from esys.escript import * |
43 |
from esys.finley import Rectangle,Brick |
44 |
from esys.escript.linearPDEs import LinearPDE |
45 |
OPTIMIZE=False |
46 |
SOLVER_VERBOSE=False |
47 |
|
48 |
try: |
49 |
FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA'] |
50 |
except KeyError: |
51 |
FINLEY_TEST_DATA='.' |
52 |
|
53 |
FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/" |
54 |
|
55 |
# number of elements in the spatial directions |
56 |
NE0=8 |
57 |
NE1=10 |
58 |
NE2=12 |
59 |
|
60 |
NE0=12 |
61 |
NE1=12 |
62 |
NE2=8 |
63 |
|
64 |
SOLVER_TOL=1.e-8 |
65 |
REL_TOL=1.e-6 |
66 |
|
67 |
FAC_DIAG=1. |
68 |
FAC_OFFDIAG=-0.4 |
69 |
|
70 |
class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase): |
71 |
def test_solve(self): |
72 |
domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE) |
73 |
x=Solution(domain).getX() |
74 |
# --- set exact solution ---- |
75 |
u_ex=Scalar(0,Solution(domain)) |
76 |
u_ex=1.+2.*x[0]+3.*x[1] |
77 |
# --- set exact gradient ----------- |
78 |
g_ex=Data(0.,(2,),Solution(domain)) |
79 |
g_ex[0]=2. |
80 |
g_ex[1]=3. |
81 |
# -------- test gradient -------------------------------- |
82 |
g=grad(u_ex) |
83 |
self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex)) |
84 |
# -------- set-up PDE ----------------------------------- |
85 |
pde=LinearPDE(domain,numEquations=1) |
86 |
mask=whereZero(x[0]) |
87 |
pde.setValue(r=u_ex,q=mask) |
88 |
pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal())) |
89 |
# -------- get the solution --------------------------- |
90 |
pde.setTolerance(SOLVER_TOL) |
91 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
92 |
pde.setSolverPackage(pde.PASO) |
93 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
94 |
# -------- test the solution --------------------------- |
95 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
96 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
97 |
class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase): |
98 |
def test_solve(self): |
99 |
domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE) |
100 |
x=Solution(domain).getX() |
101 |
# --- set exact solution ---- |
102 |
u_ex=Vector(0,Solution(domain)) |
103 |
u_ex[0]=1.+2.*x[0]+3.*x[1] |
104 |
u_ex[1]=-1.+3.*x[0]+2.*x[1] |
105 |
# --- set exact gradient ----------- |
106 |
g_ex=Data(0.,(2,2),Solution(domain)) |
107 |
g_ex[0,0]=2. |
108 |
g_ex[0,1]=3. |
109 |
g_ex[1,0]=3. |
110 |
g_ex[1,1]=2. |
111 |
# -------- test gradient -------------------------------- |
112 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
113 |
# -------- set-up PDE ----------------------------------- |
114 |
pde=LinearPDE(domain,numEquations=2) |
115 |
mask=whereZero(x[0]) |
116 |
pde.setValue(r=u_ex,q=mask*numarray.ones(2,)) |
117 |
A=Tensor4(0,Function(domain)) |
118 |
A[0,:,0,:]=kronecker(2) |
119 |
A[1,:,1,:]=kronecker(2) |
120 |
Y=Vector(0.,Function(domain)) |
121 |
Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG |
122 |
Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG |
123 |
pde.setValue(A=A, |
124 |
D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG, |
125 |
Y=Y, |
126 |
y=matrixmult(g_ex,domain.getNormal())) |
127 |
# -------- get the solution --------------------------- |
128 |
pde.setTolerance(SOLVER_TOL) |
129 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
130 |
pde.setSolverPackage(pde.PASO) |
131 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
132 |
# -------- test the solution --------------------------- |
133 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
134 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
135 |
class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase): |
136 |
def test_solve(self): |
137 |
domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE) |
138 |
x=Solution(domain).getX() |
139 |
# --- set exact solution ---- |
140 |
u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2 |
141 |
# --- set exact gradient ----------- |
142 |
g_ex=Data(0.,(2,),Solution(domain)) |
143 |
g_ex[0]=2.+8.*x[0]+5.*x[1] |
144 |
g_ex[1]=3.+5.*x[0]+12.*x[1] |
145 |
# -------- test gradient -------------------------------- |
146 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
147 |
# -------- set-up PDE ----------------------------------- |
148 |
pde=LinearPDE(domain,numEquations=1) |
149 |
mask=whereZero(x[0]) |
150 |
pde.setValue(r=u_ex,q=mask) |
151 |
pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.) |
152 |
# -------- get the solution --------------------------- |
153 |
pde.setTolerance(SOLVER_TOL) |
154 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
155 |
pde.setSolverPackage(pde.PASO) |
156 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
157 |
# -------- test the solution --------------------------- |
158 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
159 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
160 |
class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase): |
161 |
def test_solve(self): |
162 |
domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE) |
163 |
x=Solution(domain).getX() |
164 |
# --- set exact solution ---- |
165 |
u_ex=Vector(0,Solution(domain)) |
166 |
u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2 |
167 |
u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2 |
168 |
# --- set exact gradient ----------- |
169 |
g_ex=Data(0.,(2,2),Solution(domain)) |
170 |
g_ex[0,0]=2.+8.*x[0]+5.*x[1] |
171 |
g_ex[0,1]=3.+5.*x[0]+12.*x[1] |
172 |
g_ex[1,0]=4.+2.*x[0]+6.*x[1] |
173 |
g_ex[1,1]=2.+6.*x[0]+8.*x[1] |
174 |
# -------- test gradient -------------------------------- |
175 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
176 |
# -------- set-up PDE ----------------------------------- |
177 |
pde=LinearPDE(domain,numEquations=2) |
178 |
mask=whereZero(x[0]) |
179 |
pde.setValue(r=u_ex,q=mask*numarray.ones(2,)) |
180 |
A=Tensor4(0,Function(domain)) |
181 |
A[0,:,0,:]=kronecker(2) |
182 |
A[1,:,1,:]=kronecker(2) |
183 |
Y=Vector(0.,Function(domain)) |
184 |
Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG |
185 |
Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG |
186 |
pde.setValue(A=A, |
187 |
D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG, |
188 |
Y=Y-[20.,10.], |
189 |
y=matrixmult(g_ex,domain.getNormal())) |
190 |
# -------- get the solution --------------------------- |
191 |
pde.setTolerance(SOLVER_TOL) |
192 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
193 |
pde.setSolverPackage(pde.PASO) |
194 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
195 |
# -------- test the solution --------------------------- |
196 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
197 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
198 |
class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase): |
199 |
def test_solve(self): |
200 |
domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE) |
201 |
x=Solution(domain).getX() |
202 |
u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2] |
203 |
# --- set exact gradient ----------- |
204 |
g_ex=Data(0.,(3,),Solution(domain)) |
205 |
g_ex[0]=2. |
206 |
g_ex[1]=3. |
207 |
g_ex[2]=4. |
208 |
# -------- test gradient -------------------------------- |
209 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
210 |
# -------- set-up PDE ----------------------------------- |
211 |
pde=LinearPDE(domain,numEquations=1) |
212 |
mask=whereZero(x[0]) |
213 |
pde.setValue(r=u_ex,q=mask) |
214 |
pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal())) |
215 |
# -------- get the solution --------------------------- |
216 |
pde.setTolerance(SOLVER_TOL) |
217 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
218 |
pde.setSolverPackage(pde.PASO) |
219 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
220 |
# -------- test the solution --------------------------- |
221 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
222 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
223 |
class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase): |
224 |
def test_solve(self): |
225 |
domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE) |
226 |
x=Solution(domain).getX() |
227 |
# --- set exact solution ---- |
228 |
u_ex=Vector(0,Solution(domain)) |
229 |
u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2] |
230 |
u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2] |
231 |
u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2] |
232 |
# --- set exact gradient ----------- |
233 |
g_ex=Data(0.,(3,3),Solution(domain)) |
234 |
g_ex[0,0]=2. |
235 |
g_ex[0,1]=3. |
236 |
g_ex[0,2]=4. |
237 |
g_ex[1,0]=4. |
238 |
g_ex[1,1]=1. |
239 |
g_ex[1,2]=-2. |
240 |
g_ex[2,0]=8. |
241 |
g_ex[2,1]=4. |
242 |
g_ex[2,2]=5. |
243 |
# -------- test gradient -------------------------------- |
244 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
245 |
# -------- set-up PDE ----------------------------------- |
246 |
pde=LinearPDE(domain,numEquations=3) |
247 |
mask=whereZero(x[0]) |
248 |
pde.setValue(r=u_ex,q=mask*numarray.ones(3,)) |
249 |
A=Tensor4(0,Function(domain)) |
250 |
A[0,:,0,:]=kronecker(3) |
251 |
A[1,:,1,:]=kronecker(3) |
252 |
A[2,:,2,:]=kronecker(3) |
253 |
Y=Vector(0.,Function(domain)) |
254 |
Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG |
255 |
Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG |
256 |
Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG |
257 |
pde.setValue(A=A, |
258 |
D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG, |
259 |
Y=Y, |
260 |
y=matrixmult(g_ex,domain.getNormal())) |
261 |
# -------- get the solution --------------------------- |
262 |
pde.setTolerance(SOLVER_TOL) |
263 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
264 |
pde.setSolverPackage(pde.PASO) |
265 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
266 |
# -------- test the solution --------------------------- |
267 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
268 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
269 |
class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase): |
270 |
def test_solve(self): |
271 |
domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE) |
272 |
x=Solution(domain).getX() |
273 |
# --- set exact solution ---- |
274 |
u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2 |
275 |
# --- set exact gradient ----------- |
276 |
g_ex=Data(0.,(3,),Solution(domain)) |
277 |
g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0] |
278 |
g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1] |
279 |
g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2] |
280 |
# -------- test gradient -------------------------------- |
281 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
282 |
# -------- set-up PDE ----------------------------------- |
283 |
pde=LinearPDE(domain,numEquations=1) |
284 |
mask=whereZero(x[0]) |
285 |
pde.setValue(r=u_ex,q=mask) |
286 |
pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.) |
287 |
# -------- get the solution --------------------------- |
288 |
pde.setTolerance(SOLVER_TOL) |
289 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
290 |
pde.setSolverPackage(pde.PASO) |
291 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
292 |
# -------- test the solution --------------------------- |
293 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
294 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
295 |
class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase): |
296 |
def test_solve(self): |
297 |
domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE) |
298 |
x=Solution(domain).getX() |
299 |
# --- set exact solution ---- |
300 |
u_ex=Vector(0,Solution(domain)) |
301 |
u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]+6.*x[0]*x[1]+7.*x[1]*x[2]+8.*x[2]*x[0]+9.*x[0]**2+10.*x[1]**2+11.*x[2]**2 |
302 |
u_ex[1]=2.+4.*x[0]+1.*x[1]-6.*x[2]+3.*x[0]*x[1]+2.*x[1]*x[2]-8.*x[2]*x[0]-2.*x[0]**2+7.*x[1]**2+5.*x[2]**2 |
303 |
u_ex[2]=-2.+7.*x[0]+9.*x[1]+2*x[2]-6.*x[0]*x[1]+8.*x[1]*x[2]+2.*x[2]*x[0]+2.*x[0]**2+8.*x[1]**2+1.*x[2]**2 |
304 |
# --- set exact gradient ----------- |
305 |
g_ex=Data(0.,(3,3),Solution(domain)) |
306 |
g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0] |
307 |
g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1] |
308 |
g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2] |
309 |
g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0] |
310 |
g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1] |
311 |
g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2] |
312 |
g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0] |
313 |
g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1] |
314 |
g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2] |
315 |
# -------- test gradient -------------------------------- |
316 |
self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex)) |
317 |
# -------- set-up PDE ----------------------------------- |
318 |
pde=LinearPDE(domain,numEquations=3) |
319 |
mask=whereZero(x[0]) |
320 |
pde.setValue(r=u_ex,q=mask*numarray.ones(3,)) |
321 |
Y=Vector(0.,Function(domain)) |
322 |
Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG |
323 |
Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG |
324 |
Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG |
325 |
A=Tensor4(0,Function(domain)) |
326 |
A[0,:,0,:]=kronecker(3) |
327 |
A[1,:,1,:]=kronecker(3) |
328 |
A[2,:,2,:]=kronecker(3) |
329 |
pde.setValue(A=A, |
330 |
D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG, |
331 |
Y=Y-numarray.array([60.,20.,22.]), |
332 |
y=matrixmult(g_ex,domain.getNormal())) |
333 |
# -------- get the solution --------------------------- |
334 |
pde.setTolerance(SOLVER_TOL) |
335 |
pde.setSolverMethod(pde.PCG,pde.JACOBI) |
336 |
pde.setSolverPackage(pde.PASO) |
337 |
u=pde.getSolution(verbose=SOLVER_VERBOSE) |
338 |
# -------- test the solution --------------------------- |
339 |
error=Lsup(u-u_ex)/Lsup(u_ex) |
340 |
self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error) |
341 |
|
342 |
if __name__ == '__main__': |
343 |
suite = unittest.TestSuite() |
344 |
suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi)) |
345 |
suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi)) |
346 |
suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi)) |
347 |
suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi)) |
348 |
suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi)) |
349 |
suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi)) |
350 |
suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi)) |
351 |
suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi)) |
352 |
s=unittest.TextTestRunner(verbosity=2).run(suite) |