/[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 3097 - (show annotations)
Fri Aug 20 04:59:12 2010 UTC (9 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 127414 byte(s)
some modifications to the GaussSeidel
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2010 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-2010 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 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)
1577 A[2,:,2,:]=kronecker(3)
1578 pde.setValue(A=A,
1579 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1580 Y=Y-numpy.array([60.,20.,22.]),
1581 y=matrixmult(g_ex,domain.getNormal()))
1582 # -------- get the solution ---------------------------
1583 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1584 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1585 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1586 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1587 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1588 u=pde.getSolution()
1589 # -------- test the solution ---------------------------
1590 error=Lsup(u-u_ex)
1591 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1592
1593 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1594 def test_solve(self):
1595 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1596 x=Solution(domain).getX()
1597 # --- set exact solution ----
1598 u_ex=Scalar(0,Solution(domain))
1599 u_ex=1.+2.*x[0]+3.*x[1]
1600 # --- set exact gradient -----------
1601 g_ex=Data(0.,(2,),Solution(domain))
1602 g_ex[0]=2.
1603 g_ex[1]=3.
1604 # -------- test gradient --------------------------------
1605 g=grad(u_ex)
1606 self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1607 # -------- set-up PDE -----------------------------------
1608 pde=LinearPDE(domain,numEquations=1)
1609 mask=whereZero(x[0])
1610 pde.setValue(r=u_ex,q=mask)
1611 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1612 # -------- get the solution ---------------------------
1613 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1614 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1615 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1616 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1617 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1618 u=pde.getSolution()
1619 # -------- test the solution ---------------------------
1620 error=Lsup(u-u_ex)
1621 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1622
1623 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1624 def test_solve(self):
1625 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1626 x=Solution(domain).getX()
1627 # --- set exact solution ----
1628 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1629 # --- set exact gradient -----------
1630 g_ex=Data(0.,(2,),Solution(domain))
1631 g_ex[0]=2.+8.*x[0]+5.*x[1]
1632 g_ex[1]=3.+5.*x[0]+12.*x[1]
1633 # -------- test gradient --------------------------------
1634 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1635 # -------- set-up PDE -----------------------------------
1636 pde=LinearPDE(domain,numEquations=1)
1637 mask=whereZero(x[0])
1638 pde.setValue(r=u_ex,q=mask)
1639 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1640 # -------- get the solution ---------------------------
1641 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1642 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1643 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1644 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1645 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1646 u=pde.getSolution()
1647 # -------- test the solution ---------------------------
1648 error=Lsup(u-u_ex)
1649 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1650
1651 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1652 def test_solve(self):
1653 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1654 x=Solution(domain).getX()
1655 # --- set exact solution ----
1656 u_ex=Vector(0,Solution(domain))
1657 u_ex[0]=1.+2.*x[0]+3.*x[1]
1658 u_ex[1]=-1.+3.*x[0]+2.*x[1]
1659 # --- set exact gradient -----------
1660 g_ex=Data(0.,(2,2),Solution(domain))
1661 g_ex[0,0]=2.
1662 g_ex[0,1]=3.
1663 g_ex[1,0]=3.
1664 g_ex[1,1]=2.
1665 # -------- test gradient --------------------------------
1666 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1667 # -------- set-up PDE -----------------------------------
1668 pde=LinearPDE(domain,numEquations=2)
1669 mask=whereZero(x[0])
1670 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1671 A=Tensor4(0,Function(domain))
1672 A[0,:,0,:]=kronecker(2)
1673 A[1,:,1,:]=kronecker(2)
1674 Y=Vector(0.,Function(domain))
1675 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1676 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1677 pde.setValue(A=A,
1678 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1679 Y=Y,
1680 y=matrixmult(g_ex,domain.getNormal()))
1681 # -------- get the solution ---------------------------
1682 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1683 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1684 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1685 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1686 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1687 u=pde.getSolution()
1688 # -------- test the solution ---------------------------
1689 error=Lsup(u-u_ex)
1690 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1691
1692 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1693 def test_solve(self):
1694 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1695 x=Solution(domain).getX()
1696 # --- set exact solution ----
1697 u_ex=Vector(0,Solution(domain))
1698 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1699 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1700 # --- set exact gradient -----------
1701 g_ex=Data(0.,(2,2),Solution(domain))
1702 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1703 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1704 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1705 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1706 # -------- test gradient --------------------------------
1707 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1708 # -------- set-up PDE -----------------------------------
1709 pde=LinearPDE(domain,numEquations=2)
1710 mask=whereZero(x[0])
1711 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1712 A=Tensor4(0,Function(domain))
1713 A[0,:,0,:]=kronecker(2)
1714 A[1,:,1,:]=kronecker(2)
1715 Y=Vector(0.,Function(domain))
1716 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1717 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1718 pde.setValue(A=A,
1719 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1720 Y=Y-[20.,10.],
1721 y=matrixmult(g_ex,domain.getNormal()))
1722 # -------- get the solution ---------------------------
1723 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1724 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1725 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1726 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1727 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1728 u=pde.getSolution()
1729 # -------- test the solution ---------------------------
1730 error=Lsup(u-u_ex)
1731 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1732
1733 class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1734 def test_solve(self):
1735 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1736 x=Solution(domain).getX()
1737 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1738 # --- set exact gradient -----------
1739 g_ex=Data(0.,(3,),Solution(domain))
1740 g_ex[0]=2.
1741 g_ex[1]=3.
1742 g_ex[2]=4.
1743 # -------- test gradient --------------------------------
1744 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1745 # -------- set-up PDE -----------------------------------
1746 pde=LinearPDE(domain,numEquations=1)
1747 mask=whereZero(x[0])
1748 pde.setValue(r=u_ex,q=mask)
1749 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1750 # -------- get the solution ---------------------------
1751 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1752 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1753 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1754 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1755 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1756 u=pde.getSolution()
1757 # -------- test the solution ---------------------------
1758 error=Lsup(u-u_ex)
1759 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1760
1761 class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1762 def test_solve(self):
1763 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1764 x=Solution(domain).getX()
1765 # --- set exact solution ----
1766 u_ex=Vector(0,Solution(domain))
1767 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1768 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1769 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1770 # --- set exact gradient -----------
1771 g_ex=Data(0.,(3,3),Solution(domain))
1772 g_ex[0,0]=2.
1773 g_ex[0,1]=3.
1774 g_ex[0,2]=4.
1775 g_ex[1,0]=4.
1776 g_ex[1,1]=1.
1777 g_ex[1,2]=-2.
1778 g_ex[2,0]=8.
1779 g_ex[2,1]=4.
1780 g_ex[2,2]=5.
1781 # -------- test gradient --------------------------------
1782 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1783 # -------- set-up PDE -----------------------------------
1784 pde=LinearPDE(domain,numEquations=3)
1785 mask=whereZero(x[0])
1786 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1787 A=Tensor4(0,Function(domain))
1788 A[0,:,0,:]=kronecker(3)
1789 A[1,:,1,:]=kronecker(3)
1790 A[2,:,2,:]=kronecker(3)
1791 Y=Vector(0.,Function(domain))
1792 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1793 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1794 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1795 pde.setValue(A=A,
1796 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1797 Y=Y,
1798 y=matrixmult(g_ex,domain.getNormal()))
1799 # -------- get the solution ---------------------------
1800 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1801 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1802 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1803 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1804 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1805 u=pde.getSolution()
1806 # -------- test the solution ---------------------------
1807 error=Lsup(u-u_ex)
1808 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1809
1810 class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1811 def test_solve(self):
1812 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1813 x=Solution(domain).getX()
1814 # --- set exact solution ----
1815 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
1816 # --- set exact gradient -----------
1817 g_ex=Data(0.,(3,),Solution(domain))
1818 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1819 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1820 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1821 # -------- test gradient --------------------------------
1822 self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1823 # -------- set-up PDE -----------------------------------
1824 pde=LinearPDE(domain,numEquations=1)
1825 mask=whereZero(x[0])
1826 pde.setValue(r=u_ex,q=mask)
1827 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1828 # -------- get the solution ---------------------------
1829 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1830 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1831 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1832 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1833 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1834 u=pde.getSolution()
1835 # -------- test the solution ---------------------------
1836 error=Lsup(u-u_ex)
1837 self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1838
1839 class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1840 def test_solve(self):
1841 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1842 x=Solution(domain).getX()
1843 # --- set exact solution ----
1844 u_ex=Vector(0,Solution(domain))
1845 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
1846 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
1847 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
1848 # --- set exact gradient -----------
1849 g_ex=Data(0.,(3,3),Solution(domain))
1850 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1851 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1852 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1853 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1854 g_ex[1,1]=1+3.*x[0]+2.*x[2