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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3097 - (hide annotations)
Fri Aug 20 04:59:12 2010 UTC (9 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 127414 byte(s)
some modifications to the GaussSeidel
1 ksteube 1809
2     ########################################################
3 ksteube 1310 #
4 jfenwick 2881 # Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1310 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1310 #
12 ksteube 1809 ########################################################
13 ksteube 1310
14 jfenwick 2881 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15 ksteube 1809 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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 ksteube 1310 """
23     Test suite for the linearPDE and pdetools test on finley
24    
25 jfenwick 2625 :remark:
26 ksteube 1310
27 jfenwick 2625 :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 ksteube 1310 """
33    
34     __author__="Lutz Gross, l.gross@uq.edu.au"
35    
36 gross 1553 import unittest, sys
37 ksteube 1310
38     from esys.escript import *
39     from esys.finley import Rectangle,Brick
40 gross 2474 from esys.escript.linearPDEs import LinearPDE, SolverOptions
41 jfenwick 2455 import numpy
42 gross 1804 OPTIMIZE=True
43 gross 1553 SOLVER_VERBOSE=False
44 gross 1556 # setNumberOfThreads(2)
45 ksteube 1310
46     try:
47     FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
48     except KeyError:
49     FINLEY_TEST_DATA='.'
50    
51 gross 2416 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
52 ksteube 1310
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 gross 1556
69 artak 1425 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
70     def test_solve(self):
71 ksteube 1567 # Tell about how many MPI CPUs and OpenMP threads
72 artak 1425 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 gross 2474 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 artak 1425 # -------- test the solution ---------------------------
97 gross 2561 error=Lsup(u-u_ex)
98 artak 1425 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
99 ksteube 1310 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 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
126 gross 2561 error=Lsup(u-u_ex)
127 ksteube 1310 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
148 ksteube 1310 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
156 ksteube 1310 Y=Y,
157     y=matrixmult(g_ex,domain.getNormal()))
158     # -------- get the solution ---------------------------
159 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
166 gross 2561 error=Lsup(u-u_ex)
167 ksteube 1310 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 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
193 gross 2561 error=Lsup(u-u_ex)
194 ksteube 1310 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
215 ksteube 1310 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
223 ksteube 1310 Y=Y-[20.,10.],
224     y=matrixmult(g_ex,domain.getNormal()))
225     # -------- get the solution ---------------------------
226 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
233 gross 2561 error=Lsup(u-u_ex)
234 ksteube 1310 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 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
260 gross 2561 error=Lsup(u-u_ex)
261 ksteube 1310 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
288 ksteube 1310 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
298 ksteube 1310 Y=Y,
299     y=matrixmult(g_ex,domain.getNormal()))
300     # -------- get the solution ---------------------------
301 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
308 gross 2561 error=Lsup(u-u_ex)
309 ksteube 1310 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 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
336 gross 2561 error=Lsup(u-u_ex)
337 ksteube 1310 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
364 ksteube 1310 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
374     Y=Y-numpy.array([60.,20.,22.]),
375 ksteube 1310 y=matrixmult(g_ex,domain.getNormal()))
376     # -------- get the solution ---------------------------
377 gross 2474 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 ksteube 1310 # -------- test the solution ---------------------------
384 gross 2561 error=Lsup(u-u_ex)
385 ksteube 1310 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
386 artak 1703
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 gross 2474 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 artak 1703 # -------- test the solution ---------------------------
414 gross 2561 error=Lsup(u-u_ex)
415 artak 1703 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
416 artak 1787
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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
442 gross 2561 error=Lsup(u-u_ex)
443 artak 1787 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
444 ksteube 1310
445 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
465 artak 1787 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
473 artak 1787 Y=Y,
474     y=matrixmult(g_ex,domain.getNormal()))
475     # -------- get the solution ---------------------------
476 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
483 gross 2561 error=Lsup(u-u_ex)
484 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
506 artak 1787 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
514 artak 1787 Y=Y-[20.,10.],
515     y=matrixmult(g_ex,domain.getNormal()))
516     # -------- get the solution ---------------------------
517 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
524 gross 2561 error=Lsup(u-u_ex)
525 artak 1787 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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
552 gross 2561 error=Lsup(u-u_ex)
553 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
581 artak 1787 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
591 artak 1787 Y=Y,
592     y=matrixmult(g_ex,domain.getNormal()))
593     # -------- get the solution ---------------------------
594 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
601 gross 2561 error=Lsup(u-u_ex)
602 artak 1787 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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
630 gross 2561 error=Lsup(u-u_ex)
631 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
659 artak 1787 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
669     Y=Y-numpy.array([60.,20.,22.]),
670 artak 1787 y=matrixmult(g_ex,domain.getNormal()))
671     # -------- get the solution ---------------------------
672 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
679 gross 2561 error=Lsup(u-u_ex)
680 artak 1787 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 gross 2474 pde.getSolverOptions().setTolerance(SOLVER_TOL)
704     pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
705     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
706     pde.getSolverOptions().setPackage(SolverOptions.PASO)
707 artak 1787
708 gross 2474 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
709     u=pde.getSolution()
710 artak 1787 # -------- test the solution ---------------------------
711 gross 2561 error=Lsup(u-u_ex)
712 artak 1787 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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
739 gross 2561 error=Lsup(u-u_ex)
740 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
762 artak 1787 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
770 artak 1787 Y=Y,
771     y=matrixmult(g_ex,domain.getNormal()))
772     # -------- get the solution ---------------------------
773 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
780 gross 2561 error=Lsup(u-u_ex)
781 artak 1787 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
782    
783 artak 2330 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
784 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
803 artak 1787 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 jfenwick 2455 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
811 artak 1787 Y=Y-[20.,10.],
812     y=matrixmult(g_ex,domain.getNormal()))
813     # -------- get the solution ---------------------------
814 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
821 gross 2561 error=Lsup(u-u_ex)
822 artak 1787 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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
849 gross 2561 error=Lsup(u-u_ex)
850 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
878 artak 1787 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
888 artak 1787 Y=Y,
889     y=matrixmult(g_ex,domain.getNormal()))
890     # -------- get the solution ---------------------------
891 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
898 gross 2561 error=Lsup(u-u_ex)
899 artak 1787 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
900    
901 artak 2330 class SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
902 artak 1787 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 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
927 gross 2561 error=Lsup(u-u_ex)
928 artak 1787 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
929    
930 artak 2330 class SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
931 artak 1787 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 jfenwick 2455 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
956 artak 1787 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 jfenwick 2455 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
966     Y=Y-numpy.array([60.,20.,22.]),
967 artak 1787 y=matrixmult(g_ex,domain.getNormal()))
968     # -------- get the solution ---------------------------
969 gross 2474 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 artak 1787 # -------- test the solution ---------------------------
976 gross 2561 error=Lsup(u-u_ex)
977 artak 1787 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
978    
979 gross 3097 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
980     def test_solve(self):
981     domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
982     x=Solution(domain).getX()
983     # --- set exact solution ----
984     u_ex=Scalar(0,Solution(domain))
985     u_ex=1.+2.*x[0]+3.*x[1]
986     # --- set exact gradient -----------
987     g_ex=Data(0.,(2,),Solution(domain))
988     g_ex[0]=2.
989     g_ex[1]=3.
990     # -------- test gradient --------------------------------
991     g=grad(u_ex)
992     self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
993     # -------- set-up PDE -----------------------------------
994     pde=LinearPDE(domain,numEquations=1)
995     mask=whereZero(x[0])
996     pde.setValue(r=u_ex,q=mask)
997     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
998     # -------- get the solution ---------------------------
999     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1000     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1001     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1002     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1003     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1004     u=pde.getSolution()
1005     # -------- test the solution ---------------------------
1006     error=Lsup(u-u_ex)
1007     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1008    
1009     class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1010     def test_solve(self):
1011     domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1012     x=Solution(domain).getX()
1013     # --- set exact solution ----
1014     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1015     # --- set exact gradient -----------
1016     g_ex=Data(0.,(2,),Solution(domain))
1017     g_ex[0]=2.+8.*x[0]+5.*x[1]
1018     g_ex[1]=3.+5.*x[0]+12.*x[1]
1019     # -------- test gradient --------------------------------
1020     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1021     # -------- set-up PDE -----------------------------------
1022     pde=LinearPDE(domain,numEquations=1)
1023     mask=whereZero(x[0])
1024     pde.setValue(r=u_ex,q=mask)
1025     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1026     # -------- get the solution ---------------------------
1027     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1028     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1029     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1030     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1031     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1032     u=pde.getSolution()
1033     # -------- test the solution ---------------------------
1034     error=Lsup(u-u_ex)
1035     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1036    
1037     class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1038     def test_solve(self):
1039     domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1040     x=Solution(domain).getX()
1041     # --- set exact solution ----
1042     u_ex=Vector(0,Solution(domain))
1043     u_ex[0]=1.+2.*x[0]+3.*x[1]
1044     u_ex[1]=-1.+3.*x[0]+2.*x[1]
1045     # --- set exact gradient -----------
1046     g_ex=Data(0.,(2,2),Solution(domain))
1047     g_ex[0,0]=2.
1048     g_ex[0,1]=3.
1049     g_ex[1,0]=3.
1050     g_ex[1,1]=2.
1051     # -------- test gradient --------------------------------
1052     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1053     # -------- set-up PDE -----------------------------------
1054     pde=LinearPDE(domain,numEquations=2)
1055     mask=whereZero(x[0])
1056     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1057     A=Tensor4(0,Function(domain))
1058     A[0,:,0,:]=kronecker(2)
1059     A[1,:,1,:]=kronecker(2)
1060     Y=Vector(0.,Function(domain))
1061     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1062     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1063     pde.setValue(A=A,
1064     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1065     Y=Y,
1066     y=matrixmult(g_ex,domain.getNormal()))
1067     # -------- get the solution ---------------------------
1068     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1069     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1070     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1071     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1072     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1073     u=pde.getSolution()
1074     # -------- test the solution ---------------------------
1075     error=Lsup(u-u_ex)
1076     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1077    
1078     class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1079     def test_solve(self):
1080     domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1081     x=Solution(domain).getX()
1082     # --- set exact solution ----
1083     u_ex=Vector(0,Solution(domain))
1084     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1085     u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1086     # --- set exact gradient -----------
1087     g_ex=Data(0.,(2,2),Solution(domain))
1088     g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1089     g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1090     g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1091     g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1092     # -------- test gradient --------------------------------
1093     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1094     # -------- set-up PDE -----------------------------------
1095     pde=LinearPDE(domain,numEquations=2)
1096     mask=whereZero(x[0])
1097     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1098     A=Tensor4(0,Function(domain))
1099     A[0,:,0,:]=kronecker(2)
1100     A[1,:,1,:]=kronecker(2)
1101     Y=Vector(0.,Function(domain))
1102     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1103     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1104     pde.setValue(A=A,
1105     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1106     Y=Y-[20.,10.],
1107     y=matrixmult(g_ex,domain.getNormal()))
1108     # -------- get the solution ---------------------------
1109     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1110     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1111     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1112     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1113     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1114     u=pde.getSolution()
1115     # -------- test the solution ---------------------------
1116     error=Lsup(u-u_ex)
1117     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1118    
1119     class SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1120     def test_solve(self):
1121     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1122     x=Solution(domain).getX()
1123     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1124     # --- set exact gradient -----------
1125     g_ex=Data(0.,(3,),Solution(domain))
1126     g_ex[0]=2.
1127     g_ex[1]=3.
1128     g_ex[2]=4.
1129     # -------- test gradient --------------------------------
1130     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1131     # -------- set-up PDE -----------------------------------
1132     pde=LinearPDE(domain,numEquations=1)
1133     mask=whereZero(x[0])
1134     pde.setValue(r=u_ex,q=mask)
1135     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1136     # -------- get the solution ---------------------------
1137     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1138     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1139     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1140     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1141     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1142     u=pde.getSolution()
1143     # -------- test the solution ---------------------------
1144     error=Lsup(u-u_ex)
1145     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1146    
1147     class SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1148     def test_solve(self):
1149     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1150     x=Solution(domain).getX()
1151     # --- set exact solution ----
1152     u_ex=Vector(0,Solution(domain))
1153     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1154     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1155     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1156     # --- set exact gradient -----------
1157     g_ex=Data(0.,(3,3),Solution(domain))
1158     g_ex[0,0]=2.
1159     g_ex[0,1]=3.
1160     g_ex[0,2]=4.
1161     g_ex[1,0]=4.
1162     g_ex[1,1]=1.
1163     g_ex[1,2]=-2.
1164     g_ex[2,0]=8.
1165     g_ex[2,1]=4.
1166     g_ex[2,2]=5.
1167     # -------- test gradient --------------------------------
1168     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1169     # -------- set-up PDE -----------------------------------
1170     pde=LinearPDE(domain,numEquations=3)
1171     mask=whereZero(x[0])
1172     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1173     A=Tensor4(0,Function(domain))
1174     A[0,:,0,:]=kronecker(3)
1175     A[1,:,1,:]=kronecker(3)
1176     A[2,:,2,:]=kronecker(3)
1177     Y=Vector(0.,Function(domain))
1178     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1179     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1180     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1181     pde.setValue(A=A,
1182     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1183     Y=Y,
1184     y=matrixmult(g_ex,domain.getNormal()))
1185     # -------- get the solution ---------------------------
1186     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1187     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1188     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1189     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1190     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1191     u=pde.getSolution()
1192     # -------- test the solution ---------------------------
1193     error=Lsup(u-u_ex)
1194     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1195    
1196     class SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1197     def test_solve(self):
1198     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1199     x=Solution(domain).getX()
1200     # --- set exact solution ----
1201     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
1202     # --- set exact gradient -----------
1203     g_ex=Data(0.,(3,),Solution(domain))
1204     g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1205     g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1206     g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1207     # -------- test gradient --------------------------------
1208     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1209     # -------- set-up PDE -----------------------------------
1210     pde=LinearPDE(domain,numEquations=1)
1211     mask=whereZero(x[0])
1212     pde.setValue(r=u_ex,q=mask)
1213     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1214     # -------- get the solution ---------------------------
1215     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1216     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1217     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1218     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1219     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1220     u=pde.getSolution()
1221     # -------- test the solution ---------------------------
1222     error=Lsup(u-u_ex)
1223     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1224    
1225     class SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1226     def test_solve(self):
1227     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1228     x=Solution(domain).getX()
1229     # --- set exact solution ----
1230     u_ex=Vector(0,Solution(domain))
1231     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
1232     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
1233     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
1234     # --- set exact gradient -----------
1235     g_ex=Data(0.,(3,3),Solution(domain))
1236     g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1237     g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1238     g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1239     g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1240     g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1241     g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1242     g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1243     g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1244     g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1245     # -------- test gradient --------------------------------
1246     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1247     # -------- set-up PDE -----------------------------------
1248     pde=LinearPDE(domain,numEquations=3)
1249     mask=whereZero(x[0])
1250     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1251     Y=Vector(0.,Function(domain))
1252     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1253     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1254     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1255     A=Tensor4(0,Function(domain))
1256     A[0,:,0,:]=kronecker(3)
1257     A[1,:,1,:]=kronecker(3)
1258     A[2,:,2,:]=kronecker(3)
1259     pde.setValue(A=A,
1260     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1261     Y=Y-numpy.array([60.,20.,22.]),
1262     y=matrixmult(g_ex,domain.getNormal()))
1263     # -------- get the solution ---------------------------
1264     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1265     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1266     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1267     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1268     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1269     u=pde.getSolution()
1270     # -------- test the solution ---------------------------
1271     error=Lsup(u-u_ex)
1272     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1273    
1274    
1275     class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
1276     def test_solve(self):
1277     # Tell about how many MPI CPUs and OpenMP threads
1278     domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1279     x=Solution(domain).getX()
1280     # --- set exact solution ----
1281     u_ex=Scalar(0,Solution(domain))
1282     u_ex=1.+2.*x[0]+3.*x[1]
1283     # --- set exact gradient -----------
1284     g_ex=Data(0.,(2,),Solution(domain))
1285     g_ex[0]=2.
1286     g_ex[1]=3.
1287     # -------- test gradient --------------------------------
1288     g=grad(u_ex)
1289     self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1290     # -------- set-up PDE -----------------------------------
1291     pde=LinearPDE(domain,numEquations=1)
1292     mask=whereZero(x[0])
1293     pde.setValue(r=u_ex,q=mask)
1294     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1295     # -------- get the solution ---------------------------
1296     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1297     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1298     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1299     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1300     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1301     u=pde.getSolution()
1302     # -------- test the solution ---------------------------
1303     error=Lsup(u-u_ex)
1304     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1305     class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1306     def test_solve(self):
1307     domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1308     x=Solution(domain).getX()
1309     # --- set exact solution ----
1310     u_ex=Scalar(0,Solution(domain))
1311     u_ex=1.+2.*x[0]+3.*x[1]
1312     # --- set exact gradient -----------
1313     g_ex=Data(0.,(2,),Solution(domain))
1314     g_ex[0]=2.
1315     g_ex[1]=3.
1316     # -------- test gradient --------------------------------
1317     g=grad(u_ex)
1318     self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1319     # -------- set-up PDE -----------------------------------
1320     pde=LinearPDE(domain,numEquations=1)
1321     mask=whereZero(x[0])
1322     pde.setValue(r=u_ex,q=mask)
1323     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1324     # -------- get the solution ---------------------------
1325     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1326     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1327     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1328     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1329     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1330     u=pde.getSolution()
1331     # -------- test the solution ---------------------------
1332     error=Lsup(u-u_ex)
1333     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1334     class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1335     def test_solve(self):
1336     domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1337     x=Solution(domain).getX()
1338     # --- set exact solution ----
1339     u_ex=Vector(0,Solution(domain))
1340     u_ex[0]=1.+2.*x[0]+3.*x[1]
1341     u_ex[1]=-1.+3.*x[0]+2.*x[1]
1342     # --- set exact gradient -----------
1343     g_ex=Data(0.,(2,2),Solution(domain))
1344     g_ex[0,0]=2.
1345     g_ex[0,1]=3.
1346     g_ex[1,0]=3.
1347     g_ex[1,1]=2.
1348     # -------- test gradient --------------------------------
1349     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1350     # -------- set-up PDE -----------------------------------
1351     pde=LinearPDE(domain,numEquations=2)
1352     mask=whereZero(x[0])
1353     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1354     A=Tensor4(0,Function(domain))
1355     A[0,:,0,:]=kronecker(2)
1356     A[1,:,1,:]=kronecker(2)
1357     Y=Vector(0.,Function(domain))
1358     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1359     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1360     pde.setValue(A=A,
1361     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1362     Y=Y,
1363     y=matrixmult(g_ex,domain.getNormal()))
1364     # -------- get the solution ---------------------------
1365     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1366     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1367     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1368     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1369     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1370     u=pde.getSolution()
1371     # -------- test the solution ---------------------------
1372     error=Lsup(u-u_ex)
1373     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1374     class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1375     def test_solve(self):
1376     domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1377     x=Solution(domain).getX()
1378     # --- set exact solution ----
1379     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1380     # --- set exact gradient -----------
1381     g_ex=Data(0.,(2,),Solution(domain))
1382     g_ex[0]=2.+8.*x[0]+5.*x[1]
1383     g_ex[1]=3.+5.*x[0]+12.*x[1]
1384     # -------- test gradient --------------------------------
1385     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1386     # -------- set-up PDE -----------------------------------
1387     pde=LinearPDE(domain,numEquations=1)
1388     mask=whereZero(x[0])
1389     pde.setValue(r=u_ex,q=mask)
1390     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1391     # -------- get the solution ---------------------------
1392     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1393     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1394     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1395     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1396     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1397     u=pde.getSolution()
1398     # -------- test the solution ---------------------------
1399     error=Lsup(u-u_ex)
1400     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1401     class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1402     def test_solve(self):
1403     domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1404     x=Solution(domain).getX()
1405     # --- set exact solution ----
1406     u_ex=Vector(0,Solution(domain))
1407     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1408     u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1409     # --- set exact gradient -----------
1410     g_ex=Data(0.,(2,2),Solution(domain))
1411     g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1412     g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1413     g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1414     g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1415     # -------- test gradient --------------------------------
1416     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1417     # -------- set-up PDE -----------------------------------
1418     pde=LinearPDE(domain,numEquations=2)
1419     mask=whereZero(x[0])
1420     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1421     A=Tensor4(0,Function(domain))
1422     A[0,:,0,:]=kronecker(2)
1423     A[1,:,1,:]=kronecker(2)
1424     Y=Vector(0.,Function(domain))
1425     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1426     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1427     pde.setValue(A=A,
1428     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1429     Y=Y-[20.,10.],
1430     y=matrixmult(g_ex,domain.getNormal()))
1431     # -------- get the solution ---------------------------
1432     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1433     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1434     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1435     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1436     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1437     u=pde.getSolution()
1438     # -------- test the solution ---------------------------
1439     error=Lsup(u-u_ex)
1440     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1441     class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1442     def test_solve(self):
1443     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1444     x=Solution(domain).getX()
1445     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1446     # --- set exact gradient -----------
1447     g_ex=Data(0.,(3,),Solution(domain))
1448     g_ex[0]=2.
1449     g_ex[1]=3.
1450     g_ex[2]=4.
1451     # -------- test gradient --------------------------------
1452     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1453     # -------- set-up PDE -----------------------------------
1454     pde=LinearPDE(domain,numEquations=1)
1455     mask=whereZero(x[0])
1456     pde.setValue(r=u_ex,q=mask)
1457     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1458     # -------- get the solution ---------------------------
1459     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1460     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1461     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1462     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1463     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1464     u=pde.getSolution()
1465     # -------- test the solution ---------------------------
1466     error=Lsup(u-u_ex)
1467     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1468     class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1469     def test_solve(self):
1470     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1471     x=Solution(domain).getX()
1472     # --- set exact solution ----
1473     u_ex=Vector(0,Solution(domain))
1474     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1475     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1476     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1477     # --- set exact gradient -----------
1478     g_ex=Data(0.,(3,3),Solution(domain))
1479     g_ex[0,0]=2.
1480     g_ex[0,1]=3.
1481     g_ex[0,2]=4.
1482     g_ex[1,0]=4.
1483     g_ex[1,1]=1.
1484     g_ex[1,2]=-2.
1485     g_ex[2,0]=8.
1486     g_ex[2,1]=4.
1487     g_ex[2,2]=5.
1488     # -------- test gradient --------------------------------
1489     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1490     # -------- set-up PDE -----------------------------------
1491     pde=LinearPDE(domain,numEquations=3)
1492     mask=whereZero(x[0])
1493     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1494     A=Tensor4(0,Function(domain))
1495     A[0,:,0,:]=kronecker(3)
1496     A[1,:,1,:]=kronecker(3)
1497     A[2,:,2,:]=kronecker(3)
1498     Y=Vector(0.,Function(domain))
1499     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1500     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1501     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1502     pde.setValue(A=A,
1503     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1504     Y=Y,
1505     y=matrixmult(g_ex,domain.getNormal()))
1506     # -------- get the solution ---------------------------
1507     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1508     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1509     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1510     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1511     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1512     u=pde.getSolution()
1513     # -------- test the solution ---------------------------
1514     error=Lsup(u-u_ex)
1515     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1516     class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1517     def test_solve(self):
1518     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1519     x=Solution(domain).getX()
1520     # --- set exact solution ----
1521     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
1522     # --- set exact gradient -----------
1523     g_ex=Data(0.,(3,),Solution(domain))
1524     g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1525     g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1526     g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1527     # -------- test gradient --------------------------------
1528     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1529     # -------- set-up PDE -----------------------------------
1530     pde=LinearPDE(domain,numEquations=1)
1531     mask=whereZero(x[0])
1532     pde.setValue(r=u_ex,q=mask)
1533     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1534     # -------- get the solution ---------------------------
1535     pde.getSolverOptions().setTolerance(SOLVER_TOL)
1536     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1537     pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1538     pde.getSolverOptions().setPackage(SolverOptions.PASO)
1539     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1540     u=pde.getSolution()
1541     # -------- test the solution ---------------------------
1542     error=Lsup(u-u_ex)
1543     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1544     class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1545     def test_solve(self):
1546     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1547     x=Solution(domain).getX()
1548     # --- set exact solution ----
1549     u_ex=Vector(0,Solution(domain))
1550     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
1551     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
1552     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
1553     # --- set exact gradient -----------
1554     g_ex=Data(0.,(3,3),Solution(domain))
1555     g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1556     g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1557     g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1558     g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1559     g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1560     g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1561     g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1562     g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1563     g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1564     # -------- test gradient --------------------------------
1565     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1566     # -------- set-up PDE -----------------------------------
1567     pde=LinearPDE(domain,numEquations=3)
1568     mask=whereZero(x[0])
1569     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1570     Y=Vector(0.,Function(domain))
1571     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1572     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1573     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1574     A=Tensor4(0,Function(domain))
1575     A[0,:,0,:]=kronecker(3)
1576     A[1,:,1,:]=kronecker(3)