/[escript]/trunk/downunder/py_src/forwardmodels/dcresistivity.py
ViewVC logotype

Contents of /trunk/downunder/py_src/forwardmodels/dcresistivity.py

Parent Directory Parent Directory | Revision Log Revision Log


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

  ViewVC Help
Powered by ViewVC 1.1.26