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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1488 - (hide annotations)
Fri Apr 11 00:22:31 2008 UTC (11 years, 5 months ago) by artak
File MIME type: text/x-python
File size: 15864 byte(s)
restart parameter is changed in GMRES(m)
1 gross 1414 #######################################################
2     #
3     # Copyright 2003-2007 by ACceSS MNRF
4     # Copyright 2007 by University of Queensland
5     #
6     # http://esscc.uq.edu.au
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     #######################################################
12     #
13    
14     __copyright__=""" Copyright (c) 2006 by ACcESS MNRF
15     http://www.access.edu.au
16     Primary Business: Queensland, Australia"""
17     __license__="""Licensed under the Open Software License version 3.0
18     http://www.opensource.org/licenses/osl-3.0.php"""
19     import unittest
20     import tempfile
21    
22     from esys.escript import *
23     from esys.finley import Rectangle
24     import sys
25     import os
26     try:
27     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
28     except KeyError:
29     FINLEY_WORKDIR='.'
30    
31    
32     NE=6
33     TOL=1.e-5
34     VERBOSE=False
35    
36     from esys.escript import *
37     from esys.escript.models import StokesProblemCartesian
38     from esys.finley import Rectangle, Brick
39     class Test_Simple2DModels(unittest.TestCase):
40     def setUp(self):
41     self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
42     def tearDown(self):
43     del self.domain
44 gross 1471 def test_StokesProblemCartesian_PCG_P_0(self):
45 gross 1414 ETA=1.
46     P1=0.
47    
48     x=self.domain.getX()
49     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
50     mask=whereZero(x[0]) * [1.,1.] \
51     +whereZero(x[0]-1) * [1.,1.] \
52     +whereZero(x[1]) * [1.,0.] \
53     +whereZero(x[1]-1) * [1.,1.]
54    
55     sp=StokesProblemCartesian(self.domain)
56    
57     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
58     u0=(1-x[0])*x[0]*[0.,1.]
59     p0=Scalar(P1,ReducedSolution(self.domain))
60 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
61 gross 1414
62     error_v0=Lsup(u[0]-u0[0])
63     error_v1=Lsup(u[1]-u0[1])/0.25
64     zz=P1*x[0]*x[1]-p
65     error_p=Lsup(zz-integrate(zz))
66     # print error_p, error_v0,error_v1
67     self.failUnless(error_p<TOL,"pressure error too large.")
68     self.failUnless(error_v0<TOL,"0-velocity error too large.")
69     self.failUnless(error_v1<TOL,"1-velocity error too large.")
70    
71 gross 1471 def test_StokesProblemCartesian_PCG_P_small(self):
72 gross 1414 ETA=1.
73     P1=1.
74    
75     x=self.domain.getX()
76     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
77     mask=whereZero(x[0]) * [1.,1.] \
78     +whereZero(x[0]-1) * [1.,1.] \
79     +whereZero(x[1]) * [1.,0.] \
80     +whereZero(x[1]-1) * [1.,1.]
81    
82     sp=StokesProblemCartesian(self.domain)
83    
84     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
85     u0=(1-x[0])*x[0]*[0.,1.]
86     p0=Scalar(P1,ReducedSolution(self.domain))
87 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
88 gross 1414
89     error_v0=Lsup(u[0]-u0[0])
90     error_v1=Lsup(u[1]-u0[1])/0.25
91     zz=P1*x[0]*x[1]-p
92     error_p=Lsup(zz-integrate(zz))
93     # print error_p, error_v0,error_v1
94     self.failUnless(error_p<TOL,"pressure error too large.")
95     self.failUnless(error_v0<TOL,"0-velocity error too large.")
96     self.failUnless(error_v1<TOL,"1-velocity error too large.")
97    
98 gross 1471 def test_StokesProblemCartesian_PCG_P_large(self):
99 gross 1414 ETA=1.
100     P1=1000.
101    
102     x=self.domain.getX()
103     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
104     mask=whereZero(x[0]) * [1.,1.] \
105     +whereZero(x[0]-1) * [1.,1.] \
106     +whereZero(x[1]) * [1.,0.] \
107     +whereZero(x[1]-1) * [1.,1.]
108    
109     sp=StokesProblemCartesian(self.domain)
110    
111     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
112     u0=(1-x[0])*x[0]*[0.,1.]
113     p0=Scalar(P1,ReducedSolution(self.domain))
114 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
115 gross 1414
116     error_v0=Lsup(u[0]-u0[0])
117     error_v1=Lsup(u[1]-u0[1])/0.25
118     zz=P1*x[0]*x[1]-p
119     error_p=Lsup(zz-integrate(zz))/P1
120     # print error_p, error_v0,error_v1
121     self.failUnless(error_p<TOL,"pressure error too large.")
122     self.failUnless(error_v0<TOL,"0-velocity error too large.")
123     self.failUnless(error_v1<TOL,"1-velocity error too large.")
124    
125 gross 1471 def test_StokesProblemCartesian_GMRES_P_0(self):
126     ETA=1.
127     P1=0.
128 gross 1414
129 gross 1471 x=self.domain.getX()
130     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
131     mask=whereZero(x[0]) * [1.,1.] \
132     +whereZero(x[0]-1) * [1.,1.] \
133     +whereZero(x[1]) * [1.,0.] \
134     +whereZero(x[1]-1) * [1.,1.]
135    
136     sp=StokesProblemCartesian(self.domain)
137    
138     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
139     u0=(1-x[0])*x[0]*[0.,1.]
140     p0=Scalar(P1,ReducedSolution(self.domain))
141 artak 1488 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=12)
142 gross 1471
143     error_v0=Lsup(u[0]-u0[0])
144     error_v1=Lsup(u[1]-u0[1])/0.25
145     zz=P1*x[0]*x[1]-p
146     error_p=Lsup(zz-integrate(zz))
147     # print error_p, error_v0,error_v1
148     self.failUnless(error_p<TOL,"pressure error too large.")
149     self.failUnless(error_v0<TOL,"0-velocity error too large.")
150     self.failUnless(error_v1<TOL,"1-velocity error too large.")
151    
152     def test_StokesProblemCartesian_GMRES_P_small(self):
153     ETA=1.
154     P1=1.
155    
156     x=self.domain.getX()
157     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
158     mask=whereZero(x[0]) * [1.,1.] \
159     +whereZero(x[0]-1) * [1.,1.] \
160     +whereZero(x[1]) * [1.,0.] \
161     +whereZero(x[1]-1) * [1.,1.]
162    
163     sp=StokesProblemCartesian(self.domain)
164    
165     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
166     u0=(1-x[0])*x[0]*[0.,1.]
167     p0=Scalar(P1,ReducedSolution(self.domain))
168     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
169    
170     error_v0=Lsup(u[0]-u0[0])
171     error_v1=Lsup(u[1]-u0[1])/0.25
172     zz=P1*x[0]*x[1]-p
173     error_p=Lsup(zz-integrate(zz))
174     # print error_p, error_v0,error_v1
175     self.failUnless(error_p<TOL,"pressure error too large.")
176     self.failUnless(error_v0<TOL,"0-velocity error too large.")
177     self.failUnless(error_v1<TOL,"1-velocity error too large.")
178    
179     def test_StokesProblemCartesian_GMRES_P_large(self):
180     ETA=1.
181     P1=1000.
182    
183     x=self.domain.getX()
184     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
185     mask=whereZero(x[0]) * [1.,1.] \
186     +whereZero(x[0]-1) * [1.,1.] \
187     +whereZero(x[1]) * [1.,0.] \
188     +whereZero(x[1]-1) * [1.,1.]
189    
190     sp=StokesProblemCartesian(self.domain)
191    
192     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
193     u0=(1-x[0])*x[0]*[0.,1.]
194     p0=Scalar(P1,ReducedSolution(self.domain))
195     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
196    
197     error_v0=Lsup(u[0]-u0[0])
198     error_v1=Lsup(u[1]-u0[1])/0.25
199     zz=P1*x[0]*x[1]-p
200     error_p=Lsup(zz-integrate(zz))/P1
201     # print error_p, error_v0,error_v1
202     self.failUnless(error_p<TOL,"pressure error too large.")
203     self.failUnless(error_v0<TOL,"0-velocity error too large.")
204     self.failUnless(error_v1<TOL,"1-velocity error too large.")
205    
206    
207 gross 1414 class Test_Simple3DModels(unittest.TestCase):
208     def setUp(self):
209     self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
210     def tearDown(self):
211     del self.domain
212 gross 1471 def test_StokesProblemCartesian_PCG_P_0(self):
213 gross 1414 ETA=1.
214     P1=0.
215    
216     x=self.domain.getX()
217     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
218     x=self.domain.getX()
219     mask=whereZero(x[0]) * [1.,1.,1.] \
220     +whereZero(x[0]-1) * [1.,1.,1.] \
221     +whereZero(x[1]) * [1.,1.,1.] \
222     +whereZero(x[1]-1) * [1.,1.,1.] \
223     +whereZero(x[2]) * [1.,1.,0.] \
224     +whereZero(x[2]-1) * [1.,1.,1.]
225    
226    
227     sp=StokesProblemCartesian(self.domain)
228    
229     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
230     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
231     p0=Scalar(P1,ReducedSolution(self.domain))
232 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
233 gross 1414
234     error_v0=Lsup(u[0]-u0[0])
235     error_v1=Lsup(u[1]-u0[1])
236     error_v2=Lsup(u[2]-u0[2])/0.25**2
237     zz=P1*x[0]*x[1]*x[2]-p
238     error_p=Lsup(zz-integrate(zz))
239     # print error_p, error_v0,error_v1,error_v2
240     self.failUnless(error_p<TOL,"pressure error too large.")
241     self.failUnless(error_v0<TOL,"0-velocity error too large.")
242     self.failUnless(error_v1<TOL,"1-velocity error too large.")
243     self.failUnless(error_v2<TOL,"2-velocity error too large.")
244 gross 1471 def test_StokesProblemCartesian_PCG_P_small(self):
245 gross 1414 ETA=1.
246     P1=1.
247    
248     x=self.domain.getX()
249     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
250     mask=whereZero(x[0]) * [1.,1.,1.] \
251     +whereZero(x[0]-1) * [1.,1.,1.] \
252     +whereZero(x[1]) * [1.,1.,1.] \
253     +whereZero(x[1]-1) * [1.,1.,1.] \
254     +whereZero(x[2]) * [1.,1.,0.] \
255     +whereZero(x[2]-1) * [1.,1.,1.]
256    
257    
258     sp=StokesProblemCartesian(self.domain)
259    
260     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
261     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
262     p0=Scalar(P1,ReducedSolution(self.domain))
263 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
264 gross 1414
265     error_v0=Lsup(u[0]-u0[0])
266     error_v1=Lsup(u[1]-u0[1])
267     error_v2=Lsup(u[2]-u0[2])/0.25**2
268     zz=P1*x[0]*x[1]*x[2]-p
269     error_p=Lsup(zz-integrate(zz))/P1
270     # print error_p, error_v0,error_v1,error_v2
271     self.failUnless(error_p<TOL,"pressure error too large.")
272     self.failUnless(error_v0<TOL,"0-velocity error too large.")
273     self.failUnless(error_v1<TOL,"1-velocity error too large.")
274     self.failUnless(error_v2<TOL,"2-velocity error too large.")
275 gross 1471 def test_StokesProblemCartesian_PCG_P_large(self):
276 gross 1414 ETA=1.
277     P1=1000.
278    
279     x=self.domain.getX()
280     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
281     mask=whereZero(x[0]) * [1.,1.,1.] \
282     +whereZero(x[0]-1) * [1.,1.,1.] \
283     +whereZero(x[1]) * [1.,1.,1.] \
284     +whereZero(x[1]-1) * [1.,1.,1.] \
285     +whereZero(x[2]) * [1.,1.,0.] \
286     +whereZero(x[2]-1) * [1.,1.,1.]
287    
288    
289     sp=StokesProblemCartesian(self.domain)
290    
291     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
292     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
293     p0=Scalar(P1,ReducedSolution(self.domain))
294 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
295 gross 1414
296     error_v0=Lsup(u[0]-u0[0])
297     error_v1=Lsup(u[1]-u0[1])
298     error_v2=Lsup(u[2]-u0[2])/0.25**2
299     zz=P1*x[0]*x[1]*x[2]-p
300     error_p=Lsup(zz-integrate(zz))/P1
301     # print error_p, error_v0,error_v1,error_v2
302     self.failUnless(error_p<TOL,"pressure error too large.")
303     self.failUnless(error_v0<TOL,"0-velocity error too large.")
304     self.failUnless(error_v1<TOL,"1-velocity error too large.")
305     self.failUnless(error_v2<TOL,"2-velocity error too large.")
306    
307 gross 1471 def test_StokesProblemCartesian_GMRES_P_0(self):
308     ETA=1.
309     P1=0.
310    
311     x=self.domain.getX()
312     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
313     x=self.domain.getX()
314     mask=whereZero(x[0]) * [1.,1.,1.] \
315     +whereZero(x[0]-1) * [1.,1.,1.] \
316     +whereZero(x[1]) * [1.,1.,1.] \
317     +whereZero(x[1]-1) * [1.,1.,1.] \
318     +whereZero(x[2]) * [1.,1.,0.] \
319     +whereZero(x[2]-1) * [1.,1.,1.]
320    
321    
322     sp=StokesProblemCartesian(self.domain)
323    
324     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
325     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
326     p0=Scalar(P1,ReducedSolution(self.domain))
327 artak 1488 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=14)
328 gross 1471
329     error_v0=Lsup(u[0]-u0[0])
330     error_v1=Lsup(u[1]-u0[1])
331     error_v2=Lsup(u[2]-u0[2])/0.25**2
332     zz=P1*x[0]*x[1]*x[2]-p
333     error_p=Lsup(zz-integrate(zz))
334     # print error_p, error_v0,error_v1,error_v2
335     self.failUnless(error_p<TOL,"pressure error too large.")
336     self.failUnless(error_v0<TOL,"0-velocity error too large.")
337     self.failUnless(error_v1<TOL,"1-velocity error too large.")
338     self.failUnless(error_v2<TOL,"2-velocity error too large.")
339     def test_StokesProblemCartesian_GMRES_P_small(self):
340     ETA=1.
341     P1=1.
342    
343     x=self.domain.getX()
344     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
345     mask=whereZero(x[0]) * [1.,1.,1.] \
346     +whereZero(x[0]-1) * [1.,1.,1.] \
347     +whereZero(x[1]) * [1.,1.,1.] \
348     +whereZero(x[1]-1) * [1.,1.,1.] \
349     +whereZero(x[2]) * [1.,1.,0.] \
350     +whereZero(x[2]-1) * [1.,1.,1.]
351    
352    
353     sp=StokesProblemCartesian(self.domain)
354    
355     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
356     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
357     p0=Scalar(P1,ReducedSolution(self.domain))
358     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
359    
360     error_v0=Lsup(u[0]-u0[0])
361     error_v1=Lsup(u[1]-u0[1])
362     error_v2=Lsup(u[2]-u0[2])/0.25**2
363     zz=P1*x[0]*x[1]*x[2]-p
364     error_p=Lsup(zz-integrate(zz))/P1
365     # print error_p, error_v0,error_v1,error_v2
366     self.failUnless(error_p<TOL,"pressure error too large.")
367     self.failUnless(error_v0<TOL,"0-velocity error too large.")
368     self.failUnless(error_v1<TOL,"1-velocity error too large.")
369     self.failUnless(error_v2<TOL,"2-velocity error too large.")
370     def test_StokesProblemCartesian_GMRES_P_large(self):
371     ETA=1.
372     P1=1000.
373    
374     x=self.domain.getX()
375     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
376     mask=whereZero(x[0]) * [1.,1.,1.] \
377     +whereZero(x[0]-1) * [1.,1.,1.] \
378     +whereZero(x[1]) * [1.,1.,1.] \
379     +whereZero(x[1]-1) * [1.,1.,1.] \
380     +whereZero(x[2]) * [1.,1.,0.] \
381     +whereZero(x[2]-1) * [1.,1.,1.]
382    
383    
384     sp=StokesProblemCartesian(self.domain)
385    
386     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
387     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
388     p0=Scalar(P1,ReducedSolution(self.domain))
389     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
390    
391     error_v0=Lsup(u[0]-u0[0])
392     error_v1=Lsup(u[1]-u0[1])
393     error_v2=Lsup(u[2]-u0[2])/0.25**2
394     zz=P1*x[0]*x[1]*x[2]-p
395     error_p=Lsup(zz-integrate(zz))/P1
396     # print error_p, error_v0,error_v1,error_v2
397     self.failUnless(error_p<TOL,"pressure error too large.")
398     self.failUnless(error_v0<TOL,"0-velocity error too large.")
399     self.failUnless(error_v1<TOL,"1-velocity error too large.")
400     self.failUnless(error_v2<TOL,"2-velocity error too large.")
401    
402 gross 1414 if __name__ == '__main__':
403     suite = unittest.TestSuite()
404     suite.addTest(unittest.makeSuite(Test_Simple2DModels))
405     suite.addTest(unittest.makeSuite(Test_Simple3DModels))
406     s=unittest.TextTestRunner(verbosity=2).run(suite)
407     if not s.wasSuccessful(): sys.exit(1)
408    

  ViewVC Help
Powered by ViewVC 1.1.26