/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

Contents of /trunk/finley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2208 - (show annotations)
Mon Jan 12 06:37:07 2009 UTC (11 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 28894 byte(s)
more work on the dary solver 


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

  ViewVC Help
Powered by ViewVC 1.1.26