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

  ViewVC Help
Powered by ViewVC 1.1.26