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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 2564 - (hide annotations)
Tue Jul 28 04:46:38 2009 UTC (10 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 47206 byte(s)
reduced PDE solutions added to topography
1 ksteube 1809
2     ########################################################
3 gross 1414 #
4 jfenwick 2548 # Copyright (c) 2003-2009 by University of Queensland
5 ksteube 1809 # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 jfenwick 2549 __copyright__="""Copyright (c) 2003-2009 by University of Queensland
15 ksteube 1809 Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18 gross 1414 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 import unittest
23     import tempfile
24    
25     from esys.escript import *
26     from esys.finley import Rectangle
27 gross 2100 from esys.escript.models import DarcyFlow
28 gross 1414 import sys
29     import os
30     try:
31     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
32     except KeyError:
33     FINLEY_WORKDIR='.'
34    
35    
36 gross 2389 VERBOSE=False # or True
37 gross 1414
38     from esys.escript import *
39 gross 2432 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian
40 gross 1414 from esys.finley import Rectangle, Brick
41 gross 2100
42 artak 2234 from esys.escript.models import Mountains
43     from math import pi
44    
45 gross 2100 #====================================================================================================================
46     class Test_StokesProblemCartesian2D(unittest.TestCase):
47 gross 1414 def setUp(self):
48 gross 2100 NE=6
49 gross 2208 self.TOL=1e-3
50 gross 1414 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
51     def tearDown(self):
52     del self.domain
53 gross 2100 def test_PCG_P_0(self):
54 gross 1414 ETA=1.
55     P1=0.
56    
57     x=self.domain.getX()
58     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
59     mask=whereZero(x[0]) * [1.,1.] \
60     +whereZero(x[0]-1) * [1.,1.] \
61     +whereZero(x[1]) * [1.,0.] \
62     +whereZero(x[1]-1) * [1.,1.]
63    
64     sp=StokesProblemCartesian(self.domain)
65    
66     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
67     u0=(1-x[0])*x[0]*[0.,1.]
68 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
69 gross 2100 sp.setTolerance(self.TOL)
70 gross 2474 u,p=sp.solve(u0,p0,verbose=VERBOSE,max_iter=100,usePCG=True)
71 gross 1414
72     error_v0=Lsup(u[0]-u0[0])
73     error_v1=Lsup(u[1]-u0[1])/0.25
74 gross 2100 error_p=Lsup(p+P1*x[0]*x[1])
75     self.failUnless(error_p<10*self.TOL, "pressure error too large.")
76     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
77     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
78 gross 1414
79 gross 2100 def test_PCG_P_small(self):
80 gross 1414 ETA=1.
81     P1=1.
82    
83     x=self.domain.getX()
84     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
85     mask=whereZero(x[0]) * [1.,1.] \
86     +whereZero(x[0]-1) * [1.,1.] \
87     +whereZero(x[1]) * [1.,0.] \
88     +whereZero(x[1]-1) * [1.,1.]
89    
90     sp=StokesProblemCartesian(self.domain)
91    
92     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
93     u0=(1-x[0])*x[0]*[0.,1.]
94 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
95 gross 2156 sp.setTolerance(self.TOL*0.2)
96 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
97 gross 1414 error_v0=Lsup(u[0]-u0[0])
98     error_v1=Lsup(u[1]-u0[1])/0.25
99 gross 2100 error_p=Lsup(P1*x[0]*x[1]+p)
100     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 gross 2156 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
103 gross 1414
104 gross 2100 def test_PCG_P_large(self):
105 gross 1414 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
120 gross 2100 sp.setTolerance(self.TOL)
121 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=True)
122 gross 1414
123     error_v0=Lsup(u[0]-u0[0])
124     error_v1=Lsup(u[1]-u0[1])/0.25
125 gross 2100 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 gross 1414
130 gross 2100 def test_GMRES_P_0(self):
131 gross 1471 ETA=1.
132     P1=0.
133 gross 1414
134 gross 1471 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
146 gross 2100 sp.setTolerance(self.TOL)
147 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
148 gross 1471
149     error_v0=Lsup(u[0]-u0[0])
150     error_v1=Lsup(u[1]-u0[1])/0.25
151 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
152 gross 2100 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 artak 1518
156 gross 2100 def test_GMRES_P_small(self):
157 gross 1471 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
172 gross 2156 sp.setTolerance(self.TOL*0.1)
173 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=20,usePCG=False)
174 gross 1471
175     error_v0=Lsup(u[0]-u0[0])
176     error_v1=Lsup(u[1]-u0[1])/0.25
177 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
178 gross 2100 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 gross 1471
182 gross 2100 def test_GMRES_P_large(self):
183 gross 1471 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
198 gross 2100 sp.setTolerance(self.TOL)
199 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
200 gross 1471
201     error_v0=Lsup(u[0]-u0[0])
202     error_v1=Lsup(u[1]-u0[1])/0.25
203 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)/P1
204 gross 2100 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 gross 1414 def setUp(self):
210 gross 2100 NE=6
211 gross 2208 self.TOL=1e-4
212 gross 1414 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
213     def tearDown(self):
214     del self.domain
215 gross 2100 def test_PCG_P_0(self):
216 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
225 gross 1414 +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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
235 gross 2156 sp.setTolerance(self.TOL)
236 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
237 gross 1414
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
242 gross 2100 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 artak 1518
247 gross 2100 def test_PCG_P_small(self):
248 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
256 gross 1414 +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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
266 gross 2251 sp.setTolerance(self.TOL*0.1)
267 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
268 gross 1414 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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
272 gross 2100 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 gross 2445
277 gross 2100 def test_PCG_P_large(self):
278 gross 1414 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
286 gross 1414 +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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
296 gross 2156 sp.setTolerance(self.TOL)
297 gross 2474 u,p=sp.solve(u0,-p0, verbose=VERBOSE ,max_iter=100,usePCG=True)
298 gross 1414
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
303 gross 2100 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 gross 1414
308 gross 2100 def test_GMRES_P_0(self):
309 gross 1471 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
328 gross 2156 sp.setTolerance(self.TOL)
329 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
330 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
335 gross 2100 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 gross 1471 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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
358 gross 2251 sp.setTolerance(self.TOL*0.1)
359 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE,max_iter=100,usePCG=False)
360 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
365 gross 2100 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 gross 2251 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
369 gross 2100 def test_GMRES_P_large(self):
370 gross 1471 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 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
378 gross 1471 +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 gross 2445 p0=Scalar(-P1,ReducedSolution(self.domain))
388 gross 2156 sp.setTolerance(self.TOL)
389 gross 2474 u,p=sp.solve(u0,p0, verbose=VERBOSE ,max_iter=100,usePCG=False)
390 gross 1471
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 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
395 gross 2100 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 gross 1471
445 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
467     v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
468 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
483     v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
484 gross 2208 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
485 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
486 artak 1518
487 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
500     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
501 gross 2100 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 artak 1518
504 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
517     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
518 gross 2100 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 artak 1518
521 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
534     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
535 gross 2100 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 artak 1518
538 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
551     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
552 gross 2100 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 artak 1518
555 gross 2100 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 gross 2208 location_of_fixed_pressure=mp,
565 gross 2100 location_of_fixed_flux=mv,
566     permeability=Scalar(k,Function(self.dom)))
567 gross 2264 df.setTolerance(rtol=self.TOL)
568     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
569 gross 2100 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 artak 1518
572 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
585     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
586 gross 2100 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 artak 1518
589 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
602     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
603 gross 2100 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 artak 1518
606 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
619     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
620 gross 2208 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 artak 1518
623 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
636     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
637 gross 2100 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 artak 1518
640 gross 2100 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 gross 2264 df.setTolerance(rtol=self.TOL)
653     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
654 gross 2100 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 artak 1518
657 gross 2100 class Test_Darcy2D(Test_Darcy):
658 gross 2208 TOL=1e-4
659 gross 2100 WIDTH=1.
660     def setUp(self):
661 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
662 gross 2100 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 artak 1518
676 gross 2100
677 artak 2234 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 gross 2100
687 artak 2234 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 gross 2563 mts=Mountains(self.domain,eps=EPS)
700     mts.setVelocity(v)
701     Z=mts.update()
702 artak 2234
703 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[2]-1.)))
704 artak 2234 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 gross 2563 mts=Mountains(self.domain,eps=EPS)
727     mts.setVelocity(v)
728     Z=mts.update()
729 artak 2234
730 gross 2564 error_int=abs(integrate(Z*whereZero(FunctionOnBoundary(self.domain).getX()[1]-1.)))
731 artak 2234 self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
732    
733 gross 2301
734    
735     class Test_Rheologies(unittest.TestCase):
736     """
737     this is the program used to generate the powerlaw tests:
738    
739     TAU_Y=100.
740     N=10
741     M=5
742    
743     def getE(tau):
744     if tau<=TAU_Y:
745     return 1./(0.5+20*sqrt(tau))
746     else:
747     raise ValueError,"out of range."
748     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
749     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
750     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
751    
752     print tau
753     print e
754     """
755     TOL=1.e-8
756     def test_PowerLaw_Init(self):
757     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
758    
759     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
760     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
761     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
762     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
763     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
764    
765     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
766     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
767     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
768     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
769     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
770    
771     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
772     pl.setElasticShearModulus(1000)
773     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
774    
775     e=pl.getEtaN()
776     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
777     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
778     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
779     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
780     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
781     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
782    
783     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
784     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
785     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
786    
787     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
788     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
789     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
790     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
791    
792     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
793     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
794     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
795     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
796    
797     self.failUnlessRaises(ValueError,pl.getPower,-1)
798     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
799     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
800     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
801     self.failUnlessRaises(ValueError,pl.getPower,3)
802    
803     self.failUnlessRaises(ValueError,pl.getTauT,-1)
804     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
805     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
806     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
807     self.failUnlessRaises(ValueError,pl.getTauT,3)
808    
809     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
810     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
811     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
812     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
813     self.failUnlessRaises(ValueError,pl.getEtaN,3)
814    
815     def checkResult(self,id,gamma_dot_, eta, tau_ref):
816     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
817     error=abs(gamma_dot_*eta-tau_ref)
818     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))
819    
820     def test_PowerLaw_Linear(self):
821     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]
822     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]
823     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
824     pl.setDruckerPragerLaw(tau_Y=100.)
825     pl.setPowerLaw(eta_N=2.)
826     pl.setEtaTolerance(self.TOL)
827     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
828    
829     def test_PowerLaw_QuadLarge(self):
830     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]
831 gross 2438 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]
832 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
833     pl.setDruckerPragerLaw(tau_Y=100.)
834     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
835     pl.setEtaTolerance(self.TOL)
836     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
837    
838     def test_PowerLaw_QuadSmall(self):
839     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]
840 gross 2438 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]
841 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
842     pl.setDruckerPragerLaw(tau_Y=100.)
843     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
844     pl.setEtaTolerance(self.TOL)
845     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
846    
847     def test_PowerLaw_CubeLarge(self):
848     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]
849 gross 2438 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]
850 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
851     pl.setDruckerPragerLaw(tau_Y=100.)
852     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
853     pl.setEtaTolerance(self.TOL)
854     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
855    
856     def test_PowerLaw_CubeSmall(self):
857     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]
858 gross 2438 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]
859 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
860     pl.setDruckerPragerLaw(tau_Y=100.)
861     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
862     pl.setEtaTolerance(self.TOL)
863     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
864    
865     def test_PowerLaw_QuadLarge_CubeLarge(self):
866     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]
867 gross 2438 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]
868    
869 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
870     pl.setDruckerPragerLaw(tau_Y=100.)
871     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
872     pl.setEtaTolerance(self.TOL)
873     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
874    
875     def test_PowerLaw_QuadLarge_CubeSmall(self):
876     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]
877 gross 2438 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]
878    
879 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
880     pl.setDruckerPragerLaw(tau_Y=100.)
881     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
882     pl.setEtaTolerance(self.TOL)
883     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
884    
885     def test_PowerLaw_QuadSmall_CubeLarge(self):
886     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]
887 gross 2438 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]
888 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
889     pl.setDruckerPragerLaw(tau_Y=100.)
890     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
891     pl.setEtaTolerance(self.TOL)
892     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
893    
894     def test_PowerLaw_QuadSmall_CubeSmall(self):
895     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]
896 gross 2438 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]
897 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
898     pl.setDruckerPragerLaw(tau_Y=100.)
899     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
900     pl.setEtaTolerance(self.TOL)
901     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
902    
903     def test_PowerLaw_withShear(self):
904     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]
905     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]
906     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
907     pl.setDruckerPragerLaw(tau_Y=100.)
908     pl.setPowerLaw(eta_N=2.)
909     pl.setElasticShearModulus(3.)
910     dt=1./3.
911     pl.setEtaTolerance(self.TOL)
912     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
913     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
914    
915 gross 2432 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
916     TOL=1.e-6
917     VERBOSE=False
918     A=1.
919     P_max=100
920     NE=2*getMPISizeWorld()
921     tau_Y=10.
922     N_dt=10
923    
924     # material parameter:
925     tau_1=5.
926     tau_2=5.
927     eta_0=100.
928     eta_1=50.
929     eta_2=400.
930     N_1=2.
931     N_2=3.
932     def getReference(self, t):
933    
934     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
935     x=self.dom.getX()
936    
937     s_00=min(self.A*t,B)
938     tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
939 gross 2438 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.)
940 gross 2432
941     alpha=0.5*inv_eta*s_00
942     if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
943     u_ref=x*alpha
944     u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
945     sigma_ref=kronecker(self.dom)*s_00
946     sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
947    
948     p_ref=self.P_max
949     for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
950     p_ref-=integrate(p_ref)/vol(self.dom)
951     return u_ref, sigma_ref, p_ref
952    
953     def runIt(self, free=None):
954     x=self.dom.getX()
955     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
956     dt=B/int(self.N_dt/2)
957     if self.VERBOSE: print "dt =",dt
958     if self.latestart:
959     t=dt
960     else:
961     t=0
962     v,s,p=self.getReference(t)
963    
964     mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
965     mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
966     mod.setElasticShearModulus(self.mu)
967     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])
968     mod.setTolerance(self.TOL)
969     mod.setFlowTolerance(self.TOL)
970    
971     BF=Vector(self.P_max,Function(self.dom))
972     for d in range(self.dom.getDim()):
973     for d2 in range(self.dom.getDim()):
974     if d!=d2: BF[d]*=x[d2]
975     v_mask=Vector(0,Solution(self.dom))
976     if free==None:
977     for d in range(self.dom.getDim()):
978     v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
979     else:
980     for d in range(self.dom.getDim()):
981     if d == self.dom.getDim()-1:
982     v_mask[d]=whereZero(x[d]-1.)
983     else:
984     v_mask[d]=whereZero(x[d])
985     mod.setExternals(F=BF,fixed_v_mask=v_mask)
986    
987     n=self.dom.getNormal()
988     N_t=0
989     errors=[]
990     while N_t < self.N_dt:
991     t_ref=t+dt
992     v_ref, s_ref,p_ref=self.getReference(t_ref)
993     mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
994     mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
995     self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
996     t+=dt
997     N_t+=1
998    
999     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1000     p=mod.getPressure()
1001     p-=integrate(p)/vol(self.dom)
1002     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1003     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1004     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1005     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1006     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1007     self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1008     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1009     self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1010     self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1011     def tearDown(self):
1012     del self.dom
1013    
1014     def test_D2_Fixed_MuNone_LateStart(self):
1015     self.dom = Rectangle(self.NE,self.NE,order=2)
1016     self.mu=None
1017     self.latestart=True
1018     self.runIt()
1019     def test_D2_Fixed_Mu_LateStart(self):
1020     self.dom = Rectangle(self.NE,self.NE,order=2)
1021     self.mu=555.
1022     self.latestart=True
1023     self.runIt()
1024     def test_D2_Fixed_MuNone(self):
1025     self.dom = Rectangle(self.NE,self.NE,order=2)
1026     self.mu=None
1027     self.latestart=False
1028     self.runIt()
1029     def test_D2_Fixed_Mu(self):
1030     self.dom = Rectangle(self.NE,self.NE,order=2)
1031     self.mu=555.
1032     self.latestart=False
1033     self.runIt()
1034     def test_D2_Free_MuNone_LateStart(self):
1035     self.dom = Rectangle(self.NE,self.NE,order=2)
1036     self.mu=None
1037     self.latestart=True
1038     self.runIt(free=0)
1039     def test_D2_Free_Mu_LateStart(self):
1040     self.dom = Rectangle(self.NE,self.NE,order=2)
1041     self.mu=555.
1042     self.latestart=True
1043     self.runIt(free=0)
1044     def test_D2_Free_MuNone(self):
1045     self.dom = Rectangle(self.NE,self.NE,order=2)
1046     self.mu=None
1047     self.latestart=False
1048     self.runIt(free=0)
1049     def test_D2_Free_Mu(self):
1050     self.dom = Rectangle(self.NE,self.NE,order=2)
1051     self.mu=555.
1052     self.latestart=False
1053     self.runIt(free=0)
1054    
1055     def test_D3_Fixed_MuNone_LateStart(self):
1056     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1057     self.mu=None
1058     self.latestart=True
1059     self.runIt()
1060     def test_D3_Fixed_Mu_LateStart(self):
1061     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1062     self.mu=555.
1063     self.latestart=True
1064     self.runIt()
1065     def test_D3_Fixed_MuNone(self):
1066     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1067     self.mu=None
1068     self.latestart=False
1069     self.runIt()
1070     def test_D3_Fixed_Mu(self):
1071     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1072     self.mu=555.
1073     self.latestart=False
1074     self.runIt()
1075     def test_D3_Free_MuNone_LateStart(self):
1076     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1077     self.mu=None
1078     self.latestart=True
1079     self.runIt(free=0)
1080     def test_D3_Free_Mu_LateStart(self):
1081     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1082     self.mu=555.
1083     self.latestart=True
1084     self.runIt(free=0)
1085     def test_D3_Free_MuNone(self):
1086     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1087     self.mu=None
1088     self.latestart=False
1089     self.runIt(free=0)
1090     def test_D3_Free_Mu(self):
1091     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1092     self.mu=555.
1093     self.latestart=False
1094     self.runIt(free=0)
1095    
1096 gross 1414 if __name__ == '__main__':
1097     suite = unittest.TestSuite()
1098 gross 2100 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1099     suite.addTest(unittest.makeSuite(Test_Darcy3D))
1100     suite.addTest(unittest.makeSuite(Test_Darcy2D))
1101     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1102 artak 2234 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1103     suite.addTest(unittest.makeSuite(Test_Mountains2D))
1104 gross 2301 suite.addTest(unittest.makeSuite(Test_Rheologies))
1105 gross 2432 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1106 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
1107     if not s.wasSuccessful(): sys.exit(1)
1108    

  ViewVC Help
Powered by ViewVC 1.1.26