/[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 4984 - (show annotations)
Mon Jun 2 02:50:34 2014 UTC (5 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 121403 byte(s)
revamping testrunners, now uses automated discovery and allows running specific tests without modifying files (see escriptcore/py_src/testing.py for more info/examples)

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