/[escript]/trunk/finley/test/python/run_models.py
ViewVC logotype

Contents of /trunk/finley/test/python/run_models.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2301 - (show annotations)
Thu Mar 12 01:35:26 2009 UTC (10 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 42532 byte(s)
test for power law added
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2008 by University of Queensland
5 # Earth Systems Science Computational Center (ESSCC)
6 # http://www.uq.edu.au/esscc
7 #
8 # 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 #
12 ########################################################
13
14 __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 __license__="""Licensed under the Open Software License version 3.0
19 http://www.opensource.org/licenses/osl-3.0.php"""
20 __url__="http://www.uq.edu.au/esscc/escript-finley"
21
22 import unittest
23 import tempfile
24
25 from esys.escript import *
26 from esys.finley import Rectangle
27 from esys.escript.models import DarcyFlow
28 import sys
29 import os
30 try:
31 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
32 except KeyError:
33 FINLEY_WORKDIR='.'
34
35
36 VERBOSE=False # or True
37 DETAIL_VERBOSE=False
38
39 from esys.escript import *
40 from esys.escript.models import StokesProblemCartesian, PowerLaw
41 from esys.finley import Rectangle, Brick
42
43 from esys.escript.models import Mountains
44 from math import pi
45
46 #====================================================================================================================
47 class Test_StokesProblemCartesian2D(unittest.TestCase):
48 def setUp(self):
49 NE=6
50 self.TOL=1e-3
51 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
52 def tearDown(self):
53 del self.domain
54 def test_PCG_P_0(self):
55 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 sp.setTolerance(self.TOL)
71 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
72
73 error_v0=Lsup(u[0]-u0[0])
74 error_v1=Lsup(u[1]-u0[1])/0.25
75 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
80 def test_PCG_P_small(self):
81 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 sp.setTolerance(self.TOL*0.2)
97 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
98 error_v0=Lsup(u[0]-u0[0])
99 error_v1=Lsup(u[1]-u0[1])/0.25
100 error_p=Lsup(P1*x[0]*x[1]+p)
101 saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
102 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 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
105
106 def test_PCG_P_large(self):
107 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 sp.setTolerance(self.TOL)
123 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
124
125 error_v0=Lsup(u[0]-u0[0])
126 error_v1=Lsup(u[1]-u0[1])/0.25
127 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
132 def test_GMRES_P_0(self):
133 ETA=1.
134 P1=0.
135
136 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 sp.setTolerance(self.TOL)
149 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
150
151 error_v0=Lsup(u[0]-u0[0])
152 error_v1=Lsup(u[1]-u0[1])/0.25
153 error_p=Lsup(P1*x[0]*x[1]+p)
154 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
158 def test_GMRES_P_small(self):
159 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 sp.setTolerance(self.TOL*0.1)
175 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=20,usePCG=False)
176
177 error_v0=Lsup(u[0]-u0[0])
178 error_v1=Lsup(u[1]-u0[1])/0.25
179 error_p=Lsup(P1*x[0]*x[1]+p)
180 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
184 def test_GMRES_P_large(self):
185 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 sp.setTolerance(self.TOL)
201 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False)
202
203 error_v0=Lsup(u[0]-u0[0])
204 error_v1=Lsup(u[1]-u0[1])/0.25
205 error_p=Lsup(P1*x[0]*x[1]+p)/P1
206 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 def setUp(self):
212 NE=6
213 self.TOL=1e-4
214 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
215 def tearDown(self):
216 del self.domain
217 def test_PCG_P_0(self):
218 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 +whereZero(x[1]) * [1.,0.,1.] \
227 +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 sp.setTolerance(self.TOL)
238 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
239
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
244 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
249 def test_PCG_P_small(self):
250 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 +whereZero(x[1]) * [1.,0.,1.] \
258 +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 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 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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
274 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 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 +whereZero(x[1]) * [1.,0.,1.] \
287 +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 sp.setTolerance(self.TOL)
298 u,p=sp.solve(u0,-p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
299
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
304 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
309 def test_GMRES_P_0(self):
310 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 sp.setTolerance(self.TOL)
330 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
331
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
336 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 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 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
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
366 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 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
370 def test_GMRES_P_large(self):
371 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 +whereZero(x[1]) * [1.,0.,1.] \
379 +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 sp.setTolerance(self.TOL)
390 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=False)
391
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 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
396 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
446 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 df.setTolerance(rtol=self.TOL)
468 df.setSubProblemTolerance()
469 v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
470 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 df.setTolerance(rtol=self.TOL)
485 df.setSubProblemTolerance()
486 v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
487 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
488 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
489
490 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 df.setTolerance(rtol=self.TOL)
503 df.setSubProblemTolerance()
504 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
505 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
508 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 df.setTolerance(rtol=self.TOL)
521 df.setSubProblemTolerance()
522 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
523 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
526 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 df.setTolerance(rtol=self.TOL)
539 df.setSubProblemTolerance()
540 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
541 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
544 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 df.setTolerance(rtol=self.TOL)
557 df.setSubProblemTolerance()
558 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
559 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
562 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 location_of_fixed_pressure=mp,
572 location_of_fixed_flux=mv,
573 permeability=Scalar(k,Function(self.dom)))
574 df.setTolerance(rtol=self.TOL)
575 df.setSubProblemTolerance()
576 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
577 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
580 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 df.setTolerance(rtol=self.TOL)
593 df.setSubProblemTolerance()
594 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
595 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
598 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 df.setTolerance(rtol=self.TOL)
611 df.setSubProblemTolerance()
612 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
613 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
616 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 df.setTolerance(rtol=self.TOL)
629 df.setSubProblemTolerance()
630 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
631 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
634 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 df.setTolerance(rtol=self.TOL)
647 df.setSubProblemTolerance()
648 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
649 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
652 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 df.setTolerance(rtol=self.TOL)
665 df.setSubProblemTolerance()
666 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
667 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
670 class Test_Darcy2D(Test_Darcy):
671 TOL=1e-4
672 WIDTH=1.
673 def setUp(self):
674 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
675 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
689
690 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
700 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
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 if __name__ == '__main__':
928 suite = unittest.TestSuite()
929
930 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 suite.addTest(unittest.makeSuite(Test_Mountains3D))
935 suite.addTest(unittest.makeSuite(Test_Mountains2D))
936 suite.addTest(unittest.makeSuite(Test_Rheologies))
937 s=unittest.TextTestRunner(verbosity=2).run(suite)
938 if not s.wasSuccessful(): sys.exit(1)
939

  ViewVC Help
Powered by ViewVC 1.1.26