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