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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2100 - (hide 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 ksteube 1809
2     ########################################################
3 gross 1414 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # 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 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 ksteube 1809 __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 gross 1414 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 gross 1414 import unittest
23     import tempfile
24    
25     from esys.escript import *
26     from esys.finley import Rectangle
27 gross 2100 from esys.escript.models import DarcyFlow
28 gross 1414 import sys
29     import os
30     try:
31     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
32     except KeyError:
33     FINLEY_WORKDIR='.'
34    
35    
36 gross 2100 VERBOSE=False # or True
37 gross 1414
38     from esys.escript import *
39     from esys.escript.models import StokesProblemCartesian
40     from esys.finley import Rectangle, Brick
41 gross 2100
42     #====================================================================================================================
43     class Test_StokesProblemCartesian2D(unittest.TestCase):
44 gross 1414 def setUp(self):
45 gross 2100 NE=6
46     self.TOL=1.e-5
47 gross 1414 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
48     def tearDown(self):
49     del self.domain
50 gross 2100 def test_PCG_P_0(self):
51 gross 1414 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 gross 2100 sp.setTolerance(self.TOL)
67     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
68 gross 1414
69     error_v0=Lsup(u[0]-u0[0])
70     error_v1=Lsup(u[1]-u0[1])/0.25
71 gross 2100 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 gross 1414
76 gross 2100 def test_PCG_P_small(self):
77 gross 1414 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 gross 2100 sp.setTolerance(self.TOL)
93     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE or True,max_iter=100,useUzawa=True)
94 gross 1414 error_v0=Lsup(u[0]-u0[0])
95     error_v1=Lsup(u[1]-u0[1])/0.25
96 gross 2100 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 gross 1414
101 gross 2100 def test_PCG_P_large(self):
102 gross 1414 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 gross 2100 sp.setTolerance(self.TOL)
118     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
119 gross 1414
120     error_v0=Lsup(u[0]-u0[0])
121     error_v1=Lsup(u[1]-u0[1])/0.25
122 gross 2100 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 gross 1414
127 gross 2100 def test_GMRES_P_0(self):
128 gross 1471 ETA=1.
129     P1=0.
130 gross 1414
131 gross 1471 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 gross 2100 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 gross 1471
147     error_v0=Lsup(u[0]-u0[0])
148     error_v1=Lsup(u[1]-u0[1])/0.25
149 gross 2100 zz=P1*x[0]*x[1]+p
150 gross 1471 error_p=Lsup(zz-integrate(zz))
151 gross 2100 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 artak 1518
155 gross 2100 def test_GMRES_P_small(self):
156 gross 1471 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 gross 2100 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 gross 1471
175     error_v0=Lsup(u[0]-u0[0])
176     error_v1=Lsup(u[1]-u0[1])/0.25
177 gross 2100 zz=P1*x[0]*x[1]+p
178 gross 1471 error_p=Lsup(zz-integrate(zz))
179     # print error_p, error_v0,error_v1
180 gross 2100 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 gross 1471
184 gross 2100 def test_GMRES_P_large(self):
185 gross 1471 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 gross 2100 sp.setTolerance(self.TOL)
201     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
202 gross 1471
203     error_v0=Lsup(u[0]-u0[0])
204     error_v1=Lsup(u[1]-u0[1])/0.25
205 gross 2100 zz=P1*x[0]*x[1]+p
206 gross 1471 error_p=Lsup(zz-integrate(zz))/P1
207 gross 2100 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 gross 1414 def setUp(self):
213 gross 2100 NE=6
214     self.TOL=1.e-4
215 gross 1414 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
216     def tearDown(self):
217     del self.domain
218 gross 2100 def test_PCG_P_0(self):
219 gross 1414 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 gross 2100 sp.setTolerance(1.e-8)
239     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
240 gross 1414
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
245 gross 1414 error_p=Lsup(zz-integrate(zz))
246     # print error_p, error_v0,error_v1,error_v2
247 gross 2100 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 artak 1518
252 gross 2100 def test_PCG_P_small(self):
253 gross 1414 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 gross 2100 sp.setTolerance(1.e-8)
272     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
273 gross 1414
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
278 gross 1414 error_p=Lsup(zz-integrate(zz))/P1
279 gross 2100 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 gross 1414 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 gross 2100 sp.setTolerance(1.e-8)
303     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
304 gross 1414
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
309 gross 1414 error_p=Lsup(zz-integrate(zz))/P1
310     # print error_p, error_v0,error_v1,error_v2
311 gross 2100 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 gross 1414
316 gross 2100 def test_GMRES_P_0(self):
317 gross 1471 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 gross 2100 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 gross 1471
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
343 gross 1471 error_p=Lsup(zz-integrate(zz))
344     # print error_p, error_v0,error_v1,error_v2
345 gross 2100 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 gross 1471 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 gross 2100 sp.setTolerance(1.e-8)
369     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
370 gross 1471
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
375 gross 1471 error_p=Lsup(zz-integrate(zz))/P1
376     # print error_p, error_v0,error_v1,error_v2
377 gross 2100 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 gross 1471 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 gross 2100 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 gross 1471
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 gross 2100 zz=P1*x[0]*x[1]*x[2]+p
409 gross 1471 error_p=Lsup(zz-integrate(zz))/P1
410     # print error_p, error_v0,error_v1,error_v2
411 gross 2100 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 gross 1471
461 gross 2100 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 artak 1518
486 gross 2100 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 artak 1518
502 gross 2100 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 artak 1518
518 gross 2100 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 artak 1518
534 gross 2100 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 artak 1518
550 gross 2100 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 artak 1518
566 gross 2100 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 artak 1518
582 gross 2100 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 artak 1518
598 gross 2100 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 artak 1518
614 gross 2100 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 artak 1518
630 gross 2100 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 artak 1518
646 gross 2100 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 artak 1518
662 gross 2100 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 artak 1518
681 gross 2100
682    
683 gross 1414 if __name__ == '__main__':
684     suite = unittest.TestSuite()
685 gross 2100 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 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
690     if not s.wasSuccessful(): sys.exit(1)
691    

  ViewVC Help
Powered by ViewVC 1.1.26