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

Contents of /trunk/finley/test/python/run_darcy.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3569 - (show annotations)
Thu Sep 1 02:42:36 2011 UTC (8 years, 1 month ago) by gross
File MIME type: text/x-python
File size: 14181 byte(s)
Darcy flow solver and docu reviewed. Coupled solvers have been removed.

1 # -*- coding: utf-8 -*-
2
3 ########################################################
4 #
5 # Copyright (c) 2003-2010 by University of Queensland
6 # Earth Systems Science Computational Center (ESSCC)
7 # http://www.uq.edu.au/esscc
8 #
9 # Primary Business: Queensland, Australia
10 # Licensed under the Open Software License version 3.0
11 # http://www.opensource.org/licenses/osl-3.0.php
12 #
13 ########################################################
14
15 __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16 Earth Systems Science Computational Center (ESSCC)
17 http://www.uq.edu.au/esscc
18 Primary Business: Queensland, Australia"""
19 __license__="""Licensed under the Open Software License version 3.0
20 http://www.opensource.org/licenses/osl-3.0.php"""
21 __url__="https://launchpad.net/escript-finley"
22
23 import unittest
24 import tempfile
25
26
27
28 VERBOSE=False and True
29
30 from esys.escript import *
31 from esys.escript.models import DarcyFlow
32 from esys.finley import Rectangle, Brick
33
34 from math import pi
35 import numpy
36 import sys
37 import os
38 #====================================================================================================================
39 try:
40 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
41 except KeyError:
42 FINLEY_WORKDIR='.'
43
44 class Test_Darcy(unittest.TestCase):
45 # this is a simple test for the darcy flux problem
46 #
47 #
48 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
49 #
50 # with f = (fb-f0)/xb* x + f0
51 #
52 # u = f - k * p,x = ub
53 #
54 # we prescribe pressure at x=x0=0 to p0
55 #
56 # if we prescribe the pressure on the bottom x=xb we set
57 #
58 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
59 #
60 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
61 #
62 def rescaleDomain(self):
63 x=self.dom.getX().copy()
64 for i in xrange(self.dom.getDim()):
65 x_inf=inf(x[i])
66 x_sup=sup(x[i])
67 if i == self.dom.getDim()-1:
68 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
69 else:
70 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
71 self.dom.setX(x)
72 def getScalarMask(self,include_bottom=True):
73 x=self.dom.getX().copy()
74 x_inf=inf(x[self.dom.getDim()-1])
75 x_sup=sup(x[self.dom.getDim()-1])
76 out=whereZero(x[self.dom.getDim()-1]-x_sup)
77 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
78 return wherePositive(out)
79 def getVectorMask(self,include_bottom=True):
80 x=self.dom.getX().copy()
81 out=Vector(0.,Solution(self.dom))
82 for i in xrange(self.dom.getDim()):
83 x_inf=inf(x[i])
84 x_sup=sup(x[i])
85 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
86 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
87 return wherePositive(out)
88
89 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
90 d=self.dom.getDim()
91 x=self.dom.getX()[d-1]
92 xb=inf(x)
93 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
94 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
95 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
96 return u,p,f
97
98 def testConstF_FixedBottom_smallK(self):
99 k=1.e-8
100 mp=self.getScalarMask(include_bottom=True)
101 mv=self.getVectorMask(include_bottom=False)
102 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
103 p=p_ref*mp
104 u=u_ref*mv
105 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
106 df.setValue(g=f,
107 location_of_fixed_pressure=mp,
108 location_of_fixed_flux=mv,
109 permeability=Scalar(k,Function(self.dom)))
110 v,p=df.solve(u_ref,p)
111
112
113 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
114 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
115 def testConstF_FixedBottom_mediumK(self):
116 k=1.
117 mp=self.getScalarMask(include_bottom=True)
118 mv=self.getVectorMask(include_bottom=False)
119 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
120 p=p_ref*mp
121 u=u_ref*mv
122 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
123 df.setValue(g=f,
124 location_of_fixed_pressure=mp,
125 location_of_fixed_flux=mv,
126 permeability=Scalar(k,Function(self.dom)))
127 v,p=df.solve(u,p )
128
129
130 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
131 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
132
133 def testConstF_FixedBottom_largeK(self):
134 k=1.e8
135 mp=self.getScalarMask(include_bottom=True)
136 mv=self.getVectorMask(include_bottom=False)
137 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
138 p=p_ref*mp
139 u=u_ref*mv
140 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
141 df.setValue(g=f,
142 location_of_fixed_pressure=mp,
143 location_of_fixed_flux=mv,
144 permeability=Scalar(k,Function(self.dom)))
145 v,p=df.solve(u,p)
146 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
147 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
148
149 def testVarioF_FixedBottom_smallK(self):
150 k=1.e-8
151 mp=self.getScalarMask(include_bottom=True)
152 mv=self.getVectorMask(include_bottom=False)
153 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
154 p=p_ref*mp
155 u=u_ref*mv
156 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
157 #df.getSolverOptionsPressure().setVerbosityOn()
158 df.setValue(g=f,
159 location_of_fixed_pressure=mp,
160 location_of_fixed_flux=mv,
161 permeability=Scalar(k,Function(self.dom)))
162 v,p=df.solve(u,p)
163
164 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
165 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
166
167 def testVarioF_FixedBottom_mediumK(self):
168 k=1.
169 mp=self.getScalarMask(include_bottom=True)
170 mv=self.getVectorMask(include_bottom=False)
171 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
172 p=p_ref *mp
173 u=u_ref *mv
174 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
175 df.setValue(g=f,
176 location_of_fixed_pressure=mp,
177 location_of_fixed_flux=mv,
178 permeability=Scalar(k,Function(self.dom)))
179 v,p=df.solve(u,p)
180
181 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
182 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
183
184 def testVarioF_FixedBottom_largeK(self):
185 k=1.e8
186 mp=self.getScalarMask(include_bottom=True)
187 mv=self.getVectorMask(include_bottom=False)
188 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
189 p=p_ref*mp
190
191 u=u_ref*mv
192 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
193 df.setValue(g=f,
194 location_of_fixed_pressure=mp,
195 location_of_fixed_flux=mv,
196 permeability=Scalar(k,Function(self.dom)))
197 v,p=df.solve(u,p)
198
199 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
200 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
201
202 def testConstF_FreeBottom_smallK(self):
203 k=1.e-8
204 mp=self.getScalarMask(include_bottom=False)
205 mv=self.getVectorMask(include_bottom=True)
206 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
207 p=p_ref # *mp
208 u=u_ref*mv
209 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
210 df.setValue(g=f,
211 location_of_fixed_pressure=mp,
212 location_of_fixed_flux=mv,
213 permeability=Scalar(k,Function(self.dom)))
214
215 v,p=df.solve(u,p)
216
217 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
218 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
219
220 def testConstF_FreeBottom_mediumK(self):
221 k=1.
222 mp=self.getScalarMask(include_bottom=False)
223 mv=self.getVectorMask(include_bottom=True)
224 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
225 p=p_ref *mp
226 u=u_ref*mv
227 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
228 df.setValue(g=f,
229 location_of_fixed_pressure=mp,
230 location_of_fixed_flux=mv,
231 permeability=Scalar(k,Function(self.dom)))
232 v,p=df.solve(u,p)
233
234
235 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
236 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
237
238 def testConstF_FreeBottom_largeK(self):
239 k=1.e8
240 mp=self.getScalarMask(include_bottom=False)
241 mv=self.getVectorMask(include_bottom=True)
242 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
243 p=p_ref*mp
244 u=u_ref*mv
245 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
246 df.setValue(g=f,
247 location_of_fixed_pressure=mp,
248 location_of_fixed_flux=mv,
249 permeability=Scalar(k,Function(self.dom)))
250 v,p=df.solve(u,p)
251
252
253 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
254 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
255
256 def testVarioF_FreeBottom_smallK(self):
257 k=1.e-8
258 mp=self.getScalarMask(include_bottom=False)
259 mv=self.getVectorMask(include_bottom=True)
260 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
261 p=p_ref*mp
262 u=u_ref*mv
263 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
264 df.setValue(g=f,
265 location_of_fixed_pressure=mp,
266 location_of_fixed_flux=mv,
267 permeability=Scalar(k,Function(self.dom)))
268 v,p=df.solve(u,p)
269
270 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
271 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
272
273 def testVarioF_FreeBottom_mediumK(self):
274 k=1.
275 mp=self.getScalarMask(include_bottom=False)
276 mv=self.getVectorMask(include_bottom=True)
277 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
278 p=p_ref*mp
279 u=u_ref*mv
280 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
281 df.setValue(g=f,
282 location_of_fixed_pressure=mp,
283 location_of_fixed_flux=mv,
284 permeability=Scalar(k,Function(self.dom)))
285 v,p=df.solve(u,p)
286 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
287 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
288
289 def testVarioF_FreeBottom_largeK(self):
290 k=1.e8
291 mp=self.getScalarMask(include_bottom=False)
292 mv=self.getVectorMask(include_bottom=True)
293 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
294 p=p_ref*mp
295 u=u_ref*mv
296 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
297 df.setValue(g=f,
298 location_of_fixed_pressure=mp,
299 location_of_fixed_flux=mv,
300 permeability=Scalar(k,Function(self.dom)))
301 v,p=df.solve(u,p)
302 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
303 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
304
305 class Test_Darcy2D(Test_Darcy):
306 TOL=1e-4
307 WIDTH=1.
308 SOLVER=DarcyFlow.POST
309 def setUp(self):
310 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
311 self.dom = Rectangle(NE,NE)
312 self.rescaleDomain()
313 def tearDown(self):
314 del self.dom
315
316 class Test_Darcy2D_EVAL(Test_Darcy2D):
317 TEST_TOL=0.01
318 SOLVER=DarcyFlow.EVAL
319
320 class Test_Darcy2D_POST(Test_Darcy2D):
321 TEST_TOL=1.e-3
322 SOLVER=DarcyFlow.POST
323
324 class Test_Darcy2D_SMOOTH(Test_Darcy2D):
325 TEST_TOL=0.01
326 SOLVER=DarcyFlow.SMOOTH
327
328 class Test_Darcy3D(Test_Darcy):
329 WIDTH=1.
330 SOLVER=DarcyFlow.POST
331 def setUp(self):
332 NE=29 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
333 self.dom = Brick(NE,NE,NE)
334 self.rescaleDomain()
335 def tearDown(self):
336 del self.dom
337
338 class Test_Darcy3D_EVAL(Test_Darcy3D):
339 TEST_TOL=0.01
340 SOLVER=DarcyFlow.EVAL
341
342 class Test_Darcy3D_POST(Test_Darcy3D):
343 TEST_TOL=1.e-3
344 SOLVER=DarcyFlow.POST
345
346 class Test_Darcy3D_SMOOTH(Test_Darcy3D):
347 TEST_TOL=0.01
348 SOLVER=DarcyFlow.SMOOTH
349
350
351 if __name__ == '__main__':
352 suite = unittest.TestSuite()
353 suite.addTest(unittest.makeSuite(Test_Darcy2D_SMOOTH))
354 suite.addTest(unittest.makeSuite(Test_Darcy2D_POST))
355 suite.addTest(unittest.makeSuite(Test_Darcy2D_EVAL))
356
357 suite.addTest(unittest.makeSuite(Test_Darcy3D_SMOOTH))
358 suite.addTest(unittest.makeSuite(Test_Darcy3D_POST))
359 suite.addTest(unittest.makeSuite(Test_Darcy3D_EVAL))
360
361
362 #suite.addTest(Test_Darcy2D_SMOOTH("testConstF_FreeBottom_largeK"))
363 s=unittest.TextTestRunner(verbosity=2).run(suite)
364 if not s.wasSuccessful(): sys.exit(1)
365
366

  ViewVC Help
Powered by ViewVC 1.1.26