1 |
####################################################### |
2 |
# |
3 |
# Copyright 2003-2007 by ACceSS MNRF |
4 |
# Copyright 2007 by University of Queensland |
5 |
# |
6 |
# http://esscc.uq.edu.au |
7 |
# Primary Business: Queensland, Australia |
8 |
# Licensed under the Open Software License version 3.0 |
9 |
# http://www.opensource.org/licenses/osl-3.0.php |
10 |
# |
11 |
####################################################### |
12 |
# |
13 |
|
14 |
__copyright__=""" Copyright (c) 2006 by ACcESS MNRF |
15 |
http://www.access.edu.au |
16 |
Primary Business: Queensland, Australia""" |
17 |
__license__="""Licensed under the Open Software License version 3.0 |
18 |
http://www.opensource.org/licenses/osl-3.0.php""" |
19 |
import unittest |
20 |
import tempfile |
21 |
|
22 |
from esys.escript import * |
23 |
from esys.finley import Rectangle |
24 |
import sys |
25 |
import os |
26 |
try: |
27 |
FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR'] |
28 |
except KeyError: |
29 |
FINLEY_WORKDIR='.' |
30 |
|
31 |
|
32 |
NE=6 |
33 |
TOL=1.e-5 |
34 |
VERBOSE=False |
35 |
|
36 |
from esys.escript import * |
37 |
from esys.escript.models import StokesProblemCartesian |
38 |
from esys.finley import Rectangle, Brick |
39 |
class Test_Simple2DModels(unittest.TestCase): |
40 |
def setUp(self): |
41 |
self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True) |
42 |
def tearDown(self): |
43 |
del self.domain |
44 |
def test_StokesProblemCartesian_PCG_P_0(self): |
45 |
ETA=1. |
46 |
P1=0. |
47 |
|
48 |
x=self.domain.getX() |
49 |
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] |
50 |
mask=whereZero(x[0]) * [1.,1.] \ |
51 |
+whereZero(x[0]-1) * [1.,1.] \ |
52 |
+whereZero(x[1]) * [1.,0.] \ |
53 |
+whereZero(x[1]-1) * [1.,1.] |
54 |
|
55 |
sp=StokesProblemCartesian(self.domain) |
56 |
|
57 |
sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) |
58 |
u0=(1-x[0])*x[0]*[0.,1.] |
59 |
p0=Scalar(P1,ReducedSolution(self.domain)) |
60 |
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG") |
61 |
|
62 |
error_v0=Lsup(u[0]-u0[0]) |
63 |
error_v1=Lsup(u[1]-u0[1])/0.25 |
64 |
zz=P1*x[0]*x[1]-p |
65 |
error_p=Lsup(zz-integrate(zz)) |
66 |
# print error_p, error_v0,error_v1 |
67 |
self.failUnless(error_p<TOL,"pressure error too large.") |
68 |
self.failUnless(error_v0<TOL,"0-velocity error too large.") |
69 |
self.failUnless(error_v1<TOL,"1-velocity error too large.") |
70 |
|
71 |
def test_StokesProblemCartesian_PCG_P_small(self): |
72 |
ETA=1. |
73 |
P1=1. |
74 |
|
75 |
x=self.domain.getX() |
76 |
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] |
77 |
mask=whereZero(x[0]) * [1.,1.] \ |
78 |
+whereZero(x[0]-1) * [1.,1.] \ |
79 |
+whereZero(x[1]) * [1.,0.] \ |
80 |
+whereZero(x[1]-1) * [1.,1.] |
81 |
|
82 |
sp=StokesProblemCartesian(self.domain) |
83 |
|
84 |
sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) |
85 |
u0=(1-x[0])*x[0]*[0.,1.] |
86 |
p0=Scalar(P1,ReducedSolution(self.domain)) |
87 |
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG") |
88 |
|
89 |
error_v0=Lsup(u[0]-u0[0]) |
90 |
error_v1=Lsup(u[1]-u0[1])/0.25 |
91 |
zz=P1*x[0]*x[1]-p |
92 |
error_p=Lsup(zz-integrate(zz)) |
93 |
# print error_p, error_v0,error_v1 |
94 |
self.failUnless(error_p<TOL,"pressure error too large.") |
95 |
self.failUnless(error_v0<TOL,"0-velocity error too large.") |
96 |
self.failUnless(error_v1<TOL,"1-velocity error too large.") |
97 |
|
98 |
def test_StokesProblemCartesian_PCG_P_large(self): |
99 |
ETA=1. |
100 |
P1=1000. |
101 |
|
102 |
x=self.domain.getX() |
103 |
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] |
104 |
mask=whereZero(x[0]) * [1.,1.] \ |
105 |
+whereZero(x[0]-1) * [1.,1.] \ |
106 |
+whereZero(x[1]) * [1.,0.] \ |
107 |
+whereZero(x[1]-1) * [1.,1.] |
108 |
|
109 |
sp=StokesProblemCartesian(self.domain) |
110 |
|
111 |
sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) |
112 |
u0=(1-x[0])*x[0]*[0.,1.] |
113 |
p0=Scalar(P1,ReducedSolution(self.domain)) |
114 |
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="PCG") |
115 |
|
116 |
error_v0=Lsup(u[0]-u0[0]) |
117 |
error_v1=Lsup(u[1]-u0[1])/0.25 |
118 |
zz=P1*x[0]*x[1]-p |
119 |
error_p=Lsup(zz-integrate(zz))/P1 |
120 |
# print error_p, error_v0,error_v1 |
121 |
self.failUnless(error_p<TOL,"pressure error too large.") |
122 |
self.failUnless(error_v0<TOL,"0-velocity error too large.") |
123 |
self.failUnless(error_v1<TOL,"1-velocity error too large.") |
124 |
|
125 |
def test_StokesProblemCartesian_GMRES_P_0(self): |
126 |
ETA=1. |
127 |
P1=0. |
128 |
|
129 |
x=self.domain.getX() |
130 |
F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.] |
131 |
mask=whereZero(x[0]) * [1.,1.] \ |
132 |
+whereZero(x[0]-1) * [1.,1.] \ |
133 |
+whereZero(x[1]) * [1.,0.] \ |
134 |
+whereZero(x[1]-1) * [1.,1.] |
135 |
|
136 |
sp=StokesProblemCartesian(self.domain) |
137 |
|
138 |
sp.initialize(f=F,fixed_u_mask=mask,eta=ETA) |
139 |
u0=(1-x[0])*x[0]*[0.,1.] |
140 |
p0=Scalar(P1,ReducedSolution(self.domain)) |
141 |
u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,solver="GMRES",iter_restart=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 |
|