/[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 5790 - (show annotations)
Wed Sep 9 04:56:37 2015 UTC (3 years ago) by jduplessis
File MIME type: text/x-python
File size: 9188 byte(s)
wrong order for subtraction. based on ralfs notes.

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

  ViewVC Help
Powered by ViewVC 1.1.26