/[escript]/trunk/finley/test/python/run_simplesolve.py
ViewVC logotype

Contents of /trunk/finley/test/python/run_simplesolve.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2561 - (show annotations)
Mon Jul 27 07:38:45 2009 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 49084 byte(s)
small bug in the calculation of the errors fixed.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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__="https://launchpad.net/escript-finley"
21
22 """
23 Test suite for the linearPDE and pdetools test on finley
24
25 @remark:
26
27 @var __author__: name of author
28 @var __licence__: licence agreement
29 @var __url__: url entry point on documentation
30 @var __version__: version
31 @var __date__: date of the version
32 """
33
34 __author__="Lutz Gross, l.gross@uq.edu.au"
35
36 import unittest, sys
37
38 from esys.escript import *
39 from esys.finley import Rectangle,Brick
40 from esys.escript.linearPDEs import LinearPDE, SolverOptions
41 import numpy
42 OPTIMIZE=True
43 SOLVER_VERBOSE=False
44 # setNumberOfThreads(2)
45
46 try:
47 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
48 except KeyError:
49 FINLEY_TEST_DATA='.'
50
51 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
52
53 # number of elements in the spatial directions
54 NE0=8
55 NE1=10
56 NE2=12
57
58 NE0=12
59 NE1=12
60 NE2=8
61
62 SOLVER_TOL=1.e-8
63 REL_TOL=1.e-6
64
65 FAC_DIAG=1.
66 FAC_OFFDIAG=-0.4
67
68
69 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
70 def test_solve(self):
71 # Tell about how many MPI CPUs and OpenMP threads
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.getSolverOptions().setTolerance(SOLVER_TOL)
91 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
92 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
93 pde.getSolverOptions().setPackage(SolverOptions.PASO)
94 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
95 u=pde.getSolution()
96 # -------- test the solution ---------------------------
97 error=Lsup(u-u_ex)
98 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
99 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
100 def test_solve(self):
101 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
102 x=Solution(domain).getX()
103 # --- set exact solution ----
104 u_ex=Scalar(0,Solution(domain))
105 u_ex=1.+2.*x[0]+3.*x[1]
106 # --- set exact gradient -----------
107 g_ex=Data(0.,(2,),Solution(domain))
108 g_ex[0]=2.
109 g_ex[1]=3.
110 # -------- test gradient --------------------------------
111 g=grad(u_ex)
112 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
113 # -------- set-up PDE -----------------------------------
114 pde=LinearPDE(domain,numEquations=1)
115 mask=whereZero(x[0])
116 pde.setValue(r=u_ex,q=mask)
117 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
118 # -------- get the solution ---------------------------
119 pde.getSolverOptions().setTolerance(SOLVER_TOL)
120 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
121 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
122 pde.getSolverOptions().setPackage(SolverOptions.PASO)
123 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
124 u=pde.getSolution()
125 # -------- test the solution ---------------------------
126 error=Lsup(u-u_ex)
127 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
128 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
129 def test_solve(self):
130 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
131 x=Solution(domain).getX()
132 # --- set exact solution ----
133 u_ex=Vector(0,Solution(domain))
134 u_ex[0]=1.+2.*x[0]+3.*x[1]
135 u_ex[1]=-1.+3.*x[0]+2.*x[1]
136 # --- set exact gradient -----------
137 g_ex=Data(0.,(2,2),Solution(domain))
138 g_ex[0,0]=2.
139 g_ex[0,1]=3.
140 g_ex[1,0]=3.
141 g_ex[1,1]=2.
142 # -------- test gradient --------------------------------
143 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
144 # -------- set-up PDE -----------------------------------
145 pde=LinearPDE(domain,numEquations=2)
146 mask=whereZero(x[0])
147 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
148 A=Tensor4(0,Function(domain))
149 A[0,:,0,:]=kronecker(2)
150 A[1,:,1,:]=kronecker(2)
151 Y=Vector(0.,Function(domain))
152 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
153 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
154 pde.setValue(A=A,
155 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
156 Y=Y,
157 y=matrixmult(g_ex,domain.getNormal()))
158 # -------- get the solution ---------------------------
159 pde.getSolverOptions().setTolerance(SOLVER_TOL)
160 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
161 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
162 pde.getSolverOptions().setPackage(SolverOptions.PASO)
163 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
164 u=pde.getSolution()
165 # -------- test the solution ---------------------------
166 error=Lsup(u-u_ex)
167 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
168 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
169 def test_solve(self):
170 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
171 x=Solution(domain).getX()
172 # --- set exact solution ----
173 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
174 # --- set exact gradient -----------
175 g_ex=Data(0.,(2,),Solution(domain))
176 g_ex[0]=2.+8.*x[0]+5.*x[1]
177 g_ex[1]=3.+5.*x[0]+12.*x[1]
178 # -------- test gradient --------------------------------
179 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
180 # -------- set-up PDE -----------------------------------
181 pde=LinearPDE(domain,numEquations=1)
182 mask=whereZero(x[0])
183 pde.setValue(r=u_ex,q=mask)
184 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
185 # -------- get the solution ---------------------------
186 pde.getSolverOptions().setTolerance(SOLVER_TOL)
187 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
188 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
189 pde.getSolverOptions().setPackage(SolverOptions.PASO)
190 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
191 u=pde.getSolution()
192 # -------- test the solution ---------------------------
193 error=Lsup(u-u_ex)
194 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
195 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
196 def test_solve(self):
197 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
198 x=Solution(domain).getX()
199 # --- set exact solution ----
200 u_ex=Vector(0,Solution(domain))
201 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
202 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
203 # --- set exact gradient -----------
204 g_ex=Data(0.,(2,2),Solution(domain))
205 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
206 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
207 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
208 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
209 # -------- test gradient --------------------------------
210 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
211 # -------- set-up PDE -----------------------------------
212 pde=LinearPDE(domain,numEquations=2)
213 mask=whereZero(x[0])
214 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
215 A=Tensor4(0,Function(domain))
216 A[0,:,0,:]=kronecker(2)
217 A[1,:,1,:]=kronecker(2)
218 Y=Vector(0.,Function(domain))
219 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
220 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
221 pde.setValue(A=A,
222 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
223 Y=Y-[20.,10.],
224 y=matrixmult(g_ex,domain.getNormal()))
225 # -------- get the solution ---------------------------
226 pde.getSolverOptions().setTolerance(SOLVER_TOL)
227 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
228 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
229 pde.getSolverOptions().setPackage(SolverOptions.PASO)
230 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
231 u=pde.getSolution()
232 # -------- test the solution ---------------------------
233 error=Lsup(u-u_ex)
234 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
235 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
236 def test_solve(self):
237 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
238 x=Solution(domain).getX()
239 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
240 # --- set exact gradient -----------
241 g_ex=Data(0.,(3,),Solution(domain))
242 g_ex[0]=2.
243 g_ex[1]=3.
244 g_ex[2]=4.
245 # -------- test gradient --------------------------------
246 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
247 # -------- set-up PDE -----------------------------------
248 pde=LinearPDE(domain,numEquations=1)
249 mask=whereZero(x[0])
250 pde.setValue(r=u_ex,q=mask)
251 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
252 # -------- get the solution ---------------------------
253 pde.getSolverOptions().setTolerance(SOLVER_TOL)
254 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
255 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
256 pde.getSolverOptions().setPackage(SolverOptions.PASO)
257 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
258 u=pde.getSolution()
259 # -------- test the solution ---------------------------
260 error=Lsup(u-u_ex)
261 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
262 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
263 def test_solve(self):
264 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
265 x=Solution(domain).getX()
266 # --- set exact solution ----
267 u_ex=Vector(0,Solution(domain))
268 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
269 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
270 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
271 # --- set exact gradient -----------
272 g_ex=Data(0.,(3,3),Solution(domain))
273 g_ex[0,0]=2.
274 g_ex[0,1]=3.
275 g_ex[0,2]=4.
276 g_ex[1,0]=4.
277 g_ex[1,1]=1.
278 g_ex[1,2]=-2.
279 g_ex[2,0]=8.
280 g_ex[2,1]=4.
281 g_ex[2,2]=5.
282 # -------- test gradient --------------------------------
283 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
284 # -------- set-up PDE -----------------------------------
285 pde=LinearPDE(domain,numEquations=3)
286 mask=whereZero(x[0])
287 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
288 A=Tensor4(0,Function(domain))
289 A[0,:,0,:]=kronecker(3)
290 A[1,:,1,:]=kronecker(3)
291 A[2,:,2,:]=kronecker(3)
292 Y=Vector(0.,Function(domain))
293 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
294 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
295 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
296 pde.setValue(A=A,
297 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
298 Y=Y,
299 y=matrixmult(g_ex,domain.getNormal()))
300 # -------- get the solution ---------------------------
301 pde.getSolverOptions().setTolerance(SOLVER_TOL)
302 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
303 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
304 pde.getSolverOptions().setPackage(SolverOptions.PASO)
305 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
306 u=pde.getSolution()
307 # -------- test the solution ---------------------------
308 error=Lsup(u-u_ex)
309 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
310 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
311 def test_solve(self):
312 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
313 x=Solution(domain).getX()
314 # --- set exact solution ----
315 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
316 # --- set exact gradient -----------
317 g_ex=Data(0.,(3,),Solution(domain))
318 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
319 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
320 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
321 # -------- test gradient --------------------------------
322 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
323 # -------- set-up PDE -----------------------------------
324 pde=LinearPDE(domain,numEquations=1)
325 mask=whereZero(x[0])
326 pde.setValue(r=u_ex,q=mask)
327 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
328 # -------- get the solution ---------------------------
329 pde.getSolverOptions().setTolerance(SOLVER_TOL)
330 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
331 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
332 pde.getSolverOptions().setPackage(SolverOptions.PASO)
333 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
334 u=pde.getSolution()
335 # -------- test the solution ---------------------------
336 error=Lsup(u-u_ex)
337 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
338 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
339 def test_solve(self):
340 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
341 x=Solution(domain).getX()
342 # --- set exact solution ----
343 u_ex=Vector(0,Solution(domain))
344 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
345 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
346 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
347 # --- set exact gradient -----------
348 g_ex=Data(0.,(3,3),Solution(domain))
349 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
350 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
351 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
352 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
353 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
354 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
355 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
356 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
357 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
358 # -------- test gradient --------------------------------
359 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
360 # -------- set-up PDE -----------------------------------
361 pde=LinearPDE(domain,numEquations=3)
362 mask=whereZero(x[0])
363 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
364 Y=Vector(0.,Function(domain))
365 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
366 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
367 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
368 A=Tensor4(0,Function(domain))
369 A[0,:,0,:]=kronecker(3)
370 A[1,:,1,:]=kronecker(3)
371 A[2,:,2,:]=kronecker(3)
372 pde.setValue(A=A,
373 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
374 Y=Y-numpy.array([60.,20.,22.]),
375 y=matrixmult(g_ex,domain.getNormal()))
376 # -------- get the solution ---------------------------
377 pde.getSolverOptions().setTolerance(SOLVER_TOL)
378 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
379 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
380 pde.getSolverOptions().setPackage(SolverOptions.PASO)
381 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
382 u=pde.getSolution()
383 # -------- test the solution ---------------------------
384 error=Lsup(u-u_ex)
385 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
386
387 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
388 def test_solve(self):
389 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
390 x=Solution(domain).getX()
391 # --- set exact solution ----
392 u_ex=Scalar(0,Solution(domain))
393 u_ex=1.+2.*x[0]+3.*x[1]
394 # --- set exact gradient -----------
395 g_ex=Data(0.,(2,),Solution(domain))
396 g_ex[0]=2.
397 g_ex[1]=3.
398 # -------- test gradient --------------------------------
399 g=grad(u_ex)
400 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
401 # -------- set-up PDE -----------------------------------
402 pde=LinearPDE(domain,numEquations=1)
403 mask=whereZero(x[0])
404 pde.setValue(r=u_ex,q=mask)
405 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
406 # -------- get the solution ---------------------------
407 pde.getSolverOptions().setTolerance(SOLVER_TOL)
408 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
409 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
410 pde.getSolverOptions().setPackage(SolverOptions.PASO)
411 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
412 u=pde.getSolution()
413 # -------- test the solution ---------------------------
414 error=Lsup(u-u_ex)
415 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
416
417 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
418 def test_solve(self):
419 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
420 x=Solution(domain).getX()
421 # --- set exact solution ----
422 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
423 # --- set exact gradient -----------
424 g_ex=Data(0.,(2,),Solution(domain))
425 g_ex[0]=2.+8.*x[0]+5.*x[1]
426 g_ex[1]=3.+5.*x[0]+12.*x[1]
427 # -------- test gradient --------------------------------
428 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
429 # -------- set-up PDE -----------------------------------
430 pde=LinearPDE(domain,numEquations=1)
431 mask=whereZero(x[0])
432 pde.setValue(r=u_ex,q=mask)
433 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
434 # -------- get the solution ---------------------------
435 pde.getSolverOptions().setTolerance(SOLVER_TOL)
436 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
437 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
438 pde.getSolverOptions().setPackage(SolverOptions.PASO)
439 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
440 u=pde.getSolution()
441 # -------- test the solution ---------------------------
442 error=Lsup(u-u_ex)
443 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
444
445 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
446 def test_solve(self):
447 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
448 x=Solution(domain).getX()
449 # --- set exact solution ----
450 u_ex=Vector(0,Solution(domain))
451 u_ex[0]=1.+2.*x[0]+3.*x[1]
452 u_ex[1]=-1.+3.*x[0]+2.*x[1]
453 # --- set exact gradient -----------
454 g_ex=Data(0.,(2,2),Solution(domain))
455 g_ex[0,0]=2.
456 g_ex[0,1]=3.
457 g_ex[1,0]=3.
458 g_ex[1,1]=2.
459 # -------- test gradient --------------------------------
460 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
461 # -------- set-up PDE -----------------------------------
462 pde=LinearPDE(domain,numEquations=2)
463 mask=whereZero(x[0])
464 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
465 A=Tensor4(0,Function(domain))
466 A[0,:,0,:]=kronecker(2)
467 A[1,:,1,:]=kronecker(2)
468 Y=Vector(0.,Function(domain))
469 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
470 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
471 pde.setValue(A=A,
472 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
473 Y=Y,
474 y=matrixmult(g_ex,domain.getNormal()))
475 # -------- get the solution ---------------------------
476 pde.getSolverOptions().setTolerance(SOLVER_TOL)
477 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
478 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
479 pde.getSolverOptions().setPackage(SolverOptions.PASO)
480 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
481 u=pde.getSolution()
482 # -------- test the solution ---------------------------
483 error=Lsup(u-u_ex)
484 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
485
486 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
487 def test_solve(self):
488 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
489 x=Solution(domain).getX()
490 # --- set exact solution ----
491 u_ex=Vector(0,Solution(domain))
492 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
493 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
494 # --- set exact gradient -----------
495 g_ex=Data(0.,(2,2),Solution(domain))
496 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
497 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
498 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
499 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
500 # -------- test gradient --------------------------------
501 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
502 # -------- set-up PDE -----------------------------------
503 pde=LinearPDE(domain,numEquations=2)
504 mask=whereZero(x[0])
505 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
506 A=Tensor4(0,Function(domain))
507 A[0,:,0,:]=kronecker(2)
508 A[1,:,1,:]=kronecker(2)
509 Y=Vector(0.,Function(domain))
510 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
511 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
512 pde.setValue(A=A,
513 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
514 Y=Y-[20.,10.],
515 y=matrixmult(g_ex,domain.getNormal()))
516 # -------- get the solution ---------------------------
517 pde.getSolverOptions().setTolerance(SOLVER_TOL)
518 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
519 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
520 pde.getSolverOptions().setPackage(SolverOptions.PASO)
521 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
522 u=pde.getSolution()
523 # -------- test the solution ---------------------------
524 error=Lsup(u-u_ex)
525 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
526
527 class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
528 def test_solve(self):
529 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
530 x=Solution(domain).getX()
531 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
532 # --- set exact gradient -----------
533 g_ex=Data(0.,(3,),Solution(domain))
534 g_ex[0]=2.
535 g_ex[1]=3.
536 g_ex[2]=4.
537 # -------- test gradient --------------------------------
538 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
539 # -------- set-up PDE -----------------------------------
540 pde=LinearPDE(domain,numEquations=1)
541 mask=whereZero(x[0])
542 pde.setValue(r=u_ex,q=mask)
543 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
544 # -------- get the solution ---------------------------
545 pde.getSolverOptions().setTolerance(SOLVER_TOL)
546 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
547 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
548 pde.getSolverOptions().setPackage(SolverOptions.PASO)
549 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
550 u=pde.getSolution()
551 # -------- test the solution ---------------------------
552 error=Lsup(u-u_ex)
553 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
554
555 class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
556 def test_solve(self):
557 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
558 x=Solution(domain).getX()
559 # --- set exact solution ----
560 u_ex=Vector(0,Solution(domain))
561 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
562 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
563 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
564 # --- set exact gradient -----------
565 g_ex=Data(0.,(3,3),Solution(domain))
566 g_ex[0,0]=2.
567 g_ex[0,1]=3.
568 g_ex[0,2]=4.
569 g_ex[1,0]=4.
570 g_ex[1,1]=1.
571 g_ex[1,2]=-2.
572 g_ex[2,0]=8.
573 g_ex[2,1]=4.
574 g_ex[2,2]=5.
575 # -------- test gradient --------------------------------
576 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
577 # -------- set-up PDE -----------------------------------
578 pde=LinearPDE(domain,numEquations=3)
579 mask=whereZero(x[0])
580 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
581 A=Tensor4(0,Function(domain))
582 A[0,:,0,:]=kronecker(3)
583 A[1,:,1,:]=kronecker(3)
584 A[2,:,2,:]=kronecker(3)
585 Y=Vector(0.,Function(domain))
586 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
587 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
588 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
589 pde.setValue(A=A,
590 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
591 Y=Y,
592 y=matrixmult(g_ex,domain.getNormal()))
593 # -------- get the solution ---------------------------
594 pde.getSolverOptions().setTolerance(SOLVER_TOL)
595 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
596 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
597 pde.getSolverOptions().setPackage(SolverOptions.PASO)
598 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
599 u=pde.getSolution()
600 # -------- test the solution ---------------------------
601 error=Lsup(u-u_ex)
602 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
603
604 class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
605 def test_solve(self):
606 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
607 x=Solution(domain).getX()
608 # --- set exact solution ----
609 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
610 # --- set exact gradient -----------
611 g_ex=Data(0.,(3,),Solution(domain))
612 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
613 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
614 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
615 # -------- test gradient --------------------------------
616 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
617 # -------- set-up PDE -----------------------------------
618 pde=LinearPDE(domain,numEquations=1)
619 mask=whereZero(x[0])
620 pde.setValue(r=u_ex,q=mask)
621 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
622 # -------- get the solution ---------------------------
623 pde.getSolverOptions().setTolerance(SOLVER_TOL)
624 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
625 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
626 pde.getSolverOptions().setPackage(SolverOptions.PASO)
627 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
628 u=pde.getSolution()
629 # -------- test the solution ---------------------------
630 error=Lsup(u-u_ex)
631 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
632
633 class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
634 def test_solve(self):
635 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
636 x=Solution(domain).getX()
637 # --- set exact solution ----
638 u_ex=Vector(0,Solution(domain))
639 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
640 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
641 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
642 # --- set exact gradient -----------
643 g_ex=Data(0.,(3,3),Solution(domain))
644 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
645 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
646 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
647 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
648 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
649 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
650 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
651 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
652 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
653 # -------- test gradient --------------------------------
654 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
655 # -------- set-up PDE -----------------------------------
656 pde=LinearPDE(domain,numEquations=3)
657 mask=whereZero(x[0])
658 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
659 Y=Vector(0.,Function(domain))
660 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
661 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
662 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
663 A=Tensor4(0,Function(domain))
664 A[0,:,0,:]=kronecker(3)
665 A[1,:,1,:]=kronecker(3)
666 A[2,:,2,:]=kronecker(3)
667 pde.setValue(A=A,
668 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
669 Y=Y-numpy.array([60.,20.,22.]),
670 y=matrixmult(g_ex,domain.getNormal()))
671 # -------- get the solution ---------------------------
672 pde.getSolverOptions().setTolerance(SOLVER_TOL)
673 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
674 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
675 pde.getSolverOptions().setPackage(SolverOptions.PASO)
676 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
677 u=pde.getSolution()
678 # -------- test the solution ---------------------------
679 error=Lsup(u-u_ex)
680 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
681
682
683 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
684 def test_solve(self):
685 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
686 x=Solution(domain).getX()
687 # --- set exact solution ----
688 u_ex=Scalar(0,Solution(domain))
689 u_ex=1.+2.*x[0]+3.*x[1]
690 # --- set exact gradient -----------
691 g_ex=Data(0.,(2,),Solution(domain))
692 g_ex[0]=2.
693 g_ex[1]=3.
694 # -------- test gradient --------------------------------
695 g=grad(u_ex)
696 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
697 # -------- set-up PDE -----------------------------------
698 pde=LinearPDE(domain,numEquations=1)
699 mask=whereZero(x[0])
700 pde.setValue(r=u_ex,q=mask)
701 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
702 # -------- get the solution ---------------------------
703 pde.getSolverOptions().setTolerance(SOLVER_TOL)
704 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
705 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
706 pde.getSolverOptions().setPackage(SolverOptions.PASO)
707
708 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
709 u=pde.getSolution()
710 # -------- test the solution ---------------------------
711 error=Lsup(u-u_ex)
712 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
713
714 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
715 def test_solve(self):
716 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
717 x=Solution(domain).getX()
718 # --- set exact solution ----
719 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
720 # --- set exact gradient -----------
721 g_ex=Data(0.,(2,),Solution(domain))
722 g_ex[0]=2.+8.*x[0]+5.*x[1]
723 g_ex[1]=3.+5.*x[0]+12.*x[1]
724 # -------- test gradient --------------------------------
725 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
726 # -------- set-up PDE -----------------------------------
727 pde=LinearPDE(domain,numEquations=1)
728 mask=whereZero(x[0])
729 pde.setValue(r=u_ex,q=mask)
730 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
731 # -------- get the solution ---------------------------
732 pde.getSolverOptions().setTolerance(SOLVER_TOL)
733 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
734 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
735 pde.getSolverOptions().setPackage(SolverOptions.PASO)
736 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
737 u=pde.getSolution()
738 # -------- test the solution ---------------------------
739 error=Lsup(u-u_ex)
740 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
741
742 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
743 def test_solve(self):
744 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
745 x=Solution(domain).getX()
746 # --- set exact solution ----
747 u_ex=Vector(0,Solution(domain))
748 u_ex[0]=1.+2.*x[0]+3.*x[1]
749 u_ex[1]=-1.+3.*x[0]+2.*x[1]
750 # --- set exact gradient -----------
751 g_ex=Data(0.,(2,2),Solution(domain))
752 g_ex[0,0]=2.
753 g_ex[0,1]=3.
754 g_ex[1,0]=3.
755 g_ex[1,1]=2.
756 # -------- test gradient --------------------------------
757 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
758 # -------- set-up PDE -----------------------------------
759 pde=LinearPDE(domain,numEquations=2)
760 mask=whereZero(x[0])
761 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
762 A=Tensor4(0,Function(domain))
763 A[0,:,0,:]=kronecker(2)
764 A[1,:,1,:]=kronecker(2)
765 Y=Vector(0.,Function(domain))
766 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
767 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
768 pde.setValue(A=A,
769 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
770 Y=Y,
771 y=matrixmult(g_ex,domain.getNormal()))
772 # -------- get the solution ---------------------------
773 pde.getSolverOptions().setTolerance(SOLVER_TOL)
774 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
775 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
776 pde.getSolverOptions().setPackage(SolverOptions.PASO)
777 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
778 u=pde.getSolution()
779 # -------- test the solution ---------------------------
780 error=Lsup(u-u_ex)
781 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
782
783 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
784 def test_solve(self):
785 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
786 x=Solution(domain).getX()
787 # --- set exact solution ----
788 u_ex=Vector(0,Solution(domain))
789 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
790 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
791 # --- set exact gradient -----------
792 g_ex=Data(0.,(2,2),Solution(domain))
793 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
794 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
795 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
796 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
797 # -------- test gradient --------------------------------
798 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
799 # -------- set-up PDE -----------------------------------
800 pde=LinearPDE(domain,numEquations=2)
801 mask=whereZero(x[0])
802 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
803 A=Tensor4(0,Function(domain))
804 A[0,:,0,:]=kronecker(2)
805 A[1,:,1,:]=kronecker(2)
806 Y=Vector(0.,Function(domain))
807 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
808 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
809 pde.setValue(A=A,
810 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
811 Y=Y-[20.,10.],
812 y=matrixmult(g_ex,domain.getNormal()))
813 # -------- get the solution ---------------------------
814 pde.getSolverOptions().setTolerance(SOLVER_TOL)
815 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
816 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
817 pde.getSolverOptions().setPackage(SolverOptions.PASO)
818 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
819 u=pde.getSolution()
820 # -------- test the solution ---------------------------
821 error=Lsup(u-u_ex)
822 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
823
824 class SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
825 def test_solve(self):
826 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
827 x=Solution(domain).getX()
828 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
829 # --- set exact gradient -----------
830 g_ex=Data(0.,(3,),Solution(domain))
831 g_ex[0]=2.
832 g_ex[1]=3.
833 g_ex[2]=4.
834 # -------- test gradient --------------------------------
835 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
836 # -------- set-up PDE -----------------------------------
837 pde=LinearPDE(domain,numEquations=1)
838 mask=whereZero(x[0])
839 pde.setValue(r=u_ex,q=mask)
840 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
841 # -------- get the solution ---------------------------
842 pde.getSolverOptions().setTolerance(SOLVER_TOL)
843 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
844 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
845 pde.getSolverOptions().setPackage(SolverOptions.PASO)
846 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
847 u=pde.getSolution()
848 # -------- test the solution ---------------------------
849 error=Lsup(u-u_ex)
850 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
851
852 class SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
853 def test_solve(self):
854 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
855 x=Solution(domain).getX()
856 # --- set exact solution ----
857 u_ex=Vector(0,Solution(domain))
858 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
859 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
860 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
861 # --- set exact gradient -----------
862 g_ex=Data(0.,(3,3),Solution(domain))
863 g_ex[0,0]=2.
864 g_ex[0,1]=3.
865 g_ex[0,2]=4.
866 g_ex[1,0]=4.
867 g_ex[1,1]=1.
868 g_ex[1,2]=-2.
869 g_ex[2,0]=8.
870 g_ex[2,1]=4.
871 g_ex[2,2]=5.
872 # -------- test gradient --------------------------------
873 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
874 # -------- set-up PDE -----------------------------------
875 pde=LinearPDE(domain,numEquations=3)
876 mask=whereZero(x[0])
877 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
878 A=Tensor4(0,Function(domain))
879 A[0,:,0,:]=kronecker(3)
880 A[1,:,1,:]=kronecker(3)
881 A[2,:,2,:]=kronecker(3)
882 Y=Vector(0.,Function(domain))
883 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
884 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
885 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
886 pde.setValue(A=A,
887 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
888 Y=Y,
889 y=matrixmult(g_ex,domain.getNormal()))
890 # -------- get the solution ---------------------------
891 pde.getSolverOptions().setTolerance(SOLVER_TOL)
892 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
893 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
894 pde.getSolverOptions().setPackage(SolverOptions.PASO)
895 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
896 u=pde.getSolution()
897 # -------- test the solution ---------------------------
898 error=Lsup(u-u_ex)
899 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
900
901 class SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
902 def test_solve(self):
903 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
904 x=Solution(domain).getX()
905 # --- set exact solution ----
906 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
907 # --- set exact gradient -----------
908 g_ex=Data(0.,(3,),Solution(domain))
909 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
910 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
911 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
912 # -------- test gradient --------------------------------
913 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
914 # -------- set-up PDE -----------------------------------
915 pde=LinearPDE(domain,numEquations=1)
916 mask=whereZero(x[0])
917 pde.setValue(r=u_ex,q=mask)
918 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
919 # -------- get the solution ---------------------------
920 pde.getSolverOptions().setTolerance(SOLVER_TOL)
921 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
922 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
923 pde.getSolverOptions().setPackage(SolverOptions.PASO)
924 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
925 u=pde.getSolution()
926 # -------- test the solution ---------------------------
927 error=Lsup(u-u_ex)
928 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
929
930 class SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
931 def test_solve(self):
932 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
933 x=Solution(domain).getX()
934 # --- set exact solution ----
935 u_ex=Vector(0,Solution(domain))
936 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
937 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
938 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
939 # --- set exact gradient -----------
940 g_ex=Data(0.,(3,3),Solution(domain))
941 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
942 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
943 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
944 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
945 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
946 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
947 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
948 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
949 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
950 # -------- test gradient --------------------------------
951 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
952 # -------- set-up PDE -----------------------------------
953 pde=LinearPDE(domain,numEquations=3)
954 mask=whereZero(x[0])
955 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
956 Y=Vector(0.,Function(domain))
957 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
958 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
959 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
960 A=Tensor4(0,Function(domain))
961 A[0,:,0,:]=kronecker(3)
962 A[1,:,1,:]=kronecker(3)
963 A[2,:,2,:]=kronecker(3)
964 pde.setValue(A=A,
965 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
966 Y=Y-numpy.array([60.,20.,22.]),
967 y=matrixmult(g_ex,domain.getNormal()))
968 # -------- get the solution ---------------------------
969 pde.getSolverOptions().setTolerance(SOLVER_TOL)
970 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
971 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
972 pde.getSolverOptions().setPackage(SolverOptions.PASO)
973 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
974 u=pde.getSolution()
975 # -------- test the solution ---------------------------
976 error=Lsup(u-u_ex)
977 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
978
979 if __name__ == '__main__':
980 suite = unittest.TestSuite()
981 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
982 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
983 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
984 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
985 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
986 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
987 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
988 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
989 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
990 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_Jacobi))
991 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_Jacobi))
992 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_Jacobi))
993 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_Jacobi))
994 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_Jacobi))
995 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_Jacobi))
996 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_Jacobi))
997 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_Jacobi))
998 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi))
999 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_Jacobi))
1000 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi))
1001 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_Jacobi))
1002 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_Jacobi))
1003 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi))
1004 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi))
1005 suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi))
1006
1007 s=unittest.TextTestRunner(verbosity=2).run(suite)
1008 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26