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

  ViewVC Help
Powered by ViewVC 1.1.26