/[escript]/trunk/downunder/test/python/run_dcforward.py
ViewVC logotype

Contents of /trunk/downunder/test/python/run_dcforward.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5706 - (show annotations)
Mon Jun 29 03:41:36 2015 UTC (4 years, 2 months ago) by sshaw
File MIME type: text/x-python
File size: 14327 byte(s)
all python files now force use of python3 prints and division syntax to stop sneaky errors appearing in py3 environs
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16 from __future__ import print_function, division
17
18 from esys.downunder import *
19 from esys.escript import *
20 from esys.escriptcore.testing import *
21 import esys.escriptcore.utestselect as unittest
22 import numpy as np
23
24 try:
25 TEST_DATA_ROOT=os.environ['DOWNUNDER_TEST_DATA_ROOT']
26 except KeyError:
27 TEST_DATA_ROOT='ref_data'
28
29 try:
30 WORKDIR=os.environ['DOWNUNDER_WORKDIR']
31 except KeyError:
32 WORKDIR='.'
33
34 try:
35 from esys.finley import Rectangle, Brick
36 HAVE_FINLEY = True
37 except ImportError:
38 HAVE_FINLEY = False
39
40 HAVE_GMSH = getEscriptParamInt("GMSH_SUPPORT")
41
42 mpisize = getMPISizeWorld()
43 mpirank = getMPIRankWorld()
44
45 @unittest.skipIf(not HAVE_FINLEY, "Finley module not available")
46 class TestDCResistivityForward(unittest.TestCase):
47 def test_getPotential2d(self):
48 dom=Rectangle(20,20,l0=1000,l1=-1000,d1=mpisize)
49 extents=[1000,1000]
50 primaryConductivity=Data(1/100., ContinuousFunction(dom))
51 secondaryConductivity=Data(1/130., ContinuousFunction(dom))
52 current=1000.
53 a=0.05*extents[0]
54 start = [0.25*extents[0]]
55 directionVector=[1]
56 numElectrodes = 10
57
58 self.assertRaises(NotImplementedError,
59 lambda:PolePoleSurvey(dom, primaryConductivity,
60 secondaryConductivity, current, a, start,
61 directionVector, numElectrodes))
62
63 def test_getpotential3dPolePole(self):
64 structured=False
65 if structured:
66 extents=[1000,1000,1000]
67 dom=Brick(50,50,50,l0=extents[0],l1=extents[1],l2=-extents[2])
68 else:
69 if not HAVE_GMSH:
70 raise unittest.SkipTest("gmsh required for test")
71 lc=30.0
72 bufferThickness=3000
73 extents=[1000,2000,2000]
74 electrodeLst=[]
75 lcDiv=10
76 electrodeLst.append(("e1",[-0.4*extents[0], 0, 0,lc/lcDiv]))
77 electrodeLst.append(("e2",[-0.2*extents[0], 0, 0,lc/lcDiv]))
78 electrodeLst.append(("e3",[0.2*extents[0], 0, 0,lc/lcDiv]))
79 electrodeLst.append(("e4",[0.4*extents[0], 0, 0,lc/lcDiv]))
80 runName=os.path.join(WORKDIR, "dcResPolePole%d-%d"%(lc,lc/lcDiv))
81 domGen=DCResDomGenerator(extents, electrodeLst,lc=lc,tmpDir=WORKDIR,bufferThickness=bufferThickness,prism=None)
82 dom = domGen.getDom(mshName=runName+".msh")
83 if mpirank==0:
84 os.unlink(runName+".msh")
85 totalApparentRes = 130.
86 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
87 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
88 current=1000.
89 a=4*0.05*extents[0]
90 midPoint = [0,0]
91 directionVector=[1.,0.]
92 numElectrodes = 4
93
94 pps=PolePoleSurvey(dom, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)
95 delPhi=pps.getPotential()
96 totalApparentResList = pps.getApparentResistivityTotal()
97 for i in totalApparentResList:
98 res_a = abs(i-totalApparentRes)
99 res_b = 0.1 * totalApparentRes
100 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
101
102 def test_getPotential3dSchlumberger(self):
103 structured=False
104 numElectrodes = 12
105 directionVector=[1.,0.]
106 midPoint=[]
107 totalApparentResVal = 130.
108 current=0.5
109 interval_a = 5
110 interval_n = 5
111
112
113 if structured:
114 #does not work because finley does not allow the specification of domain origin
115 extents=[200,200,200]
116 dom=Brick(25,25,25,l0=(-extents[0]/2,extents[0]/2),l1=(-extents[1]/2,extents[1]/2),l2=-extents[2])
117 midPoint = [0,0]
118 else:
119 if not HAVE_GMSH:
120 raise unittest.SkipTest("gmsh required for test")
121 lc=10.0
122 bufferThickness=100
123 extents=[200,200,200]
124 midPoint = [0,0]
125 lcDiv=10.0
126 electrodes=[]
127 start=[]
128 start.append(midPoint[0] - (((numElectrodes-1)*interval_a)/2. * directionVector[0]))
129 start.append(midPoint[1] - (((numElectrodes-1)*interval_a)/2. * directionVector[1]))
130 electrodeTags=[]
131 electrodeLst=[]
132 for i in range(numElectrodes):
133 electrodes.append([start[0]+(directionVector[0]*i*interval_a), start[1]+(directionVector[1]*i*interval_a),0])
134 electrodeTags.append("e%d"%i)
135 electrodeLst.append((electrodeTags[i],[electrodes[i][0], electrodes[i][1], electrodes[i][2], lc/lcDiv]))
136 runName=os.path.join(WORKDIR, "dcResSchlum%d-%d"%(lc,lc/lcDiv))
137 domGen=DCResDomGenerator(extents, electrodeLst,lc=lc,tmpDir=WORKDIR,bufferThickness=bufferThickness,prism=None)
138 dom = domGen.getDom(mshName=runName+".msh",fieldSize=[70,100])
139 fn = domGen.getFileName()
140 if mpirank==0:
141 os.unlink(runName+".msh")
142
143 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
144 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
145 schs=SchlumbergerSurvey(dom, primaryConductivity, secondaryConductivity,
146 current, interval_a, interval_n, midPoint, directionVector,
147 numElectrodes)
148 potentials=schs.getPotentialAnalytic()
149 totalApparentRes=schs.getApparentResistivity(potentials[2])
150 for i in totalApparentRes:
151 for j in i:
152 res_a = abs(j-totalApparentResVal)
153 res_b = 0.1 * totalApparentResVal
154 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
155
156 def test_getPotentialDipDip(self):
157 structured=False
158 totalApparentRes = 130.
159 if structured:
160 extents=[100,100,100]
161 dom=Brick(25,25,25,l0=extents[0],l1=extents[1],l2=-extents[2])
162 else:
163 if not HAVE_GMSH:
164 raise unittest.SkipTest("gmsh required for test")
165 lc=10.0
166 bufferThickness=300
167 extents=[100,100,100]
168 electrodeLst=[]
169 lcDiv=10
170 electrodeLst.append(("e0" , [-22.0, 0.0, 0, lc/lcDiv]))
171 electrodeLst.append(("e1" , [-18.0, 0.0, 0, lc/lcDiv]))
172 electrodeLst.append(("e2" , [-14.0, 0.0, 0, lc/lcDiv]))
173 electrodeLst.append(("e3" , [-10.0, 0.0, 0, lc/lcDiv]))
174 electrodeLst.append(("e4" , [-6.0, 0.0, 0, lc/lcDiv]))
175 electrodeLst.append(("e5" , [-2.0, 0.0, 0, lc/lcDiv]))
176 electrodeLst.append(("e6" , [ 2.0, 0.0, 0, lc/lcDiv]))
177 electrodeLst.append(("e7" , [ 6.0, 0.0, 0, lc/lcDiv]))
178 electrodeLst.append(("e8" , [ 10.0, 0.0, 0, lc/lcDiv]))
179 electrodeLst.append(("e9" , [ 14.0, 0.0, 0, lc/lcDiv]))
180 electrodeLst.append(("e10", [ 18.0, 0.0, 0, lc/lcDiv]))
181 electrodeLst.append(("e11", [ 22.0, 0.0, 0, lc/lcDiv]))
182 runName=os.path.join(WORKDIR, "dcResdipdip%d-%d"%(lc,lc/lcDiv))
183 domGen=DCResDomGenerator(extents, electrodeLst,lc=lc,tmpDir=WORKDIR,bufferThickness=bufferThickness,prism=None)
184 dom = domGen.getDom(mshName=runName+".msh",reUse=False)
185 if mpirank==0:
186 os.unlink(runName+".msh")
187 n=5
188 totalApparentResVal = 130.
189 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
190 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
191 current=1000.
192 numElectrodes = 12
193 # a=(.8*extents[0])/numElectrodes
194 a=4
195 midPoint = [0,0]
196 directionVector=[1.,0.]
197 dipdips=DipoleDipoleSurvey(dom, primaryConductivity, secondaryConductivity, current, a,n, midPoint, directionVector, numElectrodes)
198 dipdips.getPotential()
199 primaryApparentRes=dipdips.getApparentResistivityPrimary()
200 SecondaryApparentRes=dipdips.getApparentResistivitySecondary()
201 totalApparentRes=dipdips.getApparentResistivityTotal()
202 for i in totalApparentRes:
203 for j in i:
204 res_a = abs(j-totalApparentResVal)
205 res_b = 0.1 * totalApparentResVal
206 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
207
208 def test_getPotentialWenner(self):
209 structured=True
210 totalApparentResVal = 130.
211 if structured:
212 extents=[100,100,100]
213 dom=Brick(50,50,50,l0=extents[0],l1=extents[1],l2=-extents[2])
214 else:
215 if not HAVE_GMSH:
216 raise unittest.SkipTest("gmsh required for test")
217 lc=50.0
218 bufferThickness=3000
219 extents=[1000,2000,2000]
220 electrodeLst=[]
221 lcDiv=10
222
223 electrodeLst.append(("e0" , [28.0, 48.0, 0, lc/lcDiv]))
224 electrodeLst.append(("e1" , [32.0, 48.0, 0, lc/lcDiv]))
225 electrodeLst.append(("e2" , [36.0, 48.0, 0, lc/lcDiv]))
226 electrodeLst.append(("e3" , [40.0, 48.0, 0, lc/lcDiv]))
227 electrodeLst.append(("e4" , [44.0, 48.0, 0, lc/lcDiv]))
228 electrodeLst.append(("e5" , [48.0, 48.0, 0, lc/lcDiv]))
229 electrodeLst.append(("e6" , [52.0, 48.0, 0, lc/lcDiv]))
230 electrodeLst.append(("e7" , [56.0, 48.0, 0, lc/lcDiv]))
231 electrodeLst.append(("e8" , [60.0, 48.0, 0, lc/lcDiv]))
232 electrodeLst.append(("e9" , [64.0, 48.0, 0, lc/lcDiv]))
233 electrodeLst.append(("e10" , [68.0, 48.0, 0, lc/lcDiv]))
234 electrodeLst.append(("e11" , [72.0, 48.0, 0, lc/lcDiv]))
235
236 domGen=DCResDomGenerator(extents, electrodeLst,lc=lc,tmpDir=WORKDIR,bufferThickness=bufferThickness,prism=None)
237 runName=os.path.join(WORKDIR, "wenner%d-%d"%(lc,lc/lcDiv))
238 dom = domGen.getDom(mshName=runName+".msh")
239 if mpirank==0:
240 os.unlink(runName+".msh")
241 totalApparentRes = 130.
242 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
243 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
244 current=1000.
245 numElectrodes = 8
246 # a=(.8*extents[0])/numElectrodes
247 a=2
248 midPoint = [0.5*extents[0]+1,0.5*extents[1]]
249 directionVector=[1.,0.]
250 wenSurv=WennerSurvey(dom, primaryConductivity, secondaryConductivity,
251 current, a, midPoint, directionVector, numElectrodes)
252 wenSurv.getPotential()
253 primaryApparentRes=wenSurv.getApparentResistivityPrimary()
254 SecondaryApparentRes=wenSurv.getApparentResistivitySecondary()
255 totalApparentRes=wenSurv.getApparentResistivityTotal()
256 for i in totalApparentRes:
257 res_a = abs(i-totalApparentResVal)
258 res_b = 0.1 * totalApparentResVal
259 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
260
261 def test_getPotentialPolDip(self):
262 structured=False
263 totalApparentRes = 130.
264 if structured:
265 extents=[100,100,100]
266 dom=Brick(25,25,25,l0=extents[0],l1=extents[1],l2=-extents[2])
267 else:
268 if not HAVE_GMSH:
269 raise unittest.SkipTest("gmsh required for test")
270 lc=10.0
271 bufferThickness=300
272 extents=[100,100,100]
273 electrodeLst=[]
274 lcDiv=10
275 electrodeLst.append(("e0" , [28.0, 48.0, 0, lc/lcDiv]))
276 electrodeLst.append(("e1" , [32.0, 48.0, 0, lc/lcDiv]))
277 electrodeLst.append(("e2" , [36.0, 48.0, 0, lc/lcDiv]))
278 electrodeLst.append(("e3" , [40.0, 48.0, 0, lc/lcDiv]))
279 electrodeLst.append(("e4" , [44.0, 48.0, 0, lc/lcDiv]))
280 electrodeLst.append(("e5" , [48.0, 48.0, 0, lc/lcDiv]))
281 electrodeLst.append(("e6" , [52.0, 48.0, 0, lc/lcDiv]))
282 electrodeLst.append(("e7" , [56.0, 48.0, 0, lc/lcDiv]))
283 electrodeLst.append(("e8" , [60.0, 48.0, 0, lc/lcDiv]))
284 electrodeLst.append(("e9" , [64.0, 48.0, 0, lc/lcDiv]))
285 electrodeLst.append(("e10" , [68.0, 48.0, 0, lc/lcDiv]))
286 electrodeLst.append(("e11" , [72.0, 48.0, 0, lc/lcDiv]))
287 runName=os.path.join(WORKDIR, "dcRespoldip%d-%d"%(lc,lc/lcDiv))
288 domGen=DCResDomGenerator(extents, electrodeLst,lc=lc,tmpDir=WORKDIR,bufferThickness=bufferThickness,prism=None)
289 dom = domGen.getDom(mshName=runName+".msh")
290 if mpirank==0:
291 os.unlink(runName+".msh")
292 n=5
293 totalApparentResVal = 130.
294 primaryConductivity = Scalar(1/100., ContinuousFunction(dom))
295 secondaryConductivity = Scalar(1/130., ContinuousFunction(dom))
296 current=1000.
297 numElectrodes = 12
298 # a=(.8*extents[0])/numElectrodes
299 a=4
300 midPoint = [0.5*extents[0],0.5*extents[1] - 2]
301 directionVector=[1.,0.]
302 poldips=PoleDipoleSurvey(dom, primaryConductivity,
303 secondaryConductivity, current, a,n, midPoint,
304 directionVector, numElectrodes)
305 poldips.getPotential()
306 primaryApparentRes=poldips.getApparentResistivityPrimary()
307 SecondaryApparentRes=poldips.getApparentResistivitySecondary()
308 totalApparentRes=poldips.getApparentResistivityTotal()
309 for i in totalApparentRes:
310 for j in i:
311 res_a = abs(j-totalApparentResVal)
312 res_b = 0.1 * totalApparentResVal
313 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
314
315 ################################
316 if __name__ == '__main__':
317 run_tests(__name__, exit_on_failure=True)
318
319

  ViewVC Help
Powered by ViewVC 1.1.26