/[escript]/trunk/downunder/py_src/dcresistivityforwardmodeling.py
ViewVC logotype

Contents of /trunk/downunder/py_src/dcresistivityforwardmodeling.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: 41202 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.escript import Data, kronecker, whereZero,inf,sup,ContinuousFunction,grad,Function,Lsup, Scalar, DiracDeltaFunctions
19 from esys.escript.linearPDEs import LinearPDE
20 from esys.escript.pdetools import Locator
21 from math import pi
22 from esys.weipa import saveSilo
23
24 try:
25 xrange
26 except NameError:
27 xrange = range
28
29
30
31 class DcResistivityForward(object):
32 """
33 This class allows for the solution of dc resistivity forward problems via
34 the calculation of a primary and secondary potential. Conductivity values
35 are to be provided for the primary problem which is a homogeneous half space
36 of a chosen conductivity and for the secondary problem which typically
37 varies it conductivity spatially across the domain. The primary potential
38 acts as a reference point typically based of some know reference conductivity
39 however any value will suffice. The primary potential is implemented to
40 avoid the use of dirac delta functions.
41 """
42 def __init__(self):
43 """
44 This is a skeleton class for all the other forward modeling classes.
45 """
46 pass
47
48 def getPotential(self):
49 raise NotImplementedError
50
51
52 def getApparentResistivity(self):
53 raise NotImplementedError
54
55 def checkBounds(self):
56 X = ContinuousFunction(self.domain).getX()
57 xDim=[inf(X[0]),sup(X[0])]
58 yDim=[inf(X[1]),sup(X[1])]
59 zDim=[inf(X[2]),sup(X[2])]
60 for i in range(self.numElectrodes):
61 if (self.electrodes[i][0] < xDim[0] or self.electrodes[i][0] > xDim[1]
62 or self.electrodes[i][1] < yDim[0] or self.electrodes[i][1] > yDim[1]):
63 print (self.electrodes[i])
64 raise ValueError("Electrode setup extents past domain dimentions")
65 def getElectrodes(self):
66 """
67 retuns the list of electrodes with locations
68 """
69 return self.electrodes
70
71 class SchlumbergerSurvey(DcResistivityForward):
72 """
73 Schlumberger survey forward calculation
74 """
75 def __init__(self, domain, primaryConductivity, secondaryConductivity,
76 current, a, n, midPoint, directionVector, numElectrodes):
77 super(SchlumbergerSurvey, self).__init__()
78 """
79 :param domain: Domain of the model
80 :type domain: `Domain`
81 :param primaryConductivity: preset primary conductivity data object
82 :type primaryConductivity: data
83 :param secondaryConductivity:preset secondary conductivity data object
84 :type secondaryConductivity: data
85 :param current: amount of current to be injected at the current electrode
86 :type current: float or int
87 :param a: the spacing between current and potential electrodes
88 :type a: float or int
89 :param n: multiple of spacing between electrodes. typicaly interger
90 :type n: float or int
91 :param midPoint: midPoint of the survey, as a list containing x,y coords
92 :type a: list
93 :param directionVector: two element list specifying the direction the
94 survey should extend from the midpoint
95 :type a: list
96 :param numElectrodes: the number of electrodes to be used in the survey
97 must be a multiple of 2 for polepole survey:
98 :type numElectrodes: int
99 """
100 self.domain=domain
101 self.primaryConductivity=primaryConductivity
102 self.secondaryConductivity=secondaryConductivity
103 self.a=a
104 self.n=n
105 self.current=current
106 self.numElectrodes=numElectrodes
107 self.delPhiPrimaryList=[]
108 self.delPhiSecondaryList=[]
109 self.samples=[]
110 self.sources=[]
111 self.electrodeDict={}
112 self.electrodeTags=[]
113 if (numElectrodes < 4):
114 raise ValueError("numElectrodes must atleast 4 for schlumberger surveys")
115 if n > ((numElectrodes-2)//2):
116 raise ValueError("specified n does not fit max n = %d"%((numElectrodes-2)//2))
117 if len(directionVector) == 2:
118 electrodes=[]
119 start=[]
120 start.append(midPoint[0] - (((numElectrodes-1)*a)/2. * directionVector[0]))
121 start.append(midPoint[1] - (((numElectrodes-1)*a)/2. * directionVector[1]))
122 for i in range(numElectrodes):
123 electrodes.append([start[0]+(directionVector[0]*i*a), start[1]+(directionVector[1]*i*a),0])
124 self.electrodeTags.append("e%d"%i)
125 self.electrodeDict[self.electrodeTags[i]]=electrodes[i]
126 else:
127 raise NotImplementedError("2d forward model is not yet implemented please provide a 2 component directionVector for a 3d survey")
128 self.electrodes=electrodes
129 self.checkBounds()
130
131 def getPotentialNumeric(self):
132 """
133 Returns 3 list each made up of a number of list containing primary, secondary and total
134 potentials diferences. Each of the lists contain a list for each value of n.
135 """
136 primCon=self.primaryConductivity
137 coords=self.domain.getX()
138 tol=1e-8
139 pde=LinearPDE(self.domain, numEquations=1)
140 pde.getSolverOptions().setTolerance(tol)
141 pde.setSymmetryOn()
142 primaryPde=LinearPDE(self.domain, numEquations=1)
143 primaryPde.getSolverOptions().setTolerance(tol)
144 primaryPde.setSymmetryOn()
145 DIM=self.domain.getDim()
146 x=self.domain.getX()
147 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
148 for i in xrange(DIM-1):
149 xi=x[i]
150 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
151 A = self.secondaryConductivity * kronecker(self.domain)
152 APrimary = self.primaryConductivity * kronecker(self.domain)
153 pde.setValue(A=A,q=q)
154 primaryPde.setValue(A=APrimary,q=q)
155
156 delPhiSecondaryList = []
157 delPhiPrimaryList = []
158 delPhiTotalList = []
159 for i in range(1,self.n+1): # 1 to n
160 maxR = self.numElectrodes - 1 - (2*i) #max amount of readings that will fit in the survey
161 delPhiSecondary = []
162 delPhiPrimary = []
163 delPhiTotal = []
164 for j in range(maxR):
165 y_dirac=Scalar(0,DiracDeltaFunctions(self.domain))
166 y_dirac.setTaggedValue(self.electrodeTags[j],self.current)
167 y_dirac.setTaggedValue(self.electrodeTags[j + ((2*i) + 1)],-self.current)
168 self.sources.append([self.electrodeTags[j], self.electrodeTags[j + ((2*i) + 1)]])
169 primaryPde.setValue(y_dirac=y_dirac)
170 numericPrimaryPot = primaryPde.getSolution()
171 X=(primCon-self.secondaryConductivity) * grad(numericPrimaryPot)
172 pde.setValue(X=X)
173 u=pde.getSolution()
174 loc=Locator(self.domain,[self.electrodes[j+i],self.electrodes[j+i+1]])
175 self.samples.append([self.electrodeTags[j+i],self.electrodeTags[j+i+1]])
176 valPrimary=loc.getValue(numericPrimaryPot)
177 valSecondary=loc.getValue(u)
178 delPhiPrimary.append(valPrimary[1]-valPrimary[0])
179 delPhiSecondary.append(valSecondary[1]-valSecondary[0])
180 delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
181 delPhiPrimaryList.append(delPhiPrimary)
182 delPhiSecondaryList.append(delPhiSecondary)
183 delPhiTotalList.append(delPhiTotal)
184
185 self.delPhiPrimaryList=delPhiPrimaryList
186 self.delPhiSecondaryList=delPhiSecondaryList
187 self.delPhiTotalList = delPhiTotalList
188 return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]
189
190
191
192 def getPotentialAnalytic(self):
193 """
194 Returns 3 list each made up of a number of list containing primary, secondary and total
195 potentials diferences. Each of the lists contain a list for each value of n.
196 """
197 coords=self.domain.getX()
198 pde=LinearPDE(self.domain, numEquations=1)
199 tol=1e-8
200 pde.getSolverOptions().setTolerance(tol)
201 pde.setSymmetryOn()
202 primCon=self.primaryConductivity
203 DIM=self.domain.getDim()
204 x=self.domain.getX()
205 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
206 for i in xrange(DIM-1):
207 xi=x[i]
208 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
209 A = self.secondaryConductivity * kronecker(self.domain)
210 pde.setValue(A=A,q=q)
211
212
213 delPhiSecondaryList = []
214 delPhiPrimaryList = []
215 delPhiTotalList = []
216 for i in range(1,self.n+1): # 1 to n
217 maxR = self.numElectrodes - 1 - (2*i) #max amount of readings that will fit in the survey
218 delPhiSecondary = []
219 delPhiPrimary = []
220 delPhiTotal = []
221 for j in range(maxR):
222 analyticRsOne=Data(0,(3,),ContinuousFunction(self.domain))
223 analyticRsOne[0]=(coords[0]-self.electrodes[j][0])
224 analyticRsOne[1]=(coords[1]-self.electrodes[j][1])
225 analyticRsOne[2]=(coords[2])
226 rsMagOne=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**0.5
227 analyticRsTwo=Data(0,(3,),ContinuousFunction(self.domain))
228 analyticRsTwo[0]=(coords[0]-self.electrodes[j + ((2*i) + 1)][0])
229 analyticRsTwo[1]=(coords[1]-self.electrodes[j + ((2*i) + 1)][1])
230 analyticRsTwo[2]=(coords[2])
231 rsMagTwo=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**0.5
232 self.sources.append([self.electrodeTags[j], self.electrodeTags[j + ((2*i) + 1)]])
233 rsMagOne+=(whereZero(rsMagOne)*0.0000001)
234 rsMagTwo+=(whereZero(rsMagTwo)*0.0000001)
235
236 analyticPrimaryPot=(self.current/(2*pi*primCon*rsMagOne))-(self.current/(2*pi*primCon*rsMagTwo))
237 analyticRsOnePower=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**1.5
238 analyticRsOnePower = analyticRsOnePower+(whereZero(analyticRsOnePower)*0.0001)
239 analyticRsTwoPower=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**1.5
240 analyticRsTwoPower = analyticRsTwoPower+(whereZero(analyticRsTwoPower)*0.0001)
241
242 gradAnalyticPrimaryPot = Data(0,(3,),ContinuousFunction(self.domain))
243 gradAnalyticPrimaryPot[0] =(self.current/(2*pi*primCon)) * ((-analyticRsOne[0]/analyticRsOnePower) + (analyticRsTwo[0]/analyticRsTwoPower))
244 gradAnalyticPrimaryPot[1] =(self.current/(2*pi*primCon)) * ((-analyticRsOne[1]/analyticRsOnePower) + (analyticRsTwo[1]/analyticRsTwoPower))
245 gradAnalyticPrimaryPot[2] =(self.current/(2*pi*primCon)) * ((-analyticRsOne[2]/analyticRsOnePower) + (analyticRsTwo[2]/analyticRsTwoPower))
246 X=(primCon-self.secondaryConductivity) * (gradAnalyticPrimaryPot)
247 pde.setValue(X=X)
248 u=pde.getSolution()
249 loc=Locator(self.domain,[self.electrodes[j+i],self.electrodes[j+i+1]])
250 self.samples.append([self.electrodeTags[j+i],self.electrodeTags[j+i+1]])
251 valPrimary=loc.getValue(analyticPrimaryPot)
252 valSecondary=loc.getValue(u)
253 delPhiPrimary.append(valPrimary[1]-valPrimary[0])
254 delPhiSecondary.append(valSecondary[1]-valSecondary[0])
255 delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
256 delPhiPrimaryList.append(delPhiPrimary)
257 delPhiSecondaryList.append(delPhiSecondary)
258 delPhiTotalList.append(delPhiTotal)
259
260 self.delPhiPrimaryList=delPhiPrimaryList
261 self.delPhiSecondaryList=delPhiSecondaryList
262 self.delPhiTotalList = delPhiTotalList
263 return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]
264
265 def getApparentResistivity(self, delPhiList):
266 resistivityList = []
267 n=self.n
268 if (self.delPhiPrimaryList==[]):
269 self.getPotential()
270
271 nCount=1
272 for i in delPhiList:
273 resistivity=[]
274 for j in i:
275 resistivity.append((-j/self.current)*(nCount*(nCount+1))*pi*self.a)
276 nCount=nCount+1
277 resistivityList.append(resistivity)
278 return resistivityList
279
280 def getSourcesSamples(self):
281 """
282 return a list of tuples of sample locations followed by a list of tuples
283 of source locations.
284 """
285 return self.samples, self.sources
286
287 def getElectrodeDict(self):
288 """
289 retuns the electrode dictionary
290 """
291 return self.electrodeDict
292
293 class WennerSurvey(DcResistivityForward):
294 """
295 WennerSurvey forward calculation
296 """
297 def __init__(self, domain, primaryConductivity, secondaryConductivity,
298 current, a, midPoint, directionVector, numElectrodes):
299 """
300 :param domain: domain of the model
301 :type domain: `Domain`
302 :param primaryConductivity: preset primary conductivity data object
303 :type primaryConductivity: ``data``
304 :param secondaryConductivity:preset secondary conductivity data object
305 :type secondaryConductivity: ``data``
306 :param current: amount of current to be injected at the current electrode
307 :type current: ``float`` or ``int``
308 :param a: the spacing between current and potential electrodes
309 :type a: ``float`` or ``int``
310 :param midPoint: midPoint of the survey, as a list containing x,y coords
311 :type a: ``list``
312 :param directionVector: two element list specifying the direction the
313 survey should extend from the midpoint
314 :type a: ``list``
315 :param numElectrodes: the number of electrodes to be used in the survey
316 must be a multiple of 2 for polepole survey
317 :type numElectrodes: ``int``
318 """
319 super(WennerSurvey, self).__init__()
320 self.domain=domain
321 self.primaryConductivity=primaryConductivity
322 self.secondaryConductivity=secondaryConductivity
323 self.a=a
324 self.current=current
325 self.numElectrodes=numElectrodes
326 self.delPhiSecondary=[]
327 self.delPhiPrimary=[]
328 if (numElectrodes<4 ):
329 raise ValueError("numElectrodes must be atleast 4 for Wenner surveys")
330 if len(directionVector) == 2:
331 electrodes = []
332 start=[]
333 start.append(midPoint[0] - (((numElectrodes-1)*a)/2. * directionVector[0]))
334 start.append(midPoint[1] - (((numElectrodes-1)*a)/2. * directionVector[1]))
335 for i in range(numElectrodes):
336 electrodes.append([start[0]+(directionVector[0]*i*a),
337 start[1]+(directionVector[1]*i*a),0])
338 else:
339 raise NotImplementedError("2d forward model is not yet implemented please provide a 2 component directionVector for a 3d survey")
340 self.electrodes=electrodes
341 self.checkBounds()
342
343 def getPotential(self):
344 """
345 returns a list containing 3 lists one for each the primary, secondary
346 and total potential.
347 """
348
349 primCon=self.primaryConductivity
350 coords=self.domain.getX()
351 pde=LinearPDE(self.domain, numEquations=1)
352 tol=1e-8
353 pde.getSolverOptions().setTolerance(tol)
354 pde.setSymmetryOn()
355
356 DIM=self.domain.getDim()
357 x=self.domain.getX()
358 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
359 for i in xrange(DIM-1):
360 xi=x[i]
361 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
362 A = self.secondaryConductivity * kronecker(self.domain)
363 pde.setValue(A=A,q=q)
364
365 delPhiSecondary = []
366 delPhiPrimary = []
367 delPhiTotal = []
368 if(len(self.electrodes[0])==3):
369
370 for i in range(self.numElectrodes-3):
371 analyticRsOne=Data(0,(3,),ContinuousFunction(self.domain))
372 analyticRsOne[0]=(coords[0]-self.electrodes[i][0])
373 analyticRsOne[1]=(coords[1]-self.electrodes[i][1])
374 analyticRsOne[2]=(coords[2])
375 rsMagOne=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**0.5
376 analyticRsTwo=Data(0,(3,),ContinuousFunction(self.domain))
377 analyticRsTwo[0]=(coords[0]-self.electrodes[i+3][0])
378 analyticRsTwo[1]=(coords[1]-self.electrodes[i+3][1])
379 analyticRsTwo[2]=(coords[2])
380 rsMagTwo=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**0.5
381 rsMagOne+=(whereZero(rsMagOne)*0.0000001)
382 rsMagTwo+=(whereZero(rsMagTwo)*0.0000001)
383 analyticPrimaryPot=(self.current/(2*pi*primCon*rsMagOne))-(self.current/(2*pi*primCon*rsMagTwo))
384
385 analyticRsOnePower=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**1.5
386 analyticRsOnePower = analyticRsOnePower+(whereZero(analyticRsOnePower)*0.0001)
387 analyticRsTwoPower=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**1.5
388 analyticRsTwoPower = analyticRsTwoPower+(whereZero(analyticRsTwoPower)*0.0001)
389
390 gradAnalyticPrimaryPot = Data(0,(3,),ContinuousFunction(self.domain))
391 gradAnalyticPrimaryPot[0] =(self.current/(2*pi*primCon)) \
392 * ((-analyticRsOne[0]/analyticRsOnePower) \
393 + (analyticRsTwo[0]/analyticRsTwoPower))
394 gradAnalyticPrimaryPot[1] =(self.current/(2*pi*primCon)) \
395 * ((-analyticRsOne[1]/analyticRsOnePower) \
396 + (analyticRsTwo[1]/analyticRsTwoPower))
397 gradAnalyticPrimaryPot[2] =(self.current/(2*pi*primCon)) \
398 * ((-analyticRsOne[2]/analyticRsOnePower)
399 + (analyticRsTwo[2]/analyticRsTwoPower))
400 X=(primCon-self.secondaryConductivity) * (gradAnalyticPrimaryPot)
401 pde.setValue(X=X)
402 u=pde.getSolution()
403 loc=Locator(self.domain,[self.electrodes[i+1],self.electrodes[i+2]])
404 valPrimary=loc.getValue(analyticPrimaryPot)
405 valSecondary=loc.getValue(u)
406 delPhiPrimary.append(valPrimary[1]-valPrimary[0])
407 delPhiSecondary.append(valSecondary[1]-valSecondary[0])
408 delPhiTotal.append(delPhiPrimary[i]+delPhiSecondary[i])
409 else:
410 raise NotImplementedError("2d forward model is not yet implemented")
411
412 self.delPhiSecondary = delPhiSecondary
413 self.delPhiPrimary = delPhiPrimary
414 self.delPhiTotal=delPhiTotal
415 return [delPhiPrimary, delPhiSecondary, delPhiTotal]
416
417 def getApparentResistivityPrimary(self):
418 resistivity = []
419
420 if (self.delPhiPrimary==[]):
421 self.getPotential()
422
423 for i in self.delPhiPrimary:
424 resistivity.append((-i/self.current)*2*pi*self.a)
425
426 return resistivity
427
428
429 def getApparentResistivitySecondary(self):
430 resistivity = []
431
432 if (self.delPhiSecondary==[]):
433 self.getPotential()
434
435 for i in self.delPhiSecondary:
436 resistivity.append((-i/self.current)*2*pi*self.a)
437
438 return resistivity
439
440 def getApparentResistivityTotal(self):
441 resistivity = []
442
443 if (self.delPhiSecondary==[]):
444 self.getPotential()
445
446 for i in range(len(self.delPhiSecondary)):
447 resistivity.append((-(self.delPhiSecondary[i]+self.delPhiPrimary[i])/self.current)*2*pi*self.a)
448
449 return resistivity
450
451
452 class DipoleDipoleSurvey(DcResistivityForward):
453 """
454 DipoleDipoleSurvey forward modeling
455 """
456 def __init__(self, domain, primaryConductivity, secondaryConductivity,
457 current, a, n, midPoint, directionVector, numElectrodes):
458 super(DipoleDipoleSurvey, self).__init__()
459 """
460 :param domain: domain of the model
461 :type domain: `Domain`
462 :param primaryConductivity: preset primary conductivity data object
463 :type primaryConductivity: data
464 :param secondaryConductivity:preset secondary conductivity data object
465 :type secondaryConductivity: data
466 :param current: amount of current to be injected at the current electrode
467 :type current: float or int
468 :param a: the spacing between current and potential electrodes
469 :type a: float or int
470 :param n: multiple of spacing between electrodes. typicaly interger
471 :type n: float or int
472 :param midPoint: midPoint of the survey, as a list containing x,y coords
473 :type a: list
474 :param directionVector: two element list specifying the direction the
475 survey should extend from the midpoint
476 :type a: list
477 :param numElectrodes: the number of electrodes to be used in the survey
478 must be a multiple of 2 for polepole survey:
479 :type numElectrodes: int
480 """
481 self.domain=domain
482 self.primaryConductivity=primaryConductivity
483 self.secondaryConductivity=secondaryConductivity
484 self.a=a
485 self.n=n
486 self.current=current
487 self.numElectrodes=numElectrodes
488 self.delPhiPrimaryList=[]
489 self.delPhiSecondaryList=[]
490 if (numElectrodes<4 ):
491 raise ValueError("numElectrodes must be atleast 4 for DipoleDipole surveys")
492 if n > (numElectrodes-3):
493 raise ValueError("specified n does not fit max n = %d"%((numElectrodes-2)//2))
494 if len(directionVector) == 2:
495 electrodes=[]
496 start=[]
497 start.append(midPoint[0] - (((numElectrodes-1)*a)/2. * directionVector[0]))
498 start.append(midPoint[1] - (((numElectrodes-1)*a)/2. * directionVector[1]))
499 for i in range(numElectrodes):
500 electrodes.append([start[0]+(directionVector[0]*i*a), start[1]+(directionVector[1]*i*a),0])
501 else:
502 raise NotImplementedError("2d forward model is not yet implemented please provide a 2 component directionVector for a 3d survey")
503 self.electrodes=electrodes
504 self.checkBounds()
505
506
507 def getPotential(self):
508 """
509 Returns 3 list each made up of a number of list containing primary, secondary and total
510 potentials diferences. Each of the lists contain a list for each value of n.
511 """
512 coords=self.domain.getX()
513 pde=LinearPDE(self.domain, numEquations=1)
514 tol=1e-8
515 pde.getSolverOptions().setTolerance(tol)
516 pde.setSymmetryOn()
517 primCon=self.primaryConductivity
518 DIM=self.domain.getDim()
519 x=self.domain.getX()
520 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
521 for i in xrange(DIM-1):
522 xi=x[i]
523 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
524 A = self.secondaryConductivity * kronecker(self.domain)
525 pde.setValue(A=A,q=q)
526
527
528 delPhiSecondaryList = []
529 delPhiPrimaryList = []
530 delPhiTotalList = []
531 for i in range(1,self.n+1): # 1 to n
532 maxR = self.numElectrodes - 2 - (i) #max amount of readings that will fit in the survey
533 delPhiSecondary = []
534 delPhiPrimary = []
535 delPhiTotal = []
536 for j in range(maxR):
537 analyticRsOne=Data(0,(3,),ContinuousFunction(self.domain))
538 analyticRsOne[0]=(coords[0]-self.electrodes[j][0])
539 analyticRsOne[1]=(coords[1]-self.electrodes[j][1])
540 analyticRsOne[2]=(coords[2])
541 rsMagOne=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**0.5
542 analyticRsTwo=Data(0,(3,),ContinuousFunction(self.domain))
543 analyticRsTwo[0]=(coords[0]-self.electrodes[j + 1][0])
544 analyticRsTwo[1]=(coords[1]-self.electrodes[j + 1][1])
545 analyticRsTwo[2]=(coords[2])
546 rsMagTwo=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**0.5
547 rsMagOne+=(whereZero(rsMagOne)*0.0000001)
548 rsMagTwo+=(whereZero(rsMagTwo)*0.0000001)
549 analyticPrimaryPot=(self.current/(2*pi*primCon*rsMagTwo))-(self.current/(2*pi*primCon*rsMagOne))
550
551 analyticRsOnePower=(analyticRsOne[0]**2+analyticRsOne[1]**2+analyticRsOne[2]**2)**1.5
552 analyticRsOnePower = analyticRsOnePower+(whereZero(analyticRsOnePower)*0.0001)
553 analyticRsTwoPower=(analyticRsTwo[0]**2+analyticRsTwo[1]**2+analyticRsTwo[2]**2)**1.5
554 analyticRsTwoPower = analyticRsTwoPower+(whereZero(analyticRsTwoPower)*0.0001)
555
556 gradAnalyticPrimaryPot = Data(0,(3,),ContinuousFunction(self.domain))
557 gradAnalyticPrimaryPot[0] =(self.current/(2*pi*primCon)) * ((analyticRsOne[0]/analyticRsOnePower) - (analyticRsTwo[0]/analyticRsTwoPower))
558 gradAnalyticPrimaryPot[1] =(self.current/(2*pi*primCon)) * ((analyticRsOne[1]/analyticRsOnePower) - (analyticRsTwo[1]/analyticRsTwoPower))
559 gradAnalyticPrimaryPot[2] =(self.current/(2*pi*primCon)) * ((analyticRsOne[2]/analyticRsOnePower) - (analyticRsTwo[2]/analyticRsTwoPower))
560 X=(primCon-self.secondaryConductivity) * (gradAnalyticPrimaryPot)
561 pde.setValue(X=X)
562 u=pde.getSolution()
563 loc=Locator(self.domain,[self.electrodes[1+j+i],self.electrodes[j+i+2]])
564 valPrimary=loc.getValue(analyticPrimaryPot)
565 valSecondary=loc.getValue(u)
566 delPhiPrimary.append(valPrimary[1]-valPrimary[0])
567 delPhiSecondary.append(valSecondary[1]-valSecondary[0])
568 delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
569 delPhiPrimaryList.append(delPhiPrimary)
570 delPhiSecondaryList.append(delPhiSecondary)
571 delPhiTotalList.append(delPhiTotal)
572
573 self.delPhiPrimaryList=delPhiPrimaryList
574 self.delPhiSecondaryList=delPhiSecondaryList
575 self.delPhiTotalList = delPhiTotalList
576 return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]
577
578 def getApparentResistivityPrimary(self):
579 resistivityList = []
580 n=self.n
581 if (self.delPhiPrimaryList==[]):
582 self.getPotential()
583
584 nCount=1
585 for i in self.delPhiPrimaryList:
586 resistivity=[]
587 for j in i:
588 resistivity.append((-j/self.current)*(nCount*(nCount+1)*(nCount+2))*pi*self.a)
589 nCount=nCount+1
590 resistivityList.append(resistivity)
591 return resistivityList
592
593
594 def getApparentResistivitySecondary(self):
595 resistivityList = []
596 n=self.n
597 if (self.delPhiSecondaryList==[]):
598 self.getPotential()
599
600 nCount=1
601 for i in self.delPhiSecondaryList:
602 resistivity=[]
603 for j in i:
604 resistivity.append((-j/self.current)*(nCount*(nCount+1)*(nCount+2))*pi*self.a)
605 nCount=nCount+1
606 resistivityList.append(resistivity)
607 return resistivityList
608
609 def getApparentResistivityTotal(self):
610 resistivityList = []
611 n=self.n
612 if (self.delPhiSecondaryList==[]):
613 self.getPotential()
614
615 nCount=1
616 for i in self.delPhiTotalList:
617 resistivity=[]
618 for j in i:
619 resistivity.append((-j/self.current)*(nCount*(nCount+1)*(nCount+2))*pi*self.a)
620 nCount=nCount+1
621 resistivityList.append(resistivity)
622 return resistivityList
623
624 class PoleDipoleSurvey(DcResistivityForward):
625 """
626 Forward model class for a poledipole setup
627 """
628 def __init__(self, domain, primaryConductivity, secondaryConductivity, current, a, n, midPoint, directionVector, numElectrodes):
629 """
630 :param domain: domain of the model
631 :type domain: `Domain`
632 :param primaryConductivity: preset primary conductivity data object
633 :type primaryConductivity: data
634 :param secondaryConductivity:preset secondary conductivity data object
635 :type secondaryConductivity: data
636 :param current: amount of current to be injected at the current electrode
637 :type current: float or int
638 :param a: the spacing between current and potential electrodes
639 :type a: float or int
640 :param n: multiple of spacing between electrodes. typicaly interger
641 :type n: float or int
642 :param midPoint: midPoint of the survey, as a list containing x,y coords
643 :type a: list
644 :param directionVector: two element list specifying the direction the
645 survey should extend from the midpoint
646 :type a: list
647 :param numElectrodes: the number of electrodes to be used in the survey
648 must be a multiple of 2 for polepole survey:
649 :type numElectrodes: int
650 """
651
652 super(PoleDipoleSurvey, self).__init__()
653
654 self.domain=domain
655 self.primaryConductivity=primaryConductivity
656 self.secondaryConductivity=secondaryConductivity
657 self.a=a
658 self.n=n
659 self.current=current
660 self.numElectrodes=numElectrodes
661 self.delPhiPrimaryList=[]
662 self.delPhiSecondaryList=[]
663 if ((numElectrodes%4) != 0 ):
664 raise ValueError("numElectrodes must be a multiple of 4 for schlumberger surveys")
665 if n > (numElectrodes-3):
666 raise ValueError("specified n does not fit max n = %d"%((numElectrodes-2)//2))
667 if len(directionVector) == 2:
668 electrodes=[]
669 start=[]
670 start.append(midPoint[0] - (((numElectrodes-1)*a)/2. * directionVector[0]))
671 start.append(midPoint[1] - (((numElectrodes-1)*a)/2. * directionVector[1]))
672 for i in range(numElectrodes):
673 electrodes.append([start[0]+(directionVector[0]*i*a), start[1]+(directionVector[1]*i*a),0])
674 else:
675 raise NotImplementedError("2d forward model is not yet implemented please provide a 2 component directionVector for a 3d survey")
676 self.electrodes=electrodes
677 self.checkBounds()
678
679 def getPotential(self):
680 """
681 Returns 3 list each made up of a number of list containing primary, secondary and total
682 potentials diferences. Each of the lists contain a list for each value of n.
683 """
684
685 primCon=self.primaryConductivity
686 coords=self.domain.getX()
687 pde=LinearPDE(self.domain, numEquations=1)
688 tol=1e-8
689 pde.getSolverOptions().setTolerance(tol)
690 pde.setSymmetryOn()
691
692 DIM=self.domain.getDim()
693 x=self.domain.getX()
694 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
695 for i in xrange(DIM-1):
696 xi=x[i]
697 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
698 A = self.secondaryConductivity * kronecker(self.domain)
699 pde.setValue(A=A,q=q)
700
701 delPhiSecondaryList = []
702 delPhiPrimaryList = []
703 delPhiTotalList = []
704 for i in range(1,self.n+1): # 1 to n
705 maxR = self.numElectrodes - 1 - i #max amount of readings that will fit in the survey
706 delPhiSecondary = []
707 delPhiPrimary = []
708 delPhiTotal = []
709 for j in range(maxR):
710 analyticRs=Data(0,(3,),ContinuousFunction(self.domain))
711 analyticRs[0]=(coords[0]-self.electrodes[j][0])
712 analyticRs[1]=(coords[1]-self.electrodes[j][1])
713 analyticRs[2]=(coords[2])
714 rsMag=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**0.5
715 analyticPrimaryPot=(self.current*(1./primCon))/(2*pi*(rsMag+(whereZero(rsMag)*0.0000001))) #the magic number 0.0000001 is to avoid devide by 0
716
717 analyticRsPolePower=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**1.5
718 analyticRsPolePower = analyticRsPolePower+(whereZero(analyticRsPolePower)*0.0000001)
719 gradUPrimary = Data(0,(3,),ContinuousFunction(self.domain))
720 gradUPrimary[0] =(self.current/(2*pi*primCon)) * (analyticRs[0]/analyticRsPolePower)
721 gradUPrimary[1] =(self.current/(2*pi*primCon)) * (analyticRs[1]/analyticRsPolePower)
722 gradUPrimary[2] =(self.current/(2*pi*primCon)) * (analyticRs[2]/analyticRsPolePower)
723 gradUPrimary=-gradUPrimary
724 X=(primCon-self.secondaryConductivity) * gradUPrimary
725 pde.setValue(X=X)
726 u=pde.getSolution()
727 loc=Locator(self.domain,[self.electrodes[i+j],self.electrodes[i+j+1]])
728 valPrimary=loc.getValue(analyticPrimaryPot)
729 valSecondary=loc.getValue(u)
730 delPhiPrimary.append(valPrimary[1]-valPrimary[0])
731 delPhiSecondary.append(valSecondary[1]-valSecondary[0])
732 delPhiTotal.append(delPhiPrimary[j]+delPhiSecondary[j])
733
734 delPhiPrimaryList.append(delPhiPrimary)
735 delPhiSecondaryList.append(delPhiSecondary)
736 delPhiTotalList.append(delPhiTotal)
737
738
739
740 self.delPhiPrimaryList=delPhiPrimaryList
741 self.delPhiSecondaryList=delPhiSecondaryList
742 self.delPhiTotalList = delPhiTotalList
743
744 return [delPhiPrimaryList, delPhiSecondaryList, delPhiTotalList]
745
746
747 def getApparentResistivityPrimary(self):
748 resistivityList = []
749 n=self.n
750 if (self.delPhiPrimaryList==[]):
751 self.getPotential()
752
753 nCount=1
754 for i in self.delPhiPrimaryList:
755 resistivity=[]
756 for j in i:
757 resistivity.append((-j/self.current)*(nCount*(nCount+1))*pi*self.a*2)
758 nCount=nCount+1
759 resistivityList.append(resistivity)
760 return resistivityList
761
762
763 def getApparentResistivitySecondary(self):
764 resistivityList = []
765 n=self.n
766 if (self.delPhiSecondaryList==[]):
767 self.getPotential()
768
769 nCount=1
770 for i in self.delPhiSecondaryList:
771 resistivity=[]
772 for j in i:
773 resistivity.append((-j/self.current)*(nCount*(nCount+1))*pi*self.a*2)
774 nCount=nCount+1
775 resistivityList.append(resistivity)
776 return resistivityList
777
778 def getApparentResistivityTotal(self):
779 resistivityList = []
780 n=self.n
781 if (self.delPhiSecondaryList==[]):
782 self.getPotential()
783
784 nCount=1
785 for i in self.delPhiTotalList:
786 resistivity=[]
787 for j in i:
788 resistivity.append((-j/self.current)*(nCount*(nCount+1))*pi*self.a*2)
789 nCount=nCount+1
790 resistivityList.append(resistivity)
791 return resistivityList
792
793 class PolePoleSurvey(DcResistivityForward):
794 """
795 Forward model class for a polepole setup
796 """
797 def __init__(self, domain, primaryConductivity, secondaryConductivity, current, a, midPoint, directionVector, numElectrodes):
798 """
799 :param domain: domain of the model
800 :type domain: `Domain`
801 :param primaryConductivity: preset primary conductivity data object
802 :type primaryConductivity: data
803 :param secondaryConductivity:preset secondary conductivity data object
804 :type secondaryConductivity: data
805 :param current: amount of current to be injected at the current electrode
806 :type current: float or int
807 :param a: the spacing between current and potential electrodes
808 :type a: float or int
809 :param midPoint: midPoint of the survey, as a list containing x,y coords
810 :type a: list
811 :param directionVector: two element list specifying the direction the
812 survey should extend from the midpoint
813 :type a: list
814 :param numElectrodes: the number of electrodes to be used in the survey
815 must be a multiple of 2 for polepole survey:
816 :type numElectrodes: int
817 """
818
819 super(PolePoleSurvey, self).__init__()
820 self.domain=domain
821 self.primaryConductivity=primaryConductivity
822 self.secondaryConductivity=secondaryConductivity
823 self.a=a
824 self.current=current
825 self.numElectrodes=numElectrodes
826 self.delPhiSecondary=[]
827 self.delPhiPrimary=[]
828 if ((numElectrodes<2) != 0 ):
829 raise ValueError("numElectrodes must be atleast 2 for pole-pole surveys")
830 if len(directionVector) == 2:
831 electrodes = []
832 start=[]
833 start.append(midPoint[0] - (((numElectrodes-1)*a)/2. * directionVector[0]))
834 start.append(midPoint[1] - (((numElectrodes-1)*a)/2. * directionVector[1]))
835 for i in range(numElectrodes):
836 electrodes.append([start[0]+(directionVector[0]*i*a), start[1]+(directionVector[1]*i*a),0])
837 else:
838 raise NotImplementedError("2d forward model is not yet implemented please provide a 2 component directionVector for a 3d survey")
839 self.electrodes=electrodes
840 self.checkBounds()
841
842 def getPotential(self):
843 """
844 returns a list containing 3 lists one for each the primary, secondary
845 and total potential.
846 """
847
848
849 primCon=self.primaryConductivity
850 coords=self.domain.getX()
851 pde=LinearPDE(self.domain, numEquations=1)
852 tol=1e-8
853 pde.getSolverOptions().setTolerance(tol)
854 pde.setSymmetryOn()
855
856 DIM=self.domain.getDim()
857 x=self.domain.getX()
858 q=whereZero(x[DIM-1]-inf(x[DIM-1]))
859 for i in xrange(DIM-1):
860 xi=x[i]
861 q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
862 A = self.secondaryConductivity * kronecker(self.domain)
863 pde.setValue(A=A,q=q)
864
865 delPhiSecondary = []
866 delPhiPrimary = []
867 delPhiTotal = []
868 if(len(self.electrodes[0])==3):
869
870 for i in range(self.numElectrodes-1):
871 analyticRs=Data(0,(3,),ContinuousFunction(self.domain))
872 analyticRs[0]=(coords[0]-self.electrodes[i][0])
873 analyticRs[1]=(coords[1]-self.electrodes[i][1])
874 analyticRs[2]=(coords[2])
875 rsMag=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**0.5
876 analyticPrimaryPot=(self.current*(1./primCon))/(2*pi*(rsMag+(whereZero(rsMag)*0.0000001))) #the magic number 0.0000001 is to avoid devide by 0
877 analyticRsPolePower=(analyticRs[0]**2+analyticRs[1]**2+analyticRs[2]**2)**1.5
878 analyticRsPolePower = analyticRsPolePower+(whereZero(analyticRsPolePower)*0.0000001)
879 gradUPrimary = Data(0,(3,),ContinuousFunction(self.domain))
880 gradUPrimary[0] =(self.current/(2*pi*primCon)) * (analyticRs[0]/analyticRsPolePower)
881 gradUPrimary[1] =(self.current/(2*pi*primCon)) * (analyticRs[1]/analyticRsPolePower)
882 gradUPrimary[2] =(self.current/(2*pi*primCon)) * (analyticRs[2]/analyticRsPolePower)
883 gradUPrimary=-gradUPrimary
884 X=(primCon-self.secondaryConductivity) * gradUPrimary
885 pde.setValue(X=X)
886 u=pde.getSolution()
887 loc=Locator(self.domain,self.electrodes[i+1])
888 delPhiSecondary.append(loc.getValue(u))
889 delPhiPrimary.append(loc.getValue(analyticPrimaryPot))
890 else:
891 raise NotImplementedError("2d forward model is not yet implemented")
892
893 self.delPhiSecondary = delPhiSecondary
894 self.delPhiPrimary = delPhiPrimary
895 for i in range(len(delPhiPrimary)):
896 delPhiTotal.append(delPhiPrimary[i] + delPhiSecondary[i])
897 self.delPhiTotal=delPhiTotal
898 return [delPhiPrimary, delPhiSecondary, delPhiTotal]
899
900 def getApparentResistivityPrimary(self):
901 resistivity = []
902
903 if (self.delPhiPrimary==[]):
904 self.getPotential()
905
906 for i in self.delPhiPrimary:
907 resistivity.append((i/self.current)*2*pi*self.a)
908
909 return resistivity
910
911
912 def getApparentResistivitySecondary(self):
913 resistivity = []
914
915 if (self.delPhiSecondary==[]):
916 self.getPotential()
917
918 for i in self.delPhiSecondary:
919 resistivity.append((i/self.current)*2*pi*self.a)
920
921 return resistivity
922
923 def getApparentResistivityTotal(self):
924 resistivity = []
925
926 if (self.delPhiSecondary==[]):
927 self.getPotential()
928
929 for i in range(len(self.delPhiSecondary)):
930 resistivity.append(((self.delPhiSecondary[i]+self.delPhiPrimary[i])/self.current)*2*pi*self.a)
931
932 return resistivity
933
934
935
936
937
938
939
940
941

  ViewVC Help
Powered by ViewVC 1.1.26