/[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 2344 - (hide annotations)
Mon Mar 30 02:13:58 2009 UTC (10 years, 5 months ago) by jfenwick
File MIME type: text/x-python
File size: 42527 byte(s)
Change __url__ to launchpad site

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 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
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 gross 2301 from esys.escript.models import StokesProblemCartesian, PowerLaw
41 gross 1414 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 2301
746    
747     class Test_Rheologies(unittest.TestCase):
748    
749     """
750     this is the program used to generate the powerlaw tests:
751    
752     TAU_Y=100.
753     N=10
754     M=5
755    
756     def getE(tau):
757     if tau<=TAU_Y:
758     return 1./(0.5+20*sqrt(tau))
759     else:
760     raise ValueError,"out of range."
761     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
762     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
763     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
764    
765     print tau
766     print e
767     """
768     TOL=1.e-8
769     def test_PowerLaw_Init(self):
770     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
771    
772     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
773     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
774     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
775     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
776     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
777    
778     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
779     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
780     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
781     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
782     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
783    
784     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
785     pl.setElasticShearModulus(1000)
786     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
787    
788     e=pl.getEtaN()
789     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
790     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
791     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
792     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
793     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
794     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
795    
796     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
797     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
798     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
799    
800     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
801     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
802     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
803     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
804    
805     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
806     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
807     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
808     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
809    
810     self.failUnlessRaises(ValueError,pl.getPower,-1)
811     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
812     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
813     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
814     self.failUnlessRaises(ValueError,pl.getPower,3)
815    
816     self.failUnlessRaises(ValueError,pl.getTauT,-1)
817     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
818     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
819     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
820     self.failUnlessRaises(ValueError,pl.getTauT,3)
821    
822     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
823     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
824     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
825     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
826     self.failUnlessRaises(ValueError,pl.getEtaN,3)
827    
828     def checkResult(self,id,gamma_dot_, eta, tau_ref):
829     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
830     error=abs(gamma_dot_*eta-tau_ref)
831     self.failUnless(error<=self.TOL*tau_ref,"eta is wrong: error = gamma_dot_*eta-tau_ref = %s * %s - %s = %s (test %s)"%(gamma_dot_,eta,tau_ref,error,id))
832    
833     def test_PowerLaw_Linear(self):
834     taus= [0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
835     gamma_dot_s=[0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0, 45.0, 50.0, 55.0, 60.0, 65.0, 70.0, 75.0]
836     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
837     pl.setDruckerPragerLaw(tau_Y=100.)
838     pl.setPowerLaw(eta_N=2.)
839     pl.setEtaTolerance(self.TOL)
840     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
841    
842     def test_PowerLaw_QuadLarge(self):
843     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
844     gamma_dot_s=[0.0, 637.45553203367592, 1798.8543819998317, 3301.3353450309969, 5079.6442562694074, 7096.067811865476, 9325.1600308977995, 11748.240371477059, 14350.835055998656, 17121.29936490925, 20050.0, 22055.0, 24060.0, 26065.0, 28070.0, 30075.0]
845     pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
846     pl.setDruckerPragerLaw(tau_Y=100.)
847     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
848     pl.setEtaTolerance(self.TOL)
849     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
850    
851     def test_PowerLaw_QuadSmall(self):
852     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
853     gamma_dot_s=[0.0, 5.632455532033676, 11.788854381999831, 18.286335345030995, 25.059644256269408, 32.071067811865476, 39.295160030897804, 46.713240371477056, 54.310835055998652, 62.076299364909254, 70.0, 77.0, 84.0, 91.0, 98.0, 105.0]
854     pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
855     pl.setDruckerPragerLaw(tau_Y=100.)
856     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
857     pl.setEtaTolerance(self.TOL)
858     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
859    
860     def test_PowerLaw_CubeLarge(self):
861     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
862     gamma_dot_s=[0.0, 51.415888336127786, 157.36125994561547, 304.6468153816889, 487.84283811405851, 703.60440414872664, 949.57131887226353, 1223.9494765692673, 1525.3084267560891, 1852.4689652218574, 2204.4346900318833, 2424.8781590350718, 2645.3216280382603, 2865.7650970414484, 3086.2085660446369, 3306.6520350478249]
863     pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
864     pl.setDruckerPragerLaw(tau_Y=100.)
865     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
866     pl.setEtaTolerance(self.TOL)
867     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
868    
869     def test_PowerLaw_CubeSmall(self):
870     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
871     gamma_dot_s=[0.0, 5.4641588833612778, 11.473612599456157, 17.89646815381689, 24.678428381140588, 31.786044041487269, 39.195713188722635, 46.889494765692675, 54.853084267560895, 63.074689652218574, 71.544346900318828, 78.698781590350706, 85.853216280382583, 93.007650970414474, 100.16208566044635, 107.316520350478]
872     pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
873     pl.setDruckerPragerLaw(tau_Y=100.)
874     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
875     pl.setEtaTolerance(self.TOL)
876     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
877    
878     def test_PowerLaw_QuadLarge_CubeLarge(self):
879     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
880     gamma_dot_s=[0.0, 683.87142036980367, 1946.2156419454475, 3590.982160412686, 5547.4870943834667, 7774.6722160142008, 10244.731349770063, 12937.189848046326, 15836.143482754744, 18928.768330131104, 22204.434690031885, 24424.878159035074, 26645.321628038262, 28865.765097041451, 31086.208566044639, 33306.652035047824]
881     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
882     pl.setDruckerPragerLaw(tau_Y=100.)
883     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
884     pl.setEtaTolerance(self.TOL)
885     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
886    
887     def test_PowerLaw_QuadLarge_CubeSmall(self):
888     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
889     gamma_dot_s=[0.0, 637.9196909170372, 1800.3279945992881, 3304.2318131848137, 5084.3226846505486, 7102.853855906963, 9334.3557440865225, 11760.129866242751, 14365.688140266215, 17139.374054561471, 20071.544346900318, 22078.698781590349, 24085.853216280382, 26093.007650970416, 28100.162085660446, 30107.316520350476]
890     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
891     pl.setDruckerPragerLaw(tau_Y=100.)
892     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
893     pl.setEtaTolerance(self.TOL)
894     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
895    
896     def test_PowerLaw_QuadSmall_CubeLarge(self):
897     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
898     gamma_dot_s=[0.0, 52.04834386816146, 159.15011432761528, 307.93315072671987, 492.9024823703279, 710.67547196059206, 958.86647890316135, 1235.6627169407443, 1539.6192618120876, 1869.5452645867665, 2224.4346900318833, 2446.8781590350718, 2669.3216280382603, 2891.7650970414484, 3114.2085660446369, 3336.6520350478249]
899    
900     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
901     pl.setDruckerPragerLaw(tau_Y=100.)
902     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
903     pl.setEtaTolerance(self.TOL)
904     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
905    
906     def test_PowerLaw_QuadSmall_CubeSmall(self):
907     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
908     gamma_dot_s=[0.0, 6.0966144153949529, 13.262466981455987, 21.182803498847885, 29.738072637409996, 38.857111853352741, 48.49087321962044, 58.602735137169738, 69.16391932355954, 80.150989017127827, 91.544346900318828, 100.69878159035071, 109.85321628038258, 119.00765097041447, 128.16208566044637, 137.31652035047824]
909     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
910     pl.setDruckerPragerLaw(tau_Y=100.)
911     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
912     pl.setEtaTolerance(self.TOL)
913     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
914    
915     def test_PowerLaw_withShear(self):
916     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
917     gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
918     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
919     pl.setDruckerPragerLaw(tau_Y=100.)
920     pl.setPowerLaw(eta_N=2.)
921     pl.setElasticShearModulus(3.)
922     dt=1./3.
923     pl.setEtaTolerance(self.TOL)
924     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
925     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
926    
927 gross 1414 if __name__ == '__main__':
928     suite = unittest.TestSuite()
929 artak 2234
930 gross 2100 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
931     suite.addTest(unittest.makeSuite(Test_Darcy3D))
932     suite.addTest(unittest.makeSuite(Test_Darcy2D))
933     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
934 artak 2234 suite.addTest(unittest.makeSuite(Test_Mountains3D))
935     suite.addTest(unittest.makeSuite(Test_Mountains2D))
936 gross 2301 suite.addTest(unittest.makeSuite(Test_Rheologies))
937 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
938     if not s.wasSuccessful(): sys.exit(1)
939    

  ViewVC Help
Powered by ViewVC 1.1.26