/[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 4576 - (show annotations)
Mon Dec 9 23:35:30 2013 UTC (5 years, 9 months ago) by sshaw
File MIME type: text/x-python
File size: 127490 byte(s)
python3ified things, replaced mixed whitespace and xrange calls
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 http://www.uq.edu.au
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 """
24 Test suite for the linearPDE and pdetools test on finley
25
26 :remark:
27
28 :var __author__: name of author
29 :var __licence__: licence agreement
30 :var __url__: url entry point on documentation
31 :var __version__: version
32 :var __date__: date of the version
33 """
34
35 __author__="Lutz Gross, l.gross@uq.edu.au"
36
37 import unittest, sys
38
39 from esys.escript import *
40 from esys.finley import Rectangle,Brick
41 from esys.escript.linearPDEs import LinearPDE, SolverOptions
42 import numpy
43 OPTIMIZE=True
44 SOLVER_VERBOSE=False
45 # setNumberOfThreads(2)
46
47 try:
48 FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
49 except KeyError:
50 FINLEY_TEST_DATA='.'
51
52 FINLEY_TEST_MESH_PATH=os.path.join(FINLEY_TEST_DATA,"data_meshes")
53
54 # number of elements in the spatial directions
55 NE0=8
56 NE1=10
57 NE2=12
58
59 NE0=12
60 NE1=12
61 NE2=8
62
63 SOLVER_TOL=1.e-8
64 REL_TOL=1.e-6
65
66 FAC_DIAG=1.
67 FAC_OFFDIAG=-0.4
68
69
70 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
71 def test_solve(self):
72 # Tell about how many MPI CPUs and OpenMP threads
73 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
74 x=Solution(domain).getX()
75 # --- set exact solution ----
76 u_ex=Scalar(0,Solution(domain))
77 u_ex=1.+2.*x[0]+3.*x[1]
78 # --- set exact gradient -----------
79 g_ex=Data(0.,(2,),Solution(domain))
80 g_ex[0]=2.
81 g_ex[1]=3.
82 # -------- test gradient --------------------------------
83 g=grad(u_ex)
84 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
85 # -------- set-up PDE -----------------------------------
86 pde=LinearPDE(domain,numEquations=1)
87 mask=whereZero(x[0])
88 pde.setValue(r=u_ex,q=mask)
89 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
90 # -------- get the solution ---------------------------
91 pde.getSolverOptions().setTolerance(SOLVER_TOL)
92 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
93 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
94 pde.getSolverOptions().setPackage(SolverOptions.PASO)
95 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
96 u=pde.getSolution()
97 # -------- test the solution ---------------------------
98 error=Lsup(u-u_ex)
99 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
100 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
101 def test_solve(self):
102 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
103 x=Solution(domain).getX()
104 # --- set exact solution ----
105 u_ex=Scalar(0,Solution(domain))
106 u_ex=1.+2.*x[0]+3.*x[1]
107 # --- set exact gradient -----------
108 g_ex=Data(0.,(2,),Solution(domain))
109 g_ex[0]=2.
110 g_ex[1]=3.
111 # -------- test gradient --------------------------------
112 g=grad(u_ex)
113 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
114 # -------- set-up PDE -----------------------------------
115 pde=LinearPDE(domain,numEquations=1)
116 mask=whereZero(x[0])
117 pde.setValue(r=u_ex,q=mask)
118 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
119 # -------- get the solution ---------------------------
120 pde.getSolverOptions().setTolerance(SOLVER_TOL)
121 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
122 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
123 pde.getSolverOptions().setPackage(SolverOptions.PASO)
124 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
125 u=pde.getSolution()
126 # -------- test the solution ---------------------------
127 error=Lsup(u-u_ex)
128 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
129 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
130 def test_solve(self):
131 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
132 x=Solution(domain).getX()
133 # --- set exact solution ----
134 u_ex=Vector(0,Solution(domain))
135 u_ex[0]=1.+2.*x[0]+3.*x[1]
136 u_ex[1]=-1.+3.*x[0]+2.*x[1]
137 # --- set exact gradient -----------
138 g_ex=Data(0.,(2,2),Solution(domain))
139 g_ex[0,0]=2.
140 g_ex[0,1]=3.
141 g_ex[1,0]=3.
142 g_ex[1,1]=2.
143 # -------- test gradient --------------------------------
144 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
145 # -------- set-up PDE -----------------------------------
146 pde=LinearPDE(domain,numEquations=2)
147 mask=whereZero(x[0])
148 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
149 A=Tensor4(0,Function(domain))
150 A[0,:,0,:]=kronecker(2)
151 A[1,:,1,:]=kronecker(2)
152 Y=Vector(0.,Function(domain))
153 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
154 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
155 pde.setValue(A=A,
156 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
157 Y=Y,
158 y=matrixmult(g_ex,domain.getNormal()))
159 # -------- get the solution ---------------------------
160 pde.getSolverOptions().setTolerance(SOLVER_TOL)
161 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
162 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
163 pde.getSolverOptions().setPackage(SolverOptions.PASO)
164 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
165 u=pde.getSolution()
166 # -------- test the solution ---------------------------
167 error=Lsup(u-u_ex)
168 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
169 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
170 def test_solve(self):
171 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
172 x=Solution(domain).getX()
173 # --- set exact solution ----
174 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
175 # --- set exact gradient -----------
176 g_ex=Data(0.,(2,),Solution(domain))
177 g_ex[0]=2.+8.*x[0]+5.*x[1]
178 g_ex[1]=3.+5.*x[0]+12.*x[1]
179 # -------- test gradient --------------------------------
180 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
181 # -------- set-up PDE -----------------------------------
182 pde=LinearPDE(domain,numEquations=1)
183 mask=whereZero(x[0])
184 pde.setValue(r=u_ex,q=mask)
185 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
186 # -------- get the solution ---------------------------
187 pde.getSolverOptions().setTolerance(SOLVER_TOL)
188 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
189 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
190 pde.getSolverOptions().setPackage(SolverOptions.PASO)
191 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
192 u=pde.getSolution()
193 # -------- test the solution ---------------------------
194 error=Lsup(u-u_ex)
195 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
196 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
197 def test_solve(self):
198 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
199 x=Solution(domain).getX()
200 # --- set exact solution ----
201 u_ex=Vector(0,Solution(domain))
202 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
203 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
204 # --- set exact gradient -----------
205 g_ex=Data(0.,(2,2),Solution(domain))
206 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
207 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
208 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
209 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
210 # -------- test gradient --------------------------------
211 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
212 # -------- set-up PDE -----------------------------------
213 pde=LinearPDE(domain,numEquations=2)
214 mask=whereZero(x[0])
215 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
216 A=Tensor4(0,Function(domain))
217 A[0,:,0,:]=kronecker(2)
218 A[1,:,1,:]=kronecker(2)
219 Y=Vector(0.,Function(domain))
220 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
221 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
222 pde.setValue(A=A,
223 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
224 Y=Y-[20.,10.],
225 y=matrixmult(g_ex,domain.getNormal()))
226 # -------- get the solution ---------------------------
227 pde.getSolverOptions().setTolerance(SOLVER_TOL)
228 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
229 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
230 pde.getSolverOptions().setPackage(SolverOptions.PASO)
231 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
232 u=pde.getSolution()
233 # -------- test the solution ---------------------------
234 error=Lsup(u-u_ex)
235 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
236 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
237 def test_solve(self):
238 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
239 x=Solution(domain).getX()
240 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
241 # --- set exact gradient -----------
242 g_ex=Data(0.,(3,),Solution(domain))
243 g_ex[0]=2.
244 g_ex[1]=3.
245 g_ex[2]=4.
246 # -------- test gradient --------------------------------
247 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
248 # -------- set-up PDE -----------------------------------
249 pde=LinearPDE(domain,numEquations=1)
250 mask=whereZero(x[0])
251 pde.setValue(r=u_ex,q=mask)
252 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
253 # -------- get the solution ---------------------------
254 pde.getSolverOptions().setTolerance(SOLVER_TOL)
255 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
256 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
257 pde.getSolverOptions().setPackage(SolverOptions.PASO)
258 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
259 u=pde.getSolution()
260 # -------- test the solution ---------------------------
261 error=Lsup(u-u_ex)
262 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
263 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
264 def test_solve(self):
265 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
266 x=Solution(domain).getX()
267 # --- set exact solution ----
268 u_ex=Vector(0,Solution(domain))
269 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
270 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
271 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
272 # --- set exact gradient -----------
273 g_ex=Data(0.,(3,3),Solution(domain))
274 g_ex[0,0]=2.
275 g_ex[0,1]=3.
276 g_ex[0,2]=4.
277 g_ex[1,0]=4.
278 g_ex[1,1]=1.
279 g_ex[1,2]=-2.
280 g_ex[2,0]=8.
281 g_ex[2,1]=4.
282 g_ex[2,2]=5.
283 # -------- test gradient --------------------------------
284 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
285 # -------- set-up PDE -----------------------------------
286 pde=LinearPDE(domain,numEquations=3)
287 mask=whereZero(x[0])
288 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
289 A=Tensor4(0,Function(domain))
290 A[0,:,0,:]=kronecker(3)
291 A[1,:,1,:]=kronecker(3)
292 A[2,:,2,:]=kronecker(3)
293 Y=Vector(0.,Function(domain))
294 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
295 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
296 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
297 pde.setValue(A=A,
298 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
299 Y=Y,
300 y=matrixmult(g_ex,domain.getNormal()))
301 # -------- get the solution ---------------------------
302 pde.getSolverOptions().setTolerance(SOLVER_TOL)
303 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
304 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
305 pde.getSolverOptions().setPackage(SolverOptions.PASO)
306 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
307 u=pde.getSolution()
308 # -------- test the solution ---------------------------
309 error=Lsup(u-u_ex)
310 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
311 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
312 def test_solve(self):
313 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
314 x=Solution(domain).getX()
315 # --- set exact solution ----
316 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
317 # --- set exact gradient -----------
318 g_ex=Data(0.,(3,),Solution(domain))
319 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
320 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
321 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
322 # -------- test gradient --------------------------------
323 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
324 # -------- set-up PDE -----------------------------------
325 pde=LinearPDE(domain,numEquations=1)
326 mask=whereZero(x[0])
327 pde.setValue(r=u_ex,q=mask)
328 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
329 # -------- get the solution ---------------------------
330 pde.getSolverOptions().setTolerance(SOLVER_TOL)
331 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
332 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
333 pde.getSolverOptions().setPackage(SolverOptions.PASO)
334 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
335 u=pde.getSolution()
336 # -------- test the solution ---------------------------
337 error=Lsup(u-u_ex)
338 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
339 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
340 def test_solve(self):
341 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
342 x=Solution(domain).getX()
343 # --- set exact solution ----
344 u_ex=Vector(0,Solution(domain))
345 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
346 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
347 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
348 # --- set exact gradient -----------
349 g_ex=Data(0.,(3,3),Solution(domain))
350 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
351 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
352 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
353 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
354 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
355 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
356 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
357 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
358 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
359 # -------- test gradient --------------------------------
360 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
361 # -------- set-up PDE -----------------------------------
362 pde=LinearPDE(domain,numEquations=3)
363 mask=whereZero(x[0])
364 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
365 Y=Vector(0.,Function(domain))
366 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
367 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
368 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
369 A=Tensor4(0,Function(domain))
370 A[0,:,0,:]=kronecker(3)
371 A[1,:,1,:]=kronecker(3)
372 A[2,:,2,:]=kronecker(3)
373 pde.setValue(A=A,
374 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
375 Y=Y-numpy.array([60.,20.,22.]),
376 y=matrixmult(g_ex,domain.getNormal()))
377 # -------- get the solution ---------------------------
378 pde.getSolverOptions().setTolerance(SOLVER_TOL)
379 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
380 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
381 pde.getSolverOptions().setPackage(SolverOptions.PASO)
382 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
383 u=pde.getSolution()
384 # -------- test the solution ---------------------------
385 error=Lsup(u-u_ex)
386 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
387
388 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
389 def test_solve(self):
390 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
391 x=Solution(domain).getX()
392 # --- set exact solution ----
393 u_ex=Scalar(0,Solution(domain))
394 u_ex=1.+2.*x[0]+3.*x[1]
395 # --- set exact gradient -----------
396 g_ex=Data(0.,(2,),Solution(domain))
397 g_ex[0]=2.
398 g_ex[1]=3.
399 # -------- test gradient --------------------------------
400 g=grad(u_ex)
401 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
402 # -------- set-up PDE -----------------------------------
403 pde=LinearPDE(domain,numEquations=1)
404 mask=whereZero(x[0])
405 pde.setValue(r=u_ex,q=mask)
406 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
407 # -------- get the solution ---------------------------
408 pde.getSolverOptions().setTolerance(SOLVER_TOL)
409 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
410 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
411 pde.getSolverOptions().setPackage(SolverOptions.PASO)
412 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
413 u=pde.getSolution()
414 # -------- test the solution ---------------------------
415 error=Lsup(u-u_ex)
416 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
417
418 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
419 def test_solve(self):
420 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
421 x=Solution(domain).getX()
422 # --- set exact solution ----
423 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
424 # --- set exact gradient -----------
425 g_ex=Data(0.,(2,),Solution(domain))
426 g_ex[0]=2.+8.*x[0]+5.*x[1]
427 g_ex[1]=3.+5.*x[0]+12.*x[1]
428 # -------- test gradient --------------------------------
429 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
430 # -------- set-up PDE -----------------------------------
431 pde=LinearPDE(domain,numEquations=1)
432 mask=whereZero(x[0])
433 pde.setValue(r=u_ex,q=mask)
434 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
435 # -------- get the solution ---------------------------
436 pde.getSolverOptions().setTolerance(SOLVER_TOL)
437 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
438 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
439 pde.getSolverOptions().setPackage(SolverOptions.PASO)
440 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
441 u=pde.getSolution()
442 # -------- test the solution ---------------------------
443 error=Lsup(u-u_ex)
444 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
445
446 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
447 def test_solve(self):
448 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
449 x=Solution(domain).getX()
450 # --- set exact solution ----
451 u_ex=Vector(0,Solution(domain))
452 u_ex[0]=1.+2.*x[0]+3.*x[1]
453 u_ex[1]=-1.+3.*x[0]+2.*x[1]
454 # --- set exact gradient -----------
455 g_ex=Data(0.,(2,2),Solution(domain))
456 g_ex[0,0]=2.
457 g_ex[0,1]=3.
458 g_ex[1,0]=3.
459 g_ex[1,1]=2.
460 # -------- test gradient --------------------------------
461 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
462 # -------- set-up PDE -----------------------------------
463 pde=LinearPDE(domain,numEquations=2)
464 mask=whereZero(x[0])
465 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
466 A=Tensor4(0,Function(domain))
467 A[0,:,0,:]=kronecker(2)
468 A[1,:,1,:]=kronecker(2)
469 Y=Vector(0.,Function(domain))
470 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
471 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
472 pde.setValue(A=A,
473 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
474 Y=Y,
475 y=matrixmult(g_ex,domain.getNormal()))
476 # -------- get the solution ---------------------------
477 pde.getSolverOptions().setTolerance(SOLVER_TOL)
478 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
479 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
480 pde.getSolverOptions().setPackage(SolverOptions.PASO)
481 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
482 u=pde.getSolution()
483 # -------- test the solution ---------------------------
484 error=Lsup(u-u_ex)
485 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
486
487 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
488 def test_solve(self):
489 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
490 x=Solution(domain).getX()
491 # --- set exact solution ----
492 u_ex=Vector(0,Solution(domain))
493 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
494 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
495 # --- set exact gradient -----------
496 g_ex=Data(0.,(2,2),Solution(domain))
497 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
498 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
499 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
500 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
501 # -------- test gradient --------------------------------
502 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
503 # -------- set-up PDE -----------------------------------
504 pde=LinearPDE(domain,numEquations=2)
505 mask=whereZero(x[0])
506 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
507 A=Tensor4(0,Function(domain))
508 A[0,:,0,:]=kronecker(2)
509 A[1,:,1,:]=kronecker(2)
510 Y=Vector(0.,Function(domain))
511 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
512 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
513 pde.setValue(A=A,
514 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
515 Y=Y-[20.,10.],
516 y=matrixmult(g_ex,domain.getNormal()))
517 # -------- get the solution ---------------------------
518 pde.getSolverOptions().setTolerance(SOLVER_TOL)
519 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
520 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
521 pde.getSolverOptions().setPackage(SolverOptions.PASO)
522 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
523 u=pde.getSolution()
524 # -------- test the solution ---------------------------
525 error=Lsup(u-u_ex)
526 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
527
528 class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
529 def test_solve(self):
530 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
531 x=Solution(domain).getX()
532 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
533 # --- set exact gradient -----------
534 g_ex=Data(0.,(3,),Solution(domain))
535 g_ex[0]=2.
536 g_ex[1]=3.
537 g_ex[2]=4.
538 # -------- test gradient --------------------------------
539 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
540 # -------- set-up PDE -----------------------------------
541 pde=LinearPDE(domain,numEquations=1)
542 mask=whereZero(x[0])
543 pde.setValue(r=u_ex,q=mask)
544 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
545 # -------- get the solution ---------------------------
546 pde.getSolverOptions().setTolerance(SOLVER_TOL)
547 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
548 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
549 pde.getSolverOptions().setPackage(SolverOptions.PASO)
550 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
551 u=pde.getSolution()
552 # -------- test the solution ---------------------------
553 error=Lsup(u-u_ex)
554 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
555
556 class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
557 def test_solve(self):
558 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
559 x=Solution(domain).getX()
560 # --- set exact solution ----
561 u_ex=Vector(0,Solution(domain))
562 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
563 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
564 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
565 # --- set exact gradient -----------
566 g_ex=Data(0.,(3,3),Solution(domain))
567 g_ex[0,0]=2.
568 g_ex[0,1]=3.
569 g_ex[0,2]=4.
570 g_ex[1,0]=4.
571 g_ex[1,1]=1.
572 g_ex[1,2]=-2.
573 g_ex[2,0]=8.
574 g_ex[2,1]=4.
575 g_ex[2,2]=5.
576 # -------- test gradient --------------------------------
577 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
578 # -------- set-up PDE -----------------------------------
579 pde=LinearPDE(domain,numEquations=3)
580 mask=whereZero(x[0])
581 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
582 A=Tensor4(0,Function(domain))
583 A[0,:,0,:]=kronecker(3)
584 A[1,:,1,:]=kronecker(3)
585 A[2,:,2,:]=kronecker(3)
586 Y=Vector(0.,Function(domain))
587 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
588 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
589 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
590 pde.setValue(A=A,
591 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
592 Y=Y,
593 y=matrixmult(g_ex,domain.getNormal()))
594 # -------- get the solution ---------------------------
595 pde.getSolverOptions().setTolerance(SOLVER_TOL)
596 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
597 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
598 pde.getSolverOptions().setPackage(SolverOptions.PASO)
599 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
600 u=pde.getSolution()
601 # -------- test the solution ---------------------------
602 error=Lsup(u-u_ex)
603 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
604
605 class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
606 def test_solve(self):
607 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
608 x=Solution(domain).getX()
609 # --- set exact solution ----
610 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
611 # --- set exact gradient -----------
612 g_ex=Data(0.,(3,),Solution(domain))
613 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
614 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
615 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
616 # -------- test gradient --------------------------------
617 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
618 # -------- set-up PDE -----------------------------------
619 pde=LinearPDE(domain,numEquations=1)
620 mask=whereZero(x[0])
621 pde.setValue(r=u_ex,q=mask)
622 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
623 # -------- get the solution ---------------------------
624 pde.getSolverOptions().setTolerance(SOLVER_TOL)
625 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
626 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
627 pde.getSolverOptions().setPackage(SolverOptions.PASO)
628 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
629 u=pde.getSolution()
630 # -------- test the solution ---------------------------
631 error=Lsup(u-u_ex)
632 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
633
634 class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_Jacobi(unittest.TestCase):
635 def test_solve(self):
636 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
637 x=Solution(domain).getX()
638 # --- set exact solution ----
639 u_ex=Vector(0,Solution(domain))
640 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
641 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
642 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
643 # --- set exact gradient -----------
644 g_ex=Data(0.,(3,3),Solution(domain))
645 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
646 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
647 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
648 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
649 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
650 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
651 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
652 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
653 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
654 # -------- test gradient --------------------------------
655 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
656 # -------- set-up PDE -----------------------------------
657 pde=LinearPDE(domain,numEquations=3)
658 mask=whereZero(x[0])
659 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
660 Y=Vector(0.,Function(domain))
661 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
662 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
663 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
664 A=Tensor4(0,Function(domain))
665 A[0,:,0,:]=kronecker(3)
666 A[1,:,1,:]=kronecker(3)
667 A[2,:,2,:]=kronecker(3)
668 pde.setValue(A=A,
669 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
670 Y=Y-numpy.array([60.,20.,22.]),
671 y=matrixmult(g_ex,domain.getNormal()))
672 # -------- get the solution ---------------------------
673 pde.getSolverOptions().setTolerance(SOLVER_TOL)
674 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
675 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
676 pde.getSolverOptions().setPackage(SolverOptions.PASO)
677 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
678 u=pde.getSolution()
679 # -------- test the solution ---------------------------
680 error=Lsup(u-u_ex)
681 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
682
683
684 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
685 def test_solve(self):
686 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
687 x=Solution(domain).getX()
688 # --- set exact solution ----
689 u_ex=Scalar(0,Solution(domain))
690 u_ex=1.+2.*x[0]+3.*x[1]
691 # --- set exact gradient -----------
692 g_ex=Data(0.,(2,),Solution(domain))
693 g_ex[0]=2.
694 g_ex[1]=3.
695 # -------- test gradient --------------------------------
696 g=grad(u_ex)
697 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
698 # -------- set-up PDE -----------------------------------
699 pde=LinearPDE(domain,numEquations=1)
700 mask=whereZero(x[0])
701 pde.setValue(r=u_ex,q=mask)
702 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
703 # -------- get the solution ---------------------------
704 pde.getSolverOptions().setTolerance(SOLVER_TOL)
705 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
706 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
707 pde.getSolverOptions().setPackage(SolverOptions.PASO)
708
709 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
710 u=pde.getSolution()
711 # -------- test the solution ---------------------------
712 error=Lsup(u-u_ex)
713 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
714
715 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
716 def test_solve(self):
717 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
718 x=Solution(domain).getX()
719 # --- set exact solution ----
720 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
721 # --- set exact gradient -----------
722 g_ex=Data(0.,(2,),Solution(domain))
723 g_ex[0]=2.+8.*x[0]+5.*x[1]
724 g_ex[1]=3.+5.*x[0]+12.*x[1]
725 # -------- test gradient --------------------------------
726 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
727 # -------- set-up PDE -----------------------------------
728 pde=LinearPDE(domain,numEquations=1)
729 mask=whereZero(x[0])
730 pde.setValue(r=u_ex,q=mask)
731 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
732 # -------- get the solution ---------------------------
733 pde.getSolverOptions().setTolerance(SOLVER_TOL)
734 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
735 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
736 pde.getSolverOptions().setPackage(SolverOptions.PASO)
737 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
738 u=pde.getSolution()
739 # -------- test the solution ---------------------------
740 error=Lsup(u-u_ex)
741 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
742
743 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
744 def test_solve(self):
745 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
746 x=Solution(domain).getX()
747 # --- set exact solution ----
748 u_ex=Vector(0,Solution(domain))
749 u_ex[0]=1.+2.*x[0]+3.*x[1]
750 u_ex[1]=-1.+3.*x[0]+2.*x[1]
751 # --- set exact gradient -----------
752 g_ex=Data(0.,(2,2),Solution(domain))
753 g_ex[0,0]=2.
754 g_ex[0,1]=3.
755 g_ex[1,0]=3.
756 g_ex[1,1]=2.
757 # -------- test gradient --------------------------------
758 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
759 # -------- set-up PDE -----------------------------------
760 pde=LinearPDE(domain,numEquations=2)
761 mask=whereZero(x[0])
762 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
763 A=Tensor4(0,Function(domain))
764 A[0,:,0,:]=kronecker(2)
765 A[1,:,1,:]=kronecker(2)
766 Y=Vector(0.,Function(domain))
767 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
768 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
769 pde.setValue(A=A,
770 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
771 Y=Y,
772 y=matrixmult(g_ex,domain.getNormal()))
773 # -------- get the solution ---------------------------
774 pde.getSolverOptions().setTolerance(SOLVER_TOL)
775 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
776 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
777 pde.getSolverOptions().setPackage(SolverOptions.PASO)
778 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
779 u=pde.getSolution()
780 # -------- test the solution ---------------------------
781 error=Lsup(u-u_ex)
782 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
783
784 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
785 def test_solve(self):
786 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
787 x=Solution(domain).getX()
788 # --- set exact solution ----
789 u_ex=Vector(0,Solution(domain))
790 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
791 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
792 # --- set exact gradient -----------
793 g_ex=Data(0.,(2,2),Solution(domain))
794 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
795 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
796 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
797 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
798 # -------- test gradient --------------------------------
799 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
800 # -------- set-up PDE -----------------------------------
801 pde=LinearPDE(domain,numEquations=2)
802 mask=whereZero(x[0])
803 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
804 A=Tensor4(0,Function(domain))
805 A[0,:,0,:]=kronecker(2)
806 A[1,:,1,:]=kronecker(2)
807 Y=Vector(0.,Function(domain))
808 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
809 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
810 pde.setValue(A=A,
811 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
812 Y=Y-[20.,10.],
813 y=matrixmult(g_ex,domain.getNormal()))
814 # -------- get the solution ---------------------------
815 pde.getSolverOptions().setTolerance(SOLVER_TOL)
816 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
817 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
818 pde.getSolverOptions().setPackage(SolverOptions.PASO)
819 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
820 u=pde.getSolution()
821 # -------- test the solution ---------------------------
822 error=Lsup(u-u_ex)
823 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
824
825 class SimpleSolve_Brick_Order1_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
826 def test_solve(self):
827 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
828 x=Solution(domain).getX()
829 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
830 # --- set exact gradient -----------
831 g_ex=Data(0.,(3,),Solution(domain))
832 g_ex[0]=2.
833 g_ex[1]=3.
834 g_ex[2]=4.
835 # -------- test gradient --------------------------------
836 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
837 # -------- set-up PDE -----------------------------------
838 pde=LinearPDE(domain,numEquations=1)
839 mask=whereZero(x[0])
840 pde.setValue(r=u_ex,q=mask)
841 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
842 # -------- get the solution ---------------------------
843 pde.getSolverOptions().setTolerance(SOLVER_TOL)
844 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
845 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
846 pde.getSolverOptions().setPackage(SolverOptions.PASO)
847 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
848 u=pde.getSolution()
849 # -------- test the solution ---------------------------
850 error=Lsup(u-u_ex)
851 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
852
853 class SimpleSolve_Brick_Order1_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
854 def test_solve(self):
855 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
856 x=Solution(domain).getX()
857 # --- set exact solution ----
858 u_ex=Vector(0,Solution(domain))
859 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
860 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
861 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
862 # --- set exact gradient -----------
863 g_ex=Data(0.,(3,3),Solution(domain))
864 g_ex[0,0]=2.
865 g_ex[0,1]=3.
866 g_ex[0,2]=4.
867 g_ex[1,0]=4.
868 g_ex[1,1]=1.
869 g_ex[1,2]=-2.
870 g_ex[2,0]=8.
871 g_ex[2,1]=4.
872 g_ex[2,2]=5.
873 # -------- test gradient --------------------------------
874 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
875 # -------- set-up PDE -----------------------------------
876 pde=LinearPDE(domain,numEquations=3)
877 mask=whereZero(x[0])
878 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
879 A=Tensor4(0,Function(domain))
880 A[0,:,0,:]=kronecker(3)
881 A[1,:,1,:]=kronecker(3)
882 A[2,:,2,:]=kronecker(3)
883 Y=Vector(0.,Function(domain))
884 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
885 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
886 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
887 pde.setValue(A=A,
888 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
889 Y=Y,
890 y=matrixmult(g_ex,domain.getNormal()))
891 # -------- get the solution ---------------------------
892 pde.getSolverOptions().setTolerance(SOLVER_TOL)
893 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
894 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
895 pde.getSolverOptions().setPackage(SolverOptions.PASO)
896 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
897 u=pde.getSolution()
898 # -------- test the solution ---------------------------
899 error=Lsup(u-u_ex)
900 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
901
902 class SimpleSolve_Brick_Order2_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
903 def test_solve(self):
904 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
905 x=Solution(domain).getX()
906 # --- set exact solution ----
907 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
908 # --- set exact gradient -----------
909 g_ex=Data(0.,(3,),Solution(domain))
910 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
911 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
912 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
913 # -------- test gradient --------------------------------
914 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
915 # -------- set-up PDE -----------------------------------
916 pde=LinearPDE(domain,numEquations=1)
917 mask=whereZero(x[0])
918 pde.setValue(r=u_ex,q=mask)
919 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
920 # -------- get the solution ---------------------------
921 pde.getSolverOptions().setTolerance(SOLVER_TOL)
922 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
923 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
924 pde.getSolverOptions().setPackage(SolverOptions.PASO)
925 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
926 u=pde.getSolution()
927 # -------- test the solution ---------------------------
928 error=Lsup(u-u_ex)
929 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
930
931 class SimpleSolve_Brick_Order2_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
932 def test_solve(self):
933 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
934 x=Solution(domain).getX()
935 # --- set exact solution ----
936 u_ex=Vector(0,Solution(domain))
937 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
938 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
939 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
940 # --- set exact gradient -----------
941 g_ex=Data(0.,(3,3),Solution(domain))
942 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
943 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
944 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
945 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
946 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
947 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
948 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
949 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
950 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
951 # -------- test gradient --------------------------------
952 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
953 # -------- set-up PDE -----------------------------------
954 pde=LinearPDE(domain,numEquations=3)
955 mask=whereZero(x[0])
956 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
957 Y=Vector(0.,Function(domain))
958 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
959 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
960 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
961 A=Tensor4(0,Function(domain))
962 A[0,:,0,:]=kronecker(3)
963 A[1,:,1,:]=kronecker(3)
964 A[2,:,2,:]=kronecker(3)
965 pde.setValue(A=A,
966 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
967 Y=Y-numpy.array([60.,20.,22.]),
968 y=matrixmult(g_ex,domain.getNormal()))
969 # -------- get the solution ---------------------------
970 pde.getSolverOptions().setTolerance(SOLVER_TOL)
971 pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
972 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
973 pde.getSolverOptions().setPackage(SolverOptions.PASO)
974 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
975 u=pde.getSolution()
976 # -------- test the solution ---------------------------
977 error=Lsup(u-u_ex)
978 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
979
980 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
981 def test_solve(self):
982 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
983 x=Solution(domain).getX()
984 # --- set exact solution ----
985 u_ex=Scalar(0,Solution(domain))
986 u_ex=1.+2.*x[0]+3.*x[1]
987 # --- set exact gradient -----------
988 g_ex=Data(0.,(2,),Solution(domain))
989 g_ex[0]=2.
990 g_ex[1]=3.
991 # -------- test gradient --------------------------------
992 g=grad(u_ex)
993 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
994 # -------- set-up PDE -----------------------------------
995 pde=LinearPDE(domain,numEquations=1)
996 mask=whereZero(x[0])
997 pde.setValue(r=u_ex,q=mask)
998 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
999 # -------- get the solution ---------------------------
1000 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1001 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1002 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1003 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1004 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1005 u=pde.getSolution()
1006 # -------- test the solution ---------------------------
1007 error=Lsup(u-u_ex)
1008 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1009
1010 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1011 def test_solve(self):
1012 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1013 x=Solution(domain).getX()
1014 # --- set exact solution ----
1015 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1016 # --- set exact gradient -----------
1017 g_ex=Data(0.,(2,),Solution(domain))
1018 g_ex[0]=2.+8.*x[0]+5.*x[1]
1019 g_ex[1]=3.+5.*x[0]+12.*x[1]
1020 # -------- test gradient --------------------------------
1021 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1022 # -------- set-up PDE -----------------------------------
1023 pde=LinearPDE(domain,numEquations=1)
1024 mask=whereZero(x[0])
1025 pde.setValue(r=u_ex,q=mask)
1026 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1027 # -------- get the solution ---------------------------
1028 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1029 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1030 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1031 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1032 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1033 u=pde.getSolution()
1034 # -------- test the solution ---------------------------
1035 error=Lsup(u-u_ex)
1036 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1037
1038 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1039 def test_solve(self):
1040 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1041 x=Solution(domain).getX()
1042 # --- set exact solution ----
1043 u_ex=Vector(0,Solution(domain))
1044 u_ex[0]=1.+2.*x[0]+3.*x[1]
1045 u_ex[1]=-1.+3.*x[0]+2.*x[1]
1046 # --- set exact gradient -----------
1047 g_ex=Data(0.,(2,2),Solution(domain))
1048 g_ex[0,0]=2.
1049 g_ex[0,1]=3.
1050 g_ex[1,0]=3.
1051 g_ex[1,1]=2.
1052 # -------- test gradient --------------------------------
1053 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1054 # -------- set-up PDE -----------------------------------
1055 pde=LinearPDE(domain,numEquations=2)
1056 mask=whereZero(x[0])
1057 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1058 A=Tensor4(0,Function(domain))
1059 A[0,:,0,:]=kronecker(2)
1060 A[1,:,1,:]=kronecker(2)
1061 Y=Vector(0.,Function(domain))
1062 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1063 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1064 pde.setValue(A=A,
1065 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1066 Y=Y,
1067 y=matrixmult(g_ex,domain.getNormal()))
1068 # -------- get the solution ---------------------------
1069 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1070 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1071 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1072 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1073 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1074 u=pde.getSolution()
1075 # -------- test the solution ---------------------------
1076 error=Lsup(u-u_ex)
1077 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1078
1079 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1080 def test_solve(self):
1081 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1082 x=Solution(domain).getX()
1083 # --- set exact solution ----
1084 u_ex=Vector(0,Solution(domain))
1085 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1086 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1087 # --- set exact gradient -----------
1088 g_ex=Data(0.,(2,2),Solution(domain))
1089 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1090 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1091 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1092 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1093 # -------- test gradient --------------------------------
1094 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1095 # -------- set-up PDE -----------------------------------
1096 pde=LinearPDE(domain,numEquations=2)
1097 mask=whereZero(x[0])
1098 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1099 A=Tensor4(0,Function(domain))
1100 A[0,:,0,:]=kronecker(2)
1101 A[1,:,1,:]=kronecker(2)
1102 Y=Vector(0.,Function(domain))
1103 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1104 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1105 pde.setValue(A=A,
1106 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1107 Y=Y-[20.,10.],
1108 y=matrixmult(g_ex,domain.getNormal()))
1109 # -------- get the solution ---------------------------
1110 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1111 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1112 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1113 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1114 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1115 u=pde.getSolution()
1116 # -------- test the solution ---------------------------
1117 error=Lsup(u-u_ex)
1118 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1119
1120 class SimpleSolve_Brick_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1121 def test_solve(self):
1122 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1123 x=Solution(domain).getX()
1124 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1125 # --- set exact gradient -----------
1126 g_ex=Data(0.,(3,),Solution(domain))
1127 g_ex[0]=2.
1128 g_ex[1]=3.
1129 g_ex[2]=4.
1130 # -------- test gradient --------------------------------
1131 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1132 # -------- set-up PDE -----------------------------------
1133 pde=LinearPDE(domain,numEquations=1)
1134 mask=whereZero(x[0])
1135 pde.setValue(r=u_ex,q=mask)
1136 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1137 # -------- get the solution ---------------------------
1138 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1139 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1140 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1141 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1142 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1143 u=pde.getSolution()
1144 # -------- test the solution ---------------------------
1145 error=Lsup(u-u_ex)
1146 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1147
1148 class SimpleSolve_Brick_Order1_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1149 def test_solve(self):
1150 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1151 x=Solution(domain).getX()
1152 # --- set exact solution ----
1153 u_ex=Vector(0,Solution(domain))
1154 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1155 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1156 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1157 # --- set exact gradient -----------
1158 g_ex=Data(0.,(3,3),Solution(domain))
1159 g_ex[0,0]=2.
1160 g_ex[0,1]=3.
1161 g_ex[0,2]=4.
1162 g_ex[1,0]=4.
1163 g_ex[1,1]=1.
1164 g_ex[1,2]=-2.
1165 g_ex[2,0]=8.
1166 g_ex[2,1]=4.
1167 g_ex[2,2]=5.
1168 # -------- test gradient --------------------------------
1169 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1170 # -------- set-up PDE -----------------------------------
1171 pde=LinearPDE(domain,numEquations=3)
1172 mask=whereZero(x[0])
1173 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1174 A=Tensor4(0,Function(domain))
1175 A[0,:,0,:]=kronecker(3)
1176 A[1,:,1,:]=kronecker(3)
1177 A[2,:,2,:]=kronecker(3)
1178 Y=Vector(0.,Function(domain))
1179 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1180 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1181 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1182 pde.setValue(A=A,
1183 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1184 Y=Y,
1185 y=matrixmult(g_ex,domain.getNormal()))
1186 # -------- get the solution ---------------------------
1187 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1188 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1189 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1190 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1191 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1192 u=pde.getSolution()
1193 # -------- test the solution ---------------------------
1194 error=Lsup(u-u_ex)
1195 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1196
1197 class SimpleSolve_Brick_Order2_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1198 def test_solve(self):
1199 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1200 x=Solution(domain).getX()
1201 # --- set exact solution ----
1202 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
1203 # --- set exact gradient -----------
1204 g_ex=Data(0.,(3,),Solution(domain))
1205 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1206 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1207 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1208 # -------- test gradient --------------------------------
1209 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1210 # -------- set-up PDE -----------------------------------
1211 pde=LinearPDE(domain,numEquations=1)
1212 mask=whereZero(x[0])
1213 pde.setValue(r=u_ex,q=mask)
1214 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1215 # -------- get the solution ---------------------------
1216 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1217 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1218 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1219 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1220 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1221 u=pde.getSolution()
1222 # -------- test the solution ---------------------------
1223 error=Lsup(u-u_ex)
1224 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1225
1226 class SimpleSolve_Brick_Order2_SystemPDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
1227 def test_solve(self):
1228 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1229 x=Solution(domain).getX()
1230 # --- set exact solution ----
1231 u_ex=Vector(0,Solution(domain))
1232 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
1233 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
1234 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
1235 # --- set exact gradient -----------
1236 g_ex=Data(0.,(3,3),Solution(domain))
1237 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1238 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1239 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1240 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1241 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1242 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1243 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1244 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1245 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1246 # -------- test gradient --------------------------------
1247 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1248 # -------- set-up PDE -----------------------------------
1249 pde=LinearPDE(domain,numEquations=3)
1250 mask=whereZero(x[0])
1251 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1252 Y=Vector(0.,Function(domain))
1253 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1254 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1255 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1256 A=Tensor4(0,Function(domain))
1257 A[0,:,0,:]=kronecker(3)
1258 A[1,:,1,:]=kronecker(3)
1259 A[2,:,2,:]=kronecker(3)
1260 pde.setValue(A=A,
1261 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1262 Y=Y-numpy.array([60.,20.,22.]),
1263 y=matrixmult(g_ex,domain.getNormal()))
1264 # -------- get the solution ---------------------------
1265 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1266 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1267 pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
1268 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1269 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1270 u=pde.getSolution()
1271 # -------- test the solution ---------------------------
1272 error=Lsup(u-u_ex)
1273 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1274
1275
1276 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_GaussSeidel(unittest.TestCase):
1277 def test_solve(self):
1278 # Tell about how many MPI CPUs and OpenMP threads
1279 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1280 x=Solution(domain).getX()
1281 # --- set exact solution ----
1282 u_ex=Scalar(0,Solution(domain))
1283 u_ex=1.+2.*x[0]+3.*x[1]
1284 # --- set exact gradient -----------
1285 g_ex=Data(0.,(2,),Solution(domain))
1286 g_ex[0]=2.
1287 g_ex[1]=3.
1288 # -------- test gradient --------------------------------
1289 g=grad(u_ex)
1290 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1291 # -------- set-up PDE -----------------------------------
1292 pde=LinearPDE(domain,numEquations=1)
1293 mask=whereZero(x[0])
1294 pde.setValue(r=u_ex,q=mask)
1295 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1296 # -------- get the solution ---------------------------
1297 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1298 pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
1299 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1300 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1301 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1302 u=pde.getSolution()
1303 # -------- test the solution ---------------------------
1304 error=Lsup(u-u_ex)
1305 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1306 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1307 def test_solve(self):
1308 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1309 x=Solution(domain).getX()
1310 # --- set exact solution ----
1311 u_ex=Scalar(0,Solution(domain))
1312 u_ex=1.+2.*x[0]+3.*x[1]
1313 # --- set exact gradient -----------
1314 g_ex=Data(0.,(2,),Solution(domain))
1315 g_ex[0]=2.
1316 g_ex[1]=3.
1317 # -------- test gradient --------------------------------
1318 g=grad(u_ex)
1319 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1320 # -------- set-up PDE -----------------------------------
1321 pde=LinearPDE(domain,numEquations=1)
1322 mask=whereZero(x[0])
1323 pde.setValue(r=u_ex,q=mask)
1324 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1325 # -------- get the solution ---------------------------
1326 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1327 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1328 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1329 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1330 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1331 u=pde.getSolution()
1332 # -------- test the solution ---------------------------
1333 error=Lsup(u-u_ex)
1334 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1335 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1336 def test_solve(self):
1337 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1338 x=Solution(domain).getX()
1339 # --- set exact solution ----
1340 u_ex=Vector(0,Solution(domain))
1341 u_ex[0]=1.+2.*x[0]+3.*x[1]
1342 u_ex[1]=-1.+3.*x[0]+2.*x[1]
1343 # --- set exact gradient -----------
1344 g_ex=Data(0.,(2,2),Solution(domain))
1345 g_ex[0,0]=2.
1346 g_ex[0,1]=3.
1347 g_ex[1,0]=3.
1348 g_ex[1,1]=2.
1349 # -------- test gradient --------------------------------
1350 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1351 # -------- set-up PDE -----------------------------------
1352 pde=LinearPDE(domain,numEquations=2)
1353 mask=whereZero(x[0])
1354 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1355 A=Tensor4(0,Function(domain))
1356 A[0,:,0,:]=kronecker(2)
1357 A[1,:,1,:]=kronecker(2)
1358 Y=Vector(0.,Function(domain))
1359 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1360 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1361 pde.setValue(A=A,
1362 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1363 Y=Y,
1364 y=matrixmult(g_ex,domain.getNormal()))
1365 # -------- get the solution ---------------------------
1366 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1367 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1368 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1369 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1370 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1371 u=pde.getSolution()
1372 # -------- test the solution ---------------------------
1373 error=Lsup(u-u_ex)
1374 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1375 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1376 def test_solve(self):
1377 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1378 x=Solution(domain).getX()
1379 # --- set exact solution ----
1380 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1381 # --- set exact gradient -----------
1382 g_ex=Data(0.,(2,),Solution(domain))
1383 g_ex[0]=2.+8.*x[0]+5.*x[1]
1384 g_ex[1]=3.+5.*x[0]+12.*x[1]
1385 # -------- test gradient --------------------------------
1386 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1387 # -------- set-up PDE -----------------------------------
1388 pde=LinearPDE(domain,numEquations=1)
1389 mask=whereZero(x[0])
1390 pde.setValue(r=u_ex,q=mask)
1391 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1392 # -------- get the solution ---------------------------
1393 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1394 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1395 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1396 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1397 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1398 u=pde.getSolution()
1399 # -------- test the solution ---------------------------
1400 error=Lsup(u-u_ex)
1401 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1402 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1403 def test_solve(self):
1404 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1405 x=Solution(domain).getX()
1406 # --- set exact solution ----
1407 u_ex=Vector(0,Solution(domain))
1408 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1409 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1410 # --- set exact gradient -----------
1411 g_ex=Data(0.,(2,2),Solution(domain))
1412 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1413 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1414 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1415 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1416 # -------- test gradient --------------------------------
1417 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1418 # -------- set-up PDE -----------------------------------
1419 pde=LinearPDE(domain,numEquations=2)
1420 mask=whereZero(x[0])
1421 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1422 A=Tensor4(0,Function(domain))
1423 A[0,:,0,:]=kronecker(2)
1424 A[1,:,1,:]=kronecker(2)
1425 Y=Vector(0.,Function(domain))
1426 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1427 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1428 pde.setValue(A=A,
1429 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1430 Y=Y-[20.,10.],
1431 y=matrixmult(g_ex,domain.getNormal()))
1432 # -------- get the solution ---------------------------
1433 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1434 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1435 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1436 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1437 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1438 u=pde.getSolution()
1439 # -------- test the solution ---------------------------
1440 error=Lsup(u-u_ex)
1441 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1442 class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1443 def test_solve(self):
1444 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1445 x=Solution(domain).getX()
1446 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1447 # --- set exact gradient -----------
1448 g_ex=Data(0.,(3,),Solution(domain))
1449 g_ex[0]=2.
1450 g_ex[1]=3.
1451 g_ex[2]=4.
1452 # -------- test gradient --------------------------------
1453 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1454 # -------- set-up PDE -----------------------------------
1455 pde=LinearPDE(domain,numEquations=1)
1456 mask=whereZero(x[0])
1457 pde.setValue(r=u_ex,q=mask)
1458 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1459 # -------- get the solution ---------------------------
1460 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1461 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1462 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1463 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1464 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1465 u=pde.getSolution()
1466 # -------- test the solution ---------------------------
1467 error=Lsup(u-u_ex)
1468 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1469 class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1470 def test_solve(self):
1471 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1472 x=Solution(domain).getX()
1473 # --- set exact solution ----
1474 u_ex=Vector(0,Solution(domain))
1475 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1476 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1477 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1478 # --- set exact gradient -----------
1479 g_ex=Data(0.,(3,3),Solution(domain))
1480 g_ex[0,0]=2.
1481 g_ex[0,1]=3.
1482 g_ex[0,2]=4.
1483 g_ex[1,0]=4.
1484 g_ex[1,1]=1.
1485 g_ex[1,2]=-2.
1486 g_ex[2,0]=8.
1487 g_ex[2,1]=4.
1488 g_ex[2,2]=5.
1489 # -------- test gradient --------------------------------
1490 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1491 # -------- set-up PDE -----------------------------------
1492 pde=LinearPDE(domain,numEquations=3)
1493 mask=whereZero(x[0])
1494 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1495 A=Tensor4(0,Function(domain))
1496 A[0,:,0,:]=kronecker(3)
1497 A[1,:,1,:]=kronecker(3)
1498 A[2,:,2,:]=kronecker(3)
1499 Y=Vector(0.,Function(domain))
1500 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1501 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1502 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1503 pde.setValue(A=A,
1504 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1505 Y=Y,
1506 y=matrixmult(g_ex,domain.getNormal()))
1507 # -------- get the solution ---------------------------
1508 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1509 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1510 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1511 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1512 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1513 u=pde.getSolution()
1514 # -------- test the solution ---------------------------
1515 error=Lsup(u-u_ex)
1516 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1517 class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1518 def test_solve(self):
1519 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1520 x=Solution(domain).getX()
1521 # --- set exact solution ----
1522 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
1523 # --- set exact gradient -----------
1524 g_ex=Data(0.,(3,),Solution(domain))
1525 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1526 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1527 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1528 # -------- test gradient --------------------------------
1529 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1530 # -------- set-up PDE -----------------------------------
1531 pde=LinearPDE(domain,numEquations=1)
1532 mask=whereZero(x[0])
1533 pde.setValue(r=u_ex,q=mask)
1534 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1535 # -------- get the solution ---------------------------
1536 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1537 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1538 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1539 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1540 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1541 u=pde.getSolution()
1542 # -------- test the solution ---------------------------
1543 error=Lsup(u-u_ex)
1544 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1545 class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_GaussSeidel(unittest.TestCase):
1546 def test_solve(self):
1547 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1548 x=Solution(domain).getX()
1549 # --- set exact solution ----
1550 u_ex=Vector(0,Solution(domain))
1551 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
1552 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
1553 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
1554 # --- set exact gradient -----------
1555 g_ex=Data(0.,(3,3),Solution(domain))
1556 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1557 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1558 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1559 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1560 g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
1561 g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
1562 g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
1563 g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
1564 g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
1565 # -------- test gradient --------------------------------
1566 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1567 # -------- set-up PDE -----------------------------------
1568 pde=LinearPDE(domain,numEquations=3)
1569 mask=whereZero(x[0])
1570 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1571 Y=Vector(0.,Function(domain))
1572 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1573 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1574 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1575 A=Tensor4(0,Function(domain))
1576 A[0,:,0,:]=kronecker(3)
1577 A[1,:,1,:]=kronecker(3)
1578 A[2,:,2,:]=kronecker(3)
1579 pde.setValue(A=A,
1580 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1581 Y=Y-numpy.array([60.,20.,22.]),
1582 y=matrixmult(g_ex,domain.getNormal()))
1583 # -------- get the solution ---------------------------
1584 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1585 pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
1586 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1587 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1588 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1589 u=pde.getSolution()
1590 # -------- test the solution ---------------------------
1591 error=Lsup(u-u_ex)
1592 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1593
1594 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1595 def test_solve(self):
1596 domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
1597 x=Solution(domain).getX()
1598 # --- set exact solution ----
1599 u_ex=Scalar(0,Solution(domain))
1600 u_ex=1.+2.*x[0]+3.*x[1]
1601 # --- set exact gradient -----------
1602 g_ex=Data(0.,(2,),Solution(domain))
1603 g_ex[0]=2.
1604 g_ex[1]=3.
1605 # -------- test gradient --------------------------------
1606 g=grad(u_ex)
1607 self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
1608 # -------- set-up PDE -----------------------------------
1609 pde=LinearPDE(domain,numEquations=1)
1610 mask=whereZero(x[0])
1611 pde.setValue(r=u_ex,q=mask)
1612 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
1613 # -------- get the solution ---------------------------
1614 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1615 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1616 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1617 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1618 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1619 u=pde.getSolution()
1620 # -------- test the solution ---------------------------
1621 error=Lsup(u-u_ex)
1622 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1623
1624 class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1625 def test_solve(self):
1626 domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
1627 x=Solution(domain).getX()
1628 # --- set exact solution ----
1629 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1630 # --- set exact gradient -----------
1631 g_ex=Data(0.,(2,),Solution(domain))
1632 g_ex[0]=2.+8.*x[0]+5.*x[1]
1633 g_ex[1]=3.+5.*x[0]+12.*x[1]
1634 # -------- test gradient --------------------------------
1635 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1636 # -------- set-up PDE -----------------------------------
1637 pde=LinearPDE(domain,numEquations=1)
1638 mask=whereZero(x[0])
1639 pde.setValue(r=u_ex,q=mask)
1640 pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
1641 # -------- get the solution ---------------------------
1642 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1643 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1644 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1645 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1646 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1647 u=pde.getSolution()
1648 # -------- test the solution ---------------------------
1649 error=Lsup(u-u_ex)
1650 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1651
1652 class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1653 def test_solve(self):
1654 domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
1655 x=Solution(domain).getX()
1656 # --- set exact solution ----
1657 u_ex=Vector(0,Solution(domain))
1658 u_ex[0]=1.+2.*x[0]+3.*x[1]
1659 u_ex[1]=-1.+3.*x[0]+2.*x[1]
1660 # --- set exact gradient -----------
1661 g_ex=Data(0.,(2,2),Solution(domain))
1662 g_ex[0,0]=2.
1663 g_ex[0,1]=3.
1664 g_ex[1,0]=3.
1665 g_ex[1,1]=2.
1666 # -------- test gradient --------------------------------
1667 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1668 # -------- set-up PDE -----------------------------------
1669 pde=LinearPDE(domain,numEquations=2)
1670 mask=whereZero(x[0])
1671 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1672 A=Tensor4(0,Function(domain))
1673 A[0,:,0,:]=kronecker(2)
1674 A[1,:,1,:]=kronecker(2)
1675 Y=Vector(0.,Function(domain))
1676 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1677 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1678 pde.setValue(A=A,
1679 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1680 Y=Y,
1681 y=matrixmult(g_ex,domain.getNormal()))
1682 # -------- get the solution ---------------------------
1683 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1684 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1685 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1686 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1687 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1688 u=pde.getSolution()
1689 # -------- test the solution ---------------------------
1690 error=Lsup(u-u_ex)
1691 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1692
1693 class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1694 def test_solve(self):
1695 domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
1696 x=Solution(domain).getX()
1697 # --- set exact solution ----
1698 u_ex=Vector(0,Solution(domain))
1699 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
1700 u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
1701 # --- set exact gradient -----------
1702 g_ex=Data(0.,(2,2),Solution(domain))
1703 g_ex[0,0]=2.+8.*x[0]+5.*x[1]
1704 g_ex[0,1]=3.+5.*x[0]+12.*x[1]
1705 g_ex[1,0]=4.+2.*x[0]+6.*x[1]
1706 g_ex[1,1]=2.+6.*x[0]+8.*x[1]
1707 # -------- test gradient --------------------------------
1708 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1709 # -------- set-up PDE -----------------------------------
1710 pde=LinearPDE(domain,numEquations=2)
1711 mask=whereZero(x[0])
1712 pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
1713 A=Tensor4(0,Function(domain))
1714 A[0,:,0,:]=kronecker(2)
1715 A[1,:,1,:]=kronecker(2)
1716 Y=Vector(0.,Function(domain))
1717 Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
1718 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
1719 pde.setValue(A=A,
1720 D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
1721 Y=Y-[20.,10.],
1722 y=matrixmult(g_ex,domain.getNormal()))
1723 # -------- get the solution ---------------------------
1724 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1725 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1726 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1727 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1728 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1729 u=pde.getSolution()
1730 # -------- test the solution ---------------------------
1731 error=Lsup(u-u_ex)
1732 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1733
1734 class SimpleSolve_Brick_Order1_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1735 def test_solve(self):
1736 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1737 x=Solution(domain).getX()
1738 u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
1739 # --- set exact gradient -----------
1740 g_ex=Data(0.,(3,),Solution(domain))
1741 g_ex[0]=2.
1742 g_ex[1]=3.
1743 g_ex[2]=4.
1744 # -------- test gradient --------------------------------
1745 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1746 # -------- set-up PDE -----------------------------------
1747 pde=LinearPDE(domain,numEquations=1)
1748 mask=whereZero(x[0])
1749 pde.setValue(r=u_ex,q=mask)
1750 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
1751 # -------- get the solution ---------------------------
1752 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1753 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1754 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1755 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1756 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1757 u=pde.getSolution()
1758 # -------- test the solution ---------------------------
1759 error=Lsup(u-u_ex)
1760 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1761
1762 class SimpleSolve_Brick_Order1_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1763 def test_solve(self):
1764 domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
1765 x=Solution(domain).getX()
1766 # --- set exact solution ----
1767 u_ex=Vector(0,Solution(domain))
1768 u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
1769 u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
1770 u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
1771 # --- set exact gradient -----------
1772 g_ex=Data(0.,(3,3),Solution(domain))
1773 g_ex[0,0]=2.
1774 g_ex[0,1]=3.
1775 g_ex[0,2]=4.
1776 g_ex[1,0]=4.
1777 g_ex[1,1]=1.
1778 g_ex[1,2]=-2.
1779 g_ex[2,0]=8.
1780 g_ex[2,1]=4.
1781 g_ex[2,2]=5.
1782 # -------- test gradient --------------------------------
1783 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1784 # -------- set-up PDE -----------------------------------
1785 pde=LinearPDE(domain,numEquations=3)
1786 mask=whereZero(x[0])
1787 pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
1788 A=Tensor4(0,Function(domain))
1789 A[0,:,0,:]=kronecker(3)
1790 A[1,:,1,:]=kronecker(3)
1791 A[2,:,2,:]=kronecker(3)
1792 Y=Vector(0.,Function(domain))
1793 Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
1794 Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
1795 Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
1796 pde.setValue(A=A,
1797 D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
1798 Y=Y,
1799 y=matrixmult(g_ex,domain.getNormal()))
1800 # -------- get the solution ---------------------------
1801 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1802 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1803 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1804 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1805 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1806 u=pde.getSolution()
1807 # -------- test the solution ---------------------------
1808 error=Lsup(u-u_ex)
1809 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1810
1811 class SimpleSolve_Brick_Order2_SinglePDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1812 def test_solve(self):
1813 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1814 x=Solution(domain).getX()
1815 # --- set exact solution ----
1816 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
1817 # --- set exact gradient -----------
1818 g_ex=Data(0.,(3,),Solution(domain))
1819 g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1820 g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1821 g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1822 # -------- test gradient --------------------------------
1823 self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
1824 # -------- set-up PDE -----------------------------------
1825 pde=LinearPDE(domain,numEquations=1)
1826 mask=whereZero(x[0])
1827 pde.setValue(r=u_ex,q=mask)
1828 pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
1829 # -------- get the solution ---------------------------
1830 pde.getSolverOptions().setTolerance(SOLVER_TOL)
1831 pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
1832 pde.getSolverOptions().setPreconditioner(SolverOptions.GAUSS_SEIDEL)
1833 pde.getSolverOptions().setPackage(SolverOptions.PASO)
1834 pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
1835 u=pde.getSolution()
1836 # -------- test the solution ---------------------------
1837 error=Lsup(u-u_ex)
1838 self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
1839
1840 class SimpleSolve_Brick_Order2_SystemPDE_Paso_TFQMR_GaussSeidel(unittest.TestCase):
1841 def test_solve(self):
1842 domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
1843 x=Solution(domain).getX()
1844 # --- set exact solution ----
1845 u_ex=Vector(0,Solution(domain))
1846 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
1847 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
1848 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
1849 # --- set exact gradient -----------
1850 g_ex=Data(0.,(3,3),Solution(domain))
1851 g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
1852 g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
1853 g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
1854 g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
1855 g_ex[1,1]=1+3.*x[