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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26