# Diff of /trunk/finley/test/python/run_models.py

revision 1414 by gross, Thu Feb 14 10:01:43 2008 UTC revision 1562 by gross, Wed May 21 13:04:40 2008 UTC
# Line 31  except KeyError: Line 31  except KeyError:
31
32  NE=6  NE=6
33  TOL=1.e-5  TOL=1.e-5
34  VERBOSE=False  VERBOSE=False or True
35
36  from esys.escript import *  from esys.escript import *
37  from esys.escript.models import StokesProblemCartesian  from esys.escript.models import StokesProblemCartesian
# Line 41  class Test_Simple2DModels(unittest.TestC Line 41  class Test_Simple2DModels(unittest.TestC
41         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)         self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
42     def tearDown(self):     def tearDown(self):
43         del self.domain         del self.domain
44     def test_StokesProblemCartesian_P_0(self):     def test_StokesProblemCartesian_PCG_P_0(self):
45         ETA=1.         ETA=1.
46         P1=0.         P1=0.
47
# Line 57  class Test_Simple2DModels(unittest.TestC Line 57  class Test_Simple2DModels(unittest.TestC
58         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
59         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
60         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
63         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 68  class Test_Simple2DModels(unittest.TestC Line 68  class Test_Simple2DModels(unittest.TestC
68         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<TOL,"0-velocity error too large.")
69         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<TOL,"1-velocity error too large.")
70
71     def test_StokesProblemCartesian_P_small(self):     def test_StokesProblemCartesian_PCG_P_small(self):
72         ETA=1.         ETA=1.
73         P1=1.         P1=1.
74
# Line 84  class Test_Simple2DModels(unittest.TestC Line 84  class Test_Simple2DModels(unittest.TestC
85         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
86         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
87         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
90         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 95  class Test_Simple2DModels(unittest.TestC Line 95  class Test_Simple2DModels(unittest.TestC
95         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<TOL,"0-velocity error too large.")
96         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<TOL,"1-velocity error too large.")
97
98     def test_StokesProblemCartesian_P_large(self):     def test_StokesProblemCartesian_PCG_P_large(self):
99         ETA=1.         ETA=1.
100         P1=1000.         P1=1000.
101
# Line 111  class Test_Simple2DModels(unittest.TestC Line 111  class Test_Simple2DModels(unittest.TestC
112         u0=(1-x[0])*x[0]*[0.,1.]         u0=(1-x[0])*x[0]*[0.,1.]
113         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
114         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
117         error_v1=Lsup(u[1]-u0[1])/0.25         error_v1=Lsup(u[1]-u0[1])/0.25
# Line 122  class Test_Simple2DModels(unittest.TestC Line 122  class Test_Simple2DModels(unittest.TestC
122         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<TOL,"0-velocity error too large.")
123         self.failUnless(error_v1<TOL,"1-velocity error too large.")         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
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
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
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
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
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
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
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
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
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):  class Test_Simple3DModels(unittest.TestCase):
373     def setUp(self):     def setUp(self):
374         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)         self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
375     def tearDown(self):     def tearDown(self):
376         del self.domain         del self.domain
377     def test_StokesProblemCartesian_P_0(self):     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
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
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
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.         ETA=1.
475         P1=0.         P1=0.
476
# Line 148  class Test_Simple3DModels(unittest.TestC Line 490  class Test_Simple3DModels(unittest.TestC
491         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
492         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
493         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
496         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 160  class Test_Simple3DModels(unittest.TestC Line 502  class Test_Simple3DModels(unittest.TestC
502         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<TOL,"0-velocity error too large.")
503         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<TOL,"1-velocity error too large.")
504         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<TOL,"2-velocity error too large.")
505     def test_StokesProblemCartesian_P_small(self):     def test_StokesProblemCartesian_GMRES_P_small(self):
506         ETA=1.         ETA=1.
507         P1=1.         P1=1.
508
# Line 179  class Test_Simple3DModels(unittest.TestC Line 521  class Test_Simple3DModels(unittest.TestC
522         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
523         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
524         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
527         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 191  class Test_Simple3DModels(unittest.TestC Line 533  class Test_Simple3DModels(unittest.TestC
533         self.failUnless(error_v0<TOL,"0-velocity error too large.")         self.failUnless(error_v0<TOL,"0-velocity error too large.")
534         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<TOL,"1-velocity error too large.")
535         self.failUnless(error_v2<TOL,"2-velocity error too large.")         self.failUnless(error_v2<TOL,"2-velocity error too large.")
536     def test_StokesProblemCartesian_P_large(self):     def test_StokesProblemCartesian_GMRES_P_large(self):
537         ETA=1.         ETA=1.
538         P1=1000.         P1=1000.
539
# Line 210  class Test_Simple3DModels(unittest.TestC Line 552  class Test_Simple3DModels(unittest.TestC
553         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]         u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
554         p0=Scalar(P1,ReducedSolution(self.domain))         p0=Scalar(P1,ReducedSolution(self.domain))
555         u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100)         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])         error_v0=Lsup(u[0]-u0[0])
558         error_v1=Lsup(u[1]-u0[1])         error_v1=Lsup(u[1]-u0[1])
# Line 223  class Test_Simple3DModels(unittest.TestC Line 565  class Test_Simple3DModels(unittest.TestC
565         self.failUnless(error_v1<TOL,"1-velocity error too large.")         self.failUnless(error_v1<TOL,"1-velocity error too large.")
566         self.failUnless(error_v2<TOL,"2-velocity error too large.")         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
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
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
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
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
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
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__':  if __name__ == '__main__':
764     suite = unittest.TestSuite()     suite = unittest.TestSuite()