/[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 1425 - (hide annotations)
Wed Feb 27 06:12:04 2008 UTC (12 years, 8 months ago) by artak
File MIME type: text/x-python
File size: 17100 byte(s)
Added new test for BiCGStab in run_simplesolve.py
algorithmic change in BICGSTAB: from r=r-wt to r=s-wt

1 ksteube 1310 #
2     # $Id$
3     #
4     #######################################################
5     #
6     # Copyright 2003-2007 by ACceSS MNRF
7     # Copyright 2007 by University of Queensland
8     #
9     # http://esscc.uq.edu.au
10     # Primary Business: Queensland, Australia
11     # Licensed under the Open Software License version 3.0
12     # http://www.opensource.org/licenses/osl-3.0.php
13     #
14     #######################################################
15     #
16    
17     """
18     Test suite for the linearPDE and pdetools test on finley
19    
20     @remark:
21    
22     @var __author__: name of author
23     @var __licence__: licence agreement
24     @var __url__: url entry point on documentation
25     @var __version__: version
26     @var __date__: date of the version
27     """
28    
29     __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
30     http://www.access.edu.au
31     Primary Business: Queensland, Australia"""
32     __license__="""Licensed under the Open Software License version 3.0
33     http://www.opensource.org/licenses/osl-3.0.php"""
34     __author__="Lutz Gross, l.gross@uq.edu.au"
35     __url__="http://www.iservo.edu.au/esys/escript"
36     __version__="$Revision: 859 $"
37     __date__="$Date: 2006-09-26 12:19:18 +1000 (Tue, 26 Sep 2006) $"
38    
39     import unittest
40    
41     from esys.escript import *
42     from esys.finley import Rectangle,Brick
43     from esys.escript.linearPDEs import LinearPDE
44     OPTIMIZE=False
45     SOLVER_VERBOSE=False
46    
47     try:
48     FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
49     except KeyError:
50     FINLEY_TEST_DATA='.'
51    
52     FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/"
53    
54     # number of elements in the spatial directions
55     NE0=8
56     NE1=10
57     NE2=12
58    
59     NE0=12
60     NE1=12
61     NE2=8
62    
63     SOLVER_TOL=1.e-8
64     REL_TOL=1.e-6
65    
66     FAC_DIAG=1.
67     FAC_OFFDIAG=-0.4
68    
69 artak 1425 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi(unittest.TestCase):
70     def test_solve(self):
71     domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
72     x=Solution(domain).getX()
73     # --- set exact solution ----
74     u_ex=Scalar(0,Solution(domain))
75     u_ex=1.+2.*x[0]+3.*x[1]
76     # --- set exact gradient -----------
77     g_ex=Data(0.,(2,),Solution(domain))
78     g_ex[0]=2.
79     g_ex[1]=3.
80     # -------- test gradient --------------------------------
81     g=grad(u_ex)
82     self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
83     # -------- set-up PDE -----------------------------------
84     pde=LinearPDE(domain,numEquations=1)
85     mask=whereZero(x[0])
86     pde.setValue(r=u_ex,q=mask)
87     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
88     # -------- get the solution ---------------------------
89     pde.setTolerance(SOLVER_TOL)
90     pde.setSolverMethod(pde.BICGSTAB,pde.JACOBI)
91     pde.setSolverPackage(pde.PASO)
92     u=pde.getSolution(verbose=SOLVER_VERBOSE)
93     # -------- test the solution ---------------------------
94     error=Lsup(u-u_ex)/Lsup(u_ex)
95     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
96 ksteube 1310 class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
97     def test_solve(self):
98     domain=Rectangle(NE0,NE1,1, optimize=OPTIMIZE)
99     x=Solution(domain).getX()
100     # --- set exact solution ----
101     u_ex=Scalar(0,Solution(domain))
102     u_ex=1.+2.*x[0]+3.*x[1]
103     # --- set exact gradient -----------
104     g_ex=Data(0.,(2,),Solution(domain))
105     g_ex[0]=2.
106     g_ex[1]=3.
107     # -------- test gradient --------------------------------
108     g=grad(u_ex)
109     self.failUnless(Lsup(g_ex-g)<REL_TOL*Lsup(g_ex))
110     # -------- set-up PDE -----------------------------------
111     pde=LinearPDE(domain,numEquations=1)
112     mask=whereZero(x[0])
113     pde.setValue(r=u_ex,q=mask)
114     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()))
115     # -------- get the solution ---------------------------
116     pde.setTolerance(SOLVER_TOL)
117     pde.setSolverMethod(pde.PCG,pde.JACOBI)
118     pde.setSolverPackage(pde.PASO)
119     u=pde.getSolution(verbose=SOLVER_VERBOSE)
120     # -------- test the solution ---------------------------
121     error=Lsup(u-u_ex)/Lsup(u_ex)
122     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
123     class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
124     def test_solve(self):
125     domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
126     x=Solution(domain).getX()
127     # --- set exact solution ----
128     u_ex=Vector(0,Solution(domain))
129     u_ex[0]=1.+2.*x[0]+3.*x[1]
130     u_ex[1]=-1.+3.*x[0]+2.*x[1]
131     # --- set exact gradient -----------
132     g_ex=Data(0.,(2,2),Solution(domain))
133     g_ex[0,0]=2.
134     g_ex[0,1]=3.
135     g_ex[1,0]=3.
136     g_ex[1,1]=2.
137     # -------- test gradient --------------------------------
138     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
139     # -------- set-up PDE -----------------------------------
140     pde=LinearPDE(domain,numEquations=2)
141     mask=whereZero(x[0])
142     pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
143     A=Tensor4(0,Function(domain))
144     A[0,:,0,:]=kronecker(2)
145     A[1,:,1,:]=kronecker(2)
146     Y=Vector(0.,Function(domain))
147     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
148     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
149     pde.setValue(A=A,
150     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
151     Y=Y,
152     y=matrixmult(g_ex,domain.getNormal()))
153     # -------- get the solution ---------------------------
154     pde.setTolerance(SOLVER_TOL)
155     pde.setSolverMethod(pde.PCG,pde.JACOBI)
156     pde.setSolverPackage(pde.PASO)
157     u=pde.getSolution(verbose=SOLVER_VERBOSE)
158     # -------- test the solution ---------------------------
159     error=Lsup(u-u_ex)/Lsup(u_ex)
160     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
161     class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
162     def test_solve(self):
163     domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
164     x=Solution(domain).getX()
165     # --- set exact solution ----
166     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
167     # --- set exact gradient -----------
168     g_ex=Data(0.,(2,),Solution(domain))
169     g_ex[0]=2.+8.*x[0]+5.*x[1]
170     g_ex[1]=3.+5.*x[0]+12.*x[1]
171     # -------- test gradient --------------------------------
172     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
173     # -------- set-up PDE -----------------------------------
174     pde=LinearPDE(domain,numEquations=1)
175     mask=whereZero(x[0])
176     pde.setValue(r=u_ex,q=mask)
177     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
178     # -------- get the solution ---------------------------
179     pde.setTolerance(SOLVER_TOL)
180     pde.setSolverMethod(pde.PCG,pde.JACOBI)
181     pde.setSolverPackage(pde.PASO)
182     u=pde.getSolution(verbose=SOLVER_VERBOSE)
183     # -------- test the solution ---------------------------
184     error=Lsup(u-u_ex)/Lsup(u_ex)
185     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
186     class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
187     def test_solve(self):
188     domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
189     x=Solution(domain).getX()
190     # --- set exact solution ----
191     u_ex=Vector(0,Solution(domain))
192     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
193     u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
194     # --- set exact gradient -----------
195     g_ex=Data(0.,(2,2),Solution(domain))
196     g_ex[0,0]=2.+8.*x[0]+5.*x[1]
197     g_ex[0,1]=3.+5.*x[0]+12.*x[1]
198     g_ex[1,0]=4.+2.*x[0]+6.*x[1]
199     g_ex[1,1]=2.+6.*x[0]+8.*x[1]
200     # -------- test gradient --------------------------------
201     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
202     # -------- set-up PDE -----------------------------------
203     pde=LinearPDE(domain,numEquations=2)
204     mask=whereZero(x[0])
205     pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
206     A=Tensor4(0,Function(domain))
207     A[0,:,0,:]=kronecker(2)
208     A[1,:,1,:]=kronecker(2)
209     Y=Vector(0.,Function(domain))
210     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
211     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
212     pde.setValue(A=A,
213     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
214     Y=Y-[20.,10.],
215     y=matrixmult(g_ex,domain.getNormal()))
216     # -------- get the solution ---------------------------
217     pde.setTolerance(SOLVER_TOL)
218     pde.setSolverMethod(pde.PCG,pde.JACOBI)
219     pde.setSolverPackage(pde.PASO)
220     u=pde.getSolution(verbose=SOLVER_VERBOSE)
221     # -------- test the solution ---------------------------
222     error=Lsup(u-u_ex)/Lsup(u_ex)
223     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
224     class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
225     def test_solve(self):
226     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
227     x=Solution(domain).getX()
228     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
229     # --- set exact gradient -----------
230     g_ex=Data(0.,(3,),Solution(domain))
231     g_ex[0]=2.
232     g_ex[1]=3.
233     g_ex[2]=4.
234     # -------- test gradient --------------------------------
235     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
236     # -------- set-up PDE -----------------------------------
237     pde=LinearPDE(domain,numEquations=1)
238     mask=whereZero(x[0])
239     pde.setValue(r=u_ex,q=mask)
240     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
241     # -------- get the solution ---------------------------
242     pde.setTolerance(SOLVER_TOL)
243     pde.setSolverMethod(pde.PCG,pde.JACOBI)
244     pde.setSolverPackage(pde.PASO)
245     u=pde.getSolution(verbose=SOLVER_VERBOSE)
246     # -------- test the solution ---------------------------
247     error=Lsup(u-u_ex)/Lsup(u_ex)
248     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
249     class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
250     def test_solve(self):
251     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
252     x=Solution(domain).getX()
253     # --- set exact solution ----
254     u_ex=Vector(0,Solution(domain))
255     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
256     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
257     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
258     # --- set exact gradient -----------
259     g_ex=Data(0.,(3,3),Solution(domain))
260     g_ex[0,0]=2.
261     g_ex[0,1]=3.
262     g_ex[0,2]=4.
263     g_ex[1,0]=4.
264     g_ex[1,1]=1.
265     g_ex[1,2]=-2.
266     g_ex[2,0]=8.
267     g_ex[2,1]=4.
268     g_ex[2,2]=5.
269     # -------- test gradient --------------------------------
270     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
271     # -------- set-up PDE -----------------------------------
272     pde=LinearPDE(domain,numEquations=3)
273     mask=whereZero(x[0])
274     pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
275     A=Tensor4(0,Function(domain))
276     A[0,:,0,:]=kronecker(3)
277     A[1,:,1,:]=kronecker(3)
278     A[2,:,2,:]=kronecker(3)
279     Y=Vector(0.,Function(domain))
280     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
281     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
282     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
283     pde.setValue(A=A,
284     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
285     Y=Y,
286     y=matrixmult(g_ex,domain.getNormal()))
287     # -------- get the solution ---------------------------
288     pde.setTolerance(SOLVER_TOL)
289     pde.setSolverMethod(pde.PCG,pde.JACOBI)
290     pde.setSolverPackage(pde.PASO)
291     u=pde.getSolution(verbose=SOLVER_VERBOSE)
292     # -------- test the solution ---------------------------
293     error=Lsup(u-u_ex)/Lsup(u_ex)
294     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
295     class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
296     def test_solve(self):
297     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
298     x=Solution(domain).getX()
299     # --- set exact solution ----
300     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
301     # --- set exact gradient -----------
302     g_ex=Data(0.,(3,),Solution(domain))
303     g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
304     g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
305     g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
306     # -------- test gradient --------------------------------
307     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
308     # -------- set-up PDE -----------------------------------
309     pde=LinearPDE(domain,numEquations=1)
310     mask=whereZero(x[0])
311     pde.setValue(r=u_ex,q=mask)
312     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
313     # -------- get the solution ---------------------------
314     pde.setTolerance(SOLVER_TOL)
315     pde.setSolverMethod(pde.PCG,pde.JACOBI)
316     pde.setSolverPackage(pde.PASO)
317     u=pde.getSolution(verbose=SOLVER_VERBOSE)
318     # -------- test the solution ---------------------------
319     error=Lsup(u-u_ex)/Lsup(u_ex)
320     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
321     class SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
322     def test_solve(self):
323     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
324     x=Solution(domain).getX()
325     # --- set exact solution ----
326     u_ex=Vector(0,Solution(domain))
327     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
328     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
329     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
330     # --- set exact gradient -----------
331     g_ex=Data(0.,(3,3),Solution(domain))
332     g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
333     g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
334     g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
335     g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
336     g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
337     g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
338     g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
339     g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
340     g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
341     # -------- test gradient --------------------------------
342     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
343     # -------- set-up PDE -----------------------------------
344     pde=LinearPDE(domain,numEquations=3)
345     mask=whereZero(x[0])
346     pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
347     Y=Vector(0.,Function(domain))
348     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
349     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
350     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
351     A=Tensor4(0,Function(domain))
352     A[0,:,0,:]=kronecker(3)
353     A[1,:,1,:]=kronecker(3)
354     A[2,:,2,:]=kronecker(3)
355     pde.setValue(A=A,
356     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
357     Y=Y-numarray.array([60.,20.,22.]),
358     y=matrixmult(g_ex,domain.getNormal()))
359     # -------- get the solution ---------------------------
360     pde.setTolerance(SOLVER_TOL)
361     pde.setSolverMethod(pde.PCG,pde.JACOBI)
362     pde.setSolverPackage(pde.PASO)
363     u=pde.getSolution(verbose=SOLVER_VERBOSE)
364     # -------- test the solution ---------------------------
365     error=Lsup(u-u_ex)/Lsup(u_ex)
366     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
367    
368     if __name__ == '__main__':
369     suite = unittest.TestSuite()
370 artak 1425 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_BICGSTAB_Jacobi))
371 ksteube 1310 suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
372     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
373     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
374     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
375     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
376     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
377     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
378     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
379     s=unittest.TextTestRunner(verbosity=2).run(suite)
380 gross 1376 if not s.wasSuccessful(): sys.exit(1)

  ViewVC Help
Powered by ViewVC 1.1.26