/[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 1310 - (hide annotations)
Mon Sep 24 04:00:47 2007 UTC (13 years ago) by ksteube
File MIME type: text/x-python
File size: 15697 byte(s)
A few additions to make trunk match trunk-mpi-branch

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    
40     import unittest
41    
42     from esys.escript import *
43     from esys.finley import Rectangle,Brick
44     from esys.escript.linearPDEs import LinearPDE
45     OPTIMIZE=False
46     SOLVER_VERBOSE=False
47    
48     try:
49     FINLEY_TEST_DATA=os.environ['FINLEY_TEST_DATA']
50     except KeyError:
51     FINLEY_TEST_DATA='.'
52    
53     FINLEY_TEST_MESH_PATH=FINLEY_TEST_DATA+"/data_meshes/"
54    
55     # number of elements in the spatial directions
56     NE0=8
57     NE1=10
58     NE2=12
59    
60     NE0=12
61     NE1=12
62     NE2=8
63    
64     SOLVER_TOL=1.e-8
65     REL_TOL=1.e-6
66    
67     FAC_DIAG=1.
68     FAC_OFFDIAG=-0.4
69    
70     class SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
71     def test_solve(self):
72     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     pde.setTolerance(SOLVER_TOL)
91     pde.setSolverMethod(pde.PCG,pde.JACOBI)
92     pde.setSolverPackage(pde.PASO)
93     u=pde.getSolution(verbose=SOLVER_VERBOSE)
94     # -------- test the solution ---------------------------
95     error=Lsup(u-u_ex)/Lsup(u_ex)
96     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
97     class SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
98     def test_solve(self):
99     domain=Rectangle(NE0,NE1,1,optimize=OPTIMIZE)
100     x=Solution(domain).getX()
101     # --- set exact solution ----
102     u_ex=Vector(0,Solution(domain))
103     u_ex[0]=1.+2.*x[0]+3.*x[1]
104     u_ex[1]=-1.+3.*x[0]+2.*x[1]
105     # --- set exact gradient -----------
106     g_ex=Data(0.,(2,2),Solution(domain))
107     g_ex[0,0]=2.
108     g_ex[0,1]=3.
109     g_ex[1,0]=3.
110     g_ex[1,1]=2.
111     # -------- test gradient --------------------------------
112     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
113     # -------- set-up PDE -----------------------------------
114     pde=LinearPDE(domain,numEquations=2)
115     mask=whereZero(x[0])
116     pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
117     A=Tensor4(0,Function(domain))
118     A[0,:,0,:]=kronecker(2)
119     A[1,:,1,:]=kronecker(2)
120     Y=Vector(0.,Function(domain))
121     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
122     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
123     pde.setValue(A=A,
124     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
125     Y=Y,
126     y=matrixmult(g_ex,domain.getNormal()))
127     # -------- get the solution ---------------------------
128     pde.setTolerance(SOLVER_TOL)
129     pde.setSolverMethod(pde.PCG,pde.JACOBI)
130     pde.setSolverPackage(pde.PASO)
131     u=pde.getSolution(verbose=SOLVER_VERBOSE)
132     # -------- test the solution ---------------------------
133     error=Lsup(u-u_ex)/Lsup(u_ex)
134     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
135     class SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
136     def test_solve(self):
137     domain=Rectangle(NE0,NE1,2,l0=1.,l1=1,optimize=OPTIMIZE)
138     x=Solution(domain).getX()
139     # --- set exact solution ----
140     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
141     # --- set exact gradient -----------
142     g_ex=Data(0.,(2,),Solution(domain))
143     g_ex[0]=2.+8.*x[0]+5.*x[1]
144     g_ex[1]=3.+5.*x[0]+12.*x[1]
145     # -------- test gradient --------------------------------
146     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
147     # -------- set-up PDE -----------------------------------
148     pde=LinearPDE(domain,numEquations=1)
149     mask=whereZero(x[0])
150     pde.setValue(r=u_ex,q=mask)
151     pde.setValue(A=kronecker(2),y=inner(g_ex,domain.getNormal()),Y=-20.)
152     # -------- get the solution ---------------------------
153     pde.setTolerance(SOLVER_TOL)
154     pde.setSolverMethod(pde.PCG,pde.JACOBI)
155     pde.setSolverPackage(pde.PASO)
156     u=pde.getSolution(verbose=SOLVER_VERBOSE)
157     # -------- test the solution ---------------------------
158     error=Lsup(u-u_ex)/Lsup(u_ex)
159     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
160     class SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
161     def test_solve(self):
162     domain=Rectangle(NE0,NE1,2,optimize=OPTIMIZE)
163     x=Solution(domain).getX()
164     # --- set exact solution ----
165     u_ex=Vector(0,Solution(domain))
166     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[0]**2+5.*x[1]*x[0]+6.*x[1]**2
167     u_ex[1]=-1.+4.*x[0]+2.*x[1]+1.*x[0]**2+6.*x[1]*x[0]+4.*x[1]**2
168     # --- set exact gradient -----------
169     g_ex=Data(0.,(2,2),Solution(domain))
170     g_ex[0,0]=2.+8.*x[0]+5.*x[1]
171     g_ex[0,1]=3.+5.*x[0]+12.*x[1]
172     g_ex[1,0]=4.+2.*x[0]+6.*x[1]
173     g_ex[1,1]=2.+6.*x[0]+8.*x[1]
174     # -------- test gradient --------------------------------
175     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
176     # -------- set-up PDE -----------------------------------
177     pde=LinearPDE(domain,numEquations=2)
178     mask=whereZero(x[0])
179     pde.setValue(r=u_ex,q=mask*numarray.ones(2,))
180     A=Tensor4(0,Function(domain))
181     A[0,:,0,:]=kronecker(2)
182     A[1,:,1,:]=kronecker(2)
183     Y=Vector(0.,Function(domain))
184     Y[0]=u_ex[0]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG
185     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG
186     pde.setValue(A=A,
187     D=kronecker(2)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((2,2))*FAC_OFFDIAG,
188     Y=Y-[20.,10.],
189     y=matrixmult(g_ex,domain.getNormal()))
190     # -------- get the solution ---------------------------
191     pde.setTolerance(SOLVER_TOL)
192     pde.setSolverMethod(pde.PCG,pde.JACOBI)
193     pde.setSolverPackage(pde.PASO)
194     u=pde.getSolution(verbose=SOLVER_VERBOSE)
195     # -------- test the solution ---------------------------
196     error=Lsup(u-u_ex)/Lsup(u_ex)
197     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
198     class SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
199     def test_solve(self):
200     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
201     x=Solution(domain).getX()
202     u_ex=1.+2.*x[0]+3.*x[1]+4.*x[2]
203     # --- set exact gradient -----------
204     g_ex=Data(0.,(3,),Solution(domain))
205     g_ex[0]=2.
206     g_ex[1]=3.
207     g_ex[2]=4.
208     # -------- test gradient --------------------------------
209     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
210     # -------- set-up PDE -----------------------------------
211     pde=LinearPDE(domain,numEquations=1)
212     mask=whereZero(x[0])
213     pde.setValue(r=u_ex,q=mask)
214     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()))
215     # -------- get the solution ---------------------------
216     pde.setTolerance(SOLVER_TOL)
217     pde.setSolverMethod(pde.PCG,pde.JACOBI)
218     pde.setSolverPackage(pde.PASO)
219     u=pde.getSolution(verbose=SOLVER_VERBOSE)
220     # -------- test the solution ---------------------------
221     error=Lsup(u-u_ex)/Lsup(u_ex)
222     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
223     class SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi(unittest.TestCase):
224     def test_solve(self):
225     domain=Brick(NE0,NE1,NE2,1,optimize=OPTIMIZE)
226     x=Solution(domain).getX()
227     # --- set exact solution ----
228     u_ex=Vector(0,Solution(domain))
229     u_ex[0]=1.+2.*x[0]+3.*x[1]+4.*x[2]
230     u_ex[1]=-1.+4.*x[0]+1.*x[1]-2.*x[2]
231     u_ex[2]=5.+8.*x[0]+4.*x[1]+5.*x[2]
232     # --- set exact gradient -----------
233     g_ex=Data(0.,(3,3),Solution(domain))
234     g_ex[0,0]=2.
235     g_ex[0,1]=3.
236     g_ex[0,2]=4.
237     g_ex[1,0]=4.
238     g_ex[1,1]=1.
239     g_ex[1,2]=-2.
240     g_ex[2,0]=8.
241     g_ex[2,1]=4.
242     g_ex[2,2]=5.
243     # -------- test gradient --------------------------------
244     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
245     # -------- set-up PDE -----------------------------------
246     pde=LinearPDE(domain,numEquations=3)
247     mask=whereZero(x[0])
248     pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
249     A=Tensor4(0,Function(domain))
250     A[0,:,0,:]=kronecker(3)
251     A[1,:,1,:]=kronecker(3)
252     A[2,:,2,:]=kronecker(3)
253     Y=Vector(0.,Function(domain))
254     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
255     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
256     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
257     pde.setValue(A=A,
258     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
259     Y=Y,
260     y=matrixmult(g_ex,domain.getNormal()))
261     # -------- get the solution ---------------------------
262     pde.setTolerance(SOLVER_TOL)
263     pde.setSolverMethod(pde.PCG,pde.JACOBI)
264     pde.setSolverPackage(pde.PASO)
265     u=pde.getSolution(verbose=SOLVER_VERBOSE)
266     # -------- test the solution ---------------------------
267     error=Lsup(u-u_ex)/Lsup(u_ex)
268     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
269     class SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi(unittest.TestCase):
270     def test_solve(self):
271     domain=Brick(NE0,NE1,NE2,2,optimize=OPTIMIZE)
272     x=Solution(domain).getX()
273     # --- set exact solution ----
274     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
275     # --- set exact gradient -----------
276     g_ex=Data(0.,(3,),Solution(domain))
277     g_ex[0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
278     g_ex[1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
279     g_ex[2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
280     # -------- test gradient --------------------------------
281     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
282     # -------- set-up PDE -----------------------------------
283     pde=LinearPDE(domain,numEquations=1)
284     mask=whereZero(x[0])
285     pde.setValue(r=u_ex,q=mask)
286     pde.setValue(A=kronecker(3),y=inner(g_ex,domain.getNormal()),Y=-60.)
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_SystemPDE_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=Vector(0,Solution(domain))
301     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
302     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
303     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
304     # --- set exact gradient -----------
305     g_ex=Data(0.,(3,3),Solution(domain))
306     g_ex[0,0]=2.+6.*x[1]+8.*x[2]+18.*x[0]
307     g_ex[0,1]=3.+6.*x[0]+7.*x[2]+20.*x[1]
308     g_ex[0,2]=4.+7.*x[1]+8.*x[0]+22.*x[2]
309     g_ex[1,0]=4.+3.*x[1]-8.*x[2]-4.*x[0]
310     g_ex[1,1]=1+3.*x[0]+2.*x[2]+14.*x[1]
311     g_ex[1,2]=-6.+2.*x[1]-8.*x[0]+10.*x[2]
312     g_ex[2,0]=7.-6.*x[1]+2.*x[2]+4.*x[0]
313     g_ex[2,1]=9.-6.*x[0]+8.*x[2]+16.*x[1]
314     g_ex[2,2]=2+8.*x[1]+2.*x[0]+2.*x[2]
315     # -------- test gradient --------------------------------
316     self.failUnless(Lsup(g_ex-grad(u_ex))<REL_TOL*Lsup(g_ex))
317     # -------- set-up PDE -----------------------------------
318     pde=LinearPDE(domain,numEquations=3)
319     mask=whereZero(x[0])
320     pde.setValue(r=u_ex,q=mask*numarray.ones(3,))
321     Y=Vector(0.,Function(domain))
322     Y[0]=u_ex[0]*FAC_DIAG+u_ex[2]*FAC_OFFDIAG+u_ex[1]*FAC_OFFDIAG
323     Y[1]=u_ex[1]*FAC_DIAG+u_ex[0]*FAC_OFFDIAG+u_ex[2]*FAC_OFFDIAG
324     Y[2]=u_ex[2]*FAC_DIAG+u_ex[1]*FAC_OFFDIAG+u_ex[0]*FAC_OFFDIAG
325     A=Tensor4(0,Function(domain))
326     A[0,:,0,:]=kronecker(3)
327     A[1,:,1,:]=kronecker(3)
328     A[2,:,2,:]=kronecker(3)
329     pde.setValue(A=A,
330     D=kronecker(3)*(FAC_DIAG-FAC_OFFDIAG)+numarray.ones((3,3))*FAC_OFFDIAG,
331     Y=Y-numarray.array([60.,20.,22.]),
332     y=matrixmult(g_ex,domain.getNormal()))
333     # -------- get the solution ---------------------------
334     pde.setTolerance(SOLVER_TOL)
335     pde.setSolverMethod(pde.PCG,pde.JACOBI)
336     pde.setSolverPackage(pde.PASO)
337     u=pde.getSolution(verbose=SOLVER_VERBOSE)
338     # -------- test the solution ---------------------------
339     error=Lsup(u-u_ex)/Lsup(u_ex)
340     self.failUnless(error<REL_TOL*Lsup(u_ex), "solution error %s is too big."%error)
341    
342     if __name__ == '__main__':
343     suite = unittest.TestSuite()
344     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SinglePDE_Paso_PCG_Jacobi))
345     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order1_SystemPDE_Paso_PCG_Jacobi))
346     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SinglePDE_Paso_PCG_Jacobi))
347     suite.addTest(unittest.makeSuite(SimpleSolve_Rectangle_Order2_SystemPDE_Paso_PCG_Jacobi))
348     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SinglePDE_Paso_PCG_Jacobi))
349     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order1_SystemPDE_Paso_PCG_Jacobi))
350     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SinglePDE_Paso_PCG_Jacobi))
351     suite.addTest(unittest.makeSuite(SimpleSolve_Brick_Order2_SystemPDE_Paso_PCG_Jacobi))
352     s=unittest.TextTestRunner(verbosity=2).run(suite)

  ViewVC Help
Powered by ViewVC 1.1.26