/[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 4154 - (show annotations)
Tue Jan 22 09:30:23 2013 UTC (6 years, 9 months ago) by jfenwick
File MIME type: text/x-python
File size: 14299 byte(s)
Round 1 of copyright fixes
1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2013 by University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development since 2012 by School of Earth Sciences
14 #
15 ##############################################################################
16
17 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
18 http://www.uq.edu.au
19 Primary Business: Queensland, Australia"""
20 __license__="""Licensed under the Open Software License version 3.0
21 http://www.opensource.org/licenses/osl-3.0.php"""
22 __url__="https://launchpad.net/escript-finley"
23
24 import unittest
25 import tempfile
26
27
28
29 VERBOSE=False and True
30
31 from esys.escript import *
32 from esys.escript.models import DarcyFlow
33 from esys.finley import Rectangle, Brick
34
35 from math import pi
36 import numpy
37 import sys
38 import os
39 #====================================================================================================================
40 try:
41 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
42 except KeyError:
43 FINLEY_WORKDIR='.'
44
45 class Test_Darcy(unittest.TestCase):
46 # this is a simple test for the darcy flux problem
47 #
48 #
49 # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
50 #
51 # with f = (fb-f0)/xb* x + f0
52 #
53 # u = f - k * p,x = ub
54 #
55 # we prescribe pressure at x=x0=0 to p0
56 #
57 # if we prescribe the pressure on the bottom x=xb we set
58 #
59 # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
60 #
61 # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
62 #
63 def rescaleDomain(self):
64 x=self.dom.getX().copy()
65 for i in range(self.dom.getDim()):
66 x_inf=inf(x[i])
67 x_sup=sup(x[i])
68 if i == self.dom.getDim()-1:
69 x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
70 else:
71 x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
72 self.dom.setX(x)
73 def getScalarMask(self,include_bottom=True):
74 x=self.dom.getX().copy()
75 x_inf=inf(x[self.dom.getDim()-1])
76 x_sup=sup(x[self.dom.getDim()-1])
77 out=whereZero(x[self.dom.getDim()-1]-x_sup)
78 if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
79 return wherePositive(out)
80 def getVectorMask(self,include_bottom=True):
81 x=self.dom.getX().copy()
82 out=Vector(0.,Solution(self.dom))
83 for i in range(self.dom.getDim()):
84 x_inf=inf(x[i])
85 x_sup=sup(x[i])
86 if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
87 if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
88 return wherePositive(out)
89
90 def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
91 d=self.dom.getDim()
92 x=self.dom.getX()[d-1]
93 xb=inf(x)
94 u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
95 p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
96 f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
97 return u,p,f
98
99 def testConstF_FixedBottom_smallK(self):
100 k=1.e-8
101 mp=self.getScalarMask(include_bottom=True)
102 mv=self.getVectorMask(include_bottom=False)
103 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
104 p=p_ref*mp
105 u=u_ref*mv
106 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
107 df.setValue(g=f,
108 location_of_fixed_pressure=mp,
109 location_of_fixed_flux=mv,
110 permeability=Scalar(k,Function(self.dom)))
111 v,p=df.solve(u_ref,p)
112
113
114 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
115 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
116 def testConstF_FixedBottom_mediumK(self):
117 k=1.
118 mp=self.getScalarMask(include_bottom=True)
119 mv=self.getVectorMask(include_bottom=False)
120 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
121 p=p_ref*mp
122 u=u_ref*mv
123 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
124 df.setValue(g=f,
125 location_of_fixed_pressure=mp,
126 location_of_fixed_flux=mv,
127 permeability=Scalar(k,Function(self.dom)))
128 v,p=df.solve(u,p )
129
130
131 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
132 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
133
134 def testConstF_FixedBottom_largeK(self):
135 k=1.e8
136 mp=self.getScalarMask(include_bottom=True)
137 mv=self.getVectorMask(include_bottom=False)
138 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
139 p=p_ref*mp
140 u=u_ref*mv
141 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
142 df.setValue(g=f,
143 location_of_fixed_pressure=mp,
144 location_of_fixed_flux=mv,
145 permeability=Scalar(k,Function(self.dom)))
146 v,p=df.solve(u,p)
147 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
148 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
149
150 def testVarioF_FixedBottom_smallK(self):
151 k=1.e-8
152 mp=self.getScalarMask(include_bottom=True)
153 mv=self.getVectorMask(include_bottom=False)
154 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
155 p=p_ref*mp
156 u=u_ref*mv
157 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
158 #df.getSolverOptionsPressure().setVerbosityOn()
159 df.setValue(g=f,
160 location_of_fixed_pressure=mp,
161 location_of_fixed_flux=mv,
162 permeability=Scalar(k,Function(self.dom)))
163 v,p=df.solve(u,p)
164
165 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
166 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
167
168 def testVarioF_FixedBottom_mediumK(self):
169 k=1.
170 mp=self.getScalarMask(include_bottom=True)
171 mv=self.getVectorMask(include_bottom=False)
172 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
173 p=p_ref *mp
174 u=u_ref *mv
175 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
176 df.setValue(g=f,
177 location_of_fixed_pressure=mp,
178 location_of_fixed_flux=mv,
179 permeability=Scalar(k,Function(self.dom)))
180 v,p=df.solve(u,p)
181
182 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
183 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
184
185 def testVarioF_FixedBottom_largeK(self):
186 k=1.e8
187 mp=self.getScalarMask(include_bottom=True)
188 mv=self.getVectorMask(include_bottom=False)
189 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
190 p=p_ref*mp
191
192 u=u_ref*mv
193 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
194 df.setValue(g=f,
195 location_of_fixed_pressure=mp,
196 location_of_fixed_flux=mv,
197 permeability=Scalar(k,Function(self.dom)))
198 v,p=df.solve(u,p)
199
200 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
201 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
202
203 def testConstF_FreeBottom_smallK(self):
204 k=1.e-8
205 mp=self.getScalarMask(include_bottom=False)
206 mv=self.getVectorMask(include_bottom=True)
207 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
208 p=p_ref # *mp
209 u=u_ref*mv
210 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
211 df.setValue(g=f,
212 location_of_fixed_pressure=mp,
213 location_of_fixed_flux=mv,
214 permeability=Scalar(k,Function(self.dom)))
215
216 v,p=df.solve(u,p)
217
218 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
219 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
220
221 def testConstF_FreeBottom_mediumK(self):
222 k=1.
223 mp=self.getScalarMask(include_bottom=False)
224 mv=self.getVectorMask(include_bottom=True)
225 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
226 p=p_ref *mp
227 u=u_ref*mv
228 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
229 df.setValue(g=f,
230 location_of_fixed_pressure=mp,
231 location_of_fixed_flux=mv,
232 permeability=Scalar(k,Function(self.dom)))
233 v,p=df.solve(u,p)
234
235
236 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
237 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
238
239 def testConstF_FreeBottom_largeK(self):
240 k=1.e8
241 mp=self.getScalarMask(include_bottom=False)
242 mv=self.getVectorMask(include_bottom=True)
243 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
244 p=p_ref*mp
245 u=u_ref*mv
246 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
247 #df.getSolverOptionsPressure().setVerbosityOn()
248 df.setValue(g=f,
249 location_of_fixed_pressure=mp,
250 location_of_fixed_flux=mv,
251 permeability=Scalar(k,Function(self.dom)))
252 v,p=df.solve(u,p)
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 #suite.addTest(Test_Darcy2D_POST("testConstF_FixedBottom_largeK"))
361 s=unittest.TextTestRunner(verbosity=2).run(suite)
362 if not s.wasSuccessful(): sys.exit(1)
363
364

  ViewVC Help
Powered by ViewVC 1.1.26