/[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 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 29054 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



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

  ViewVC Help
Powered by ViewVC 1.1.26