/[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 2438 - (hide annotations)
Tue May 26 02:26:52 2009 UTC (10 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 48012 byte(s)
new version of power law (blame Hans)
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 2251 DETAIL_VERBOSE=False
38 gross 1414
39     from esys.escript import *
40 gross 2432 from esys.escript.models import StokesProblemCartesian, PowerLaw, IncompressibleIsotropicFlowCartesian
41 gross 1414 from esys.finley import Rectangle, Brick
42 gross 2100
43 artak 2234 from esys.escript.models import Mountains
44     from math import pi
45    
46 gross 2100 #====================================================================================================================
47     class Test_StokesProblemCartesian2D(unittest.TestCase):
48 gross 1414 def setUp(self):
49 gross 2100 NE=6
50 gross 2208 self.TOL=1e-3
51 gross 1414 self.domain=Rectangle(NE,NE,order=2,useFullElementOrder=True)
52     def tearDown(self):
53     del self.domain
54 gross 2100 def test_PCG_P_0(self):
55 gross 1414 ETA=1.
56     P1=0.
57    
58     x=self.domain.getX()
59     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
60     mask=whereZero(x[0]) * [1.,1.] \
61     +whereZero(x[0]-1) * [1.,1.] \
62     +whereZero(x[1]) * [1.,0.] \
63     +whereZero(x[1]-1) * [1.,1.]
64    
65     sp=StokesProblemCartesian(self.domain)
66    
67     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
68     u0=(1-x[0])*x[0]*[0.,1.]
69     p0=Scalar(P1,ReducedSolution(self.domain))
70 gross 2100 sp.setTolerance(self.TOL)
71 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
72 gross 1414
73     error_v0=Lsup(u[0]-u0[0])
74     error_v1=Lsup(u[1]-u0[1])/0.25
75 gross 2100 error_p=Lsup(p+P1*x[0]*x[1])
76     self.failUnless(error_p<10*self.TOL, "pressure error too large.")
77     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
78     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
79 gross 1414
80 gross 2100 def test_PCG_P_small(self):
81 gross 1414 ETA=1.
82     P1=1.
83    
84     x=self.domain.getX()
85     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
86     mask=whereZero(x[0]) * [1.,1.] \
87     +whereZero(x[0]-1) * [1.,1.] \
88     +whereZero(x[1]) * [1.,0.] \
89     +whereZero(x[1]-1) * [1.,1.]
90    
91     sp=StokesProblemCartesian(self.domain)
92    
93     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
94     u0=(1-x[0])*x[0]*[0.,1.]
95     p0=Scalar(P1,ReducedSolution(self.domain))
96 gross 2156 sp.setTolerance(self.TOL*0.2)
97 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
98 gross 1414 error_v0=Lsup(u[0]-u0[0])
99     error_v1=Lsup(u[1]-u0[1])/0.25
100 gross 2100 error_p=Lsup(P1*x[0]*x[1]+p)
101 gross 2156 saveVTK("d.vtu",p=p, e=P1*x[0]*x[1]+p, p_ref=P1*x[0]*x[1])
102 gross 2100 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
103     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
104 gross 2156 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
105 gross 1414
106 gross 2100 def test_PCG_P_large(self):
107 gross 1414 ETA=1.
108     P1=1000.
109    
110     x=self.domain.getX()
111     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
112     mask=whereZero(x[0]) * [1.,1.] \
113     +whereZero(x[0]-1) * [1.,1.] \
114     +whereZero(x[1]) * [1.,0.] \
115     +whereZero(x[1]-1) * [1.,1.]
116    
117     sp=StokesProblemCartesian(self.domain)
118    
119     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
120     u0=(1-x[0])*x[0]*[0.,1.]
121     p0=Scalar(P1,ReducedSolution(self.domain))
122 gross 2100 sp.setTolerance(self.TOL)
123 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=True)
124 gross 1414
125     error_v0=Lsup(u[0]-u0[0])
126     error_v1=Lsup(u[1]-u0[1])/0.25
127 gross 2100 error_p=Lsup(P1*x[0]*x[1]+p)/P1
128     self.failUnless(error_p<10*self.TOL, "pressure error too large.")
129     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
130     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
131 gross 1414
132 gross 2100 def test_GMRES_P_0(self):
133 gross 1471 ETA=1.
134     P1=0.
135 gross 1414
136 gross 1471 x=self.domain.getX()
137     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
138     mask=whereZero(x[0]) * [1.,1.] \
139     +whereZero(x[0]-1) * [1.,1.] \
140     +whereZero(x[1]) * [1.,0.] \
141     +whereZero(x[1]-1) * [1.,1.]
142    
143     sp=StokesProblemCartesian(self.domain)
144    
145     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
146     u0=(1-x[0])*x[0]*[0.,1.]
147     p0=Scalar(P1,ReducedSolution(self.domain))
148 gross 2100 sp.setTolerance(self.TOL)
149 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=50,usePCG=False,iter_restart=18)
150 gross 1471
151     error_v0=Lsup(u[0]-u0[0])
152     error_v1=Lsup(u[1]-u0[1])/0.25
153 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
154 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
155     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
156     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
157 artak 1518
158 gross 2100 def test_GMRES_P_small(self):
159 gross 1471 ETA=1.
160     P1=1.
161    
162     x=self.domain.getX()
163     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
164     mask=whereZero(x[0]) * [1.,1.] \
165     +whereZero(x[0]-1) * [1.,1.] \
166     +whereZero(x[1]) * [1.,0.] \
167     +whereZero(x[1]-1) * [1.,1.]
168    
169     sp=StokesProblemCartesian(self.domain)
170    
171     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
172     u0=(1-x[0])*x[0]*[0.,1.]
173     p0=Scalar(P1,ReducedSolution(self.domain))
174 gross 2156 sp.setTolerance(self.TOL*0.1)
175 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=20,usePCG=False)
176 gross 1471
177     error_v0=Lsup(u[0]-u0[0])
178     error_v1=Lsup(u[1]-u0[1])/0.25
179 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)
180 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
181     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
182     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
183 gross 1471
184 gross 2100 def test_GMRES_P_large(self):
185 gross 1471 ETA=1.
186     P1=1000.
187    
188     x=self.domain.getX()
189     F=-P1*x[1]*[1.,0]+(2*ETA-P1*x[0])*[0.,1.]
190     mask=whereZero(x[0]) * [1.,1.] \
191     +whereZero(x[0]-1) * [1.,1.] \
192     +whereZero(x[1]) * [1.,0.] \
193     +whereZero(x[1]-1) * [1.,1.]
194    
195     sp=StokesProblemCartesian(self.domain)
196    
197     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
198     u0=(1-x[0])*x[0]*[0.,1.]
199     p0=Scalar(P1,ReducedSolution(self.domain))
200 gross 2100 sp.setTolerance(self.TOL)
201 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False)
202 gross 1471
203     error_v0=Lsup(u[0]-u0[0])
204     error_v1=Lsup(u[1]-u0[1])/0.25
205 gross 2156 error_p=Lsup(P1*x[0]*x[1]+p)/P1
206 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
207     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
208     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
209     #====================================================================================================================
210     class Test_StokesProblemCartesian3D(unittest.TestCase):
211 gross 1414 def setUp(self):
212 gross 2100 NE=6
213 gross 2208 self.TOL=1e-4
214 gross 1414 self.domain=Brick(NE,NE,NE,order=2,useFullElementOrder=True)
215     def tearDown(self):
216     del self.domain
217 gross 2100 def test_PCG_P_0(self):
218 gross 1414 ETA=1.
219     P1=0.
220    
221     x=self.domain.getX()
222     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
223     x=self.domain.getX()
224     mask=whereZero(x[0]) * [1.,1.,1.] \
225     +whereZero(x[0]-1) * [1.,1.,1.] \
226 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
227 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
228     +whereZero(x[2]) * [1.,1.,0.] \
229     +whereZero(x[2]-1) * [1.,1.,1.]
230    
231    
232     sp=StokesProblemCartesian(self.domain)
233    
234     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
235     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
236     p0=Scalar(P1,ReducedSolution(self.domain))
237 gross 2156 sp.setTolerance(self.TOL)
238 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
239 gross 1414
240     error_v0=Lsup(u[0]-u0[0])
241     error_v1=Lsup(u[1]-u0[1])
242     error_v2=Lsup(u[2]-u0[2])/0.25**2
243 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
244 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
245     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
246     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
247     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
248 artak 1518
249 gross 2100 def test_PCG_P_small(self):
250 gross 1414 ETA=1.
251     P1=1.
252    
253     x=self.domain.getX()
254     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
255     mask=whereZero(x[0]) * [1.,1.,1.] \
256     +whereZero(x[0]-1) * [1.,1.,1.] \
257 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
258 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
259     +whereZero(x[2]) * [1.,1.,0.] \
260     +whereZero(x[2]-1) * [1.,1.,1.]
261    
262    
263     sp=StokesProblemCartesian(self.domain)
264    
265     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
266     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
267     p0=Scalar(P1,ReducedSolution(self.domain))
268 gross 2251 sp.setTolerance(self.TOL*0.1)
269     u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
270 gross 1414 error_v0=Lsup(u[0]-u0[0])
271     error_v1=Lsup(u[1]-u0[1])
272     error_v2=Lsup(u[2]-u0[2])/0.25**2
273 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
274 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
275     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
276     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
277     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
278     def test_PCG_P_large(self):
279 gross 1414 ETA=1.
280     P1=1000.
281    
282     x=self.domain.getX()
283     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
284     mask=whereZero(x[0]) * [1.,1.,1.] \
285     +whereZero(x[0]-1) * [1.,1.,1.] \
286 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
287 gross 1414 +whereZero(x[1]-1) * [1.,1.,1.] \
288     +whereZero(x[2]) * [1.,1.,0.] \
289     +whereZero(x[2]-1) * [1.,1.,1.]
290    
291    
292     sp=StokesProblemCartesian(self.domain)
293    
294     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
295     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
296     p0=Scalar(P1,ReducedSolution(self.domain))
297 gross 2156 sp.setTolerance(self.TOL)
298 gross 2251 u,p=sp.solve(u0,-p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=True)
299 gross 1414
300     error_v0=Lsup(u[0]-u0[0])
301     error_v1=Lsup(u[1]-u0[1])
302     error_v2=Lsup(u[2]-u0[2])/0.25**2
303 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
304 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
305     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
306     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
307     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
308 gross 1414
309 gross 2100 def test_GMRES_P_0(self):
310 gross 1471 ETA=1.
311     P1=0.
312    
313     x=self.domain.getX()
314     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
315     x=self.domain.getX()
316     mask=whereZero(x[0]) * [1.,1.,1.] \
317     +whereZero(x[0]-1) * [1.,1.,1.] \
318     +whereZero(x[1]) * [1.,1.,1.] \
319     +whereZero(x[1]-1) * [1.,1.,1.] \
320     +whereZero(x[2]) * [1.,1.,0.] \
321     +whereZero(x[2]-1) * [1.,1.,1.]
322    
323    
324     sp=StokesProblemCartesian(self.domain)
325    
326     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
327     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
328     p0=Scalar(P1,ReducedSolution(self.domain))
329 gross 2156 sp.setTolerance(self.TOL)
330 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE,max_iter=100,usePCG=False,iter_restart=20)
331 gross 1471
332     error_v0=Lsup(u[0]-u0[0])
333     error_v1=Lsup(u[1]-u0[1])
334     error_v2=Lsup(u[2]-u0[2])/0.25**2
335 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)
336 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
337     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
338     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
339     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
340     def test_GMRES_P_small(self):
341 gross 1471 ETA=1.
342     P1=1.
343    
344     x=self.domain.getX()
345     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
346     mask=whereZero(x[0]) * [1.,1.,1.] \
347     +whereZero(x[0]-1) * [1.,1.,1.] \
348     +whereZero(x[1]) * [1.,1.,1.] \
349     +whereZero(x[1]-1) * [1.,1.,1.] \
350     +whereZero(x[2]) * [1.,1.,0.] \
351     +whereZero(x[2]-1) * [1.,1.,1.]
352    
353    
354     sp=StokesProblemCartesian(self.domain)
355    
356     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
357     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
358     p0=Scalar(P1,ReducedSolution(self.domain))
359 gross 2251 sp.setTolerance(self.TOL*0.1)
360     u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE and False, verbose=VERBOSE,max_iter=100,usePCG=False)
361 gross 1471
362     error_v0=Lsup(u[0]-u0[0])
363     error_v1=Lsup(u[1]-u0[1])
364     error_v2=Lsup(u[2]-u0[2])/0.25**2
365 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
366 gross 2100 self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
367     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
368     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
369 gross 2251 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
370 gross 2100 def test_GMRES_P_large(self):
371 gross 1471 ETA=1.
372     P1=1000.
373    
374     x=self.domain.getX()
375     F=-P1*x[1]*x[2]*[1.,0.,0.]-P1*x[0]*x[2]*[0.,1.,0.]+(2*ETA*((1-x[0])*x[0]+(1-x[1])*x[1])-P1*x[0]*x[1])*[0.,0.,1.]
376     mask=whereZero(x[0]) * [1.,1.,1.] \
377     +whereZero(x[0]-1) * [1.,1.,1.] \
378 gross 2156 +whereZero(x[1]) * [1.,0.,1.] \
379 gross 1471 +whereZero(x[1]-1) * [1.,1.,1.] \
380     +whereZero(x[2]) * [1.,1.,0.] \
381     +whereZero(x[2]-1) * [1.,1.,1.]
382    
383    
384     sp=StokesProblemCartesian(self.domain)
385    
386     sp.initialize(f=F,fixed_u_mask=mask,eta=ETA)
387     u0=(1-x[0])*x[0]*(1-x[1])*x[1]*[0.,0.,1.]
388     p0=Scalar(P1,ReducedSolution(self.domain))
389 gross 2156 sp.setTolerance(self.TOL)
390 gross 2251 u,p=sp.solve(u0,p0,show_details=DETAIL_VERBOSE, verbose=VERBOSE ,max_iter=100,usePCG=False)
391 gross 1471
392     error_v0=Lsup(u[0]-u0[0])
393     error_v1=Lsup(u[1]-u0[1])
394     error_v2=Lsup(u[2]-u0[2])/0.25**2
395 gross 2156 error_p=Lsup(P1*x[0]*x[1]*x[2]+p)/P1
396 gross 2100 self.failUnless(error_p<10*self.TOL, "pressure error too large.")
397     self.failUnless(error_v0<10*self.TOL, "0-velocity error too large.")
398     self.failUnless(error_v1<10*self.TOL, "1-velocity error too large.")
399     self.failUnless(error_v2<10*self.TOL, "2-velocity error too large.")
400     #====================================================================================================================
401     class Test_Darcy(unittest.TestCase):
402     # this is a simple test for the darcy flux problem
403     #
404     #
405     # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
406     #
407     # with f = (fb-f0)/xb* x + f0
408     #
409     # u = f - k * p,x = ub
410     #
411     # we prescribe pressure at x=x0=0 to p0
412     #
413     # if we prescribe the pressure on the bottom x=xb we set
414     #
415     # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
416     #
417     # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
418     #
419     def rescaleDomain(self):
420     x=self.dom.getX().copy()
421     for i in xrange(self.dom.getDim()):
422     x_inf=inf(x[i])
423     x_sup=sup(x[i])
424     if i == self.dom.getDim()-1:
425     x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
426     else:
427     x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
428     self.dom.setX(x)
429     def getScalarMask(self,include_bottom=True):
430     x=self.dom.getX().copy()
431     x_inf=inf(x[self.dom.getDim()-1])
432     x_sup=sup(x[self.dom.getDim()-1])
433     out=whereZero(x[self.dom.getDim()-1]-x_sup)
434     if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
435     return wherePositive(out)
436     def getVectorMask(self,include_bottom=True):
437     x=self.dom.getX().copy()
438     out=Vector(0.,Solution(self.dom))
439     for i in xrange(self.dom.getDim()):
440     x_inf=inf(x[i])
441     x_sup=sup(x[i])
442     if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
443     if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
444     return wherePositive(out)
445 gross 1471
446 gross 2100 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
447     d=self.dom.getDim()
448     x=self.dom.getX()[d-1]
449     xb=inf(x)
450     u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
451     p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
452     f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
453     return u,p,f
454    
455     def testConstF_FixedBottom_smallK(self):
456     k=1.e-10
457     mp=self.getScalarMask(include_bottom=True)
458     mv=self.getVectorMask(include_bottom=False)
459     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
460     p=p_ref*mp
461     u=u_ref*mv
462     df=DarcyFlow(self.dom)
463     df.setValue(g=f,
464     location_of_fixed_pressure=mp,
465     location_of_fixed_flux=mv,
466     permeability=Scalar(k,Function(self.dom)))
467 gross 2264 df.setTolerance(rtol=self.TOL)
468     df.setSubProblemTolerance()
469     v,p=df.solve(u_ref,p, max_iter=100, verbose=VERBOSE)
470 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
471     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
472     def testConstF_FixedBottom_mediumK(self):
473     k=1.
474     mp=self.getScalarMask(include_bottom=True)
475     mv=self.getVectorMask(include_bottom=False)
476     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
477     p=p_ref*mp
478     u=u_ref*mv
479     df=DarcyFlow(self.dom)
480     df.setValue(g=f,
481     location_of_fixed_pressure=mp,
482     location_of_fixed_flux=mv,
483     permeability=Scalar(k,Function(self.dom)))
484 gross 2264 df.setTolerance(rtol=self.TOL)
485     df.setSubProblemTolerance()
486     v,p=df.solve(u,p,max_iter=100, verbose=VERBOSE )
487 gross 2208 self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
488 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
489 artak 1518
490 gross 2100 def testConstF_FixedBottom_largeK(self):
491     k=1.e10
492     mp=self.getScalarMask(include_bottom=True)
493     mv=self.getVectorMask(include_bottom=False)
494     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
495     p=p_ref*mp
496     u=u_ref*mv
497     df=DarcyFlow(self.dom)
498     df.setValue(g=f,
499     location_of_fixed_pressure=mp,
500     location_of_fixed_flux=mv,
501     permeability=Scalar(k,Function(self.dom)))
502 gross 2264 df.setTolerance(rtol=self.TOL)
503     df.setSubProblemTolerance()
504     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
505 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
506     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
507 artak 1518
508 gross 2100 def testVarioF_FixedBottom_smallK(self):
509     k=1.e-10
510     mp=self.getScalarMask(include_bottom=True)
511     mv=self.getVectorMask(include_bottom=False)
512     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
513     p=p_ref*mp
514     u=u_ref*mv
515     df=DarcyFlow(self.dom)
516     df.setValue(g=f,
517     location_of_fixed_pressure=mp,
518     location_of_fixed_flux=mv,
519     permeability=Scalar(k,Function(self.dom)))
520 gross 2264 df.setTolerance(rtol=self.TOL)
521     df.setSubProblemTolerance()
522     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
523 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
524     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
525 artak 1518
526 gross 2100 def testVarioF_FixedBottom_mediumK(self):
527     k=1.
528     mp=self.getScalarMask(include_bottom=True)
529     mv=self.getVectorMask(include_bottom=False)
530     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
531     p=p_ref*mp
532     u=u_ref*mv
533     df=DarcyFlow(self.dom)
534     df.setValue(g=f,
535     location_of_fixed_pressure=mp,
536     location_of_fixed_flux=mv,
537     permeability=Scalar(k,Function(self.dom)))
538 gross 2264 df.setTolerance(rtol=self.TOL)
539     df.setSubProblemTolerance()
540     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
541 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
542     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
543 artak 1518
544 gross 2100 def testVarioF_FixedBottom_largeK(self):
545     k=1.e10
546     mp=self.getScalarMask(include_bottom=True)
547     mv=self.getVectorMask(include_bottom=False)
548     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
549     p=p_ref*mp
550     u=u_ref*mv
551     df=DarcyFlow(self.dom)
552     df.setValue(g=f,
553     location_of_fixed_pressure=mp,
554     location_of_fixed_flux=mv,
555     permeability=Scalar(k,Function(self.dom)))
556 gross 2264 df.setTolerance(rtol=self.TOL)
557     df.setSubProblemTolerance()
558     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
559 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
560     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
561 artak 1518
562 gross 2100 def testConstF_FreeBottom_smallK(self):
563     k=1.e-10
564     mp=self.getScalarMask(include_bottom=False)
565     mv=self.getVectorMask(include_bottom=True)
566     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
567     p=p_ref*mp
568     u=u_ref*mv
569     df=DarcyFlow(self.dom)
570     df.setValue(g=f,
571 gross 2208 location_of_fixed_pressure=mp,
572 gross 2100 location_of_fixed_flux=mv,
573     permeability=Scalar(k,Function(self.dom)))
574 gross 2264 df.setTolerance(rtol=self.TOL)
575     df.setSubProblemTolerance()
576     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
577 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
578     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
579 artak 1518
580 gross 2100 def testConstF_FreeBottom_mediumK(self):
581     k=1.
582     mp=self.getScalarMask(include_bottom=False)
583     mv=self.getVectorMask(include_bottom=True)
584     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
585     p=p_ref*mp
586     u=u_ref*mv
587     df=DarcyFlow(self.dom)
588     df.setValue(g=f,
589     location_of_fixed_pressure=mp,
590     location_of_fixed_flux=mv,
591     permeability=Scalar(k,Function(self.dom)))
592 gross 2264 df.setTolerance(rtol=self.TOL)
593     df.setSubProblemTolerance()
594     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
595 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
596     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
597 artak 1518
598 gross 2100 def testConstF_FreeBottom_largeK(self):
599     k=1.e10
600     mp=self.getScalarMask(include_bottom=False)
601     mv=self.getVectorMask(include_bottom=True)
602     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
603     p=p_ref*mp
604     u=u_ref*mv
605     df=DarcyFlow(self.dom)
606     df.setValue(g=f,
607     location_of_fixed_pressure=mp,
608     location_of_fixed_flux=mv,
609     permeability=Scalar(k,Function(self.dom)))
610 gross 2264 df.setTolerance(rtol=self.TOL)
611     df.setSubProblemTolerance()
612     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
613 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
614     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
615 artak 1518
616 gross 2100 def testVarioF_FreeBottom_smallK(self):
617     k=1.e-10
618     mp=self.getScalarMask(include_bottom=False)
619     mv=self.getVectorMask(include_bottom=True)
620     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
621     p=p_ref*mp
622     u=u_ref*mv
623     df=DarcyFlow(self.dom)
624     df.setValue(g=f,
625     location_of_fixed_pressure=mp,
626     location_of_fixed_flux=mv,
627     permeability=Scalar(k,Function(self.dom)))
628 gross 2264 df.setTolerance(rtol=self.TOL)
629     df.setSubProblemTolerance()
630     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
631 gross 2208 self.failUnless(Lsup(v-u_ref)<self.TOL*25.*Lsup(u_ref), "flux error too big.") # 25 because of disc. error.
632     self.failUnless(Lsup(p-p_ref)<self.TOL*25.*Lsup(p_ref), "pressure error too big.")
633 artak 1518
634 gross 2100 def testVarioF_FreeBottom_mediumK(self):
635     k=1.
636     mp=self.getScalarMask(include_bottom=False)
637     mv=self.getVectorMask(include_bottom=True)
638     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
639     p=p_ref*mp
640     u=u_ref*mv
641     df=DarcyFlow(self.dom)
642     df.setValue(g=f,
643     location_of_fixed_pressure=mp,
644     location_of_fixed_flux=mv,
645     permeability=Scalar(k,Function(self.dom)))
646 gross 2264 df.setTolerance(rtol=self.TOL)
647     df.setSubProblemTolerance()
648     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
649 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
650     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
651 artak 1518
652 gross 2100 def testVarioF_FreeBottom_largeK(self):
653     k=1.e10
654     mp=self.getScalarMask(include_bottom=False)
655     mv=self.getVectorMask(include_bottom=True)
656     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
657     p=p_ref*mp
658     u=u_ref*mv
659     df=DarcyFlow(self.dom)
660     df.setValue(g=f,
661     location_of_fixed_pressure=mp,
662     location_of_fixed_flux=mv,
663     permeability=Scalar(k,Function(self.dom)))
664 gross 2264 df.setTolerance(rtol=self.TOL)
665     df.setSubProblemTolerance()
666     v,p=df.solve(u,p, max_iter=100, verbose=VERBOSE)
667 gross 2100 self.failUnless(Lsup(v-u_ref)<self.TOL*10.*Lsup(u_ref), "flux error too big.")
668     self.failUnless(Lsup(p-p_ref)<self.TOL*10.*Lsup(p_ref), "pressure error too big.")
669 artak 1518
670 gross 2100 class Test_Darcy2D(Test_Darcy):
671 gross 2208 TOL=1e-4
672 gross 2100 WIDTH=1.
673     def setUp(self):
674 gross 2208 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
675 gross 2100 self.dom = Rectangle(NE,NE)
676     self.rescaleDomain()
677     def tearDown(self):
678     del self.dom
679     class Test_Darcy3D(Test_Darcy):
680     TOL=1e-4
681     WIDTH=1.
682     def setUp(self):
683     NE=25 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
684     self.dom = Brick(NE,NE,NE)
685     self.rescaleDomain()
686     def tearDown(self):
687     del self.dom
688 artak 1518
689 gross 2100
690 artak 2234 class Test_Mountains3D(unittest.TestCase):
691     def setUp(self):
692     NE=16
693     self.TOL=1e-4
694     self.domain=Brick(NE,NE,NE,order=1,useFullElementOrder=True)
695     def tearDown(self):
696     del self.domain
697     def test_periodic(self):
698     EPS=0.01
699 gross 2100
700 artak 2234 x=self.domain.getX()
701     v = Vector(0.0, Solution(self.domain))
702     a0=1
703     a1=1
704     n0=2
705     n1=2
706     n2=0.5
707     a2=-(a0*n0+a1*n1)/n2
708     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])* cos(pi*n2*x[2])
709     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])* cos(pi*n2*x[2])
710     v[2]=a2*cos(pi*n0*x[0])* cos(pi*n1*x[1])* sin(pi*n2*x[2])
711    
712     H_t=Scalar(0.0, Solution(self.domain))
713     mts=Mountains(self.domain,v,eps=EPS,z=1)
714     u,Z=mts.update(u=v,H_t=H_t)
715    
716     error_int=integrate(Z)
717     self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
718    
719     class Test_Mountains2D(unittest.TestCase):
720     def setUp(self):
721     NE=16
722     self.TOL=1e-4
723     self.domain=Rectangle(NE,NE,order=1,useFullElementOrder=True)
724     def tearDown(self):
725     del self.domain
726     def test_periodic(self):
727     EPS=0.01
728    
729     x=self.domain.getX()
730     v = Vector(0.0, Solution(self.domain))
731     a0=1
732     n0=1
733     n1=0.5
734     a1=-(a0*n0)/n1
735     v[0]=a0*sin(pi*n0*x[0])* cos(pi*n1*x[1])
736     v[1]=a1*cos(pi*n0*x[0])* sin(pi*n1*x[1])
737    
738     H_t=Scalar(0.0, Solution(self.domain))
739     mts=Mountains(self.domain,v,eps=EPS,z=1)
740     u,Z=mts.update(u=v,H_t=H_t)
741    
742     error_int=integrate(Z)
743     self.failUnless(error_int<self.TOL, "Boundary intergral is too large.")
744    
745 gross 2301
746    
747     class Test_Rheologies(unittest.TestCase):
748     """
749     this is the program used to generate the powerlaw tests:
750    
751     TAU_Y=100.
752     N=10
753     M=5
754    
755     def getE(tau):
756     if tau<=TAU_Y:
757     return 1./(0.5+20*sqrt(tau))
758     else:
759     raise ValueError,"out of range."
760     tau=[ i* (TAU_Y/N) for i in range(N+1)] + [TAU_Y for i in range(M)]
761     e=[ tau[i]/getE(tau[i]) for i in range(N+1)]
762     e+= [ (i+1+j)* (max(e)/N) for j in range(M) ]
763    
764     print tau
765     print e
766     """
767     TOL=1.e-8
768     def test_PowerLaw_Init(self):
769     pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
770    
771     self.failUnlessEqual(pl.getNumMaterials(),3,"num materials is wrong")
772     self.failUnlessEqual(pl.validMaterialId(0),True,"material id 0 not found")
773     self.failUnlessEqual(pl.validMaterialId(1),True,"material id 1 not found")
774     self.failUnlessEqual(pl.validMaterialId(2),True,"material id 2 not found")
775     self.failUnlessEqual(pl.validMaterialId(3),False,"material id 3 not found")
776    
777     self.failUnlessEqual(pl.getFriction(),None,"initial friction wrong.")
778     self.failUnlessEqual(pl.getTauY(),None,"initial tau y wrong.")
779     pl.setDruckerPragerLaw(tau_Y=10,friction=3)
780     self.failUnlessEqual(pl.getFriction(),3,"friction wrong.")
781     self.failUnlessEqual(pl.getTauY(),10,"tau y wrong.")
782    
783     self.failUnlessEqual(pl.getElasticShearModulus(),None,"initial shear modulus is wrong")
784     pl.setElasticShearModulus(1000)
785     self.failUnlessEqual(pl.getElasticShearModulus(),1000.,"shear modulus is wrong")
786    
787     e=pl.getEtaN()
788     self.failUnlessEqual(len(e),3,"initial length of etaN is wrong.")
789     self.failUnlessEqual(e,[None, None, None],"initial etaN are wrong.")
790     self.failUnlessEqual(pl.getEtaN(0),None,"initial material id 0 is wrong")
791     self.failUnlessEqual(pl.getEtaN(1),None,"initial material id 1 is wrong")
792     self.failUnlessEqual(pl.getEtaN(2),None,"initial material id 2 is wrong")
793     self.failUnlessRaises(ValueError, pl.getEtaN, 3)
794    
795     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30,40],tau_t=[2], power=[5,6,7])
796     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,4,5,5], power=[5,6,7])
797     self.failUnlessRaises(ValueError,pl.setPowerLaws,eta_N=[10,20,30],tau_t=[2,1,2], power=[5,6,7,7])
798    
799     pl.setPowerLaws(eta_N=[10,20,30],tau_t=[100,200,300], power=[1,2,3])
800     self.failUnlessEqual(pl.getPower(),[1,2,3],"powers are wrong.")
801     self.failUnlessEqual(pl.getTauT(),[100,200,300],"tau t are wrong.")
802     self.failUnlessEqual(pl.getEtaN(),[10,20,30],"etaN are wrong.")
803    
804     pl.setPowerLaw(id=0,eta_N=40,tau_t=400, power=4)
805     self.failUnlessEqual(pl.getPower(),[4,2,3],"powers are wrong.")
806     self.failUnlessEqual(pl.getTauT(),[400,200,300],"tau t are wrong.")
807     self.failUnlessEqual(pl.getEtaN(),[40,20,30],"etaN are wrong.")
808    
809     self.failUnlessRaises(ValueError,pl.getPower,-1)
810     self.failUnlessEqual(pl.getPower(0),4,"power 0 is wrong.")
811     self.failUnlessEqual(pl.getPower(1),2,"power 1 is wrong.")
812     self.failUnlessEqual(pl.getPower(2),3,"power 2 is wrong.")
813     self.failUnlessRaises(ValueError,pl.getPower,3)
814    
815     self.failUnlessRaises(ValueError,pl.getTauT,-1)
816     self.failUnlessEqual(pl.getTauT(0),400,"tau t 0 is wrong.")
817     self.failUnlessEqual(pl.getTauT(1),200,"tau t 1 is wrong.")
818     self.failUnlessEqual(pl.getTauT(2),300,"tau t 2 is wrong.")
819     self.failUnlessRaises(ValueError,pl.getTauT,3)
820    
821     self.failUnlessRaises(ValueError,pl.getEtaN,-1)
822     self.failUnlessEqual(pl.getEtaN(0),40,"eta n 0 is wrong.")
823     self.failUnlessEqual(pl.getEtaN(1),20,"eta n 1 is wrong.")
824     self.failUnlessEqual(pl.getEtaN(2),30,"eta n 2 is wrong.")
825     self.failUnlessRaises(ValueError,pl.getEtaN,3)
826    
827     def checkResult(self,id,gamma_dot_, eta, tau_ref):
828     self.failUnless(eta>=0,"eta needs to be positive (test %s)"%id)
829     error=abs(gamma_dot_*eta-tau_ref)
830     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))
831    
832     def test_PowerLaw_Linear(self):
833     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]
834     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]
835     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
836     pl.setDruckerPragerLaw(tau_Y=100.)
837     pl.setPowerLaw(eta_N=2.)
838     pl.setEtaTolerance(self.TOL)
839     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
840    
841     def test_PowerLaw_QuadLarge(self):
842     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]
843 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]
844 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
845     pl.setDruckerPragerLaw(tau_Y=100.)
846     pl.setPowerLaws(eta_N=[2.,0.01],tau_t=[1, 25.], power=[1,2])
847     pl.setEtaTolerance(self.TOL)
848     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
849    
850     def test_PowerLaw_QuadSmall(self):
851     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]
852 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]
853 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
854     pl.setDruckerPragerLaw(tau_Y=100.)
855     pl.setPowerLaws(eta_N=[2.,10.],tau_t=[1, 25.], power=[1,2])
856     pl.setEtaTolerance(self.TOL)
857     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
858    
859     def test_PowerLaw_CubeLarge(self):
860     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]
861 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]
862 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
863     pl.setDruckerPragerLaw(tau_Y=100.)
864     pl.setPowerLaws(eta_N=[2.,1./16.],tau_t=[1, 64.], power=[1,3])
865     pl.setEtaTolerance(self.TOL)
866     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
867    
868     def test_PowerLaw_CubeSmall(self):
869     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]
870 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]
871 gross 2301 pl=PowerLaw(numMaterials=2,verbose=VERBOSE)
872     pl.setDruckerPragerLaw(tau_Y=100.)
873     pl.setPowerLaws(eta_N=[2.,25./4.],tau_t=[1, 64.], power=[1,3])
874     pl.setEtaTolerance(self.TOL)
875     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
876    
877     def test_PowerLaw_QuadLarge_CubeLarge(self):
878     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]
879 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]
880    
881 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
882     pl.setDruckerPragerLaw(tau_Y=100.)
883     pl.setPowerLaws(eta_N=[2.,0.01,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
884     pl.setEtaTolerance(self.TOL)
885     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
886    
887     def test_PowerLaw_QuadLarge_CubeSmall(self):
888     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
889 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]
890    
891 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
892     pl.setDruckerPragerLaw(tau_Y=100.)
893     pl.setPowerLaws(eta_N=[2.,0.01,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
894     pl.setEtaTolerance(self.TOL)
895     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
896    
897     def test_PowerLaw_QuadSmall_CubeLarge(self):
898     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]
899 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]
900 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
901     pl.setDruckerPragerLaw(tau_Y=100.)
902     pl.setPowerLaws(eta_N=[2.,10.,1./16.],tau_t=[1, 25.,64.], power=[1,2,3])
903     pl.setEtaTolerance(self.TOL)
904     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
905    
906     def test_PowerLaw_QuadSmall_CubeSmall(self):
907     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
908 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]
909 gross 2301 pl=PowerLaw(numMaterials=3,verbose=VERBOSE)
910     pl.setDruckerPragerLaw(tau_Y=100.)
911     pl.setPowerLaws(eta_N=[2.,10.,25./4.],tau_t=[1, 25.,64.], power=[1,2,3])
912     pl.setEtaTolerance(self.TOL)
913     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i]),taus[i])
914    
915     def test_PowerLaw_withShear(self):
916     taus=[0.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0]
917     gamma_dot_s=[0.0, 15.0, 30.0, 45.0, 60.0, 75.0, 90.0, 105.0, 120.0, 135.0, 150.0, 165.0, 180.0, 195.0, 210.0, 225.0]
918     pl=PowerLaw(numMaterials=1,verbose=VERBOSE)
919     pl.setDruckerPragerLaw(tau_Y=100.)
920     pl.setPowerLaw(eta_N=2.)
921     pl.setElasticShearModulus(3.)
922     dt=1./3.
923     pl.setEtaTolerance(self.TOL)
924     self.failUnlessRaises(ValueError, pl.getEtaEff,gamma_dot_s[0])
925     for i in xrange(len(taus)): self.checkResult(i,gamma_dot_s[i], pl.getEtaEff(gamma_dot_s[i],dt=dt),taus[i])
926    
927 gross 2432 class Test_IncompressibleIsotropicFlowCartesian(unittest.TestCase):
928     TOL=1.e-6
929     VERBOSE=False
930     A=1.
931     P_max=100
932     NE=2*getMPISizeWorld()
933     tau_Y=10.
934     N_dt=10
935    
936     # material parameter:
937     tau_1=5.
938     tau_2=5.
939     eta_0=100.
940     eta_1=50.
941     eta_2=400.
942     N_1=2.
943     N_2=3.
944     def getReference(self, t):
945    
946     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
947     x=self.dom.getX()
948    
949     s_00=min(self.A*t,B)
950     tau=sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)*abs(s_00)
951 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.)
952 gross 2432
953     alpha=0.5*inv_eta*s_00
954     if s_00 <= B and self.mu !=None: alpha+=1./(2*self.mu)*self.A
955     u_ref=x*alpha
956     u_ref[self.dom.getDim()-1]=(1.-x[self.dom.getDim()-1])*alpha*(self.dom.getDim()-1)
957     sigma_ref=kronecker(self.dom)*s_00
958     sigma_ref[self.dom.getDim()-1,self.dom.getDim()-1]=-s_00*(self.dom.getDim()-1)
959    
960     p_ref=self.P_max
961     for d in range(self.dom.getDim()): p_ref=p_ref*x[d]
962     p_ref-=integrate(p_ref)/vol(self.dom)
963     return u_ref, sigma_ref, p_ref
964    
965     def runIt(self, free=None):
966     x=self.dom.getX()
967     B=self.tau_Y/sqrt((self.dom.getDim()-1)*self.dom.getDim()*0.5)
968     dt=B/int(self.N_dt/2)
969     if self.VERBOSE: print "dt =",dt
970     if self.latestart:
971     t=dt
972     else:
973     t=0
974     v,s,p=self.getReference(t)
975    
976     mod=IncompressibleIsotropicFlowCartesian(self.dom, stress=s, v=v, p=p, t=t, numMaterials=3, verbose=self.VERBOSE)
977     mod.setDruckerPragerLaw(tau_Y=self.tau_Y,friction=None)
978     mod.setElasticShearModulus(self.mu)
979     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])
980     mod.setTolerance(self.TOL)
981     mod.setFlowSubTolerance(self.TOL**2)
982     mod.setFlowTolerance(self.TOL)
983    
984     BF=Vector(self.P_max,Function(self.dom))
985     for d in range(self.dom.getDim()):
986     for d2 in range(self.dom.getDim()):
987     if d!=d2: BF[d]*=x[d2]
988     v_mask=Vector(0,Solution(self.dom))
989     if free==None:
990     for d in range(self.dom.getDim()):
991     v_mask[d]=whereZero(x[d])+whereZero(x[d]-1.)
992     else:
993     for d in range(self.dom.getDim()):
994     if d == self.dom.getDim()-1:
995     v_mask[d]=whereZero(x[d]-1.)
996     else:
997     v_mask[d]=whereZero(x[d])
998     mod.setExternals(F=BF,fixed_v_mask=v_mask)
999    
1000     n=self.dom.getNormal()
1001     N_t=0
1002     errors=[]
1003     while N_t < self.N_dt:
1004     t_ref=t+dt
1005     v_ref, s_ref,p_ref=self.getReference(t_ref)
1006     mod.setExternals(f=matrixmult(s_ref,n)-p_ref*n, v_boundary=v_ref)
1007     mod.update(dt, iter_max=100, inner_iter_max=20, verbose=self.VERBOSE, usePCG=True)
1008     self.check(N_t,mod,t_ref,v_ref, s_ref,p_ref)
1009     t+=dt
1010     N_t+=1
1011    
1012     def check(self,N_t,mod,t_ref,v_ref, s_ref,p_ref):
1013     p=mod.getPressure()
1014     p-=integrate(p)/vol(self.dom)
1015     error_p=Lsup(mod.getPressure()-p_ref)/Lsup(p_ref)
1016     error_s=Lsup(mod.getDeviatoricStress()-s_ref)/Lsup(s_ref)
1017     error_v=Lsup(mod.getVelocity()-v_ref)/Lsup(v_ref)
1018     error_t=abs(mod.getTime()-t_ref)/abs(t_ref)
1019     if self.VERBOSE: print "time step ",N_t,"time = ",mod.getTime(),"errors s,p,v = ",error_s, error_p, error_v
1020     self.failUnless( error_p <= 10*self.TOL, "time step %s: pressure error %s too high."%(N_t,error_p) )
1021     self.failUnless( error_s <= 10*self.TOL, "time step %s: stress error %s too high."%(N_t,error_s) )
1022     self.failUnless( error_v <= 10*self.TOL, "time step %s: velocity error %s too high."%(N_t,error_v) )
1023     self.failUnless( error_t <= 10*self.TOL, "time step %s: time marker error %s too high."%(N_t,error_t) )
1024     def tearDown(self):
1025     del self.dom
1026    
1027     def test_D2_Fixed_MuNone_LateStart(self):
1028     self.dom = Rectangle(self.NE,self.NE,order=2)
1029     self.mu=None
1030     self.latestart=True
1031     self.runIt()
1032     def test_D2_Fixed_Mu_LateStart(self):
1033     self.dom = Rectangle(self.NE,self.NE,order=2)
1034     self.mu=555.
1035     self.latestart=True
1036     self.runIt()
1037     def test_D2_Fixed_MuNone(self):
1038     self.dom = Rectangle(self.NE,self.NE,order=2)
1039     self.mu=None
1040     self.latestart=False
1041     self.runIt()
1042     def test_D2_Fixed_Mu(self):
1043     self.dom = Rectangle(self.NE,self.NE,order=2)
1044     self.mu=555.
1045     self.latestart=False
1046     self.runIt()
1047     def test_D2_Free_MuNone_LateStart(self):
1048     self.dom = Rectangle(self.NE,self.NE,order=2)
1049     self.mu=None
1050     self.latestart=True
1051     self.runIt(free=0)
1052     def test_D2_Free_Mu_LateStart(self):
1053     self.dom = Rectangle(self.NE,self.NE,order=2)
1054     self.mu=555.
1055     self.latestart=True
1056     self.runIt(free=0)
1057     def test_D2_Free_MuNone(self):
1058     self.dom = Rectangle(self.NE,self.NE,order=2)
1059     self.mu=None
1060     self.latestart=False
1061     self.runIt(free=0)
1062     def test_D2_Free_Mu(self):
1063     self.dom = Rectangle(self.NE,self.NE,order=2)
1064     self.mu=555.
1065     self.latestart=False
1066     self.runIt(free=0)
1067    
1068     def test_D3_Fixed_MuNone_LateStart(self):
1069     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1070     self.mu=None
1071     self.latestart=True
1072     self.runIt()
1073     def test_D3_Fixed_Mu_LateStart(self):
1074     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1075     self.mu=555.
1076     self.latestart=True
1077     self.runIt()
1078     def test_D3_Fixed_MuNone(self):
1079     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1080     self.mu=None
1081     self.latestart=False
1082     self.runIt()
1083     def test_D3_Fixed_Mu(self):
1084     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1085     self.mu=555.
1086     self.latestart=False
1087     self.runIt()
1088     def test_D3_Free_MuNone_LateStart(self):
1089     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1090     self.mu=None
1091     self.latestart=True
1092     self.runIt(free=0)
1093     def test_D3_Free_Mu_LateStart(self):
1094     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1095     self.mu=555.
1096     self.latestart=True
1097     self.runIt(free=0)
1098     def test_D3_Free_MuNone(self):
1099     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1100     self.mu=None
1101     self.latestart=False
1102     self.runIt(free=0)
1103     def test_D3_Free_Mu(self):
1104     self.dom = Brick(self.NE,self.NE,self.NE,order=2)
1105     self.mu=555.
1106     self.latestart=False
1107     self.runIt(free=0)
1108    
1109 gross 1414 if __name__ == '__main__':
1110     suite = unittest.TestSuite()
1111 gross 2100 suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian2D))
1112     suite.addTest(unittest.makeSuite(Test_Darcy3D))
1113     suite.addTest(unittest.makeSuite(Test_Darcy2D))
1114     suite.addTest(unittest.makeSuite(Test_StokesProblemCartesian3D))
1115 artak 2234 suite.addTest(unittest.makeSuite(Test_Mountains3D))
1116     suite.addTest(unittest.makeSuite(Test_Mountains2D))
1117 gross 2301 suite.addTest(unittest.makeSuite(Test_Rheologies))
1118 gross 2432 suite.addTest(unittest.makeSuite(Test_IncompressibleIsotropicFlowCartesian))
1119 gross 1414 s=unittest.TextTestRunner(verbosity=2).run(suite)
1120     if not s.wasSuccessful(): sys.exit(1)
1121    

  ViewVC Help
Powered by ViewVC 1.1.26