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