/[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 3793 - (show annotations)
Wed Feb 1 07:39:43 2012 UTC (7 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 14231 byte(s)
new implementation of FCT solver with some modifications to the python interface
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 range(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 range(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.getSolverOptionsPressure().setVerbosityOn()
247 df.setValue(g=f,
248 location_of_fixed_pressure=mp,
249 location_of_fixed_flux=mv,
250 permeability=Scalar(k,Function(self.dom)))
251 v,p=df.solve(u,p)
252 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
253 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
254
255 def testVarioF_FreeBottom_smallK(self):
256 k=1.e-8
257 mp=self.getScalarMask(include_bottom=False)
258 mv=self.getVectorMask(include_bottom=True)
259 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
260 p=p_ref*mp
261 u=u_ref*mv
262 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
263 df.setValue(g=f,
264 location_of_fixed_pressure=mp,
265 location_of_fixed_flux=mv,
266 permeability=Scalar(k,Function(self.dom)))
267 v,p=df.solve(u,p)
268
269 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
270 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
271
272 def testVarioF_FreeBottom_mediumK(self):
273 k=1.
274 mp=self.getScalarMask(include_bottom=False)
275 mv=self.getVectorMask(include_bottom=True)
276 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
277 p=p_ref*mp
278 u=u_ref*mv
279 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
280 df.setValue(g=f,
281 location_of_fixed_pressure=mp,
282 location_of_fixed_flux=mv,
283 permeability=Scalar(k,Function(self.dom)))
284 v,p=df.solve(u,p)
285 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
286 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
287
288 def testVarioF_FreeBottom_largeK(self):
289 k=1.e8
290 mp=self.getScalarMask(include_bottom=False)
291 mv=self.getVectorMask(include_bottom=True)
292 u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=30.,k=k)
293 p=p_ref*mp
294 u=u_ref*mv
295 df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
296 df.setValue(g=f,
297 location_of_fixed_pressure=mp,
298 location_of_fixed_flux=mv,
299 permeability=Scalar(k,Function(self.dom)))
300 v,p=df.solve(u,p)
301 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
302 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
303
304 class Test_Darcy2D(Test_Darcy):
305 TOL=1e-4
306 WIDTH=1.
307 SOLVER=DarcyFlow.POST
308 def setUp(self):
309 NE=40 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
310 self.dom = Rectangle(NE,NE)
311 self.rescaleDomain()
312 def tearDown(self):
313 del self.dom
314
315 class Test_Darcy2D_EVAL(Test_Darcy2D):
316 TEST_TOL=0.01
317 SOLVER=DarcyFlow.EVAL
318
319 class Test_Darcy2D_POST(Test_Darcy2D):
320 TEST_TOL=1.e-3
321 SOLVER=DarcyFlow.POST
322
323 class Test_Darcy2D_SMOOTH(Test_Darcy2D):
324 TEST_TOL=0.01
325 SOLVER=DarcyFlow.SMOOTH
326
327 class Test_Darcy3D(Test_Darcy):
328 WIDTH=1.
329 SOLVER=DarcyFlow.POST
330 def setUp(self):
331 NE=29 # wrning smaller NE may case a failure for VarioF tests due to discretization errors.
332 self.dom = Brick(NE,NE,NE)
333 self.rescaleDomain()
334 def tearDown(self):
335 del self.dom
336
337 class Test_Darcy3D_EVAL(Test_Darcy3D):
338 TEST_TOL=0.01
339 SOLVER=DarcyFlow.EVAL
340
341 class Test_Darcy3D_POST(Test_Darcy3D):
342 TEST_TOL=1.e-3
343 SOLVER=DarcyFlow.POST
344
345 class Test_Darcy3D_SMOOTH(Test_Darcy3D):
346 TEST_TOL=0.01
347 SOLVER=DarcyFlow.SMOOTH
348
349
350 if __name__ == '__main__':
351 suite = unittest.TestSuite()
352 suite.addTest(unittest.makeSuite(Test_Darcy2D_SMOOTH))
353 suite.addTest(unittest.makeSuite(Test_Darcy2D_POST))
354 suite.addTest(unittest.makeSuite(Test_Darcy2D_EVAL))
355
356 suite.addTest(unittest.makeSuite(Test_Darcy3D_SMOOTH))
357 suite.addTest(unittest.makeSuite(Test_Darcy3D_POST))
358 suite.addTest(unittest.makeSuite(Test_Darcy3D_EVAL))
359 #suite.addTest(Test_Darcy2D_SMOOTH("testConstF_FreeBottom_largeK"))
360 s=unittest.TextTestRunner(verbosity=2).run(suite)
361 if not s.wasSuccessful(): sys.exit(1)
362
363

  ViewVC Help
Powered by ViewVC 1.1.26