/[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 5775 - (show annotations)
Thu Jul 30 08:01:06 2015 UTC (3 years, 10 months ago) by sshaw
File MIME type: text/x-python
File size: 12508 byte(s)
pushing release to trunk
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, ReadGmsh
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 @unittest.skipIf(not HAVE_GMSH, "gmsh not available")
64 def test_getpotential3dPolePole(self):
65 structured=False
66 if structured:
67 extents=[1000,1000,1000]
68 dom=Brick(50,50,50,l0=extents[0],l1=extents[1],l2=-extents[2])
69 else:
70 extents=[1000,2000,2000]
71 tags=[]
72 points=[]
73 tags.append("e1")
74 tags.append("e2")
75 tags.append("e3")
76 tags.append("e4")
77 points.append([-0.4*extents[0], 0, 0])
78 points.append([-0.2*extents[0], 0, 0])
79 points.append([0.2*extents[0], 0, 0])
80 points.append([0.4*extents[0], 0, 0])
81 verbosity = 3
82 filename = os.path.join(TEST_DATA_ROOT, "pole.geo")
83 meshname = os.path.join(TEST_DATA_ROOT, "dcResPolePole50-5.msh")
84 gmshGeo2Msh(filename, meshname, 3, 1, verbosity)
85 dom = ReadGmsh(meshname, 3, diracTags=tags, diracPoints=points)
86 totalApparentRes = 130.
87 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
88 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
89 current=1000.
90 a=4*0.05*extents[0]
91 midPoint = [0.5*extents[0],0.5*extents[1]]
92 directionVector=[1.,0.]
93 numElectrodes = 4
94
95 pps=PolePoleSurvey(dom, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes)
96 delPhi=pps.getPotential()
97 totalApparentResList = pps.getApparentResistivityTotal()
98 for i in totalApparentResList:
99 res_a = abs(i-totalApparentRes)
100 res_b = 0.05 * totalApparentRes
101 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
102
103 def test_getPotential3dSchlumberger(self):
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 if not HAVE_GMSH:
112 raise unittest.SkipTest("gmsh required for test")
113 lc=10.0
114 bufferThickness=100
115 extents=[200,200,200]
116 midPoint = [0,0]
117 lcDiv=10.0
118 electrodes=[]
119 start=[]
120 start.append(midPoint[0] - (((numElectrodes-1)*interval_a)/2. * directionVector[0]))
121 start.append(midPoint[1] - (((numElectrodes-1)*interval_a)/2. * directionVector[1]))
122 electrodeTags=[]
123 electrodeLst=[]
124 for i in range(numElectrodes):
125 electrodes.append([start[0]+(directionVector[0]*i*interval_a), start[1]+(directionVector[1]*i*interval_a),0])
126 electrodeTags.append("e%d"%i)
127 runName=os.path.join(WORKDIR, "dcResSchlum%d-%d"%(lc,lc/lcDiv))
128 filename = os.path.join(TEST_DATA_ROOT, "schlum.geo")
129 meshname = os.path.join(TEST_DATA_ROOT, "dcResSchlum10-1.msh")
130 verbosity=3
131 gmshGeo2Msh(filename, meshname, 3, 1, verbosity)
132 dom = ReadGmsh(meshname, 3, diracTags=electrodeTags, diracPoints=electrodes)
133 if mpirank==0:
134 os.unlink(meshname)
135
136 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
137 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
138 schs=SchlumbergerSurvey(dom, primaryConductivity, secondaryConductivity,
139 current, interval_a, interval_n, midPoint, directionVector,
140 numElectrodes)
141 potentials=schs.getPotentialAnalytic()
142 totalApparentRes=schs.getApparentResistivity(potentials[2])
143 for i in totalApparentRes:
144 for j in i:
145 res_a = abs(j-totalApparentResVal)
146 res_b = 0.1 * totalApparentResVal
147 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
148
149 def test_getPotentialDipDip(self):
150 structured=False
151 totalApparentRes = 130.
152 if structured:
153 extents=[100,100,100]
154 dom=Brick(25,25,25,l0=extents[0],l1=extents[1],l2=-extents[2])
155 else:
156 if not HAVE_GMSH:
157 raise unittest.SkipTest("gmsh required for test")
158 lc=10.0
159 bufferThickness=300
160 extents=[100,100,100]
161 electrodes=[]
162 tags=[]
163 lcDiv=10
164 tags.append("e0" )
165 tags.append("e1" )
166 tags.append("e2" )
167 tags.append("e3" )
168 tags.append("e4" )
169 tags.append("e5" )
170 tags.append("e6" )
171 tags.append("e7" )
172 tags.append("e8" )
173 tags.append("e9" )
174 tags.append("e10")
175 tags.append("e11")
176 electrodes.append([-22.0, 0.0, 0])
177 electrodes.append([-18.0, 0.0, 0])
178 electrodes.append([-14.0, 0.0, 0])
179 electrodes.append([-10.0, 0.0, 0])
180 electrodes.append([-6.0, 0.0, 0])
181 electrodes.append([-2.0, 0.0, 0])
182 electrodes.append([ 2.0, 0.0, 0])
183 electrodes.append([ 6.0, 0.0, 0])
184 electrodes.append([ 10.0, 0.0, 0])
185 electrodes.append([ 14.0, 0.0, 0])
186 electrodes.append([ 18.0, 0.0, 0])
187 electrodes.append([ 22.0, 0.0, 0])
188 filename = os.path.join(TEST_DATA_ROOT, "dip.geo")
189 meshname = os.path.join(TEST_DATA_ROOT, "dcResdipdip-1.msh")
190 verbosity=3
191 gmshGeo2Msh(filename, meshname, 3, 1, verbosity)
192 dom = ReadGmsh(meshname, 3, diracTags=tags, diracPoints=electrodes)
193 if mpirank==0:
194 os.unlink(meshname)
195 n=5
196 totalApparentResVal = 130.
197 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
198 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
199 current=1000.
200 numElectrodes = 12
201 # a=(.8*extents[0])/numElectrodes
202 a=4
203 midPoint = [0,0]
204 directionVector=[1.,0.]
205 dipdips=DipoleDipoleSurvey(dom, primaryConductivity, secondaryConductivity, current, a,n, midPoint, directionVector, numElectrodes)
206 dipdips.getPotential()
207 primaryApparentRes=dipdips.getApparentResistivityPrimary()
208 SecondaryApparentRes=dipdips.getApparentResistivitySecondary()
209 totalApparentRes=dipdips.getApparentResistivityTotal()
210 for i in totalApparentRes:
211 for j in i:
212 res_a = abs(j-totalApparentResVal)
213 res_b = 0.1 * totalApparentResVal
214 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
215
216 def test_getPotentialWenner(self):
217 totalApparentResVal = 130.
218 extents=[100,100,100]
219 dom=Brick(50,50,50,l0=extents[0],l1=extents[1],l2=-extents[2])
220 primaryConductivity=Scalar(1/100., ContinuousFunction(dom))
221 secondaryConductivity=Scalar(1/130., ContinuousFunction(dom))
222 current=1000.
223 numElectrodes = 8
224 a=2
225 midPoint = [0.5*extents[0]+1,0.5*extents[1]]
226 directionVector=[1.,0.]
227 wenSurv=WennerSurvey(dom, primaryConductivity, secondaryConductivity,
228 current, a, midPoint, directionVector, numElectrodes)
229 wenSurv.getPotential()
230 primaryApparentRes=wenSurv.getApparentResistivityPrimary()
231 SecondaryApparentRes=wenSurv.getApparentResistivitySecondary()
232 totalApparentRes=wenSurv.getApparentResistivityTotal()
233 for i in totalApparentRes:
234 res_a = abs(i-totalApparentResVal)
235 res_b = 0.1 * totalApparentResVal
236 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
237
238 def test_getPotentialPolDip(self):
239 structured=False
240 totalApparentRes = 130.
241 if structured:
242 extents=[100,100,100]
243 dom=Brick(25,25,25,l0=extents[0],l1=extents[1],l2=-extents[2])
244 else:
245 if not HAVE_GMSH:
246 raise unittest.SkipTest("gmsh required for test")
247
248 extents=[100,100,100]
249 electrodes=[]
250 tags=[]
251 tags.append("e0" )
252 tags.append("e1" )
253 tags.append("e2" )
254 tags.append("e3" )
255 tags.append("e4" )
256 tags.append("e5" )
257 tags.append("e6" )
258 tags.append("e7" )
259 tags.append("e8" )
260 tags.append("e9" )
261 tags.append("e10")
262 tags.append("e11")
263 electrodes.append([-22.0, 0.0, 0])
264 electrodes.append([-18.0, 0.0, 0])
265 electrodes.append([-14.0, 0.0, 0])
266 electrodes.append([-10.0, 0.0, 0])
267 electrodes.append([-6.0, 0.0, 0])
268 electrodes.append([-2.0, 0.0, 0])
269 electrodes.append([ 2.0, 0.0, 0])
270 electrodes.append([ 6.0, 0.0, 0])
271 electrodes.append([ 10.0, 0.0, 0])
272 electrodes.append([ 14.0, 0.0, 0])
273 electrodes.append([ 18.0, 0.0, 0])
274 electrodes.append([ 22.0, 0.0, 0])
275 filename = os.path.join(TEST_DATA_ROOT, "dip.geo")
276 meshname = os.path.join(TEST_DATA_ROOT, "dcRespoldip10-1.msh")
277 verbosity=3
278 gmshGeo2Msh(filename, meshname, 3, 1, verbosity)
279 dom = ReadGmsh(meshname, 3, diracTags=tags, diracPoints=electrodes)
280
281 if mpirank==0:
282 os.unlink(meshname)
283 n=5
284 totalApparentResVal = 130.
285 primaryConductivity = Scalar(1/100., ContinuousFunction(dom))
286 secondaryConductivity = Scalar(1/130., ContinuousFunction(dom))
287 current=1000.
288 numElectrodes = 12
289 # a=(.8*extents[0])/numElectrodes
290 a=4
291 midPoint = [0,0]
292 directionVector=[1.,0.]
293 poldips=PoleDipoleSurvey(dom, primaryConductivity,
294 secondaryConductivity, current, a,n, midPoint,
295 directionVector, numElectrodes)
296 poldips.getPotential()
297 primaryApparentRes=poldips.getApparentResistivityPrimary()
298 SecondaryApparentRes=poldips.getApparentResistivitySecondary()
299 totalApparentRes=poldips.getApparentResistivityTotal()
300 for i in totalApparentRes:
301 for j in i:
302 res_a = abs(j-totalApparentResVal)
303 res_b = 0.1 * totalApparentResVal
304 self.assertLess(res_a, res_b, "result of %g greater than tolerance of %g"%(res_a, res_b))
305
306 ################################
307 if __name__ == '__main__':
308 run_tests(__name__, exit_on_failure=True)
309
310

  ViewVC Help
Powered by ViewVC 1.1.26