/[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 2156 - (hide annotations)
Mon Dec 15 05:09:02 2008 UTC (10 years, 9 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 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 2156 VERBOSE=False #
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 2156 sp.setTolerance(self.TOL*0.2)
93     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,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 gross 2156 saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
98 gross 2100 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 gross 2156 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
101 gross 1414
102 gross 2100 def test_PCG_P_large(self):
103 gross 1414 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 gross 2100 sp.setTolerance(self.TOL)
119     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=True)
120 gross 1414
121     error_v0=Lsup(u[0]-u0[0])
122     error_v1=Lsup(u[1]-u0[1])/0.25
123 gross 2100 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 gross 1414
128 gross 2100 def test_GMRES_P_0(self):
129 gross 1471 ETA=1.
130     P1=0.
131 gross 1414
132 gross 1471 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 gross 2100 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 gross 1471
147     error_v0=Lsup(u[0]-u0[0])
148     error_v1=Lsup(u[1]-u0[1])/0.25
149 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
150 gross 2100 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 artak 1518
154 gross 2100 def test_GMRES_P_small(self):
155 gross 1471 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 gross 2156 sp.setTolerance(self.TOL*0.1)
171 gross 2100 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=20,useUzawa=False)
172 gross 1471
173     error_v0=Lsup(u[0]-u0[0])
174     error_v1=Lsup(u[1]-u0[1])/0.25
175 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
176 gross 2100 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 gross 1471
180 gross 2100 def test_GMRES_P_large(self):
181 gross 1471 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 gross 2100 sp.setTolerance(self.TOL)
197     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False)
198 gross 1471
199     error_v0=Lsup(u[0]-u0[0])
200     error_v1=Lsup(u[1]-u0[1])/0.25
201 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)/P1
202 gross 2100 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 gross 1414 def setUp(self):
208 gross 2100 NE=6
209     self.TOL=1.e-4
210 gross 1414 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
211     def tearDown(self):
212     del self.domain
213 gross 2100 def test_PCG_P_0(self):
214 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
223 gross 1414 +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 gross 2156 sp.setTolerance(self.TOL)
234     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
235 gross 1414
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
240 gross 2100 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 artak 1518
245 gross 2100 def test_PCG_P_small(self):
246 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
254 gross 1414 +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 gross 2156 sp.setTolerance(self.TOL)
265     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
266 gross 1414
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
271 gross 2100 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 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
284 gross 1414 +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 gross 2156 sp.setTolerance(self.TOL)
295     u,p=sp.solve(u0,-p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=True)
296 gross 1414
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
301 gross 2100 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 gross 1414
306 gross 2100 def test_GMRES_P_0(self):
307 gross 1471 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 gross 2156 sp.setTolerance(self.TOL)
327 gross 2100 u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE,max_iter=100,useUzawa=False,iter_restart=20)
328 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
333 gross 1471 # print error_p, error_v0,error_v1,error_v2
334 gross 2100 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 gross 1471 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 gross 2156 sp.setTolerance(self.TOL)
358     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)
359 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
364 gross 2100 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 gross 1471 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
377 gross 1471 +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 gross 2156 sp.setTolerance(self.TOL)
388     u,p=sp.solve(u0,p0,show_details=VERBOSE, verbose=VERBOSE ,max_iter=100,useUzawa=False)
389 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
394 gross 2100 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 gross 1471
444 gross 2100 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 artak 1518
469 gross 2100 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 artak 1518
485 gross 2100 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 gross 2156 v,p=df.solve(u,p,atol=0,rtol=self.TOL, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL*self.TOL)
498 gross 2100 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 artak 1518
501 gross 2100 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 artak 1518
517 gross 2100 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 artak 1518
533 gross 2100 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 artak 1518
549 gross 2100 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 artak 1518
565 gross 2100 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 artak 1518
581 gross 2100 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 artak 1518
597 gross 2100 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 artak 1518
613 gross 2100 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 artak 1518
629 gross 2100 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 artak 1518
645 gross 2100 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 artak 1518
664 gross 2100
665    
666 gross 1414 if __name__ == '__main__':
667     suite = unittest.TestSuite()
668 gross 2100 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 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
673     if not s.wasSuccessful(): sys.exit(1)
674    

  ViewVC Help
Powered by ViewVC 1.1.26