1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
|
##############################################################################
#
# Copyright (c) 2003-2018 by The University of Queensland
# http://www.uq.edu.au
#
# Primary Business: Queensland, Australia
# Licensed under the Apache License, version 2.0
# http://www.apache.org/licenses/LICENSE-2.0
#
# Development until 2012 by Earth Systems Science Computational Center (ESSCC)
# Development 2012-2013 by School of Earth Sciences
# Development from 2014 by Centre for Geoscience Computing (GeoComp)
#
##############################################################################
"""Forward model for DC Resistivity"""
from __future__ import division, print_function
__copyright__="""Copyright (c) 2003-2018 by The University of Queensland
http://www.uq.edu.au
Primary Business: Queensland, Australia"""
__license__="""Licensed under the Apache License, version 2.0
http://www.apache.org/licenses/LICENSE-2.0"""
__url__="https://launchpad.net/escript-finley"
__all__ = ['DcRes']
from .base import ForwardModel
from esys.downunder.coordinates import makeTransformation
from esys.escript import Scalar, DiracDeltaFunctions
from esys.escript.pdetools import Locator
from esys.escript.linearPDEs import LinearPDE
from esys.escript.util import *
class DcRes(ForwardModel):
"""
Forward Model for DC resistivity, with a given source pair.
The cost function is defined as:
:math: defect = 1/2 (sum_s sum_pq w_pqs * ((phi_sp-phi_sq)-v_pqs)**2
"""
def __init__(self, domain, locator, delphiIn, sampleTags, phiPrimary,
sigmaPrimary, w=1., coordinates=None, tol=1e-8,
saveMemory=True, b=None):
"""
setup new forward model
:param domain: the domain of the model
:type: escript domain
:param locator: contains locator to the measurement pairs
:type: `list` of ``Locator``
:param: delphiIn: this is v_pq, the potential difference for the
current source and a set of measurement pairs.
A list of measured potential differences is expected.
Note this should be the secondary potential only.
:type delphiIn: tuple
:param sampleTags: tags of measurement points from which potential
differences will be calculated.
:type sampleTags: list of tuples
:param phiPrimary: primary potential.
:type phiPrimary: `Scalar`
"""
super(DcRes, self).__init__()
if not isinstance(sampleTags, list):
raise ValueError("sampleTags must be a list.")
if not len(sampleTags) == len(delphiIn):
raise ValueError("sampleTags and delphiIn must have the same length.")
if not len(sampleTags)>0:
raise ValueError("sampleTags list is empty.")
if not isinstance(sampleTags[0], tuple) and not isinstance(sampleTags[0], list):
raise ValueError("sampleTags must be a list of tuples or a list of lists.")
if isinstance(w, float) or isinstance(w, int):
w = [float(w) for z in delphiIn]
self.__w = w
else:
self.__w = w
if not len(w) == len(delphiIn):
raise ValueError("Number of confidence factors and number of potential input values don't match.")
self.__trafo = makeTransformation(domain, coordinates)
if not self.getCoordinateTransformation().isCartesian():
raise ValueError("Non-Cartesian Coordinates are not supported yet.")
if not isinstance(locator, Locator):
raise ValueError("locator must be an escript locator object")
self.__domain = domain
self.__tol = tol
self.__locator = locator
self.delphiIn = delphiIn
self.__sampleTags = sampleTags
self.__sigmaPrimary = sigmaPrimary
self.__phiPrimary = phiPrimary
self.__pde = None
if not saveMemory:
self.__pde=self.setUpPDE()
def getDomain(self):
"""
Returns the domain of the forward model.
:rtype: `Domain`
"""
return self.__domain
def getCoordinateTransformation(self):
"""
returns the coordinate transformation being used
:rtype: ``CoordinateTransformation``
"""
return self.__trafo
def getPrimaryPotential(self):
"""
returns the primary potential
:rtype: `Data`
"""
return self.__phiPrimary
def setUpPDE(self):
"""
Return the underlying PDE.
:rtype: `LinearPDE`
"""
if self.__pde is None:
dom=self.__domain
x = dom.getX()
DIM=dom.getDim()
q=whereZero(x[DIM-1]-inf(x[DIM-1]))
for i in range(DIM-1):
xi=x[i]
q+=whereZero(xi-inf(xi))+whereZero(xi-sup(xi))
pde=LinearPDE(dom, numEquations=1)
pde.getSolverOptions().setTolerance(self.__tol)
pde.setSymmetryOn()
A=pde.createCoefficient('A')
X=pde.createCoefficient('X')
pde.setValue(A=A, X=X, q=q)
else:
pde=self.__pde
pde.resetRightHandSideCoefficients()
return pde
def getArguments(self, sigma):
"""
Returns precomputed values shared by `getDefect()` and `getGradient()`.
:param sigma: conductivity
:type sigma: ``Data`` of shape (1,)
:return: phi
:rtype: ``Data`` of shape (1,)
"""
pde = self.setUpPDE()
X = (self.__sigmaPrimary - sigma) * grad(self.__phiPrimary)
pde.setValue(A=sigma * kronecker(self.__domain), X=X)
phi = pde.getSolution()
loc = self.__locator
loc_phi = loc.getValue(phi)
return phi, loc_phi
def getDefect(self, sigma, phi, loc_phi):
"""
Returns the defect value.
:param sigma: a suggestion for conductivity
:type sigma: ``Data`` of shape (1,)
:param phi: potential field
:type phi: ``Data`` of shape (1,)
:rtype: ``float``
"""
val=loc_phi
if not isinstance(val,list):
tmp=val
val=[tmp]
# print "val=",val
length=len(val)
# print self.__sampleTags[0]
if ((self.__sampleTags[0][1]!="-" and (length%2) != 0) or \
(self.__sampleTags[0][1]!="-" and length/2 != len(self.delphiIn))):
raise ValueError("length of locator is wrong")
delphi_calc=[]
if self.__sampleTags[0][1] != "-":
for i in range(0,length,2):
delphi_calc.append(val[i+1]-val[i])
else:
for i in range(length):
delphi_calc.append(val[i])
A=0
if (self.__sampleTags[0][1] != "-"):
for i in range(length//2):
A += (self.__w[i]*(delphi_calc[i]-self.delphiIn[i])**2)
else:
for i in range(length):
A += (self.__w[i]*(delphi_calc[i]-self.delphiIn[i])**2)
return A/2
def getGradient(self, sigma, phi, loc_phi):
"""
Returns the gradient of the defect with respect to density.
:param sigma: a suggestison for conductivity
:type sigma: ``Data`` of shape (1,)
:param phi: potential field
:type phi: ``Data`` of shape (1,)
"""
val=loc_phi
if not isinstance(val,list):
tmp=val
val=[tmp]
sampleTags=self.__sampleTags
jointSamples={}
for i in range(0,2*len(sampleTags),2): #2*len because sample tags is a list of tuples
half = i//2
if sampleTags[half][1]!="-":
tmp=(val[i+1]-val[i]-self.delphiIn[half])*self.__w[i]
else:
tmp=(val[i]-self.delphiIn[i//2]) *self.__w[i]
sample = sampleTags[half]
if sample[1]!="-":
if sample[0] in jointSamples.keys():
jointSamples[sample[0]].append((sample[1], -tmp))
else:
jointSamples[sampleTags[half][0]]=[(sample[1],-tmp)]
if sample[1] in jointSamples.keys():
jointSamples[sample[1]].append((sample[0], tmp))
else:
jointSamples[sample[1]]=[(sample[0], tmp)]
else:
if sample[0] in jointSamples.keys():
jointSamples[sample[0]].append((sample[1], tmp))
else:
jointSamples[sampleTags[half][0]]=[(sample[1],tmp)]
pde =self.setUpPDE()
dom=self.__domain
# conPrimary=self.__sigmaPrimary
# APrimary = conPrimary * kronecker(dom)
y_dirac = Scalar(0,DiracDeltaFunctions(dom))
for i in jointSamples:
total=0
for j in jointSamples[i]:
total+=j[1]
y_dirac.setTaggedValue(i,total)
pde.setValue(A=sigma*kronecker(dom), y_dirac=y_dirac)
u=pde.getSolution()
retVal=-inner(grad(u),grad(phi+self.__phiPrimary))
return retVal
|
