/[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 1562 - (hide annotations)
Wed May 21 13:04:40 2008 UTC (11 years, 5 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 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 gross 1562 VERBOSE=False or True
35 gross 1414
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 1514 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)
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 artak 1518
148     # print error_p, error_v0,error_v1
149 gross 1471 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 artak 1518 def test_StokesProblemCartesian_MINRES_P_0(self):
208     ETA=1.
209     P1=0.
210 gross 1471
211 artak 1518 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 gross 1414 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 gross 1471 def test_StokesProblemCartesian_PCG_P_0(self):
378 gross 1414 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 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
398 gross 1414
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 artak 1518
410 gross 1471 def test_StokesProblemCartesian_PCG_P_small(self):
411 gross 1414 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 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
430 gross 1414
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 gross 1471 def test_StokesProblemCartesian_PCG_P_large(self):
442 gross 1414 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 gross 1471 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG")
461 gross 1414
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 gross 1471 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 artak 1514 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=20)
494 gross 1471
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 artak 1518 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 gross 1414 if __name__ == '__main__':
764     suite = unittest.TestSuite()
765     suite.addTest(unittest.makeSuite(Test_Simple2DModels))
766 gross 1562 # suite.addTest(unittest.makeSuite(Test_Simple3DModels))
767 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
768     if not s.wasSuccessful(): sys.exit(1)
769    

  ViewVC Help
Powered by ViewVC 1.1.26