/[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 2251 - (hide annotations)
Fri Feb 6 06:50:39 2009 UTC (10 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 30729 byte(s)
new compressible flow 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 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
468     v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
469 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
470     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
471 artak 1518
472 gross 2100 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 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
485     v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE ,sub_rtol=self.TOL/200)
486     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
487 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
488 artak 1518
489 gross 2100 def testConstF_FixedBottom_largeK(self):
490     k=1.e10
491     mp=self.getScalarMask(include_bottom=True)
492     mv=self.getVectorMask(include_bottom=False)
493     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
494     p=p_ref*mp
495     u=u_ref*mv
496     df=DarcyFlow(self.dom)
497     df.setValue(g=f,
498     location_of_fixed_pressure=mp,
499     location_of_fixed_flux=mv,
500     permeability=Scalar(k,Function(self.dom)))
501 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
502     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
503 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
504     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
505 artak 1518
506 gross 2100 def testVarioF_FixedBottom_smallK(self):
507     k=1.e-10
508     mp=self.getScalarMask(include_bottom=True)
509     mv=self.getVectorMask(include_bottom=False)
510     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
511     p=p_ref*mp
512     u=u_ref*mv
513     df=DarcyFlow(self.dom)
514     df.setValue(g=f,
515     location_of_fixed_pressure=mp,
516     location_of_fixed_flux=mv,
517     permeability=Scalar(k,Function(self.dom)))
518 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
519     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
520 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
521     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
522 artak 1518
523 gross 2100 def testVarioF_FixedBottom_mediumK(self):
524     k=1.
525     mp=self.getScalarMask(include_bottom=True)
526     mv=self.getVectorMask(include_bottom=False)
527     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
528     p=p_ref*mp
529     u=u_ref*mv
530     df=DarcyFlow(self.dom)
531     df.setValue(g=f,
532     location_of_fixed_pressure=mp,
533     location_of_fixed_flux=mv,
534     permeability=Scalar(k,Function(self.dom)))
535 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
536     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
537 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
538     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
539 artak 1518
540 gross 2100 def testVarioF_FixedBottom_largeK(self):
541     k=1.e10
542     mp=self.getScalarMask(include_bottom=True)
543     mv=self.getVectorMask(include_bottom=False)
544     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
545     p=p_ref*mp
546     u=u_ref*mv
547     df=DarcyFlow(self.dom)
548     df.setValue(g=f,
549     location_of_fixed_pressure=mp,
550     location_of_fixed_flux=mv,
551     permeability=Scalar(k,Function(self.dom)))
552 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
553     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
554 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
555     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
556 artak 1518
557 gross 2100 def testConstF_FreeBottom_smallK(self):
558     k=1.e-10
559     mp=self.getScalarMask(include_bottom=False)
560     mv=self.getVectorMask(include_bottom=True)
561     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
562     p=p_ref*mp
563     u=u_ref*mv
564     df=DarcyFlow(self.dom)
565     df.setValue(g=f,
566 gross 2208 location_of_fixed_pressure=mp,
567 gross 2100 location_of_fixed_flux=mv,
568     permeability=Scalar(k,Function(self.dom)))
569 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
570     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
571 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
572     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
573 artak 1518
574 gross 2100 def testConstF_FreeBottom_mediumK(self):
575     k=1.
576     mp=self.getScalarMask(include_bottom=False)
577     mv=self.getVectorMask(include_bottom=True)
578     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
579     p=p_ref*mp
580     u=u_ref*mv
581     df=DarcyFlow(self.dom)
582     df.setValue(g=f,
583     location_of_fixed_pressure=mp,
584     location_of_fixed_flux=mv,
585     permeability=Scalar(k,Function(self.dom)))
586 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
587     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
588 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
589     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
590 artak 1518
591 gross 2100 def testConstF_FreeBottom_largeK(self):
592     k=1.e10
593     mp=self.getScalarMask(include_bottom=False)
594     mv=self.getVectorMask(include_bottom=True)
595     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
596     p=p_ref*mp
597     u=u_ref*mv
598     df=DarcyFlow(self.dom)
599     df.setValue(g=f,
600     location_of_fixed_pressure=mp,
601     location_of_fixed_flux=mv,
602     permeability=Scalar(k,Function(self.dom)))
603 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
604     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
605 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
606     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
607 artak 1518
608 gross 2100 def testVarioF_FreeBottom_smallK(self):
609     k=1.e-10
610     mp=self.getScalarMask(include_bottom=False)
611     mv=self.getVectorMask(include_bottom=True)
612     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
613     p=p_ref*mp
614     u=u_ref*mv
615     df=DarcyFlow(self.dom)
616     df.setValue(g=f,
617     location_of_fixed_pressure=mp,
618     location_of_fixed_flux=mv,
619     permeability=Scalar(k,Function(self.dom)))
620 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
621     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
622     self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.") # 25 because of disc. error.
623     self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
624 artak 1518
625 gross 2100 def testVarioF_FreeBottom_mediumK(self):
626     k=1.
627     mp=self.getScalarMask(include_bottom=False)
628     mv=self.getVectorMask(include_bottom=True)
629     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
630     p=p_ref*mp
631     u=u_ref*mv
632     df=DarcyFlow(self.dom)
633     df.setValue(g=f,
634     location_of_fixed_pressure=mp,
635     location_of_fixed_flux=mv,
636     permeability=Scalar(k,Function(self.dom)))
637 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
638     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
639 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
640     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
641 artak 1518
642 gross 2100 def testVarioF_FreeBottom_largeK(self):
643     k=1.e10
644     mp=self.getScalarMask(include_bottom=False)
645     mv=self.getVectorMask(include_bottom=True)
646     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
647     p=p_ref*mp
648     u=u_ref*mv
649     df=DarcyFlow(self.dom)
650     df.setValue(g=f,
651     location_of_fixed_pressure=mp,
652     location_of_fixed_flux=mv,
653     permeability=Scalar(k,Function(self.dom)))
654 gross 2208 df.setTolerance(rtol=self.TOL,p_ref=p_ref,v_ref=u_ref)
655     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE,sub_rtol=self.TOL/200)
656 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
657     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
658 artak 1518
659 gross 2100 class Test_Darcy2D(Test_Darcy):
660 gross 2208 TOL=1e-4
661 gross 2100 WIDTH=1.
662     def setUp(self):
663 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
664 gross 2100 self.dom = Rectangle(NE,NE)
665     self.rescaleDomain()
666     def tearDown(self):
667     del self.dom
668     class Test_Darcy3D(Test_Darcy):
669     TOL=1e-4
670     WIDTH=1.
671     def setUp(self):
672     NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
673     self.dom = Brick(NE,NE,NE)
674     self.rescaleDomain()
675     def tearDown(self):
676     del self.dom
677 artak 1518
678 gross 2100
679 artak 2234 class Test_Mountains3D(unittest.TestCase):
680     def setUp(self):
681     NE=16
682     self.TOL=1e-4
683     self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
684     def tearDown(self):
685     del self.domain
686     def test_periodic(self):
687     EPS=0.01
688 gross 2100
689 artak 2234 x=self.domain.getX()
690     v = Vector(0.0, Solution(self.domain))
691     a0=1
692     a1=1
693     n0=2
694     n1=2
695     n2=0.5
696     a2=-(a0*n0+a1*n1)/n2
697     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
698     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
699     v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
700    
701     H_t=Scalar(0.0, Solution(self.domain))
702     mts=Mountains(self.domain,v,eps=EPS,z=1)
703     u,Z=mts.update(u=v,H_t=H_t)
704    
705     error_int=integrate(Z)
706     self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
707    
708     class Test_Mountains2D(unittest.TestCase):
709     def setUp(self):
710     NE=16
711     self.TOL=1e-4
712     self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
713     def tearDown(self):
714     del self.domain
715     def test_periodic(self):
716     EPS=0.01
717    
718     x=self.domain.getX()
719     v = Vector(0.0, Solution(self.domain))
720     a0=1
721     n0=1
722     n1=0.5
723     a1=-(a0*n0)/n1
724     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
725     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
726    
727     H_t=Scalar(0.0, Solution(self.domain))
728     mts=Mountains(self.domain,v,eps=EPS,z=1)
729     u,Z=mts.update(u=v,H_t=H_t)
730    
731     error_int=integrate(Z)
732     self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
733    
734 gross 1414 if __name__ == '__main__':
735     suite = unittest.TestSuite()
736 artak 2234
737 gross 2100 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
738     suite.addTest(unittest.makeSuite(Test_Darcy3D))
739     suite.addTest(unittest.makeSuite(Test_Darcy2D))
740     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
741 artak 2234 suite.addTest(unittest.makeSuite(Test_Mountains3D))
742     suite.addTest(unittest.makeSuite(Test_Mountains2D))
743 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
744     if not s.wasSuccessful(): sys.exit(1)
745    

  ViewVC Help
Powered by ViewVC 1.1.26