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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3909 - (hide annotations)
Wed Jun 13 07:03:06 2012 UTC (7 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 14237 byte(s)
darcy solver removes hydrostatic pressure from problem now
1 gross 3503 # -*- 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 jfenwick 3772 for i in range(self.dom.getDim()):
65 gross 3503 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 jfenwick 3772 for i in range(self.dom.getDim()):
83 gross 3503 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 gross 3569 v,p=df.solve(u_ref,p)
111 gross 3503
112    
113 jfenwick 3551 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 gross 3503 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 gross 3569 v,p=df.solve(u,p )
128 gross 3503
129    
130 jfenwick 3551 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 gross 3503
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 gross 3569 v,p=df.solve(u,p)
146     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
147 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
148 gross 3503
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 gross 3569 v,p=df.solve(u,p)
163 gross 3503
164 gross 3569 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
165 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
166 gross 3503
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 gross 3569 v,p=df.solve(u,p)
180 gross 3503
181 jfenwick 3551 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 gross 3503
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 gross 3569 v,p=df.solve(u,p)
198 gross 3503
199 jfenwick 3551 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 gross 3503
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 gross 3569 v,p=df.solve(u,p)
216 gross 3503
217 jfenwick 3551 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 gross 3503
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 gross 3569 v,p=df.solve(u,p)
233 gross 3503
234    
235 jfenwick 3551 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 gross 3503
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 gross 3793 #df.getSolverOptionsPressure().setVerbosityOn()
247 gross 3503 df.setValue(g=f,
248     location_of_fixed_pressure=mp,
249     location_of_fixed_flux=mv,
250     permeability=Scalar(k,Function(self.dom)))
251 gross 3569 v,p=df.solve(u,p)
252     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
253 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
254 gross 3503
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 gross 3569 v,p=df.solve(u,p)
268 gross 3503
269 gross 3569 self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
270 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
271 gross 3503
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 gross 3569 v,p=df.solve(u,p)
285 jfenwick 3551 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 gross 3503
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 gross 3569 v,p=df.solve(u,p)
301 jfenwick 3551 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 gross 3503
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 jfenwick 3892 self.dom = Rectangle(NE,NE)
311 gross 3503 self.rescaleDomain()
312     def tearDown(self):
313     del self.dom
314    
315 gross 3569 class Test_Darcy2D_EVAL(Test_Darcy2D):
316 gross 3503 TEST_TOL=0.01
317 gross 3569 SOLVER=DarcyFlow.EVAL
318 gross 3503
319     class Test_Darcy2D_POST(Test_Darcy2D):
320     TEST_TOL=1.e-3
321     SOLVER=DarcyFlow.POST
322    
323 gross 3569 class Test_Darcy2D_SMOOTH(Test_Darcy2D):
324     TEST_TOL=0.01
325     SOLVER=DarcyFlow.SMOOTH
326 gross 3503
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 gross 3569 class Test_Darcy3D_EVAL(Test_Darcy3D):
338 gross 3503 TEST_TOL=0.01
339 gross 3569 SOLVER=DarcyFlow.EVAL
340 gross 3503
341     class Test_Darcy3D_POST(Test_Darcy3D):
342     TEST_TOL=1.e-3
343     SOLVER=DarcyFlow.POST
344    
345 gross 3569 class Test_Darcy3D_SMOOTH(Test_Darcy3D):
346     TEST_TOL=0.01
347     SOLVER=DarcyFlow.SMOOTH
348 gross 3503
349    
350     if __name__ == '__main__':
351     suite = unittest.TestSuite()
352 gross 3569 suite.addTest(unittest.makeSuite(Test_Darcy2D_SMOOTH))
353 gross 3503 suite.addTest(unittest.makeSuite(Test_Darcy2D_POST))
354 gross 3569 suite.addTest(unittest.makeSuite(Test_Darcy2D_EVAL))
355    
356     suite.addTest(unittest.makeSuite(Test_Darcy3D_SMOOTH))
357 gross 3503 suite.addTest(unittest.makeSuite(Test_Darcy3D_POST))
358 gross 3569 suite.addTest(unittest.makeSuite(Test_Darcy3D_EVAL))
359 gross 3909 #suite.addTest(Test_Darcy2D_POST("testConstF_FixedBottom_largeK"))
360 gross 3503 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