/[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 6651 - (hide annotations)
Wed Feb 7 02:12:08 2018 UTC (21 months, 1 week ago) by jfenwick
File MIME type: text/x-python
File size: 13954 byte(s)
Make everyone sad by touching all the files

Copyright dates update

1 gross 3503 # -*- coding: utf-8 -*-
2    
3 jfenwick 3981 ##############################################################################
4 gross 3503 #
5 jfenwick 6651 # Copyright (c) 2003-2018 by The University of Queensland
6 jfenwick 3981 # http://www.uq.edu.au
7 gross 3503 #
8     # Primary Business: Queensland, Australia
9 jfenwick 6112 # Licensed under the Apache License, version 2.0
10     # http://www.apache.org/licenses/LICENSE-2.0
11 gross 3503 #
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 sshaw 5706 from __future__ import print_function, division
19    
20 jfenwick 6651 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
21 jfenwick 3981 http://www.uq.edu.au
22 gross 3503 Primary Business: Queensland, Australia"""
23 jfenwick 6112 __license__="""Licensed under the Apache License, version 2.0
24     http://www.apache.org/licenses/LICENSE-2.0"""
25 gross 3503 __url__="https://launchpad.net/escript-finley"
26    
27 jfenwick 4938 import esys.escriptcore.utestselect as unittest
28 sshaw 4984 from esys.escriptcore.testing import *
29 gross 3503
30     from esys.escript import *
31 sshaw 4984 from esys.escript.models import DarcyFlow
32 gross 3503 from esys.finley import Rectangle, Brick
33    
34     import os
35 caltinay 6282
36     VERBOSE=False
37    
38 gross 3503 try:
39     FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
40     except KeyError:
41     FINLEY_WORKDIR='.'
42    
43 sshaw 5059 class Darcy(unittest.TestCase): #subclassing required
44 gross 3503 # this is a simple test for the darcy flux problem
45     #
46     #
47     # p = 1/k * ( 1/2* (fb-f0)/xb* x **2 + f0 * x - ub*x ) + p0
48     #
49     # with f = (fb-f0)/xb* x + f0
50     #
51     # u = f - k * p,x = ub
52     #
53     # we prescribe pressure at x=x0=0 to p0
54     #
55     # if we prescribe the pressure on the bottom x=xb we set
56     #
57     # pb= 1/k * ( 1/2* (fb-f0)* xb + f0 * xb - ub*xb ) + p0 = 1/k * ((fb+f0)/2 - ub ) * xb + p0
58     #
59     # which leads to ub = (fb+f0)/2-k*(pb-p0)/xb
60     #
61     def rescaleDomain(self):
62     x=self.dom.getX().copy()
63 jfenwick 3772 for i in range(self.dom.getDim()):
64 gross 3503 x_inf=inf(x[i])
65     x_sup=sup(x[i])
66     if i == self.dom.getDim()-1:
67     x[i]=-self.WIDTH*(x[i]-x_sup)/(x_inf-x_sup)
68     else:
69     x[i]=self.WIDTH*(x[i]-x_inf)/(x_sup-x_inf)
70     self.dom.setX(x)
71     def getScalarMask(self,include_bottom=True):
72     x=self.dom.getX().copy()
73     x_inf=inf(x[self.dom.getDim()-1])
74     x_sup=sup(x[self.dom.getDim()-1])
75     out=whereZero(x[self.dom.getDim()-1]-x_sup)
76     if include_bottom: out+=whereZero(x[self.dom.getDim()-1]-x_inf)
77     return wherePositive(out)
78     def getVectorMask(self,include_bottom=True):
79     x=self.dom.getX().copy()
80     out=Vector(0.,Solution(self.dom))
81 jfenwick 3772 for i in range(self.dom.getDim()):
82 gross 3503 x_inf=inf(x[i])
83     x_sup=sup(x[i])
84     if i != self.dom.getDim()-1: out[i]+=whereZero(x[i]-x_sup)
85     if i != self.dom.getDim()-1 or include_bottom: out[i]+=whereZero(x[i]-x_inf)
86     return wherePositive(out)
87    
88     def setSolutionFixedBottom(self, p0, pb, f0, fb, k):
89     d=self.dom.getDim()
90     x=self.dom.getX()[d-1]
91     xb=inf(x)
92     u=Vector(0.,Solution(self.dom))+kronecker(d)[d-1]*((f0+fb)/2.-k*(pb-p0)/xb)
93     p=1./k*((fb-f0)/(xb*2.)* x**2 - (fb-f0)/2.*x)+(pb-p0)/xb*x + p0
94     f= ((fb-f0)/xb* x + f0)*kronecker(Function(self.dom))[d-1]
95     return u,p,f
96    
97     def testConstF_FixedBottom_smallK(self):
98     k=1.e-8
99     mp=self.getScalarMask(include_bottom=True)
100     mv=self.getVectorMask(include_bottom=False)
101     u_ref,p_ref,f=self.setSolutionFixedBottom(p0=2500,pb=4000.,f0=10.,fb=10.,k=k)
102     p=p_ref*mp
103     u=u_ref*mv
104     df=DarcyFlow(self.dom, verbose=VERBOSE, solver=self.SOLVER)
105     df.setValue(g=f,
106     location_of_fixed_pressure=mp,
107     location_of_fixed_flux=mv,
108     permeability=Scalar(k,Function(self.dom)))
109 gross 3569 v,p=df.solve(u_ref,p)
110 gross 3503
111    
112 jfenwick 3551 self.assertTrue(Lsup(v-u_ref)<self.TEST_TOL*Lsup(u_ref), "flux error too big.")
113     self.assertTrue(Lsup(p-p_ref)<self.TEST_TOL*Lsup(p_ref), "pressure error too big.")
114 sshaw 5059
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 sshaw 5059 class Darcy2D(Darcy): #subclassing required
305 gross 3503 TOL=1e-4
306     WIDTH=1.
307     SOLVER=DarcyFlow.POST
308     def setUp(self):
309 caltinay 6282 NE=40 # warning: 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 sshaw 4984
315     class Test_Darcy2D_EVAL(Darcy2D):
316 gross 3503 TEST_TOL=0.01
317 gross 3569 SOLVER=DarcyFlow.EVAL
318 gross 3503
319 sshaw 4984 class Test_Darcy2D_POST(Darcy2D):
320 gross 3503 TEST_TOL=1.e-3
321     SOLVER=DarcyFlow.POST
322    
323 sshaw 4984 class Test_Darcy2D_SMOOTH(Darcy2D):
324 gross 3569 TEST_TOL=0.01
325     SOLVER=DarcyFlow.SMOOTH
326 gross 3503
327 sshaw 5059 class Darcy3D(Darcy): #subclassing required
328 gross 3503 WIDTH=1.
329     SOLVER=DarcyFlow.POST
330     def setUp(self):
331 caltinay 6282 NE=29 # warning: smaller NE may case a failure for VarioF tests due to discretization errors.
332 gross 3503 self.dom = Brick(NE,NE,NE)
333     self.rescaleDomain()
334     def tearDown(self):
335     del self.dom
336    
337 sshaw 4984 class Test_Darcy3D_EVAL(Darcy3D):
338 gross 3503 TEST_TOL=0.01
339 gross 3569 SOLVER=DarcyFlow.EVAL
340 gross 3503
341 sshaw 4984 class Test_Darcy3D_POST(Darcy3D):
342 gross 3503 TEST_TOL=1.e-3
343     SOLVER=DarcyFlow.POST
344    
345 sshaw 4984 class Test_Darcy3D_SMOOTH(Darcy3D):
346 gross 3569 TEST_TOL=0.01
347     SOLVER=DarcyFlow.SMOOTH
348 gross 3503
349    
350     if __name__ == '__main__':
351 sshaw 4984 run_tests(__name__, classes=[
352     Test_Darcy3D_EVAL, Test_Darcy3D_POST, Test_Darcy3D_SMOOTH,
353     Test_Darcy2D_EVAL, Test_Darcy2D_POST, Test_Darcy2D_SMOOTH
354     ],exit_on_failure=True)
355 gross 3503
356    

  ViewVC Help
Powered by ViewVC 1.1.26