/[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 2647 - (show annotations)
Fri Sep 4 05:25:25 2009 UTC (10 years ago) by gross
File MIME type: text/x-python
File size: 56555 byte(s)
fault system add. There is still an example for the usage missing.
1
2 ########################################################
3 #
4 # Copyright (c) 2003-2009 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-2009 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__="https://launchpad.net/escript-finley"
21
22 import unittest
23 import tempfile
24
25
26
27 VERBOSE=False # or True
28
29 from esys.escript import *
30 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian, FaultSystem, DarcyFlow
31 from esys.escript.models import Mountains
32 from esys.finley import Rectangle, Brick
33
34 from math import pi
35 import numpy
36 import sys
37 import os
38 #====================================================================================================================
39 try:
40 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
41 except KeyError:
42 FINLEY_WORKDIR='.'
43
44 #====================================================================================================================
45 class Test_StokesProblemCartesian2D(unittest.TestCase):
46 def setUp(self):
47 NE=6
48 self.TOL=1e-3
49 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
50 def tearDown(self):
51 del self.domain
52 def test_PCG_P_0(self):
53 ETA=1.
54 P1=0.
55
56 x=self.domain.getX()
57 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
58 mask=whereZero(x[0]) * [1.,1.] \
59 +whereZero(x[0]-1) * [1.,1.] \
60 +whereZero(x[1]) * [1.,0.] \
61 +whereZero(x[1]-1) * [1.,1.]
62
63 sp=StokesProblemCartesian(self.domain)
64
65 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
66 u0=(1-x[0])*x[0]*[0.,1.]
67 p0=Scalar(-P1,ReducedSolution(self.domain))
68 sp.setTolerance(self.TOL)
69 u,p=sp.solve(u0,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
70
71 error_v0=Lsup(u[0]-u0[0])
72 error_v1=Lsup(u[1]-u0[1])/0.25
73 error_p=Lsup(p+P1*x[0]*x[1])
74 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
75 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
76 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
77
78 def test_PCG_P_small(self):
79 ETA=1.
80 P1=1.
81
82 x=self.domain.getX()
83 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
84 mask=whereZero(x[0]) * [1.,1.] \
85 +whereZero(x[0]-1) * [1.,1.] \
86 +whereZero(x[1]) * [1.,0.] \
87 +whereZero(x[1]-1) * [1.,1.]
88
89 sp=StokesProblemCartesian(self.domain)
90
91 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
92 u0=(1-x[0])*x[0]*[0.,1.]
93 p0=Scalar(-P1,ReducedSolution(self.domain))
94 sp.setTolerance(self.TOL*0.2)
95 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
96 error_v0=Lsup(u[0]-u0[0])
97 error_v1=Lsup(u[1]-u0[1])/0.25
98 error_p=Lsup(P1*x[0]*x[1]+p)
99 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
100 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
101 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
102
103 def test_PCG_P_large(self):
104 ETA=1.
105 P1=1000.
106
107 x=self.domain.getX()
108 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
109 mask=whereZero(x[0]) * [1.,1.] \
110 +whereZero(x[0]-1) * [1.,1.] \
111 +whereZero(x[1]) * [1.,0.] \
112 +whereZero(x[1]-1) * [1.,1.]
113
114 sp=StokesProblemCartesian(self.domain)
115
116 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
117 u0=(1-x[0])*x[0]*[0.,1.]
118 p0=Scalar(-P1,ReducedSolution(self.domain))
119 sp.setTolerance(self.TOL)
120 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
121
122 error_v0=Lsup(u[0]-u0[0])
123 error_v1=Lsup(u[1]-u0[1])/0.25
124 error_p=Lsup(P1*x[0]*x[1]+p)/P1
125 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
126 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
127 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
128
129 def test_GMRES_P_0(self):
130 ETA=1.
131 P1=0.
132
133 x=self.domain.getX()
134 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
135 mask=whereZero(x[0]) * [1.,1.] \
136 +whereZero(x[0]-1) * [1.,1.] \
137 +whereZero(x[1]) * [1.,0.] \
138 +whereZero(x[1]-1) * [1.,1.]
139
140 sp=StokesProblemCartesian(self.domain)
141
142 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
143 u0=(1-x[0])*x[0]*[0.,1.]
144 p0=Scalar(-P1,ReducedSolution(self.domain))
145 sp.setTolerance(self.TOL)
146 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
147
148 error_v0=Lsup(u[0]-u0[0])
149 error_v1=Lsup(u[1]-u0[1])/0.25
150 error_p=Lsup(P1*x[0]*x[1]+p)
151 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
152 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
153 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
154
155 def test_GMRES_P_small(self):
156 ETA=1.
157 P1=1.
158
159 x=self.domain.getX()
160 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
161 mask=whereZero(x[0]) * [1.,1.] \
162 +whereZero(x[0]-1) * [1.,1.] \
163 +whereZero(x[1]) * [1.,0.] \
164 +whereZero(x[1]-1) * [1.,1.]
165
166 sp=StokesProblemCartesian(self.domain)
167
168 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
169 u0=(1-x[0])*x[0]*[0.,1.]
170 p0=Scalar(-P1,ReducedSolution(self.domain))
171 sp.setTolerance(self.TOL*0.1)
172 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
173
174 error_v0=Lsup(u[0]-u0[0])
175 error_v1=Lsup(u[1]-u0[1])/0.25
176 error_p=Lsup(P1*x[0]*x[1]+p)
177 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
178 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
179 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
180
181 def test_GMRES_P_large(self):
182 ETA=1.
183 P1=1000.
184
185 x=self.domain.getX()
186 F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
187 mask=whereZero(x[0]) * [1.,1.] \
188 +whereZero(x[0]-1) * [1.,1.] \
189 +whereZero(x[1]) * [1.,0.] \
190 +whereZero(x[1]-1) * [1.,1.]
191
192 sp=StokesProblemCartesian(self.domain)
193
194 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
195 u0=(1-x[0])*x[0]*[0.,1.]
196 p0=Scalar(-P1,ReducedSolution(self.domain))
197 sp.setTolerance(self.TOL)
198 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
199
200 error_v0=Lsup(u[0]-u0[0])
201 error_v1=Lsup(u[1]-u0[1])/0.25
202 error_p=Lsup(P1*x[0]*x[1]+p)/P1
203 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
204 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
205 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
206 #====================================================================================================================
207 class Test_StokesProblemCartesian3D(unittest.TestCase):
208 def setUp(self):
209 NE=6
210 self.TOL=1e-4
211 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
212 def tearDown(self):
213 del self.domain
214 def test_PCG_P_0(self):
215 ETA=1.
216 P1=0.
217
218 x=self.domain.getX()
219 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.]
220 x=self.domain.getX()
221 mask=whereZero(x[0]) * [1.,1.,1.] \
222 +whereZero(x[0]-1) * [1.,1.,1.] \
223 +whereZero(x[1]) * [1.,0.,1.] \
224 +whereZero(x[1]-1) * [1.,1.,1.] \
225 +whereZero(x[2]) * [1.,1.,0.] \
226 +whereZero(x[2]-1) * [1.,1.,1.]
227
228
229 sp=StokesProblemCartesian(self.domain)
230
231 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
232 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
233 p0=Scalar(-P1,ReducedSolution(self.domain))
234 sp.setTolerance(self.TOL)
235 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
236
237 error_v0=Lsup(u[0]-u0[0])
238 error_v1=Lsup(u[1]-u0[1])
239 error_v2=Lsup(u[2]-u0[2])/0.25**2
240 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
241 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
242 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
243 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
244 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
245
246 def test_PCG_P_small(self):
247 ETA=1.
248 P1=1.
249
250 x=self.domain.getX()
251 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.]
252 mask=whereZero(x[0]) * [1.,1.,1.] \
253 +whereZero(x[0]-1) * [1.,1.,1.] \
254 +whereZero(x[1]) * [1.,0.,1.] \
255 +whereZero(x[1]-1) * [1.,1.,1.] \
256 +whereZero(x[2]) * [1.,1.,0.] \
257 +whereZero(x[2]-1) * [1.,1.,1.]
258
259
260 sp=StokesProblemCartesian(self.domain)
261
262 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
263 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
264 p0=Scalar(-P1,ReducedSolution(self.domain))
265 sp.setTolerance(self.TOL*0.1)
266 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
267 error_v0=Lsup(u[0]-u0[0])
268 error_v1=Lsup(u[1]-u0[1])
269 error_v2=Lsup(u[2]-u0[2])/0.25**2
270 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
271 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
272 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
273 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
274 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
275
276 def test_PCG_P_large(self):
277 ETA=1.
278 P1=1000.
279
280 x=self.domain.getX()
281 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.]
282 mask=whereZero(x[0]) * [1.,1.,1.] \
283 +whereZero(x[0]-1) * [1.,1.,1.] \
284 +whereZero(x[1]) * [1.,0.,1.] \
285 +whereZero(x[1]-1) * [1.,1.,1.] \
286 +whereZero(x[2]) * [1.,1.,0.] \
287 +whereZero(x[2]-1) * [1.,1.,1.]
288
289
290 sp=StokesProblemCartesian(self.domain)
291
292 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
293 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
294 p0=Scalar(-P1,ReducedSolution(self.domain))
295 sp.setTolerance(self.TOL)
296 u,p=sp.solve(u0,-p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
297
298 error_v0=Lsup(u[0]-u0[0])
299 error_v1=Lsup(u[1]-u0[1])
300 error_v2=Lsup(u[2]-u0[2])/0.25**2
301 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
302 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
303 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
304 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
305 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
306
307 def test_GMRES_P_0(self):
308 ETA=1.
309 P1=0.
310
311 x=self.domain.getX()
312 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.]
313 x=self.domain.getX()
314 mask=whereZero(x[0]) * [1.,1.,1.] \
315 +whereZero(x[0]-1) * [1.,1.,1.] \
316 +whereZero(x[1]) * [1.,1.,1.] \
317 +whereZero(x[1]-1) * [1.,1.,1.] \
318 +whereZero(x[2]) * [1.,1.,0.] \
319 +whereZero(x[2]-1) * [1.,1.,1.]
320
321
322 sp=StokesProblemCartesian(self.domain)
323
324 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
325 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
326 p0=Scalar(-P1,ReducedSolution(self.domain))
327 sp.setTolerance(self.TOL)
328 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
329
330 error_v0=Lsup(u[0]-u0[0])
331 error_v1=Lsup(u[1]-u0[1])
332 error_v2=Lsup(u[2]-u0[2])/0.25**2
333 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
334 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
335 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
336 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
337 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
338 def test_GMRES_P_small(self):
339 ETA=1.
340 P1=1.
341
342 x=self.domain.getX()
343 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.]
344 mask=whereZero(x[0]) * [1.,1.,1.] \
345 +whereZero(x[0]-1) * [1.,1.,1.] \
346 +whereZero(x[1]) * [1.,1.,1.] \
347 +whereZero(x[1]-1) * [1.,1.,1.] \
348 +whereZero(x[2]) * [1.,1.,0.] \
349 +whereZero(x[2]-1) * [1.,1.,1.]
350
351
352 sp=StokesProblemCartesian(self.domain)
353
354 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
355 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
356 p0=Scalar(-P1,ReducedSolution(self.domain))
357 sp.setTolerance(self.TOL*0.1)
358 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
359
360 error_v0=Lsup(u[0]-u0[0])
361 error_v1=Lsup(u[1]-u0[1])
362 error_v2=Lsup(u[2]-u0[2])/0.25**2
363 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
364 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
365 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
366 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
367 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
368 def test_GMRES_P_large(self):
369 ETA=1.
370 P1=1000.
371
372 x=self.domain.getX()
373 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.]
374 mask=whereZero(x[0]) * [1.,1.,1.] \
375 +whereZero(x[0]-1) * [1.,1.,1.] \
376 +whereZero(x[1]) * [1.,0.,1.] \
377 +whereZero(x[1]-1) * [1.,1.,1.] \
378 +whereZero(x[2]) * [1.,1.,0.] \
379 +whereZero(x[2]-1) * [1.,1.,1.]
380
381
382 sp=StokesProblemCartesian(self.domain)
383
384 sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
385 u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
386 p0=Scalar(-P1,ReducedSolution(self.domain))
387 sp.setTolerance(self.TOL)
388 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
389
390 error_v0=Lsup(u[0]-u0[0])
391 error_v1=Lsup(u[1]-u0[1])
392 error_v2=Lsup(u[2]-u0[2])/0.25**2
393 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
394 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
395 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
396 self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
397 self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
398 #====================================================================================================================
399 class Test_Darcy(unittest.TestCase):
400 # this is a simple test for the darcy flux problem
401 #
402 #
403 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
404 #
405 # with f = (fb-f0)/xb* x + f0
406 #
407 # u = f - k * p,x = ub
408 #
409 # we prescribe pressure at x=x0=0 to p0
410 #
411 # if we prescribe the pressure on the bottom x=xb we set
412 #
413 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
414 #
415 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
416 #
417 def rescaleDomain(self):
418 x=self.dom.getX().copy()
419 for i in xrange(self.dom.getDim()):
420 x_inf=inf(x[i])
421 x_sup=sup(x[i])
422 if i == self.dom.getDim()-1:
423 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
424 else:
425 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
426 self.dom.setX(x)
427 def getScalarMask(self,include_bottom=True):
428 x=self.dom.getX().copy()
429 x_inf=inf(x[self.dom.getDim()-1])
430 x_sup=sup(x[self.dom.getDim()-1])
431 out=whereZero(x[self.dom.getDim()-1]-x_sup)
432 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
433 return wherePositive(out)
434 def getVectorMask(self,include_bottom=True):
435 x=self.dom.getX().copy()
436 out=Vector(0.,Solution(self.dom))
437 for i in xrange(self.dom.getDim()):
438 x_inf=inf(x[i])
439 x_sup=sup(x[i])
440 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
441 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
442 return wherePositive(out)
443
444 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
445 d=self.dom.getDim()
446 x=self.dom.getX()[d-1]
447 xb=inf(x)
448 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
449 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
450 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
451 return u,p,f
452
453 def testConstF_FixedBottom_smallK(self):
454 k=1.e-10
455 mp=self.getScalarMask(include_bottom=True)
456 mv=self.getVectorMask(include_bottom=False)
457 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
458 p=p_ref*mp
459 u=u_ref*mv
460 df=DarcyFlow(self.dom)
461 df.setValue(g=f,
462 location_of_fixed_pressure=mp,
463 location_of_fixed_flux=mv,
464 permeability=Scalar(k,Function(self.dom)))
465 df.setTolerance(rtol=self.TOL)
466 v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
467 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
468 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
469 def testConstF_FixedBottom_mediumK(self):
470 k=1.
471 mp=self.getScalarMask(include_bottom=True)
472 mv=self.getVectorMask(include_bottom=False)
473 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
474 p=p_ref*mp
475 u=u_ref*mv
476 df=DarcyFlow(self.dom)
477 df.setValue(g=f,
478 location_of_fixed_pressure=mp,
479 location_of_fixed_flux=mv,
480 permeability=Scalar(k,Function(self.dom)))
481 df.setTolerance(rtol=self.TOL)
482 v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
483 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
484 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
485
486 def testConstF_FixedBottom_largeK(self):
487 k=1.e10
488 mp=self.getScalarMask(include_bottom=True)
489 mv=self.getVectorMask(include_bottom=False)
490 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
491 p=p_ref*mp
492 u=u_ref*mv
493 df=DarcyFlow(self.dom)
494 df.setValue(g=f,
495 location_of_fixed_pressure=mp,
496 location_of_fixed_flux=mv,
497 permeability=Scalar(k,Function(self.dom)))
498 df.setTolerance(rtol=self.TOL)
499 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
500 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
501 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
502
503 def testVarioF_FixedBottom_smallK(self):
504 k=1.e-10
505 mp=self.getScalarMask(include_bottom=True)
506 mv=self.getVectorMask(include_bottom=False)
507 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
508 p=p_ref*mp
509 u=u_ref*mv
510 df=DarcyFlow(self.dom)
511 df.setValue(g=f,
512 location_of_fixed_pressure=mp,
513 location_of_fixed_flux=mv,
514 permeability=Scalar(k,Function(self.dom)))
515 df.setTolerance(rtol=self.TOL)
516 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
517 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
518 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
519
520 def testVarioF_FixedBottom_mediumK(self):
521 k=1.
522 mp=self.getScalarMask(include_bottom=True)
523 mv=self.getVectorMask(include_bottom=False)
524 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
525 p=p_ref*mp
526 u=u_ref*mv
527 df=DarcyFlow(self.dom)
528 df.setValue(g=f,
529 location_of_fixed_pressure=mp,
530 location_of_fixed_flux=mv,
531 permeability=Scalar(k,Function(self.dom)))
532 df.setTolerance(rtol=self.TOL)
533 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
534 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
535 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
536
537 def testVarioF_FixedBottom_largeK(self):
538 k=1.e10
539 mp=self.getScalarMask(include_bottom=True)
540 mv=self.getVectorMask(include_bottom=False)
541 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
542 p=p_ref*mp
543 u=u_ref*mv
544 df=DarcyFlow(self.dom)
545 df.setValue(g=f,
546 location_of_fixed_pressure=mp,
547 location_of_fixed_flux=mv,
548 permeability=Scalar(k,Function(self.dom)))
549 df.setTolerance(rtol=self.TOL)
550 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
551 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
552 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
553
554 def testConstF_FreeBottom_smallK(self):
555 k=1.e-10
556 mp=self.getScalarMask(include_bottom=False)
557 mv=self.getVectorMask(include_bottom=True)
558 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
559 p=p_ref*mp
560 u=u_ref*mv
561 df=DarcyFlow(self.dom)
562 df.setValue(g=f,
563 location_of_fixed_pressure=mp,
564 location_of_fixed_flux=mv,
565 permeability=Scalar(k,Function(self.dom)))
566 df.setTolerance(rtol=self.TOL)
567 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
568 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
569 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
570
571 def testConstF_FreeBottom_mediumK(self):
572 k=1.
573 mp=self.getScalarMask(include_bottom=False)
574 mv=self.getVectorMask(include_bottom=True)
575 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
576 p=p_ref*mp
577 u=u_ref*mv
578 df=DarcyFlow(self.dom)
579 df.setValue(g=f,
580 location_of_fixed_pressure=mp,
581 location_of_fixed_flux=mv,
582 permeability=Scalar(k,Function(self.dom)))
583 df.setTolerance(rtol=self.TOL)
584 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
585 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
586 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
587
588 def testConstF_FreeBottom_largeK(self):
589 k=1.e10
590 mp=self.getScalarMask(include_bottom=False)
591 mv=self.getVectorMask(include_bottom=True)
592 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
593 p=p_ref*mp
594 u=u_ref*mv
595 df=DarcyFlow(self.dom)
596 df.setValue(g=f,
597 location_of_fixed_pressure=mp,
598 location_of_fixed_flux=mv,
599 permeability=Scalar(k,Function(self.dom)))
600 df.setTolerance(rtol=self.TOL)
601 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
602 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
603 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
604
605 def testVarioF_FreeBottom_smallK(self):
606 k=1.e-10
607 mp=self.getScalarMask(include_bottom=False)
608 mv=self.getVectorMask(include_bottom=True)
609 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
610 p=p_ref*mp
611 u=u_ref*mv
612 df=DarcyFlow(self.dom)
613 df.setValue(g=f,
614 location_of_fixed_pressure=mp,
615 location_of_fixed_flux=mv,
616 permeability=Scalar(k,Function(self.dom)))
617 df.setTolerance(rtol=self.TOL)
618 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
619 self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.") # 25 because of disc. error.
620 self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
621
622 def testVarioF_FreeBottom_mediumK(self):
623 k=1.
624 mp=self.getScalarMask(include_bottom=False)
625 mv=self.getVectorMask(include_bottom=True)
626 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
627 p=p_ref*mp
628 u=u_ref*mv
629 df=DarcyFlow(self.dom)
630 df.setValue(g=f,
631 location_of_fixed_pressure=mp,
632 location_of_fixed_flux=mv,
633 permeability=Scalar(k,Function(self.dom)))
634 df.setTolerance(rtol=self.TOL)
635 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
636 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
637 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
638
639 def testVarioF_FreeBottom_largeK(self):
640 k=1.e10
641 mp=self.getScalarMask(include_bottom=False)
642 mv=self.getVectorMask(include_bottom=True)
643 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
644 p=p_ref*mp
645 u=u_ref*mv
646 df=DarcyFlow(self.dom)
647 df.setValue(g=f,
648 location_of_fixed_pressure=mp,
649 location_of_fixed_flux=mv,
650 permeability=Scalar(k,Function(self.dom)))
651 df.setTolerance(rtol=self.TOL)
652 v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
653 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
654 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
655
656 class Test_Darcy2D(Test_Darcy):
657 TOL=1e-4
658 WIDTH=1.
659 def setUp(self):
660 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
661 self.dom = Rectangle(NE,NE)
662 self.rescaleDomain()
663 def tearDown(self):
664 del self.dom
665 class Test_Darcy3D(Test_Darcy):
666 TOL=1e-4
667 WIDTH=1.
668 def setUp(self):
669 NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
670 self.dom = Brick(NE,NE,NE)
671 self.rescaleDomain()
672 def tearDown(self):
673 del self.dom
674
675
676 class Test_Mountains3D(unittest.TestCase):
677 def setUp(self):
678 NE=16
679 self.TOL=1e-4
680 self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
681 def tearDown(self):
682 del self.domain
683 def test_periodic(self):
684 EPS=0.01
685
686 x=self.domain.getX()
687 v = Vector(0.0, Solution(self.domain))
688 a0=1
689 a1=1
690 n0=2
691 n1=2
692 n2=0.5
693 a2=-(a0*n0+a1*n1)/n2
694 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
695 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
696 v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
697
698 mts=Mountains(self.domain,eps=EPS)
699 mts.setVelocity(v)
700 Z=mts.update()
701
702 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
703 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
704
705 class Test_Mountains2D(unittest.TestCase):
706 def setUp(self):
707 NE=16
708 self.TOL=1e-4
709 self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
710 def tearDown(self):
711 del self.domain
712 def test_periodic(self):
713 EPS=0.01
714
715 x=self.domain.getX()
716 v = Vector(0.0, Solution(self.domain))
717 a0=1
718 n0=1
719 n1=0.5
720 a1=-(a0*n0)/n1
721 v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
722 v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
723
724 H_t=Scalar(0.0, Solution(self.domain))
725 mts=Mountains(self.domain,eps=EPS)
726 mts.setVelocity(v)
727 Z=mts.update()
728
729 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
730 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
731
732
733
734 class Test_Rheologies(unittest.TestCase):
735 """
736 this is the program used to generate the powerlaw tests:
737
738 TAU_Y=100.
739 N=10
740 M=5
741
742 def getE(tau):
743 if tau<=TAU_Y:
744 return 1./(0.5+20*sqrt(tau))
745 else:
746 raise ValueError,"out of range."
747 tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
748 e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
749 e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
750
751 print tau
752 print e
753 """
754 TOL=1.e-8
755 def test_PowerLaw_Init(self):
756 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
757
758 self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
759 self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
760 self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
761 self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
762 self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
763
764 self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
765 self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
766 pl.setDruckerPragerLaw(tau_Y=10,friction=3)
767 self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
768 self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
769
770 self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
771 pl.setElasticShearModulus(1000)
772 self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
773
774 e=pl.getEtaN()
775 self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
776 self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
777 self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
778 self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
779 self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
780 self.failUnlessRaises(ValueError, pl.getEtaN, 3)
781
782 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
783 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
784 self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
785
786 pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
787 self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
788 self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
789 self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
790
791 pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
792 self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
793 self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
794 self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
795
796 self.failUnlessRaises(ValueError,pl.getPower,-1)
797 self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
798 self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
799 self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
800 self.failUnlessRaises(ValueError,pl.getPower,3)
801
802 self.failUnlessRaises(ValueError,pl.getTauT,-1)
803 self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
804 self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
805 self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
806 self.failUnlessRaises(ValueError,pl.getTauT,3)
807
808 self.failUnlessRaises(ValueError,pl.getEtaN,-1)
809 self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
810 self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
811 self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
812 self.failUnlessRaises(ValueError,pl.getEtaN,3)
813
814 def checkResult(self,id,gamma_dot_, eta, tau_ref):
815 self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
816 error=abs(gamma_dot_*eta-tau_ref)
817 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))
818
819 def test_PowerLaw_Linear(self):
820 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]
821 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]
822 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
823 pl.setDruckerPragerLaw(tau_Y=100.)
824 pl.setPowerLaw(eta_N=2.)
825 pl.setEtaTolerance(self.TOL)
826 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
827
828 def test_PowerLaw_QuadLarge(self):
829 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]
830 gamma_dot_s=[0.0, 405.0, 1610.0, 3615.0, 6420.0, 10025.0, 14430.0, 19635.0, 25640.0, 32445.0, 40050.0, 44055.0, 48060.0, 52065.0, 56070.0, 60075.0]
831 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
832 pl.setDruckerPragerLaw(tau_Y=100.)
833 pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
834 pl.setEtaTolerance(self.TOL)
835 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
836
837 def test_PowerLaw_QuadSmall(self):
838 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]
839 gamma_dot_s=[0.0, 5.4, 11.6, 18.6, 26.4, 35.0, 44.4, 54.6, 65.6, 77.4, 90.0, 99.0, 108.0, 117.0, 126.0, 135.0]
840 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
841 pl.setDruckerPragerLaw(tau_Y=100.)
842 pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
843 pl.setEtaTolerance(self.TOL)
844 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
845
846 def test_PowerLaw_CubeLarge(self):
847 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]
848 gamma_dot_s=[0.0, 8.90625, 41.25, 120.46875, 270.0, 513.28125, 873.75, 1374.84375, 2040.0, 2892.65625, 3956.25, 4351.875, 4747.5, 5143.125, 5538.75, 5934.375]
849 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
850 pl.setDruckerPragerLaw(tau_Y=100.)
851 pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
852 pl.setEtaTolerance(self.TOL)
853 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
854
855 def test_PowerLaw_CubeSmall(self):
856 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]
857 gamma_dot_s=[0.0, 5.0390625, 10.3125, 16.0546875, 22.5, 29.8828125, 38.4375, 48.3984375, 60.0, 73.4765625, 89.0625, 97.96875, 106.875, 115.78125, 124.6875, 133.59375]
858 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
859 pl.setDruckerPragerLaw(tau_Y=100.)
860 pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
861 pl.setEtaTolerance(self.TOL)
862 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
863
864 def test_PowerLaw_QuadLarge_CubeLarge(self):
865 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]
866 gamma_dot_s=[0.0, 408.90625, 1641.25, 3720.46875, 6670.0, 10513.28125, 15273.75, 20974.84375, 27640.000000000004, 35292.65625, 43956.25, 48351.875, 52747.5, 57143.125, 61538.75, 65934.375]
867
868 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
869 pl.setDruckerPragerLaw(tau_Y=100.)
870 pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
871 pl.setEtaTolerance(self.TOL)
872 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
873
874 def test_PowerLaw_QuadLarge_CubeSmall(self):
875 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]
876 gamma_dot_s=[0.0, 405.0390625, 1610.3125, 3616.0546875, 6422.5, 10029.8828125, 14438.4375, 19648.3984375, 25660.0, 32473.4765625, 40089.0625, 44097.96875, 48106.875, 52115.78125, 56124.6875, 60133.59375]
877
878 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
879 pl.setDruckerPragerLaw(tau_Y=100.)
880 pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
881 pl.setEtaTolerance(self.TOL)
882 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
883
884 def test_PowerLaw_QuadSmall_CubeLarge(self):
885 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]
886 gamma_dot_s=[0.0, 9.30625, 42.85, 124.06875, 276.4, 523.28125, 888.15, 1394.44375, 2065.6, 2925.05625, 3996.25, 4395.875, 4795.5, 5195.125, 5594.75, 5994.375]
887 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
888 pl.setDruckerPragerLaw(tau_Y=100.)
889 pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
890 pl.setEtaTolerance(self.TOL)
891 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
892
893 def test_PowerLaw_QuadSmall_CubeSmall(self):
894 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]
895 gamma_dot_s=[0.0, 5.4390625, 11.9125, 19.6546875, 28.9, 39.8828125, 52.8375, 67.9984375, 85.6, 105.8765625, 129.0625, 141.96875, 154.875, 167.78125, 180.6875, 193.59375]
896 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
897 pl.setDruckerPragerLaw(tau_Y=100.)
898 pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
899 pl.setEtaTolerance(self.TOL)
900 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
901
902 def test_PowerLaw_withShear(self):
903 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]
904 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]
905 pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
906 pl.setDruckerPragerLaw(tau_Y=100.)
907 pl.setPowerLaw(eta_N=2.)
908 pl.setElasticShearModulus(3.)
909 dt=1./3.
910 pl.setEtaTolerance(self.TOL)
911 self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
912 for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
913
914 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
915 TOL=1.e-6
916 VERBOSE=False
917 A=1.
918 P_max=100
919 NE=2*getMPISizeWorld()
920 tau_Y=10.
921 N_dt=10
922
923 # material parameter:
924 tau_1=5.
925 tau_2=5.
926 eta_0=100.
927 eta_1=50.
928 eta_2=400.
929 N_1=2.
930 N_2=3.
931 def getReference(self, t):
932
933 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
934 x=self.dom.getX()
935
936 s_00=min(self.A*t,B)
937 tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
938 inv_eta= 1./self.eta_0 + 1./self.eta_1*(tau/self.tau_1)**(self.N_1-1.) + 1./self.eta_2*(tau/self.tau_2)**(self.N_2-1.)
939
940 alpha=0.5*inv_eta*s_00
941 if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
942 u_ref=x*alpha
943 u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
944 sigma_ref=kronecker(self.dom)*s_00
945 sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
946
947 p_ref=self.P_max
948 for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
949 p_ref-=integrate(p_ref)/vol(self.dom)
950 return u_ref, sigma_ref, p_ref
951
952 def runIt(self, free=None):
953 x=self.dom.getX()
954 B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
955 dt=B/int(self.N_dt/2)
956 if self.VERBOSE: print "dt =",dt
957 if self.latestart:
958 t=dt
959 else:
960 t=0
961 v,s,p=self.getReference(t)
962
963 mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
964 mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
965 mod.setElasticShearModulus(self.mu)
966 mod.setPowerLaws([self.eta_0, self.eta_1, self.eta_2], [ 1., self.tau_1, self.tau_2], [1.,self.N_1,self.N_2])
967 mod.setTolerance(self.TOL)
968 mod.setFlowTolerance(self.TOL)
969
970 BF=Vector(self.P_max,Function(self.dom))
971 for d in range(self.dom.getDim()):
972 for d2 in range(self.dom.getDim()):
973 if d!=d2: BF[d]*=x[d2]
974 v_mask=Vector(0,Solution(self.dom))
975 if free==None:
976 for d in range(self.dom.getDim()):
977 v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
978 else:
979 for d in range(self.dom.getDim()):
980 if d == self.dom.getDim()-1:
981 v_mask[d]=whereZero(x[d]-1.)
982 else:
983 v_mask[d]=whereZero(x[d])
984 mod.setExternals(F=BF,fixed_v_mask=v_mask)
985
986 n=self.dom.getNormal()
987 N_t=0
988 errors=[]
989 while N_t < self.N_dt:
990 t_ref=t+dt
991 v_ref, s_ref,p_ref=self.getReference(t_ref)
992 mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
993 mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
994 self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
995 t+=dt
996 N_t+=1
997
998 def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
999 p=mod.getPressure()
1000 p-=integrate(p)/vol(self.dom)
1001 error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1002 error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1003 error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1004 error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1005 if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1006 self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1007 self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1008 self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1009 self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1010 def tearDown(self):
1011 del self.dom
1012
1013 def test_D2_Fixed_MuNone_LateStart(self):
1014 self.dom = Rectangle(self.NE,self.NE,order=2)
1015 self.mu=None
1016 self.latestart=True
1017 self.runIt()
1018 def test_D2_Fixed_Mu_LateStart(self):
1019 self.dom = Rectangle(self.NE,self.NE,order=2)
1020 self.mu=555.
1021 self.latestart=True
1022 self.runIt()
1023 def test_D2_Fixed_MuNone(self):
1024 self.dom = Rectangle(self.NE,self.NE,order=2)
1025 self.mu=None
1026 self.latestart=False
1027 self.runIt()
1028 def test_D2_Fixed_Mu(self):
1029 self.dom = Rectangle(self.NE,self.NE,order=2)
1030 self.mu=555.
1031 self.latestart=False
1032 self.runIt()
1033 def test_D2_Free_MuNone_LateStart(self):
1034 self.dom = Rectangle(self.NE,self.NE,order=2)
1035 self.mu=None
1036 self.latestart=True
1037 self.runIt(free=0)
1038 def test_D2_Free_Mu_LateStart(self):
1039 self.dom = Rectangle(self.NE,self.NE,order=2)
1040 self.mu=555.
1041 self.latestart=True
1042 self.runIt(free=0)
1043 def test_D2_Free_MuNone(self):
1044 self.dom = Rectangle(self.NE,self.NE,order=2)
1045 self.mu=None
1046 self.latestart=False
1047 self.runIt(free=0)
1048 def test_D2_Free_Mu(self):
1049 self.dom = Rectangle(self.NE,self.NE,order=2)
1050 self.mu=555.
1051 self.latestart=False
1052 self.runIt(free=0)
1053
1054 def test_D3_Fixed_MuNone_LateStart(self):
1055 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1056 self.mu=None
1057 self.latestart=True
1058 self.runIt()
1059 def test_D3_Fixed_Mu_LateStart(self):
1060 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1061 self.mu=555.
1062 self.latestart=True
1063 self.runIt()
1064 def test_D3_Fixed_MuNone(self):
1065 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1066 self.mu=None
1067 self.latestart=False
1068 self.runIt()
1069 def test_D3_Fixed_Mu(self):
1070 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1071 self.mu=555.
1072 self.latestart=False
1073 self.runIt()
1074 def test_D3_Free_MuNone_LateStart(self):
1075 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1076 self.mu=None
1077 self.latestart=True
1078 self.runIt(free=0)
1079 def test_D3_Free_Mu_LateStart(self):
1080 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1081 self.mu=555.
1082 self.latestart=True
1083 self.runIt(free=0)
1084 def test_D3_Free_MuNone(self):
1085 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1086 self.mu=None
1087 self.latestart=False
1088 self.runIt(free=0)
1089 def test_D3_Free_Mu(self):
1090 self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1091 self.mu=555.
1092 self.latestart=False
1093 self.runIt(free=0)
1094
1095
1096 class Test_FaultSystem(unittest.TestCase):
1097 EPS=1.e-8
1098 NE=10
1099 def test_Fault2D_MaxValue(self):
1100 dom=Rectangle(2*self.NE,2*self.NE)
1101 x=dom.getX()
1102 f=FaultSystem(dim=2)
1103 f.addFault([[0.5,0.],[0.,0.5]], tag=1)
1104 f.addFault([[1.,0.5],[0.5,0.5],[0.5,1.]], tag=2)
1105
1106 m, t, l=f.getMaxValue(x[0]*(1.-x[0])*(1-x[1]))
1107 self.failUnless( m == 0.25, "wrong max value")
1108 self.failUnless( t == 1, "wrong max tag")
1109 self.failUnless( l == 0., "wrong max location")
1110 m, t, l=f.getMaxValue(x[1]*(1.-x[1])*(1-x[0])*x[0])
1111 self.failUnless( m == 0.0625, "wrong max value")
1112 self.failUnless( t == 2, "wrong max tag")
1113 self.failUnless( l == 0.5, "wrong max location")
1114 m, t, l=f.getMaxValue(x[0]*(1.-x[0])*x[1])
1115 self.failUnless( m == 0.25, "wrong max value")
1116 self.failUnless( t == 2, "wrong max tag")
1117 self.failUnless( l == 1.0, "wrong max location")
1118 m, t, l= f.getMaxValue(x[1]*(1.-x[1])*x[0])
1119 self.failUnless( m == 0.25, "wrong max value")
1120 self.failUnless( t == 2, "wrong max tag")
1121 self.failUnless( l == 0., "wrong max location")
1122 m, t, l= f.getMaxValue(x[1]*(1.-x[1])*(1.-x[0]))
1123 self.failUnless( m == 0.25, "wrong max value")
1124 self.failUnless( t == 1, "wrong max tag")
1125 self.failUnless( abs(l-0.70710678118654) <= self.EPS, "wrong max location")
1126
1127
1128 def xtest_Fault2D_TwoFaults(self):
1129 f=FaultSystem(dim=2)
1130 top1=[ [1.,0.], [1.,1.], [0.,1.] ]
1131 self.failUnlessRaises(ValueError,f.addFault,top=top1,tag=1,bottom=top1)
1132 f.addFault(top=top1,tag=1)
1133 self.failUnless(f.getDim() == 2, "wrong dimension")
1134 self.failUnless( [ 1 ] == f.getTags(), "tags wrong")
1135 self.failUnless( 2. == f.getLength(1), "length wrong")
1136 self.failUnless( 0. == f.getDepth(1), "depth wrong")
1137 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 range")
1138 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 range")
1139 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsets")
1140 segs=f.getFaultSegments(1)[0]
1141 self.failUnless( len(segs) == 3, "wrong number of segments")
1142 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1143 self.failUnless( segs[0].size == 2, "seg 0 has wrong size.")
1144 self.failUnless( numpy.linalg.norm(segs[0]-[1.,0.]) < self.EPS, "wrong vertex. 0 ")
1145 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1146 self.failUnless( segs[1].size == 2, "seg 1 has wrong size.")
1147 self.failUnless( numpy.linalg.norm(segs[1]-[1.,1.]) < self.EPS, "wrong vertex. 1 ")
1148 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex 2")
1149 self.failUnless( segs[2].size == 2, "seg 2 has wrong size.")
1150 self.failUnless( numpy.linalg.norm(segs[2]-[0.,1.]) < self.EPS, "wrong vertex. 2 ")
1151 c=f.getCenterOnSurface()
1152 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1153 self.failUnless( c.size == 2, "center size is wrong")
1154 self.failUnless( numpy.linalg.norm(c-[2./3.,2./3.]) < self.EPS, "center has wrong coordinates.")
1155 o=f.getOrientationOnSurface()/pi*180.
1156 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1157
1158 top2=[ [10.,0.], [0.,10.] ]
1159 f.addFault(top2, tag=2, w0_offsets=[0,20], w1_max=20)
1160 self.failUnless( [ 1, 2 ] == f.getTags(), "tags wrong")
1161 self.failUnless( abs(f.getLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length")
1162 self.failUnless( 0. == f.getDepth(2), "depth wrong")
1163 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range")
1164 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range")
1165 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets")
1166 segs=f.getFaultSegments(2)[0]
1167 self.failUnless( len(segs) == 2, "wrong number of segments")
1168 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0")
1169 self.failUnless( numpy.linalg.norm(segs[0]-[10.,0.]) < self.EPS, "wrong vertex. 0 ")
1170 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1")
1171 self.failUnless( numpy.linalg.norm(segs[1]-[0.,10.]) < self.EPS, "wrong vertex. 1 ")
1172 c=f.getCenterOnSurface()
1173 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1174 self.failUnless( c.size == 2, "center size is wrong")
1175 self.failUnless( numpy.linalg.norm(c-[12./5.,12./5.]) < self.EPS, "center has wrong coordinates.")
1176 o=f.getOrientationOnSurface()/pi*180.
1177 self.failUnless( abs(o-45.) < self.EPS, "wrong orientation.")
1178
1179 f.transform(rot=-pi/2., shift=[-1.,-1.])
1180 self.failUnless( [ 1, 2 ] == f.getTags(), "tags after transformation wrong")
1181 self.failUnless( 2. == f.getLength(1), "length after transformation wrong")
1182 self.failUnless( 0. == f.getDepth(1), "depth after transformation wrong")
1183 self.failUnless( (0., 2.) == f.getW0Range(1)," wrong W0 after transformation range")
1184 self.failUnless( (0., 0.) == f.getW1Range(1)," wrong W1 rangeafter transformation ")
1185 self.failUnless( [0., 1., 2.] == f.getW0Offsets(1)," wrong W0 offsetsafter transformation ")
1186 segs=f.getFaultSegments(1)[0]
1187 self.failUnless( len(segs) == 3, "wrong number of segmentsafter transformation ")
1188 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1189 self.failUnless( segs[0].size == 2, "seg 0 has wrong size after transformation.")
1190 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,0.]) < self.EPS, "wrong vertex. 0 after transformation")
1191 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex after transformation1")
1192 self.failUnless( segs[1].size == 2, "seg 1 has wrong size after transformation.")
1193 self.failUnless( numpy.linalg.norm(segs[1]-[0.,0.]) < self.EPS, "wrong vertex. after transformation1 ")
1194 self.failUnless( isinstance(segs[2], numpy.ndarray), "wrong class of vertex after transformation2")
1195 self.failUnless( segs[2].size == 2, "seg 2 has wrong size after transformation.")
1196 self.failUnless( numpy.linalg.norm(segs[2]-[0., 1.]) < self.EPS, "wrong vertex after transformation. 2 ")
1197 self.failUnless( abs(f.getLength(2)-14.1421356237) < self.EPS * 14.1421356237, "wrong length after transformation")
1198 self.failUnless( 0. == f.getDepth(2), "depth wrong after transformation")
1199 self.failUnless( (0., 20.) == f.getW0Range(2)," wrong W0 range after transformation")
1200 self.failUnless( (0., 0.) == f.getW1Range(2)," wrong W1 range after transformation")
1201 self.failUnless( [0., 20.] == f.getW0Offsets(2)," wrong W0 offsets after transformation")
1202 segs=f.getFaultSegments(2)[0]
1203 self.failUnless( len(segs) == 2, "wrong number of segments after transformation")
1204 self.failUnless( isinstance(segs[0], numpy.ndarray), "wrong class of vertex 0 after transformation")
1205 self.failUnless( numpy.linalg.norm(segs[0]-[-1.,-9]) < self.EPS, "wrong vertex. 0 after transformation")
1206 self.failUnless( isinstance(segs[1], numpy.ndarray), "wrong class of vertex 1 after transformation")
1207 self.failUnless( numpy.linalg.norm(segs[1]-[9.,1.]) < self.EPS, "wrong vertex. 1 after transformation")
1208
1209 c=f.getCenterOnSurface()
1210 self.failUnless( isinstance(c, numpy.ndarray), "center has wrong class")
1211 self.failUnless( c.size == 2, "center size is wrong")
1212 self.failUnless( numpy.linalg.norm(c-[7./5.,-7./5.]) < self.EPS, "center has wrong coordinates.")
1213 o=f.getOrientationOnSurface()/pi*180.
1214 self.failUnless( abs(o+45.) < self.EPS, "wrong orientation.")
1215
1216 p=f.getParametrization([-1.,0.],1)
1217 self.failUnless(p[1]==1., "wrong value.")
1218 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1219 p=f.getParametrization([-0.5,0.],1)
1220 self.failUnless(p[1]==1., "wrong value.")
1221 self.failUnless(abs(p[0]-0.5)<self.EPS* 0.5, "wrong value.")
1222 p=f.getParametrization([0.,0.],1)
1223 self.failUnless(p[1]==1., "wrong value.")
1224 self.failUnless(abs(p[0]-1.)<self.EPS, "wrong value.")
1225 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-8)
1226 self.failUnless(p[1]==0., "wrong value.")
1227 p=f.getParametrization([0.0000001,0.0000001],1, tol=1.e-6)
1228 self.failUnless(p[1]==1., "wrong value.")
1229 self.failUnless(abs(p[0]-1.0000001)<self.EPS, "wrong value.")
1230 p=f.getParametrization([0.,0.5],1)
1231 self.failUnless(p[1]==1., "wrong value.")
1232 self.failUnless(abs(p[0]-1.5)<self.EPS, "wrong value.")
1233 p=f.getParametrization([0,1.],1)
1234 self.failUnless(p[1]==1., "wrong value.")
1235 self.failUnless(abs(p[0]-2.)<self.EPS, "wrong value.")
1236 p=f.getParametrization([1.,1.],1)
1237 self.failUnless(p[1]==0., "wrong value.")
1238 p=f.getParametrization([0,1.11],1)
1239 self.failUnless(p[1]==0., "wrong value.")
1240 p=f.getParametrization([-1,-9.],2)
1241 self.failUnless(p[1]==1., "wrong value.")
1242 self.failUnless(abs(p[0])<self.EPS, "wrong value.")
1243 p=f.getParametrization([9,1],2)
1244 self.failUnless(p[1]==1., "wrong value.")
1245 self.failUnless(abs(p[0]-20.)<self.EPS, "wrong value.")
1246
1247
1248
1249
1250 if __name__ == '__main__':
1251 suite = unittest.TestSuite()
1252 suite.addTest(unittest.makeSuite(Test_FaultSystem))
1253 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1254 suite.addTest(unittest.makeSuite(Test_Darcy3D))
1255 suite.addTest(unittest.makeSuite(Test_Darcy2D))
1256 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1257 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1258 suite.addTest(unittest.makeSuite(Test_Mountains2D))
1259 suite.addTest(unittest.makeSuite(Test_Rheologies))
1260 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1261 s=unittest.TextTestRunner(verbosity=2).run(suite)
1262 if not s.wasSuccessful(): sys.exit(1)
1263

  ViewVC Help
Powered by ViewVC 1.1.26