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

  ViewVC Help
Powered by ViewVC 1.1.26