/[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 6651 - (show annotations)
Wed Feb 7 02:12:08 2018 UTC (20 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 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2018 by The University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Apache License, version 2.0
10 # http://www.apache.org/licenses/LICENSE-2.0
11 #
12 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 from __future__ import print_function, division
19
20 __copyright__="""Copyright (c) 2003-2018 by The University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Apache License, version 2.0
24 http://www.apache.org/licenses/LICENSE-2.0"""
25 __url__="https://launchpad.net/escript-finley"
26
27 import esys.escriptcore.utestselect as unittest
28 from esys.escriptcore.testing import *
29
30 from esys.escript import *
31 from esys.escript.models import DarcyFlow
32 from esys.finley import Rectangle, Brick
33
34 import os
35
36 VERBOSE=False
37
38 try:
39 FINLEY_WORKDIR=os.environ['FINLEY_WORKDIR']
40 except KeyError:
41 FINLEY_WORKDIR='.'
42
43 class Darcy(unittest.TestCase): #subclassing required
44 # 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 for i in range(self.dom.getDim()):
64 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 for i in range(self.dom.getDim()):
82 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 v,p=df.solve(u_ref,p)
110
111
112 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
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 Darcy2D(Darcy): #subclassing required
305 TOL=1e-4
306 WIDTH=1.
307 SOLVER=DarcyFlow.POST
308 def setUp(self):
309 NE=40 # warning: 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(Darcy2D):
316 TEST_TOL=0.01
317 SOLVER=DarcyFlow.EVAL
318
319 class Test_Darcy2D_POST(Darcy2D):
320 TEST_TOL=1.e-3
321 SOLVER=DarcyFlow.POST
322
323 class Test_Darcy2D_SMOOTH(Darcy2D):
324 TEST_TOL=0.01
325 SOLVER=DarcyFlow.SMOOTH
326
327 class Darcy3D(Darcy): #subclassing required
328 WIDTH=1.
329 SOLVER=DarcyFlow.POST
330 def setUp(self):
331 NE=29 # warning: 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(Darcy3D):
338 TEST_TOL=0.01
339 SOLVER=DarcyFlow.EVAL
340
341 class Test_Darcy3D_POST(Darcy3D):
342 TEST_TOL=1.e-3
343 SOLVER=DarcyFlow.POST
344
345 class Test_Darcy3D_SMOOTH(Darcy3D):
346 TEST_TOL=0.01
347 SOLVER=DarcyFlow.SMOOTH
348
349
350 if __name__ == '__main__':
351 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
356

  ViewVC Help
Powered by ViewVC 1.1.26