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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26