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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1488 - (show annotations)
Fri Apr 11 00:22:31 2008 UTC (11 years, 6 months ago) by artak
File MIME type: text/x-python
File size: 15864 byte(s)
restart parameter is changed in GMRES(m)
1 #######################################################
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 def test_StokesProblemCartesian_PCG_P_0(self):
45 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
61
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 def test_StokesProblemCartesian_PCG_P_small(self):
72 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
88
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 def test_StokesProblemCartesian_PCG_P_large(self):
99 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
115
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 def test_StokesProblemCartesian_GMRES_P_0(self):
126 ETA=1.
127 P1=0.
128
129 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=12)
142
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 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 def test_StokesProblemCartesian_PCG_P_0(self):
213 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
233
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 def test_StokesProblemCartesian_PCG_P_small(self):
245 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
264
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 def test_StokesProblemCartesian_PCG_P_large(self):
276 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
295
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 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 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=14)
328
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 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