/[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 1562 - (show annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 30142 byte(s)
The algebraic upwinding with MPI. The case of boundary constraint needs still some attention. 


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 or True
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=20)
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
148 # print error_p, error_v0,error_v1
149 self.failUnless(error_p<TOL,"pressure error too large.")
150 self.failUnless(error_v0<TOL,"0-velocity error too large.")
151 self.failUnless(error_v1<TOL,"1-velocity error too large.")
152
153 def test_StokesProblemCartesian_GMRES_P_small(self):
154 ETA=1.
155 P1=1.
156
157 x=self.domain.getX()
158 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
159 mask=whereZero(x[0]) * [1.,1.] \
160 +whereZero(x[0]-1) * [1.,1.] \
161 +whereZero(x[1]) * [1.,0.] \
162 +whereZero(x[1]-1) * [1.,1.]
163
164 sp=StokesProblemCartesian(self.domain)
165
166 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
167 u0=(1-x[0])*x[0]*[0.,1.]
168 p0=Scalar(P1,ReducedSolution(self.domain))
169 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
170
171 error_v0=Lsup(u[0]-u0[0])
172 error_v1=Lsup(u[1]-u0[1])/0.25
173 zz=P1*x[0]*x[1]-p
174 error_p=Lsup(zz-integrate(zz))
175 # print error_p, error_v0,error_v1
176 self.failUnless(error_p<TOL,"pressure error too large.")
177 self.failUnless(error_v0<TOL,"0-velocity error too large.")
178 self.failUnless(error_v1<TOL,"1-velocity error too large.")
179
180 def test_StokesProblemCartesian_GMRES_P_large(self):
181 ETA=1.
182 P1=1000.
183
184 x=self.domain.getX()
185 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
186 mask=whereZero(x[0]) * [1.,1.] \
187 +whereZero(x[0]-1) * [1.,1.] \
188 +whereZero(x[1]) * [1.,0.] \
189 +whereZero(x[1]-1) * [1.,1.]
190
191 sp=StokesProblemCartesian(self.domain)
192
193 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
194 u0=(1-x[0])*x[0]*[0.,1.]
195 p0=Scalar(P1,ReducedSolution(self.domain))
196 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
197
198 error_v0=Lsup(u[0]-u0[0])
199 error_v1=Lsup(u[1]-u0[1])/0.25
200 zz=P1*x[0]*x[1]-p
201 error_p=Lsup(zz-integrate(zz))/P1
202 # print error_p, error_v0,error_v1
203 self.failUnless(error_p<TOL,"pressure error too large.")
204 self.failUnless(error_v0<TOL,"0-velocity error too large.")
205 self.failUnless(error_v1<TOL,"1-velocity error too large.")
206
207 def test_StokesProblemCartesian_MINRES_P_0(self):
208 ETA=1.
209 P1=0.
210
211 x=self.domain.getX()
212 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
213 mask=whereZero(x[0]) * [1.,1.] \
214 +whereZero(x[0]-1) * [1.,1.] \
215 +whereZero(x[1]) * [1.,0.] \
216 +whereZero(x[1]-1) * [1.,1.]
217
218 sp=StokesProblemCartesian(self.domain)
219
220 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
221 u0=(1-x[0])*x[0]*[0.,1.]
222 p0=Scalar(P1,ReducedSolution(self.domain))
223 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
224
225 error_v0=Lsup(u[0]-u0[0])
226 error_v1=Lsup(u[1]-u0[1])/0.25
227 zz=P1*x[0]*x[1]-p
228 error_p=Lsup(zz-integrate(zz))
229 # print error_p, error_v0,error_v1
230 self.failUnless(error_p<TOL,"pressure error too large.")
231 self.failUnless(error_v0<TOL,"0-velocity error too large.")
232 self.failUnless(error_v1<TOL,"1-velocity error too large.")
233
234 def test_StokesProblemCartesian_MINRES_P_small(self):
235 ETA=1.
236 P1=1.
237
238 x=self.domain.getX()
239 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
240 mask=whereZero(x[0]) * [1.,1.] \
241 +whereZero(x[0]-1) * [1.,1.] \
242 +whereZero(x[1]) * [1.,0.] \
243 +whereZero(x[1]-1) * [1.,1.]
244
245 sp=StokesProblemCartesian(self.domain)
246
247 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
248 u0=(1-x[0])*x[0]*[0.,1.]
249 p0=Scalar(P1,ReducedSolution(self.domain))
250 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
251
252 error_v0=Lsup(u[0]-u0[0])
253 error_v1=Lsup(u[1]-u0[1])/0.25
254 zz=P1*x[0]*x[1]-p
255 error_p=Lsup(zz-integrate(zz))/P1
256 # print error_p, error_v0,error_v1
257 self.failUnless(error_p<TOL,"pressure error too large.")
258 self.failUnless(error_v0<TOL,"0-velocity error too large.")
259 self.failUnless(error_v1<TOL,"1-velocity error too large.")
260
261 # def test_StokesProblemCartesian_MINRES_P_large(self):
262 # ETA=1.
263 # P1=1000.
264 #
265 # x=self.domain.getX()
266 # F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
267 # mask=whereZero(x[0]) * [1.,1.] \
268 # +whereZero(x[0]-1) * [1.,1.] \
269 # +whereZero(x[1]) * [1.,0.] \
270 # +whereZero(x[1]-1) * [1.,1.]
271
272 # sp=StokesProblemCartesian(self.domain)
273
274 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
275 # u0=(1-x[0])*x[0]*[0.,1.]
276 # p0=Scalar(P1,ReducedSolution(self.domain))
277 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
278
279 # error_v0=Lsup(u[0]-u0[0])
280 # error_v1=Lsup(u[1]-u0[1])/0.25
281 # zz=P1*x[0]*x[1]-p
282 # error_p=Lsup(zz-integrate(zz))/P1
283 # print error_p, error_v0,error_v1
284 # self.failUnless(error_p<TOL,"pressure error too large.")
285 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
286 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
287
288
289 # def test_StokesProblemCartesian_TFQMR_P_0(self):
290 # ETA=1.
291 # P1=0.
292
293 # x=self.domain.getX()
294 # F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
295 # mask=whereZero(x[0]) * [1.,1.] \
296 # +whereZero(x[0]-1) * [1.,1.] \
297 # +whereZero(x[1]) * [1.,0.] \
298 # +whereZero(x[1]-1) * [1.,1.]
299
300 # sp=StokesProblemCartesian(self.domain)
301
302 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
303 # u0=(1-x[0])*x[0]*[0.,1.]
304 # p0=Scalar(P1,ReducedSolution(self.domain))
305 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
306
307 # error_v0=Lsup(u[0]-u0[0])
308 # error_v1=Lsup(u[1]-u0[1])/0.25
309 # zz=P1*x[0]*x[1]-p
310 # error_p=Lsup(zz-integrate(zz))
311 # print error_p, error_v0,error_v1
312 # self.failUnless(error_p<TOL,"pressure error too large.")
313 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
314 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
315
316 def test_StokesProblemCartesian_TFQMR_P_small(self):
317 ETA=1.
318 P1=1.
319
320 x=self.domain.getX()
321 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
322 mask=whereZero(x[0]) * [1.,1.] \
323 +whereZero(x[0]-1) * [1.,1.] \
324 +whereZero(x[1]) * [1.,0.] \
325 +whereZero(x[1]-1) * [1.,1.]
326
327 sp=StokesProblemCartesian(self.domain)
328
329 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
330 u0=(1-x[0])*x[0]*[0.,1.]
331 p0=Scalar(P1,ReducedSolution(self.domain))
332 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
333
334 error_v0=Lsup(u[0]-u0[0])
335 error_v1=Lsup(u[1]-u0[1])/0.25
336 zz=P1*x[0]*x[1]-p
337 error_p=Lsup(zz-integrate(zz))/P1
338 # print error_p, error_v0,error_v1
339 self.failUnless(error_p<TOL,"pressure error too large.")
340 self.failUnless(error_v0<TOL,"0-velocity error too large.")
341 self.failUnless(error_v1<TOL,"1-velocity error too large.")
342
343 def test_StokesProblemCartesian_TFQMR_P_large(self):
344 ETA=1.
345 P1=1000.
346
347 x=self.domain.getX()
348 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
349 mask=whereZero(x[0]) * [1.,1.] \
350 +whereZero(x[0]-1) * [1.,1.] \
351 +whereZero(x[1]) * [1.,0.] \
352 +whereZero(x[1]-1) * [1.,1.]
353
354 sp=StokesProblemCartesian(self.domain)
355
356 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
357 u0=(1-x[0])*x[0]*[0.,1.]
358 p0=Scalar(P1,ReducedSolution(self.domain))
359 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
360
361 error_v0=Lsup(u[0]-u0[0])
362 error_v1=Lsup(u[1]-u0[1])/0.25
363 zz=P1*x[0]*x[1]-p
364 error_p=Lsup(zz-integrate(zz))/P1
365 # print error_p, error_v0,error_v1
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
370
371
372 class Test_Simple3DModels(unittest.TestCase):
373 def setUp(self):
374 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
375 def tearDown(self):
376 del self.domain
377 def test_StokesProblemCartesian_PCG_P_0(self):
378 ETA=1.
379 P1=0.
380
381 x=self.domain.getX()
382 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.]
383 x=self.domain.getX()
384 mask=whereZero(x[0]) * [1.,1.,1.] \
385 +whereZero(x[0]-1) * [1.,1.,1.] \
386 +whereZero(x[1]) * [1.,1.,1.] \
387 +whereZero(x[1]-1) * [1.,1.,1.] \
388 +whereZero(x[2]) * [1.,1.,0.] \
389 +whereZero(x[2]-1) * [1.,1.,1.]
390
391
392 sp=StokesProblemCartesian(self.domain)
393
394 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
395 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
396 p0=Scalar(P1,ReducedSolution(self.domain))
397 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
398
399 error_v0=Lsup(u[0]-u0[0])
400 error_v1=Lsup(u[1]-u0[1])
401 error_v2=Lsup(u[2]-u0[2])/0.25**2
402 zz=P1*x[0]*x[1]*x[2]-p
403 error_p=Lsup(zz-integrate(zz))
404 # print error_p, error_v0,error_v1,error_v2
405 self.failUnless(error_p<TOL,"pressure error too large.")
406 self.failUnless(error_v0<TOL,"0-velocity error too large.")
407 self.failUnless(error_v1<TOL,"1-velocity error too large.")
408 self.failUnless(error_v2<TOL,"2-velocity error too large.")
409
410 def test_StokesProblemCartesian_PCG_P_small(self):
411 ETA=1.
412 P1=1.
413
414 x=self.domain.getX()
415 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.]
416 mask=whereZero(x[0]) * [1.,1.,1.] \
417 +whereZero(x[0]-1) * [1.,1.,1.] \
418 +whereZero(x[1]) * [1.,1.,1.] \
419 +whereZero(x[1]-1) * [1.,1.,1.] \
420 +whereZero(x[2]) * [1.,1.,0.] \
421 +whereZero(x[2]-1) * [1.,1.,1.]
422
423
424 sp=StokesProblemCartesian(self.domain)
425
426 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
427 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
428 p0=Scalar(P1,ReducedSolution(self.domain))
429 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
430
431 error_v0=Lsup(u[0]-u0[0])
432 error_v1=Lsup(u[1]-u0[1])
433 error_v2=Lsup(u[2]-u0[2])/0.25**2
434 zz=P1*x[0]*x[1]*x[2]-p
435 error_p=Lsup(zz-integrate(zz))/P1
436 # print error_p, error_v0,error_v1,error_v2
437 self.failUnless(error_p<TOL,"pressure error too large.")
438 self.failUnless(error_v0<TOL,"0-velocity error too large.")
439 self.failUnless(error_v1<TOL,"1-velocity error too large.")
440 self.failUnless(error_v2<TOL,"2-velocity error too large.")
441 def test_StokesProblemCartesian_PCG_P_large(self):
442 ETA=1.
443 P1=1000.
444
445 x=self.domain.getX()
446 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.]
447 mask=whereZero(x[0]) * [1.,1.,1.] \
448 +whereZero(x[0]-1) * [1.,1.,1.] \
449 +whereZero(x[1]) * [1.,1.,1.] \
450 +whereZero(x[1]-1) * [1.,1.,1.] \
451 +whereZero(x[2]) * [1.,1.,0.] \
452 +whereZero(x[2]-1) * [1.,1.,1.]
453
454
455 sp=StokesProblemCartesian(self.domain)
456
457 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
458 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
459 p0=Scalar(P1,ReducedSolution(self.domain))
460 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
461
462 error_v0=Lsup(u[0]-u0[0])
463 error_v1=Lsup(u[1]-u0[1])
464 error_v2=Lsup(u[2]-u0[2])/0.25**2
465 zz=P1*x[0]*x[1]*x[2]-p
466 error_p=Lsup(zz-integrate(zz))/P1
467 # print error_p, error_v0,error_v1,error_v2
468 self.failUnless(error_p<TOL,"pressure error too large.")
469 self.failUnless(error_v0<TOL,"0-velocity error too large.")
470 self.failUnless(error_v1<TOL,"1-velocity error too large.")
471 self.failUnless(error_v2<TOL,"2-velocity error too large.")
472
473 def test_StokesProblemCartesian_GMRES_P_0(self):
474 ETA=1.
475 P1=0.
476
477 x=self.domain.getX()
478 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.]
479 x=self.domain.getX()
480 mask=whereZero(x[0]) * [1.,1.,1.] \
481 +whereZero(x[0]-1) * [1.,1.,1.] \
482 +whereZero(x[1]) * [1.,1.,1.] \
483 +whereZero(x[1]-1) * [1.,1.,1.] \
484 +whereZero(x[2]) * [1.,1.,0.] \
485 +whereZero(x[2]-1) * [1.,1.,1.]
486
487
488 sp=StokesProblemCartesian(self.domain)
489
490 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
491 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
492 p0=Scalar(P1,ReducedSolution(self.domain))
493 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)
494
495 error_v0=Lsup(u[0]-u0[0])
496 error_v1=Lsup(u[1]-u0[1])
497 error_v2=Lsup(u[2]-u0[2])/0.25**2
498 zz=P1*x[0]*x[1]*x[2]-p
499 error_p=Lsup(zz-integrate(zz))
500 # print error_p, error_v0,error_v1,error_v2
501 self.failUnless(error_p<TOL,"pressure error too large.")
502 self.failUnless(error_v0<TOL,"0-velocity error too large.")
503 self.failUnless(error_v1<TOL,"1-velocity error too large.")
504 self.failUnless(error_v2<TOL,"2-velocity error too large.")
505 def test_StokesProblemCartesian_GMRES_P_small(self):
506 ETA=1.
507 P1=1.
508
509 x=self.domain.getX()
510 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.]
511 mask=whereZero(x[0]) * [1.,1.,1.] \
512 +whereZero(x[0]-1) * [1.,1.,1.] \
513 +whereZero(x[1]) * [1.,1.,1.] \
514 +whereZero(x[1]-1) * [1.,1.,1.] \
515 +whereZero(x[2]) * [1.,1.,0.] \
516 +whereZero(x[2]-1) * [1.,1.,1.]
517
518
519 sp=StokesProblemCartesian(self.domain)
520
521 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
522 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
523 p0=Scalar(P1,ReducedSolution(self.domain))
524 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
525
526 error_v0=Lsup(u[0]-u0[0])
527 error_v1=Lsup(u[1]-u0[1])
528 error_v2=Lsup(u[2]-u0[2])/0.25**2
529 zz=P1*x[0]*x[1]*x[2]-p
530 error_p=Lsup(zz-integrate(zz))/P1
531 # print error_p, error_v0,error_v1,error_v2
532 self.failUnless(error_p<TOL,"pressure error too large.")
533 self.failUnless(error_v0<TOL,"0-velocity error too large.")
534 self.failUnless(error_v1<TOL,"1-velocity error too large.")
535 self.failUnless(error_v2<TOL,"2-velocity error too large.")
536 def test_StokesProblemCartesian_GMRES_P_large(self):
537 ETA=1.
538 P1=1000.
539
540 x=self.domain.getX()
541 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.]
542 mask=whereZero(x[0]) * [1.,1.,1.] \
543 +whereZero(x[0]-1) * [1.,1.,1.] \
544 +whereZero(x[1]) * [1.,1.,1.] \
545 +whereZero(x[1]-1) * [1.,1.,1.] \
546 +whereZero(x[2]) * [1.,1.,0.] \
547 +whereZero(x[2]-1) * [1.,1.,1.]
548
549
550 sp=StokesProblemCartesian(self.domain)
551
552 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
553 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
554 p0=Scalar(P1,ReducedSolution(self.domain))
555 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES")
556
557 error_v0=Lsup(u[0]-u0[0])
558 error_v1=Lsup(u[1]-u0[1])
559 error_v2=Lsup(u[2]-u0[2])/0.25**2
560 zz=P1*x[0]*x[1]*x[2]-p
561 error_p=Lsup(zz-integrate(zz))/P1
562 # print error_p, error_v0,error_v1,error_v2
563 self.failUnless(error_p<TOL,"pressure error too large.")
564 self.failUnless(error_v0<TOL,"0-velocity error too large.")
565 self.failUnless(error_v1<TOL,"1-velocity error too large.")
566 self.failUnless(error_v2<TOL,"2-velocity error too large.")
567
568 def test_StokesProblemCartesian_MINRES_P_0(self):
569 ETA=1.
570 P1=0.
571
572 x=self.domain.getX()
573 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.]
574 x=self.domain.getX()
575 mask=whereZero(x[0]) * [1.,1.,1.] \
576 +whereZero(x[0]-1) * [1.,1.,1.] \
577 +whereZero(x[1]) * [1.,1.,1.] \
578 +whereZero(x[1]-1) * [1.,1.,1.] \
579 +whereZero(x[2]) * [1.,1.,0.] \
580 +whereZero(x[2]-1) * [1.,1.,1.]
581
582
583 sp=StokesProblemCartesian(self.domain)
584
585 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
586 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
587 p0=Scalar(P1,ReducedSolution(self.domain))
588 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
589
590 error_v0=Lsup(u[0]-u0[0])
591 error_v1=Lsup(u[1]-u0[1])
592 error_v2=Lsup(u[2]-u0[2])/0.25**2
593 zz=P1*x[0]*x[1]*x[2]-p
594 error_p=Lsup(zz-integrate(zz))
595 # print error_p, error_v0,error_v1,error_v2
596 self.failUnless(error_p<TOL,"pressure error too large.")
597 self.failUnless(error_v0<TOL,"0-velocity error too large.")
598 self.failUnless(error_v1<TOL,"1-velocity error too large.")
599 self.failUnless(error_v2<TOL,"2-velocity error too large.")
600
601 def test_StokesProblemCartesian_MINRES_P_small(self):
602 ETA=1.
603 P1=1.
604
605 x=self.domain.getX()
606 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.]
607 mask=whereZero(x[0]) * [1.,1.,1.] \
608 +whereZero(x[0]-1) * [1.,1.,1.] \
609 +whereZero(x[1]) * [1.,1.,1.] \
610 +whereZero(x[1]-1) * [1.,1.,1.] \
611 +whereZero(x[2]) * [1.,1.,0.] \
612 +whereZero(x[2]-1) * [1.,1.,1.]
613
614
615 sp=StokesProblemCartesian(self.domain)
616
617 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
618 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
619 p0=Scalar(P1,ReducedSolution(self.domain))
620 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
621
622 error_v0=Lsup(u[0]-u0[0])
623 error_v1=Lsup(u[1]-u0[1])
624 error_v2=Lsup(u[2]-u0[2])/0.25**2
625 zz=P1*x[0]*x[1]*x[2]-p
626 error_p=Lsup(zz-integrate(zz))/P1
627 # print error_p, error_v0,error_v1,error_v2
628 self.failUnless(error_p<TOL,"pressure error too large.")
629 self.failUnless(error_v0<TOL,"0-velocity error too large.")
630 self.failUnless(error_v1<TOL,"1-velocity error too large.")
631 self.failUnless(error_v2<TOL,"2-velocity error too large.")
632
633 # def test_StokesProblemCartesian_MINRES_P_large(self):
634 # ETA=1.
635 # P1=1000.
636
637 # x=self.domain.getX()
638 # 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.]
639 # mask=whereZero(x[0]) * [1.,1.,1.] \
640 # +whereZero(x[0]-1) * [1.,1.,1.] \
641 # +whereZero(x[1]) * [1.,1.,1.] \
642 # +whereZero(x[1]-1) * [1.,1.,1.] \
643 # +whereZero(x[2]) * [1.,1.,0.] \
644 # +whereZero(x[2]-1) * [1.,1.,1.]
645
646
647 # sp=StokesProblemCartesian(self.domain)
648
649 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
650 # u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
651 # p0=Scalar(P1,ReducedSolution(self.domain))
652 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="MINRES")
653
654 # error_v0=Lsup(u[0]-u0[0])
655 # error_v1=Lsup(u[1]-u0[1])
656 # error_v2=Lsup(u[2]-u0[2])/0.25**2
657 # zz=P1*x[0]*x[1]*x[2]-p
658 # error_p=Lsup(zz-integrate(zz))/P1
659 # print error_p, error_v0,error_v1,error_v2
660 # self.failUnless(error_p<TOL,"pressure error too large.")
661 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
662 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
663 # self.failUnless(error_v2<TOL,"2-velocity error too large.")
664
665 # def test_StokesProblemCartesian_TFQMR_P_0(self):
666 # ETA=1.
667 # P1=0.
668
669 # x=self.domain.getX()
670 # 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.]
671 # x=self.domain.getX()
672 # mask=whereZero(x[0]) * [1.,1.,1.] \
673 # +whereZero(x[0]-1) * [1.,1.,1.] \
674 # +whereZero(x[1]) * [1.,1.,1.] \
675 # +whereZero(x[1]-1) * [1.,1.,1.] \
676 # +whereZero(x[2]) * [1.,1.,0.] \
677 # +whereZero(x[2]-1) * [1.,1.,1.]
678
679
680 # sp=StokesProblemCartesian(self.domain)
681
682 # sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
683 # u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
684 # p0=Scalar(P1,ReducedSolution(self.domain))
685 # u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
686
687 # error_v0=Lsup(u[0]-u0[0])
688 # error_v1=Lsup(u[1]-u0[1])
689 # error_v2=Lsup(u[2]-u0[2])/0.25**2
690 # zz=P1*x[0]*x[1]*x[2]-p
691 # error_p=Lsup(zz-integrate(zz))
692 # print error_p, error_v0,error_v1,error_v2
693 # self.failUnless(error_p<TOL,"pressure error too large.")
694 # self.failUnless(error_v0<TOL,"0-velocity error too large.")
695 # self.failUnless(error_v1<TOL,"1-velocity error too large.")
696 # self.failUnless(error_v2<TOL,"2-velocity error too large.")
697
698 def test_StokesProblemCartesian_TFQMR_P_small(self):
699 ETA=1.
700 P1=1.
701
702 x=self.domain.getX()
703 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.]
704 mask=whereZero(x[0]) * [1.,1.,1.] \
705 +whereZero(x[0]-1) * [1.,1.,1.] \
706 +whereZero(x[1]) * [1.,1.,1.] \
707 +whereZero(x[1]-1) * [1.,1.,1.] \
708 +whereZero(x[2]) * [1.,1.,0.] \
709 +whereZero(x[2]-1) * [1.,1.,1.]
710
711
712 sp=StokesProblemCartesian(self.domain)
713
714 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
715 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
716 p0=Scalar(P1,ReducedSolution(self.domain))
717 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
718
719 error_v0=Lsup(u[0]-u0[0])
720 error_v1=Lsup(u[1]-u0[1])
721 error_v2=Lsup(u[2]-u0[2])/0.25**2
722 zz=P1*x[0]*x[1]*x[2]-p
723 error_p=Lsup(zz-integrate(zz))/P1
724 # print error_p, error_v0,error_v1,error_v2
725 self.failUnless(error_p<TOL,"pressure error too large.")
726 self.failUnless(error_v0<TOL,"0-velocity error too large.")
727 self.failUnless(error_v1<TOL,"1-velocity error too large.")
728 self.failUnless(error_v2<TOL,"2-velocity error too large.")
729
730 def test_StokesProblemCartesian_TFQMR_P_large(self):
731 ETA=1.
732 P1=1000.
733
734 x=self.domain.getX()
735 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.]
736 mask=whereZero(x[0]) * [1.,1.,1.] \
737 +whereZero(x[0]-1) * [1.,1.,1.] \
738 +whereZero(x[1]) * [1.,1.,1.] \
739 +whereZero(x[1]-1) * [1.,1.,1.] \
740 +whereZero(x[2]) * [1.,1.,0.] \
741 +whereZero(x[2]-1) * [1.,1.,1.]
742
743
744 sp=StokesProblemCartesian(self.domain)
745
746 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
747 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
748 p0=Scalar(P1,ReducedSolution(self.domain))
749 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="TFQMR")
750
751 error_v0=Lsup(u[0]-u0[0])
752 error_v1=Lsup(u[1]-u0[1])
753 error_v2=Lsup(u[2]-u0[2])/0.25**2
754 zz=P1*x[0]*x[1]*x[2]-p
755 error_p=Lsup(zz-integrate(zz))/P1
756 # print error_p, error_v0,error_v1,error_v2
757 self.failUnless(error_p<TOL,"pressure error too large.")
758 self.failUnless(error_v0<TOL,"0-velocity error too large.")
759 self.failUnless(error_v1<TOL,"1-velocity error too large.")
760 self.failUnless(error_v2<TOL,"2-velocity error too large.")
761
762
763 if __name__ == '__main__':
764 suite = unittest.TestSuite()
765 suite.addTest(unittest.makeSuite(Test_Simple2DModels))
766 # suite.addTest(unittest.makeSuite(Test_Simple3DModels))
767 s=unittest.TextTestRunner(verbosity=2).run(suite)
768 if not s.wasSuccessful(): sys.exit(1)
769

  ViewVC Help
Powered by ViewVC 1.1.26