/[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 2474 - (hide annotations)
Tue Jun 16 06:32:15 2009 UTC (10 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 47116 byte(s)
linearPDEs has is now using the SolverOptions class to talk to PASO
1 ksteube 1809
2     ########################################################
3 gross 1414 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 gross 1414 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 gross 1414 #
12 ksteube 1809 ########################################################
13 gross 1414
14 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18 gross 1414 __license__="""Licensed under the Open Software License version 3.0
19 ksteube 1809 http://www.opensource.org/licenses/osl-3.0.php"""
20 jfenwick 2344 __url__="https://launchpad.net/escript-finley"
21 ksteube 1809
22 gross 1414 import unittest
23     import tempfile
24    
25     from esys.escript import *
26     from esys.finley import Rectangle
27 gross 2100 from esys.escript.models import DarcyFlow
28 gross 1414 import sys
29     import os
30     try:
31     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
32     except KeyError:
33     FINLEY_WORKDIR='.'
34    
35    
36 gross 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     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 gross 2301
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 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]
831 gross 2301 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 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]
840 gross 2301 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 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]
849 gross 2301 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 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]
858 gross 2301 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 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]
867    
868 gross 2301 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 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]
877    
878 gross 2301 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 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]
887 gross 2301 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 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]
896 gross 2301 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 gross 2432 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 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.)
939 gross 2432
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 gross 1414 if __name__ == '__main__':
1096     suite = unittest.TestSuite()
1097 gross 2100 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 artak 2234 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1102     suite.addTest(unittest.makeSuite(Test_Mountains2D))
1103 gross 2301 suite.addTest(unittest.makeSuite(Test_Rheologies))
1104 gross 2432 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1105 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
1106     if not s.wasSuccessful(): sys.exit(1)
1107    

  ViewVC Help
Powered by ViewVC 1.1.26