/[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 2156 - (show annotations)
Mon Dec 15 05:09:02 2008 UTC (11 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 28365 byte(s)
some modifications to the iterative solver to make them easier to use. 
There are also improved versions of the Darcy flux solver and the incompressible 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=1.e-5
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=1.e-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 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
466 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
467 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
468
469 def testConstF_FixedBottom_mediumK(self):
470 k=1.
471 mp=self.getScalarMask(include_bottom=True)
472 mv=self.getVectorMask(include_bottom=False)
473 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
474 p=p_ref*mp
475 u=u_ref*mv
476 df=DarcyFlow(self.dom)
477 df.setValue(g=f,
478 location_of_fixed_pressure=mp,
479 location_of_fixed_flux=mv,
480 permeability=Scalar(k,Function(self.dom)))
481 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
482 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
483 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
484
485 def testConstF_FixedBottom_largeK(self):
486 k=1.e10
487 mp=self.getScalarMask(include_bottom=True)
488 mv=self.getVectorMask(include_bottom=False)
489 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
490 p=p_ref*mp
491 u=u_ref*mv
492 df=DarcyFlow(self.dom)
493 df.setValue(g=f,
494 location_of_fixed_pressure=mp,
495 location_of_fixed_flux=mv,
496 permeability=Scalar(k,Function(self.dom)))
497 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL*self.TOL)
498 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
499 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
500
501 def testVarioF_FixedBottom_smallK(self):
502 k=1.e-10
503 mp=self.getScalarMask(include_bottom=True)
504 mv=self.getVectorMask(include_bottom=False)
505 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
506 p=p_ref*mp
507 u=u_ref*mv
508 df=DarcyFlow(self.dom)
509 df.setValue(g=f,
510 location_of_fixed_pressure=mp,
511 location_of_fixed_flux=mv,
512 permeability=Scalar(k,Function(self.dom)))
513 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
514 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
515 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
516
517 def testVarioF_FixedBottom_mediumK(self):
518 k=1.
519 mp=self.getScalarMask(include_bottom=True)
520 mv=self.getVectorMask(include_bottom=False)
521 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
522 p=p_ref*mp
523 u=u_ref*mv
524 df=DarcyFlow(self.dom)
525 df.setValue(g=f,
526 location_of_fixed_pressure=mp,
527 location_of_fixed_flux=mv,
528 permeability=Scalar(k,Function(self.dom)))
529 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
530 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
531 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
532
533 def testVarioF_FixedBottom_largeK(self):
534 k=1.e10
535 mp=self.getScalarMask(include_bottom=True)
536 mv=self.getVectorMask(include_bottom=False)
537 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
538 p=p_ref*mp
539 u=u_ref*mv
540 df=DarcyFlow(self.dom)
541 df.setValue(g=f,
542 location_of_fixed_pressure=mp,
543 location_of_fixed_flux=mv,
544 permeability=Scalar(k,Function(self.dom)))
545 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
546 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
547 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
548
549 def testConstF_FreeBottom_smallK(self):
550 k=1.e-10
551 mp=self.getScalarMask(include_bottom=False)
552 mv=self.getVectorMask(include_bottom=True)
553 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
554 p=p_ref*mp
555 u=u_ref*mv
556 df=DarcyFlow(self.dom)
557 df.setValue(g=f,
558 location_of_fixed_pressure=mp,
559 location_of_fixed_flux=mv,
560 permeability=Scalar(k,Function(self.dom)))
561 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
562 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
563 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
564
565 def testConstF_FreeBottom_mediumK(self):
566 k=1.
567 mp=self.getScalarMask(include_bottom=False)
568 mv=self.getVectorMask(include_bottom=True)
569 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
570 p=p_ref*mp
571 u=u_ref*mv
572 df=DarcyFlow(self.dom)
573 df.setValue(g=f,
574 location_of_fixed_pressure=mp,
575 location_of_fixed_flux=mv,
576 permeability=Scalar(k,Function(self.dom)))
577 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
578 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
579 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
580
581 def testConstF_FreeBottom_largeK(self):
582 k=1.e10
583 mp=self.getScalarMask(include_bottom=False)
584 mv=self.getVectorMask(include_bottom=True)
585 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
586 p=p_ref*mp
587 u=u_ref*mv
588 df=DarcyFlow(self.dom)
589 df.setValue(g=f,
590 location_of_fixed_pressure=mp,
591 location_of_fixed_flux=mv,
592 permeability=Scalar(k,Function(self.dom)))
593 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
594 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
595 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
596
597 def testVarioF_FreeBottom_smallK(self):
598 k=1.e-10
599 mp=self.getScalarMask(include_bottom=False)
600 mv=self.getVectorMask(include_bottom=True)
601 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
602 p=p_ref*mp
603 u=u_ref*mv
604 df=DarcyFlow(self.dom)
605 df.setValue(g=f,
606 location_of_fixed_pressure=mp,
607 location_of_fixed_flux=mv,
608 permeability=Scalar(k,Function(self.dom)))
609 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
610 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
611 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
612
613 def testVarioF_FreeBottom_mediumK(self):
614 k=1.
615 mp=self.getScalarMask(include_bottom=False)
616 mv=self.getVectorMask(include_bottom=True)
617 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
618 p=p_ref*mp
619 u=u_ref*mv
620 df=DarcyFlow(self.dom)
621 df.setValue(g=f,
622 location_of_fixed_pressure=mp,
623 location_of_fixed_flux=mv,
624 permeability=Scalar(k,Function(self.dom)))
625 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
626 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
627 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
628
629 def testVarioF_FreeBottom_largeK(self):
630 k=1.e10
631 mp=self.getScalarMask(include_bottom=False)
632 mv=self.getVectorMask(include_bottom=True)
633 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
634 p=p_ref*mp
635 u=u_ref*mv
636 df=DarcyFlow(self.dom)
637 df.setValue(g=f,
638 location_of_fixed_pressure=mp,
639 location_of_fixed_flux=mv,
640 permeability=Scalar(k,Function(self.dom)))
641 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
642 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
643 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
644
645 class Test_Darcy2D(Test_Darcy):
646 TOL=1e-5
647 WIDTH=1.
648 def setUp(self):
649 NE=60 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
650 self.dom = Rectangle(NE,NE)
651 self.rescaleDomain()
652 def tearDown(self):
653 del self.dom
654 class Test_Darcy3D(Test_Darcy):
655 TOL=1e-4
656 WIDTH=1.
657 def setUp(self):
658 NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
659 self.dom = Brick(NE,NE,NE)
660 self.rescaleDomain()
661 def tearDown(self):
662 del self.dom
663
664
665
666 if __name__ == '__main__':
667 suite = unittest.TestSuite()
668 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
669 suite.addTest(unittest.makeSuite(Test_Darcy3D))
670 suite.addTest(unittest.makeSuite(Test_Darcy2D))
671 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
672 s=unittest.TextTestRunner(verbosity=2).run(suite)
673 if not s.wasSuccessful(): sys.exit(1)
674

  ViewVC Help
Powered by ViewVC 1.1.26