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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4154 - (hide annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 7 months ago) by jfenwick
File MIME type: text/x-python
File size: 21028 byte(s)
Round 1 of copyright fixes
1 caltinay 3776
2 jfenwick 3981 ##############################################################################
3 caltinay 3776 #
4 jfenwick 4154 # Copyright (c) 2003-2013 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3776 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3776
16 jfenwick 4154 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
17 jfenwick 3981 http://www.uq.edu.au
18 caltinay 3776 Primary Business: Queensland, Australia"""
19     __license__="""Licensed under the Open Software License version 3.0
20     http://www.opensource.org/licenses/osl-3.0.php"""
21     __url__="https://launchpad.net/escript-finley"
22    
23     """
24     Test suite for the linearPDE and pdetools test on ripley
25    
26     :remark:
27    
28     :var __author__: name of author
29     :var __licence__: licence agreement
30     :var __url__: url entry point on documentation
31     :var __version__: version
32     :var __date__: date of the version
33     """
34    
35     __author__="Lutz Gross, l.gross@uq.edu.au"
36    
37     import unittest, sys
38     from esys.escript import *
39     from esys.ripley import Rectangle,Brick
40     from esys.escript.linearPDEs import LinearPDE, SolverOptions
41     import numpy
42     SOLVER_VERBOSE=False
43    
44     try:
45     RIPLEY_TEST_DATA=os.environ['RIPLEY_TEST_DATA']
46     except KeyError:
47     RIPLEY_TEST_DATA='.'
48    
49     SOLVER_TOL=1.e-8
50     REL_TOL=1.e-6
51    
52     FAC_DIAG=1.
53     FAC_OFFDIAG=-0.4
54    
55 caltinay 3803 # number of elements in the spatial directions
56     NE0=12
57     NE1=12
58     NE2=8
59     mpiSize=getMPISizeWorld()
60     for x in [int(sqrt(mpiSize)),2,3,5,7,1]:
61     NX=x
62 jfenwick 3892 NY=mpiSize//x
63 caltinay 3803 if NX*NY == mpiSize:
64     break
65 caltinay 3776
66 caltinay 3803 for x in [(int(mpiSize**(1/3.)),int(mpiSize**(1/3.))),(2,3),(2,2),(1,2),(1,1)]:
67     NXb=x[0]
68     NYb=x[1]
69 jfenwick 3892 NZb=mpiSize//(x[0]*x[1])
70 caltinay 3803 if NXb*NYb*NZb == mpiSize:
71     break
72    
73 caltinay 3776 class SimpleSolve_Rectangle_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
74     def test_solve(self):
75     # Tell about how many MPI CPUs and OpenMP threads
76 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
77 caltinay 3776 x=Solution(domain).getX()
78     # --- set exact solution ----
79     u_ex=Scalar(0,Solution(domain))
80     u_ex=1.+2.*x[0]+3.*x[1]
81     # --- set exact gradient -----------
82     g_ex=Data(0.,(2,),Solution(domain))
83     g_ex[0]=2.
84     g_ex[1]=3.
85     # -------- test gradient --------------------------------
86     g=grad(u_ex)
87     self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
88     # -------- set-up PDE -----------------------------------
89     pde=LinearPDE(domain,numEquations=1)
90     mask=whereZero(x[0])
91     pde.setValue(r=u_ex,q=mask)
92     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
93     # -------- get the solution ---------------------------
94     pde.getSolverOptions().setTolerance(SOLVER_TOL)
95     pde.getSolverOptions().setSolverMethod(SolverOptions.BICGSTAB)
96     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
97     pde.getSolverOptions().setPackage(SolverOptions.PASO)
98     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
99     u=pde.getSolution()
100     # -------- test the solution ---------------------------
101     error=Lsup(u-u_ex)
102     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
103     class SimpleSolve_Rectangle_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
104     def test_solve(self):
105 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
106 caltinay 3776 x=Solution(domain).getX()
107     # --- set exact solution ----
108     u_ex=Scalar(0,Solution(domain))
109     u_ex=1.+2.*x[0]+3.*x[1]
110     # --- set exact gradient -----------
111     g_ex=Data(0.,(2,),Solution(domain))
112     g_ex[0]=2.
113     g_ex[1]=3.
114     # -------- test gradient --------------------------------
115     g=grad(u_ex)
116     self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
117     # -------- set-up PDE -----------------------------------
118     pde=LinearPDE(domain,numEquations=1)
119     mask=whereZero(x[0])
120     pde.setValue(r=u_ex,q=mask)
121     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
122     # -------- get the solution ---------------------------
123     pde.getSolverOptions().setTolerance(SOLVER_TOL)
124     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
125     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
126     pde.getSolverOptions().setPackage(SolverOptions.PASO)
127     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
128     u=pde.getSolution()
129     # -------- test the solution ---------------------------
130     error=Lsup(u-u_ex)
131     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
132     class SimpleSolve_Rectangle_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
133     def test_solve(self):
134 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
135 caltinay 3776 x=Solution(domain).getX()
136     # --- set exact solution ----
137     u_ex=Vector(0,Solution(domain))
138     u_ex[0]=1.+2.*x[0]+3.*x[1]
139     u_ex[1]=-1.+3.*x[0]+2.*x[1]
140     # --- set exact gradient -----------
141     g_ex=Data(0.,(2,2),Solution(domain))
142     g_ex[0,0]=2.
143     g_ex[0,1]=3.
144     g_ex[1,0]=3.
145     g_ex[1,1]=2.
146     # -------- test gradient --------------------------------
147     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
148     # -------- set-up PDE -----------------------------------
149     pde=LinearPDE(domain,numEquations=2)
150     mask=whereZero(x[0])
151     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
152     A=Tensor4(0,Function(domain))
153     A[0,:,0,:]=kronecker(2)
154     A[1,:,1,:]=kronecker(2)
155     Y=Vector(0.,Function(domain))
156     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
157     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
158     pde.setValue(A=A,
159     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
160     Y=Y,
161     y=matrixmult(g_ex,domain.getNormal()))
162     # -------- get the solution ---------------------------
163     pde.getSolverOptions().setTolerance(SOLVER_TOL)
164     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
165     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
166     pde.getSolverOptions().setPackage(SolverOptions.PASO)
167     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
168     u=pde.getSolution()
169     # -------- test the solution ---------------------------
170     error=Lsup(u-u_ex)
171     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
172     class SimpleSolve_Brick_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
173     def test_solve(self):
174 caltinay 3803 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
175 caltinay 3776 x=Solution(domain).getX()
176     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
177     # --- set exact gradient -----------
178     g_ex=Data(0.,(3,),Solution(domain))
179     g_ex[0]=2.
180     g_ex[1]=3.
181     g_ex[2]=4.
182     # -------- test gradient --------------------------------
183     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
184     # -------- set-up PDE -----------------------------------
185     pde=LinearPDE(domain,numEquations=1)
186     mask=whereZero(x[0])
187     pde.setValue(r=u_ex,q=mask)
188     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
189     # -------- get the solution ---------------------------
190     pde.getSolverOptions().setTolerance(SOLVER_TOL)
191     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
192     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
193     pde.getSolverOptions().setPackage(SolverOptions.PASO)
194     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
195     u=pde.getSolution()
196     # -------- test the solution ---------------------------
197     error=Lsup(u-u_ex)
198     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
199     class SimpleSolve_Brick_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
200     def test_solve(self):
201 caltinay 3803 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
202 caltinay 3776 x=Solution(domain).getX()
203     # --- set exact solution ----
204     u_ex=Vector(0,Solution(domain))
205     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
206     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
207     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
208     # --- set exact gradient -----------
209     g_ex=Data(0.,(3,3),Solution(domain))
210     g_ex[0,0]=2.
211     g_ex[0,1]=3.
212     g_ex[0,2]=4.
213     g_ex[1,0]=4.
214     g_ex[1,1]=1.
215     g_ex[1,2]=-2.
216     g_ex[2,0]=8.
217     g_ex[2,1]=4.
218     g_ex[2,2]=5.
219     # -------- test gradient --------------------------------
220     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
221     # -------- set-up PDE -----------------------------------
222     pde=LinearPDE(domain,numEquations=3)
223     mask=whereZero(x[0])
224     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
225     A=Tensor4(0,Function(domain))
226     A[0,:,0,:]=kronecker(3)
227     A[1,:,1,:]=kronecker(3)
228     A[2,:,2,:]=kronecker(3)
229     Y=Vector(0.,Function(domain))
230     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
231     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
232     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
233     pde.setValue(A=A,
234     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
235     Y=Y,
236     y=matrixmult(g_ex,domain.getNormal()))
237     # -------- get the solution ---------------------------
238     pde.getSolverOptions().setTolerance(SOLVER_TOL)
239     pde.getSolverOptions().setSolverMethod(SolverOptions.PCG)
240     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
241     pde.getSolverOptions().setPackage(SolverOptions.PASO)
242     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
243     u=pde.getSolution()
244     # -------- test the solution ---------------------------
245     error=Lsup(u-u_ex)
246     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
247    
248     class SimpleSolve_Rectangle_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
249     def test_solve(self):
250 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
251 caltinay 3776 x=Solution(domain).getX()
252     # --- set exact solution ----
253     u_ex=Scalar(0,Solution(domain))
254     u_ex=1.+2.*x[0]+3.*x[1]
255     # --- set exact gradient -----------
256     g_ex=Data(0.,(2,),Solution(domain))
257     g_ex[0]=2.
258     g_ex[1]=3.
259     # -------- test gradient --------------------------------
260     g=grad(u_ex)
261     self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
262     # -------- set-up PDE -----------------------------------
263     pde=LinearPDE(domain,numEquations=1)
264     mask=whereZero(x[0])
265     pde.setValue(r=u_ex,q=mask)
266     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
267     # -------- get the solution ---------------------------
268     pde.getSolverOptions().setTolerance(SOLVER_TOL)
269     pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
270     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
271     pde.getSolverOptions().setPackage(SolverOptions.PASO)
272     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
273     u=pde.getSolution()
274     # -------- test the solution ---------------------------
275     error=Lsup(u-u_ex)
276     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
277    
278     class SimpleSolve_Brick_SinglePDE_Paso_TFQMR_Jacobi(unittest.TestCase):
279     def test_solve(self):
280 caltinay 3803 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
281 caltinay 3776 x=Solution(domain).getX()
282     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
283     # --- set exact gradient -----------
284     g_ex=Data(0.,(3,),Solution(domain))
285     g_ex[0]=2.
286     g_ex[1]=3.
287     g_ex[2]=4.
288     # -------- test gradient --------------------------------
289     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
290     # -------- set-up PDE -----------------------------------
291     pde=LinearPDE(domain,numEquations=1)
292     mask=whereZero(x[0])
293     pde.setValue(r=u_ex,q=mask)
294     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
295     # -------- get the solution ---------------------------
296     pde.getSolverOptions().setTolerance(SOLVER_TOL)
297     pde.getSolverOptions().setSolverMethod(SolverOptions.TFQMR)
298     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
299     pde.getSolverOptions().setPackage(SolverOptions.PASO)
300     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
301     u=pde.getSolution()
302     # -------- test the solution ---------------------------
303     error=Lsup(u-u_ex)
304     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
305    
306     class SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
307     def test_solve(self):
308 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
309 caltinay 3776 x=Solution(domain).getX()
310     # --- set exact solution ----
311     u_ex=Scalar(0,Solution(domain))
312     u_ex=1.+2.*x[0]+3.*x[1]
313     # --- set exact gradient -----------
314     g_ex=Data(0.,(2,),Solution(domain))
315     g_ex[0]=2.
316     g_ex[1]=3.
317     # -------- test gradient --------------------------------
318     g=grad(u_ex)
319     self.assertTrue(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
320     # -------- set-up PDE -----------------------------------
321     pde=LinearPDE(domain,numEquations=1)
322     mask=whereZero(x[0])
323     pde.setValue(r=u_ex,q=mask)
324     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
325     # -------- get the solution ---------------------------
326     pde.getSolverOptions().setTolerance(SOLVER_TOL)
327     pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
328     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
329     pde.getSolverOptions().setPackage(SolverOptions.PASO)
330    
331     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
332     u=pde.getSolution()
333     # -------- test the solution ---------------------------
334     error=Lsup(u-u_ex)
335     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
336    
337     class SimpleSolve_Rectangle_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
338     def test_solve(self):
339 caltinay 3803 domain=Rectangle(n0=NE0*NX-1, n1=NE1*NY-1, d0=NX, d1=NY)
340 caltinay 3776 x=Solution(domain).getX()
341     # --- set exact solution ----
342     u_ex=Vector(0,Solution(domain))
343     u_ex[0]=1.+2.*x[0]+3.*x[1]
344     u_ex[1]=-1.+3.*x[0]+2.*x[1]
345     # --- set exact gradient -----------
346     g_ex=Data(0.,(2,2),Solution(domain))
347     g_ex[0,0]=2.
348     g_ex[0,1]=3.
349     g_ex[1,0]=3.
350     g_ex[1,1]=2.
351     # -------- test gradient --------------------------------
352     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
353     # -------- set-up PDE -----------------------------------
354     pde=LinearPDE(domain,numEquations=2)
355     mask=whereZero(x[0])
356     pde.setValue(r=u_ex,q=mask*numpy.ones(2,))
357     A=Tensor4(0,Function(domain))
358     A[0,:,0,:]=kronecker(2)
359     A[1,:,1,:]=kronecker(2)
360     Y=Vector(0.,Function(domain))
361     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
362     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
363     pde.setValue(A=A,
364     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((2,2))*FAC_OFFDIAG,
365     Y=Y,
366     y=matrixmult(g_ex,domain.getNormal()))
367     # -------- get the solution ---------------------------
368     pde.getSolverOptions().setTolerance(SOLVER_TOL)
369     pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
370     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
371     pde.getSolverOptions().setPackage(SolverOptions.PASO)
372     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
373     u=pde.getSolution()
374     # -------- test the solution ---------------------------
375     error=Lsup(u-u_ex)
376     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
377     class SimpleSolve_Brick_SinglePDE_Paso_MINRES_Jacobi(unittest.TestCase):
378     def test_solve(self):
379 caltinay 3803 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
380 caltinay 3776 x=Solution(domain).getX()
381     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
382     # --- set exact gradient -----------
383     g_ex=Data(0.,(3,),Solution(domain))
384     g_ex[0]=2.
385     g_ex[1]=3.
386     g_ex[2]=4.
387     # -------- test gradient --------------------------------
388     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
389     # -------- set-up PDE -----------------------------------
390     pde=LinearPDE(domain,numEquations=1)
391     mask=whereZero(x[0])
392     pde.setValue(r=u_ex,q=mask)
393     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
394     # -------- get the solution ---------------------------
395     pde.getSolverOptions().setTolerance(SOLVER_TOL)
396     pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
397     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
398     pde.getSolverOptions().setPackage(SolverOptions.PASO)
399     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
400     u=pde.getSolution()
401     # -------- test the solution ---------------------------
402     error=Lsup(u-u_ex)
403     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
404    
405     class SimpleSolve_Brick_SystemPDE_Paso_MINRES_Jacobi(unittest.TestCase):
406     def test_solve(self):
407 caltinay 3803 domain=Brick(n0=NE0*NXb-1, n1=NE1*NYb-1, n2=NE2*NZb-1, d0=NXb, d1=NYb, d2=NZb)
408 caltinay 3776 x=Solution(domain).getX()
409     # --- set exact solution ----
410     u_ex=Vector(0,Solution(domain))
411     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
412     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
413     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
414     # --- set exact gradient -----------
415     g_ex=Data(0.,(3,3),Solution(domain))
416     g_ex[0,0]=2.
417     g_ex[0,1]=3.
418     g_ex[0,2]=4.
419     g_ex[1,0]=4.
420     g_ex[1,1]=1.
421     g_ex[1,2]=-2.
422     g_ex[2,0]=8.
423     g_ex[2,1]=4.
424     g_ex[2,2]=5.
425     # -------- test gradient --------------------------------
426     self.assertTrue(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
427     # -------- set-up PDE -----------------------------------
428     pde=LinearPDE(domain,numEquations=3)
429     mask=whereZero(x[0])
430     pde.setValue(r=u_ex,q=mask*numpy.ones(3,))
431     A=Tensor4(0,Function(domain))
432     A[0,:,0,:]=kronecker(3)
433     A[1,:,1,:]=kronecker(3)
434     A[2,:,2,:]=kronecker(3)
435     Y=Vector(0.,Function(domain))
436     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
437     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
438     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
439     pde.setValue(A=A,
440     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numpy.ones((3,3))*FAC_OFFDIAG,
441     Y=Y,
442     y=matrixmult(g_ex,domain.getNormal()))
443     # -------- get the solution ---------------------------
444     pde.getSolverOptions().setTolerance(SOLVER_TOL)
445     pde.getSolverOptions().setSolverMethod(SolverOptions.MINRES)
446     pde.getSolverOptions().setPreconditioner(SolverOptions.JACOBI)
447     pde.getSolverOptions().setPackage(SolverOptions.PASO)
448     pde.getSolverOptions().setVerbosity(SOLVER_VERBOSE)
449     u=pde.getSolution()
450     # -------- test the solution ---------------------------
451     error=Lsup(u-u_ex)
452     self.assertTrue(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
453    
454    
455     if __name__ == '__main__':
456     suite = unittest.TestSuite()
457     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_BICGSTAB_Jacobi))
458     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_PCG_Jacobi))
459     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SystemPDE_Paso_PCG_Jacobi))
460     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_PCG_Jacobi))
461     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SystemPDE_Paso_PCG_Jacobi))
462     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_TFQMR_Jacobi))
463     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_TFQMR_Jacobi))
464     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi))
465     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_SinglePDE_Paso_MINRES_Jacobi))
466     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SinglePDE_Paso_MINRES_Jacobi))
467     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_SystemPDE_Paso_MINRES_Jacobi))
468    
469     s=unittest.TextTestRunner(verbosity=2).run(suite)
470     if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26