/[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 1809 - (show annotations)
Thu Sep 25 06:43:44 2008 UTC (10 years, 11 months ago) by ksteube
File MIME type: text/x-python
File size: 30202 byte(s)
Copyright updated in all python files

1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 ########################################################
13
14 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15 Earth Systems Science Computational Center (ESSCC)
16 http://www.uq.edu.au/esscc
17 Primary Business: Queensland, Australia"""
18 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 import unittest
23 import tempfile
24
25 from esys.escript import *
26 from esys.finley import Rectangle
27 import sys
28 import os
29 try:
30 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
31 except KeyError:
32 FINLEY_WORKDIR='.'
33
34
35 NE=6
36 TOL=1.e-5
37 VERBOSE=False or True
38
39 from esys.escript import *
40 from esys.escript.models import StokesProblemCartesian
41 from esys.finley import Rectangle, Brick
42 class Test_Simple2DModels(unittest.TestCase):
43 def setUp(self):
44 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
45 def tearDown(self):
46 del self.domain
47 def test_StokesProblemCartesian_PCG_P_0(self):
48 ETA=1.
49 P1=0.
50
51 x=self.domain.getX()
52 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
53 mask=whereZero(x[0]) * [1.,1.] \
54 +whereZero(x[0]-1) * [1.,1.] \
55 +whereZero(x[1]) * [1.,0.] \
56 +whereZero(x[1]-1) * [1.,1.]
57
58 sp=StokesProblemCartesian(self.domain)
59
60 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
61 u0=(1-x[0])*x[0]*[0.,1.]
62 p0=Scalar(P1,ReducedSolution(self.domain))
63 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
64
65 error_v0=Lsup(u[0]-u0[0])
66 error_v1=Lsup(u[1]-u0[1])/0.25
67 zz=P1*x[0]*x[1]-p
68 error_p=Lsup(zz-integrate(zz))
69 # print error_p, error_v0,error_v1
70 self.failUnless(error_p<TOL,"pressure error too large.")
71 self.failUnless(error_v0<TOL,"0-velocity error too large.")
72 self.failUnless(error_v1<TOL,"1-velocity error too large.")
73
74 def test_StokesProblemCartesian_PCG_P_small(self):
75 ETA=1.
76 P1=1.
77
78 x=self.domain.getX()
79 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
80 mask=whereZero(x[0]) * [1.,1.] \
81 +whereZero(x[0]-1) * [1.,1.] \
82 +whereZero(x[1]) * [1.,0.] \
83 +whereZero(x[1]-1) * [1.,1.]
84
85 sp=StokesProblemCartesian(self.domain)
86
87 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
88 u0=(1-x[0])*x[0]*[0.,1.]
89 p0=Scalar(P1,ReducedSolution(self.domain))
90 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
91
92 error_v0=Lsup(u[0]-u0[0])
93 error_v1=Lsup(u[1]-u0[1])/0.25
94 zz=P1*x[0]*x[1]-p
95 error_p=Lsup(zz-integrate(zz))
96 # print error_p, error_v0,error_v1
97 self.failUnless(error_p<TOL,"pressure error too large.")
98 self.failUnless(error_v0<TOL,"0-velocity error too large.")
99 self.failUnless(error_v1<TOL,"1-velocity error too large.")
100
101 def test_StokesProblemCartesian_PCG_P_large(self):
102 ETA=1.
103 P1=1000.
104
105 x=self.domain.getX()
106 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
107 mask=whereZero(x[0]) * [1.,1.] \
108 +whereZero(x[0]-1) * [1.,1.] \
109 +whereZero(x[1]) * [1.,0.] \
110 +whereZero(x[1]-1) * [1.,1.]
111
112 sp=StokesProblemCartesian(self.domain)
113
114 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
115 u0=(1-x[0])*x[0]*[0.,1.]
116 p0=Scalar(P1,ReducedSolution(self.domain))
117 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
118
119 error_v0=Lsup(u[0]-u0[0])
120 error_v1=Lsup(u[1]-u0[1])/0.25
121 zz=P1*x[0]*x[1]-p
122 error_p=Lsup(zz-integrate(zz))/P1
123 # print error_p, error_v0,error_v1
124 self.failUnless(error_p<TOL,"pressure error too large.")
125 self.failUnless(error_v0<TOL,"0-velocity error too large.")
126 self.failUnless(error_v1<TOL,"1-velocity error too large.")
127
128 def test_StokesProblemCartesian_GMRES_P_0(self):
129 ETA=1.
130 P1=0.
131
132 x=self.domain.getX()
133 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
134 mask=whereZero(x[0]) * [1.,1.] \
135 +whereZero(x[0]-1) * [1.,1.] \
136 +whereZero(x[1]) * [1.,0.] \
137 +whereZero(x[1]-1) * [1.,1.]
138
139 sp=StokesProblemCartesian(self.domain)
140
141 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
142 u0=(1-x[0])*x[0]*[0.,1.]
143 p0=Scalar(P1,ReducedSolution(self.domain))
144 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)
145
146 error_v0=Lsup(u[0]-u0[0])
147 error_v1=Lsup(u[1]-u0[1])/0.25
148 zz=P1*x[0]*x[1]-p
149 error_p=Lsup(zz-integrate(zz))
150
151 # print error_p, error_v0,error_v1
152 self.failUnless(error_p<TOL,"pressure error too large.")
153 self.failUnless(error_v0<TOL,"0-velocity error too large.")
154 self.failUnless(error_v1<TOL,"1-velocity error too large.")
155
156 def test_StokesProblemCartesian_GMRES_P_small(self):
157 ETA=1.
158 P1=1.
159
160 x=self.domain.getX()
161 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
162 mask=whereZero(x[0]) * [1.,1.] \
163 +whereZero(x[0]-1) * [1.,1.] \
164 +whereZero(x[1]) * [1.,0.] \
165 +whereZero(x[1]-1) * [1.,1.]
166
167 sp=StokesProblemCartesian(self.domain)
168
169 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
170 u0=(1-x[0])*x[0]*[0.,1.]
171 p0=Scalar(P1,ReducedSolution(self.domain))
172 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
173
174 error_v0=Lsup(u[0]-u0[0])
175 error_v1=Lsup(u[1]-u0[1])/0.25
176 zz=P1*x[0]*x[1]-p
177 error_p=Lsup(zz-integrate(zz))
178 # print error_p, error_v0,error_v1
179 self.failUnless(error_p<TOL,"pressure error too large.")
180 self.failUnless(error_v0<TOL,"0-velocity error too large.")
181 self.failUnless(error_v1<TOL,"1-velocity error too large.")
182
183 def test_StokesProblemCartesian_GMRES_P_large(self):
184 ETA=1.
185 P1=1000.
186
187 x=self.domain.getX()
188 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
189 mask=whereZero(x[0]) * [1.,1.] \
190 +whereZero(x[0]-1) * [1.,1.] \
191 +whereZero(x[1]) * [1.,0.] \
192 +whereZero(x[1]-1) * [1.,1.]
193
194 sp=StokesProblemCartesian(self.domain)
195
196 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
197 u0=(1-x[0])*x[0]*[0.,1.]
198 p0=Scalar(P1,ReducedSolution(self.domain))
199 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
200
201 error_v0=Lsup(u[0]-u0[0])
202 error_v1=Lsup(u[1]-u0[1])/0.25
203 zz=P1*x[0]*x[1]-p
204 error_p=Lsup(zz-integrate(zz))/P1
205 # print error_p, error_v0,error_v1
206 self.failUnless(error_p<TOL,"pressure error too large.")
207 self.failUnless(error_v0<TOL,"0-velocity error too large.")
208 self.failUnless(error_v1<TOL,"1-velocity error too large.")
209
210 def test_StokesProblemCartesian_MINRES_P_0(self):
211 ETA=1.
212 P1=0.
213
214 x=self.domain.getX()
215 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
216 mask=whereZero(x[0]) * [1.,1.] \
217 +whereZero(x[0]-1) * [1.,1.] \
218 +whereZero(x[1]) * [1.,0.] \
219 +whereZero(x[1]-1) * [1.,1.]
220
221 sp=StokesProblemCartesian(self.domain)
222
223 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
224 u0=(1-x[0])*x[0]*[0.,1.]
225 p0=Scalar(P1,ReducedSolution(self.domain))
226 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
227
228 error_v0=Lsup(u[0]-u0[0])
229 error_v1=Lsup(u[1]-u0[1])/0.25
230 zz=P1*x[0]*x[1]-p
231 error_p=Lsup(zz-integrate(zz))
232 # print error_p, error_v0,error_v1
233 self.failUnless(error_p<TOL,"pressure error too large.")
234 self.failUnless(error_v0<TOL,"0-velocity error too large.")
235 self.failUnless(error_v1<TOL,"1-velocity error too large.")
236
237 def test_StokesProblemCartesian_MINRES_P_small(self):
238 ETA=1.
239 P1=1.
240
241 x=self.domain.getX()
242 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
243 mask=whereZero(x[0]) * [1.,1.] \
244 +whereZero(x[0]-1) * [1.,1.] \
245 +whereZero(x[1]) * [1.,0.] \
246 +whereZero(x[1]-1) * [1.,1.]
247
248 sp=StokesProblemCartesian(self.domain)
249
250 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
251 u0=(1-x[0])*x[0]*[0.,1.]
252 p0=Scalar(P1,ReducedSolution(self.domain))
253 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
254
255 error_v0=Lsup(u[0]-u0[0])
256 error_v1=Lsup(u[1]-u0[1])/0.25
257 zz=P1*x[0]*x[1]-p
258 error_p=Lsup(zz-integrate(zz))/P1
259 # print error_p, error_v0,error_v1
260 self.failUnless(error_p<TOL,"pressure error too large.")
261 self.failUnless(error_v0<TOL,"0-velocity error too large.")
262 self.failUnless(error_v1<TOL,"1-velocity error too large.")
263
264 # def test_StokesProblemCartesian_MINRES_P_large(self):
265 # ETA=1.
266 # P1=1000.
267 #
268 # x=self.domain.getX()
269 # F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
270 # mask=whereZero(x[0]) * [1.,1.] \
271 # +whereZero(x[0]-1) * [1.,1.] \
272 # +whereZero(x[1]) * [1.,0.] \
273 # +whereZero(x[1]-1) * [1.,1.]
274
275 # sp=StokesProblemCartesian(self.domain)
276
277 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
278 # u0=(1-x[0])*x[0]*[0.,1.]
279 # p0=Scalar(P1,ReducedSolution(self.domain))
280 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
281
282 # error_v0=Lsup(u[0]-u0[0])
283 # error_v1=Lsup(u[1]-u0[1])/0.25
284 # zz=P1*x[0]*x[1]-p
285 # error_p=Lsup(zz-integrate(zz))/P1
286 # print error_p, error_v0,error_v1
287 # self.failUnless(error_p<TOL,"pressure error too large.")
288 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
289 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
290
291
292 # def test_StokesProblemCartesian_TFQMR_P_0(self):
293 # ETA=1.
294 # P1=0.
295
296 # x=self.domain.getX()
297 # F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
298 # mask=whereZero(x[0]) * [1.,1.] \
299 # +whereZero(x[0]-1) * [1.,1.] \
300 # +whereZero(x[1]) * [1.,0.] \
301 # +whereZero(x[1]-1) * [1.,1.]
302
303 # sp=StokesProblemCartesian(self.domain)
304
305 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
306 # u0=(1-x[0])*x[0]*[0.,1.]
307 # p0=Scalar(P1,ReducedSolution(self.domain))
308 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
309
310 # error_v0=Lsup(u[0]-u0[0])
311 # error_v1=Lsup(u[1]-u0[1])/0.25
312 # zz=P1*x[0]*x[1]-p
313 # error_p=Lsup(zz-integrate(zz))
314 # print error_p, error_v0,error_v1
315 # self.failUnless(error_p<TOL,"pressure error too large.")
316 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
317 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
318
319 def test_StokesProblemCartesian_TFQMR_P_small(self):
320 ETA=1.
321 P1=1.
322
323 x=self.domain.getX()
324 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
325 mask=whereZero(x[0]) * [1.,1.] \
326 +whereZero(x[0]-1) * [1.,1.] \
327 +whereZero(x[1]) * [1.,0.] \
328 +whereZero(x[1]-1) * [1.,1.]
329
330 sp=StokesProblemCartesian(self.domain)
331
332 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
333 u0=(1-x[0])*x[0]*[0.,1.]
334 p0=Scalar(P1,ReducedSolution(self.domain))
335 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
336
337 error_v0=Lsup(u[0]-u0[0])
338 error_v1=Lsup(u[1]-u0[1])/0.25
339 zz=P1*x[0]*x[1]-p
340 error_p=Lsup(zz-integrate(zz))/P1
341 # print error_p, error_v0,error_v1
342 self.failUnless(error_p<TOL,"pressure error too large.")
343 self.failUnless(error_v0<TOL,"0-velocity error too large.")
344 self.failUnless(error_v1<TOL,"1-velocity error too large.")
345
346 def test_StokesProblemCartesian_TFQMR_P_large(self):
347 ETA=1.
348 P1=1000.
349
350 x=self.domain.getX()
351 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
352 mask=whereZero(x[0]) * [1.,1.] \
353 +whereZero(x[0]-1) * [1.,1.] \
354 +whereZero(x[1]) * [1.,0.] \
355 +whereZero(x[1]-1) * [1.,1.]
356
357 sp=StokesProblemCartesian(self.domain)
358
359 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
360 u0=(1-x[0])*x[0]*[0.,1.]
361 p0=Scalar(P1,ReducedSolution(self.domain))
362 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
363
364 error_v0=Lsup(u[0]-u0[0])
365 error_v1=Lsup(u[1]-u0[1])/0.25
366 zz=P1*x[0]*x[1]-p
367 error_p=Lsup(zz-integrate(zz))/P1
368 # print error_p, error_v0,error_v1
369 self.failUnless(error_p<TOL,"pressure error too large.")
370 self.failUnless(error_v0<TOL,"0-velocity error too large.")
371 self.failUnless(error_v1<TOL,"1-velocity error too large.")
372
373
374
375 class Test_Simple3DModels(unittest.TestCase):
376 def setUp(self):
377 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
378 def tearDown(self):
379 del self.domain
380 def test_StokesProblemCartesian_PCG_P_0(self):
381 ETA=1.
382 P1=0.
383
384 x=self.domain.getX()
385 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.]
386 x=self.domain.getX()
387 mask=whereZero(x[0]) * [1.,1.,1.] \
388 +whereZero(x[0]-1) * [1.,1.,1.] \
389 +whereZero(x[1]) * [1.,1.,1.] \
390 +whereZero(x[1]-1) * [1.,1.,1.] \
391 +whereZero(x[2]) * [1.,1.,0.] \
392 +whereZero(x[2]-1) * [1.,1.,1.]
393
394
395 sp=StokesProblemCartesian(self.domain)
396
397 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
398 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
399 p0=Scalar(P1,ReducedSolution(self.domain))
400 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
401
402 error_v0=Lsup(u[0]-u0[0])
403 error_v1=Lsup(u[1]-u0[1])
404 error_v2=Lsup(u[2]-u0[2])/0.25**2
405 zz=P1*x[0]*x[1]*x[2]-p
406 error_p=Lsup(zz-integrate(zz))
407 # print error_p, error_v0,error_v1,error_v2
408 self.failUnless(error_p<TOL,"pressure error too large.")
409 self.failUnless(error_v0<TOL,"0-velocity error too large.")
410 self.failUnless(error_v1<TOL,"1-velocity error too large.")
411 self.failUnless(error_v2<TOL,"2-velocity error too large.")
412
413 def test_StokesProblemCartesian_PCG_P_small(self):
414 ETA=1.
415 P1=1.
416
417 x=self.domain.getX()
418 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.]
419 mask=whereZero(x[0]) * [1.,1.,1.] \
420 +whereZero(x[0]-1) * [1.,1.,1.] \
421 +whereZero(x[1]) * [1.,1.,1.] \
422 +whereZero(x[1]-1) * [1.,1.,1.] \
423 +whereZero(x[2]) * [1.,1.,0.] \
424 +whereZero(x[2]-1) * [1.,1.,1.]
425
426
427 sp=StokesProblemCartesian(self.domain)
428
429 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
430 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
431 p0=Scalar(P1,ReducedSolution(self.domain))
432 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
433
434 error_v0=Lsup(u[0]-u0[0])
435 error_v1=Lsup(u[1]-u0[1])
436 error_v2=Lsup(u[2]-u0[2])/0.25**2
437 zz=P1*x[0]*x[1]*x[2]-p
438 error_p=Lsup(zz-integrate(zz))/P1
439 # print error_p, error_v0,error_v1,error_v2
440 self.failUnless(error_p<TOL,"pressure error too large.")
441 self.failUnless(error_v0<TOL,"0-velocity error too large.")
442 self.failUnless(error_v1<TOL,"1-velocity error too large.")
443 self.failUnless(error_v2<TOL,"2-velocity error too large.")
444 def test_StokesProblemCartesian_PCG_P_large(self):
445 ETA=1.
446 P1=1000.
447
448 x=self.domain.getX()
449 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.]
450 mask=whereZero(x[0]) * [1.,1.,1.] \
451 +whereZero(x[0]-1) * [1.,1.,1.] \
452 +whereZero(x[1]) * [1.,1.,1.] \
453 +whereZero(x[1]-1) * [1.,1.,1.] \
454 +whereZero(x[2]) * [1.,1.,0.] \
455 +whereZero(x[2]-1) * [1.,1.,1.]
456
457
458 sp=StokesProblemCartesian(self.domain)
459
460 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
461 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
462 p0=Scalar(P1,ReducedSolution(self.domain))
463 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
464
465 error_v0=Lsup(u[0]-u0[0])
466 error_v1=Lsup(u[1]-u0[1])
467 error_v2=Lsup(u[2]-u0[2])/0.25**2
468 zz=P1*x[0]*x[1]*x[2]-p
469 error_p=Lsup(zz-integrate(zz))/P1
470 # print error_p, error_v0,error_v1,error_v2
471 self.failUnless(error_p<TOL,"pressure error too large.")
472 self.failUnless(error_v0<TOL,"0-velocity error too large.")
473 self.failUnless(error_v1<TOL,"1-velocity error too large.")
474 self.failUnless(error_v2<TOL,"2-velocity error too large.")
475
476 def test_StokesProblemCartesian_GMRES_P_0(self):
477 ETA=1.
478 P1=0.
479
480 x=self.domain.getX()
481 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.]
482 x=self.domain.getX()
483 mask=whereZero(x[0]) * [1.,1.,1.] \
484 +whereZero(x[0]-1) * [1.,1.,1.] \
485 +whereZero(x[1]) * [1.,1.,1.] \
486 +whereZero(x[1]-1) * [1.,1.,1.] \
487 +whereZero(x[2]) * [1.,1.,0.] \
488 +whereZero(x[2]-1) * [1.,1.,1.]
489
490
491 sp=StokesProblemCartesian(self.domain)
492
493 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
494 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
495 p0=Scalar(P1,ReducedSolution(self.domain))
496 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)
497
498 error_v0=Lsup(u[0]-u0[0])
499 error_v1=Lsup(u[1]-u0[1])
500 error_v2=Lsup(u[2]-u0[2])/0.25**2
501 zz=P1*x[0]*x[1]*x[2]-p
502 error_p=Lsup(zz-integrate(zz))
503 # print error_p, error_v0,error_v1,error_v2
504 self.failUnless(error_p<TOL,"pressure error too large.")
505 self.failUnless(error_v0<TOL,"0-velocity error too large.")
506 self.failUnless(error_v1<TOL,"1-velocity error too large.")
507 self.failUnless(error_v2<TOL,"2-velocity error too large.")
508 def test_StokesProblemCartesian_GMRES_P_small(self):
509 ETA=1.
510 P1=1.
511
512 x=self.domain.getX()
513 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.]
514 mask=whereZero(x[0]) * [1.,1.,1.] \
515 +whereZero(x[0]-1) * [1.,1.,1.] \
516 +whereZero(x[1]) * [1.,1.,1.] \
517 +whereZero(x[1]-1) * [1.,1.,1.] \
518 +whereZero(x[2]) * [1.,1.,0.] \
519 +whereZero(x[2]-1) * [1.,1.,1.]
520
521
522 sp=StokesProblemCartesian(self.domain)
523
524 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
525 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
526 p0=Scalar(P1,ReducedSolution(self.domain))
527 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
528
529 error_v0=Lsup(u[0]-u0[0])
530 error_v1=Lsup(u[1]-u0[1])
531 error_v2=Lsup(u[2]-u0[2])/0.25**2
532 zz=P1*x[0]*x[1]*x[2]-p
533 error_p=Lsup(zz-integrate(zz))/P1
534 # print error_p, error_v0,error_v1,error_v2
535 self.failUnless(error_p<TOL,"pressure error too large.")
536 self.failUnless(error_v0<TOL,"0-velocity error too large.")
537 self.failUnless(error_v1<TOL,"1-velocity error too large.")
538 self.failUnless(error_v2<TOL,"2-velocity error too large.")
539 def test_StokesProblemCartesian_GMRES_P_large(self):
540 ETA=1.
541 P1=1000.
542
543 x=self.domain.getX()
544 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.]
545 mask=whereZero(x[0]) * [1.,1.,1.] \
546 +whereZero(x[0]-1) * [1.,1.,1.] \
547 +whereZero(x[1]) * [1.,1.,1.] \
548 +whereZero(x[1]-1) * [1.,1.,1.] \
549 +whereZero(x[2]) * [1.,1.,0.] \
550 +whereZero(x[2]-1) * [1.,1.,1.]
551
552
553 sp=StokesProblemCartesian(self.domain)
554
555 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
556 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
557 p0=Scalar(P1,ReducedSolution(self.domain))
558 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
559
560 error_v0=Lsup(u[0]-u0[0])
561 error_v1=Lsup(u[1]-u0[1])
562 error_v2=Lsup(u[2]-u0[2])/0.25**2
563 zz=P1*x[0]*x[1]*x[2]-p
564 error_p=Lsup(zz-integrate(zz))/P1
565 # print error_p, error_v0,error_v1,error_v2
566 self.failUnless(error_p<TOL,"pressure error too large.")
567 self.failUnless(error_v0<TOL,"0-velocity error too large.")
568 self.failUnless(error_v1<TOL,"1-velocity error too large.")
569 self.failUnless(error_v2<TOL,"2-velocity error too large.")
570
571 def test_StokesProblemCartesian_MINRES_P_0(self):
572 ETA=1.
573 P1=0.
574
575 x=self.domain.getX()
576 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.]
577 x=self.domain.getX()
578 mask=whereZero(x[0]) * [1.,1.,1.] \
579 +whereZero(x[0]-1) * [1.,1.,1.] \
580 +whereZero(x[1]) * [1.,1.,1.] \
581 +whereZero(x[1]-1) * [1.,1.,1.] \
582 +whereZero(x[2]) * [1.,1.,0.] \
583 +whereZero(x[2]-1) * [1.,1.,1.]
584
585
586 sp=StokesProblemCartesian(self.domain)
587
588 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
589 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
590 p0=Scalar(P1,ReducedSolution(self.domain))
591 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
592
593 error_v0=Lsup(u[0]-u0[0])
594 error_v1=Lsup(u[1]-u0[1])
595 error_v2=Lsup(u[2]-u0[2])/0.25**2
596 zz=P1*x[0]*x[1]*x[2]-p
597 error_p=Lsup(zz-integrate(zz))
598 # print error_p, error_v0,error_v1,error_v2
599 self.failUnless(error_p<TOL,"pressure error too large.")
600 self.failUnless(error_v0<TOL,"0-velocity error too large.")
601 self.failUnless(error_v1<TOL,"1-velocity error too large.")
602 self.failUnless(error_v2<TOL,"2-velocity error too large.")
603
604 def test_StokesProblemCartesian_MINRES_P_small(self):
605 ETA=1.
606 P1=1.
607
608 x=self.domain.getX()
609 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.]
610 mask=whereZero(x[0]) * [1.,1.,1.] \
611 +whereZero(x[0]-1) * [1.,1.,1.] \
612 +whereZero(x[1]) * [1.,1.,1.] \
613 +whereZero(x[1]-1) * [1.,1.,1.] \
614 +whereZero(x[2]) * [1.,1.,0.] \
615 +whereZero(x[2]-1) * [1.,1.,1.]
616
617
618 sp=StokesProblemCartesian(self.domain)
619
620 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
621 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
622 p0=Scalar(P1,ReducedSolution(self.domain))
623 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
624
625 error_v0=Lsup(u[0]-u0[0])
626 error_v1=Lsup(u[1]-u0[1])
627 error_v2=Lsup(u[2]-u0[2])/0.25**2
628 zz=P1*x[0]*x[1]*x[2]-p
629 error_p=Lsup(zz-integrate(zz))/P1
630 # print error_p, error_v0,error_v1,error_v2
631 self.failUnless(error_p<TOL,"pressure error too large.")
632 self.failUnless(error_v0<TOL,"0-velocity error too large.")
633 self.failUnless(error_v1<TOL,"1-velocity error too large.")
634 self.failUnless(error_v2<TOL,"2-velocity error too large.")
635
636 # def test_StokesProblemCartesian_MINRES_P_large(self):
637 # ETA=1.
638 # P1=1000.
639
640 # x=self.domain.getX()
641 # 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.]
642 # mask=whereZero(x[0]) * [1.,1.,1.] \
643 # +whereZero(x[0]-1) * [1.,1.,1.] \
644 # +whereZero(x[1]) * [1.,1.,1.] \
645 # +whereZero(x[1]-1) * [1.,1.,1.] \
646 # +whereZero(x[2]) * [1.,1.,0.] \
647 # +whereZero(x[2]-1) * [1.,1.,1.]
648
649
650 # sp=StokesProblemCartesian(self.domain)
651
652 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
653 # u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
654 # p0=Scalar(P1,ReducedSolution(self.domain))
655 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
656
657 # error_v0=Lsup(u[0]-u0[0])
658 # error_v1=Lsup(u[1]-u0[1])
659 # error_v2=Lsup(u[2]-u0[2])/0.25**2
660 # zz=P1*x[0]*x[1]*x[2]-p
661 # error_p=Lsup(zz-integrate(zz))/P1
662 # print error_p, error_v0,error_v1,error_v2
663 # self.failUnless(error_p<TOL,"pressure error too large.")
664 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
665 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
666 # self.failUnless(error_v2<TOL,"2-velocity error too large.")
667
668 # def test_StokesProblemCartesian_TFQMR_P_0(self):
669 # ETA=1.
670 # P1=0.
671
672 # x=self.domain.getX()
673 # 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.]
674 # x=self.domain.getX()
675 # mask=whereZero(x[0]) * [1.,1.,1.] \
676 # +whereZero(x[0]-1) * [1.,1.,1.] \
677 # +whereZero(x[1]) * [1.,1.,1.] \
678 # +whereZero(x[1]-1) * [1.,1.,1.] \
679 # +whereZero(x[2]) * [1.,1.,0.] \
680 # +whereZero(x[2]-1) * [1.,1.,1.]
681
682
683 # sp=StokesProblemCartesian(self.domain)
684
685 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
686 # u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
687 # p0=Scalar(P1,ReducedSolution(self.domain))
688 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
689
690 # error_v0=Lsup(u[0]-u0[0])
691 # error_v1=Lsup(u[1]-u0[1])
692 # error_v2=Lsup(u[2]-u0[2])/0.25**2
693 # zz=P1*x[0]*x[1]*x[2]-p
694 # error_p=Lsup(zz-integrate(zz))
695 # print error_p, error_v0,error_v1,error_v2
696 # self.failUnless(error_p<TOL,"pressure error too large.")
697 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
698 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
699 # self.failUnless(error_v2<TOL,"2-velocity error too large.")
700
701 def test_StokesProblemCartesian_TFQMR_P_small(self):
702 ETA=1.
703 P1=1.
704
705 x=self.domain.getX()
706 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.]
707 mask=whereZero(x[0]) * [1.,1.,1.] \
708 +whereZero(x[0]-1) * [1.,1.,1.] \
709 +whereZero(x[1]) * [1.,1.,1.] \
710 +whereZero(x[1]-1) * [1.,1.,1.] \
711 +whereZero(x[2]) * [1.,1.,0.] \
712 +whereZero(x[2]-1) * [1.,1.,1.]
713
714
715 sp=StokesProblemCartesian(self.domain)
716
717 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
718 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
719 p0=Scalar(P1,ReducedSolution(self.domain))
720 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
721
722 error_v0=Lsup(u[0]-u0[0])
723 error_v1=Lsup(u[1]-u0[1])
724 error_v2=Lsup(u[2]-u0[2])/0.25**2
725 zz=P1*x[0]*x[1]*x[2]-p
726 error_p=Lsup(zz-integrate(zz))/P1
727 # print error_p, error_v0,error_v1,error_v2
728 self.failUnless(error_p<TOL,"pressure error too large.")
729 self.failUnless(error_v0<TOL,"0-velocity error too large.")
730 self.failUnless(error_v1<TOL,"1-velocity error too large.")
731 self.failUnless(error_v2<TOL,"2-velocity error too large.")
732
733 def test_StokesProblemCartesian_TFQMR_P_large(self):
734 ETA=1.
735 P1=1000.
736
737 x=self.domain.getX()
738 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.]
739 mask=whereZero(x[0]) * [1.,1.,1.] \
740 +whereZero(x[0]-1) * [1.,1.,1.] \
741 +whereZero(x[1]) * [1.,1.,1.] \
742 +whereZero(x[1]-1) * [1.,1.,1.] \
743 +whereZero(x[2]) * [1.,1.,0.] \
744 +whereZero(x[2]-1) * [1.,1.,1.]
745
746
747 sp=StokesProblemCartesian(self.domain)
748
749 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
750 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
751 p0=Scalar(P1,ReducedSolution(self.domain))
752 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
753
754 error_v0=Lsup(u[0]-u0[0])
755 error_v1=Lsup(u[1]-u0[1])
756 error_v2=Lsup(u[2]-u0[2])/0.25**2
757 zz=P1*x[0]*x[1]*x[2]-p
758 error_p=Lsup(zz-integrate(zz))/P1
759 # print error_p, error_v0,error_v1,error_v2
760 self.failUnless(error_p<TOL,"pressure error too large.")
761 self.failUnless(error_v0<TOL,"0-velocity error too large.")
762 self.failUnless(error_v1<TOL,"1-velocity error too large.")
763 self.failUnless(error_v2<TOL,"2-velocity error too large.")
764
765
766 if __name__ == '__main__':
767 suite = unittest.TestSuite()
768 suite.addTest(unittest.makeSuite(Test_Simple2DModels))
769 # suite.addTest(unittest.makeSuite(Test_Simple3DModels))
770 s=unittest.TextTestRunner(verbosity=2).run(suite)
771 if not s.wasSuccessful(): sys.exit(1)
772

  ViewVC Help
Powered by ViewVC 1.1.26