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

  ViewVC Help
Powered by ViewVC 1.1.26