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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4063 - (show annotations)
Fri Nov 9 06:15:00 2012 UTC (6 years, 9 months ago) by caltinay
File MIME type: text/x-python
File size: 11442 byte(s)
Fixed a typo in method description.

1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 by University of Queensland
5 # http://www.uq.edu.au
6 #
7 # Primary Business: Queensland, Australia
8 # Licensed under the Open Software License version 3.0
9 # http://www.opensource.org/licenses/osl-3.0.php
10 #
11 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 """Collection of forward models that define the inversion problem"""
17
18 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = ['ForwardModel','ForwardModelWithPotential','GravityModel','MagneticModel']
26
27 from esys.escript import unitsSI as U
28 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
29 from esys.escript.util import *
30 from esys.escript import Data, Vector, Scalar, Function
31 from math import pi as PI
32
33 class ForwardModel(object):
34 """
35 An abstract forward model that can be plugged into a cost function.
36 Subclasses need to implement `getValue()`, `getGradient()`, and possibly
37 `getArguments()`.
38 """
39 def __init__(self):
40 pass
41
42 def getArguments(self, x):
43 return ()
44
45 def getValue(self, x, *args):
46 raise NotImplementedError
47
48 def getGradient(self, x, *args):
49 raise NotImplementedError
50
51 class ForwardModelWithPotential(ForwardModel):
52 """
53 Base class for a forward model using a potential such as magnetic or
54 gravity. It defines a cost function
55
56 defect = 1/2 sum_s integrate( weight_i[s] * ( r_i - data_i[s] )**2 )
57
58 where s runs over the survey, weight_i are weighting factors, data_i are
59 the data, and r_i are the results produced by the forward model.
60 It is assumed that the forward model is produced through postprocessing
61 of the solution of a potential PDE.
62 """
63 def __init__(self, domain, weight, data, fix_all_faces=False, tol=1e-8):
64 """
65 initialization.
66
67 :param domain: domain of the model
68 :type domain: `esys.escript.Domain`
69 :param weight: data weighting factors
70 :type weight: `Vector` or list of `Vector`
71 :param data: data
72 :type data: `Vector` or list of `Vector`
73 :param fix_all_faces: if ``true`` all faces of the domain are fixed
74 otherwise only the top surface
75 :type fix_all_faces: ``bool``
76 :param tol: tolerance of underlying PDE
77 :type tol: positive ``float``
78
79 """
80 super(ForwardModelWithPotential, self).__init__()
81 self.__domain = domain
82
83 try:
84 n=len(weight)
85 m=len(data)
86 if m != n:
87 raise ValueError("Length of weight and g must be the same.")
88 self.__weight = weight
89 self.__data = data
90 except TypeError:
91 self.__weight = [weight]
92 self.__data = [data]
93
94 A=0
95 for s in xrange(len(self.__weight)):
96 A2 = integrate(inner(self.__weight[s], self.__data[s]**2))
97 if A2 < 0:
98 raise ValueError("Negative weighting factor for survey %s"%s)
99 A=max(A2, A)
100 if not A > 0:
101 raise ValueError("No reference data set.")
102
103 BX = boundingBox(domain)
104 DIM = domain.getDim()
105 x = domain.getX()
106 self.__pde=LinearSinglePDE(domain)
107 self.__pde.getSolverOptions().setTolerance(tol)
108 self.__pde.setSymmetryOn()
109
110 if fix_all_faces:
111 constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])
112 for i in xrange(DIM-1):
113 constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])
114 else:
115 constraint=whereZero(x[DIM-1]-BX[DIM-1][1])
116 self.__pde.setValue(q=constraint)
117
118 def getDomain(self):
119 """
120 Returns the domain of the forward model.
121 """
122 return self.__domain
123
124 def getPDE(self):
125 """
126 Return the underlying PDE.
127
128 :rtype: `LinearPDE`
129 """
130 return self.__pde
131
132 def getDefect(self, result):
133 """
134 Returns the defect value.
135
136 :param result: a result vector
137 :type result: `Vector`
138 :rtype: ``float``
139 """
140 A=0.
141 for s in xrange(len(self.__weight)):
142 A = inner(self.__weight[s], (result-self.__data[s])**2) + A
143 return 0.5*integrate(A)
144
145 def getSurvey(self, index=None):
146 """
147 Returns the pair (data_index, weight_index), where data_i is the data
148 of survey i, weight_i is the weighting factor for survey i.
149 If index is None, all surveys will be returned in a pair of lists.
150 """
151 if index is None:
152 return self.__data, self.__weight
153 if index>=len(self.__data):
154 raise IndexError("Forward model only has %d surveys"%len(self.__data))
155 return self.__data[index], self.__weight[index]
156
157
158 class GravityModel(ForwardModelWithPotential):
159 """
160 Forward Model for gravity inversion as described in the inversion
161 cookbook.
162 """
163 def __init__(self, domain, chi, g, fix_all_faces=False, gravity_constant=U.Gravitational_Constant, tol=1e-8):
164 """
165 Creates a new gravity model on the given domain with one or more
166 surveys (chi, g).
167
168 :param domain: domain of the model
169 :type domain: `esys.escript.Domain`
170 :param chi: data weighting factors
171 :type chi: `Vector` or list of `Vector`
172 :param g: gravity anomaly data
173 :type g: `Vector` or list of `Vector`
174 :param fix_all_faces: if ``true`` all faces of the domain are fixed
175 otherwise only the top surface
176 :type fix_all_faces: ``bool``
177 :param tol: tolerance of underlying PDE
178 :type tol: positive ``float``
179 """
180 super(GravityModel, self).__init__(domain, chi, g, fix_all_faces, tol)
181
182 self.__G = gravity_constant
183 self.getPDE().setValue(A=kronecker(self.getDomain()))
184
185 def getArguments(self, rho):
186 """
187 Returns precomputed values shared by getValue() and getGradient().
188
189 :param rho: a suggestion for the density distribution
190 :type rho: `Scalar`
191 :return: gravity potential and corresponding gravity field.
192 :rtype: ``Scalar``, ``Vector``
193 """
194 phi = self.getPotential(rho)
195 gravity_force = -grad(phi)
196 return phi, gravity_force
197
198 def getPotential(self, rho):
199 """
200 Calculates the gravity potential for a given density distribution.
201
202 :param rho: a suggestion for the density distribution
203 :type rho: `Scalar`
204 :return: gravity potential
205 :rtype: ``Scalar``
206 """
207 pde=self.getPDE()
208
209 pde.resetRightHandSideCoefficients()
210 pde.setValue(Y=-4.*PI*self.__G*rho)
211 phi=pde.getSolution()
212
213 return phi
214
215 def getValue(self, rho, phi, gravity_force):
216 """
217 Returns the value of the defect
218
219 :param rho: density distribution
220 :type rho: `Scalar`
221 :param phi: corresponding potential
222 :type phi: `Scalar`
223 :param gravity_force: gravity force
224 :type gravity_force: `Vector`
225 :rtype: ``float``
226 """
227 return self.getDefect(gravity_force)
228
229 def getGradient(self, rho, phi, gravity_force):
230 """
231 Returns the gradient of the defect with respect to density.
232
233 :param rho: density distribution
234 :type rho: `Scalar`
235 :param phi: corresponding potential
236 :type phi: `Scalar`
237 :param gravity_force: gravity force
238 :type gravity_force: `Vector`
239 :rtype: `Scalar`
240 """
241 pde=self.getPDE()
242 g, chi = self.getSurvey()
243
244 Z=0.
245 for s in xrange(len(chi)):
246 Z = chi[s] * (g[s]-gravity_force) + Z
247
248 pde.resetRightHandSideCoefficients()
249 pde.setValue(X=Z)
250 ZT=pde.getSolution()
251 return ZT*(-4*PI*self.__G)
252
253
254 class MagneticModel(ForwardModelWithPotential):
255 """
256 Forward Model for magnetic inversion as described in the inversion
257 cookbook.
258 """
259 def __init__(self, domain, chi, B, background_field, fix_all_faces=False, tol=1e-8):
260 """
261 Creates a new magnetic model on the given domain with one or more
262 surveys (chi, B).
263
264 :param domain: domain of the model
265 :type domain: `Domain`
266 :param chi: data weighting factors
267 :type chi: `Vector` or list of `Vector`
268 :param B: magnetic field data
269 :type B: `Vector` or list of `Vector`
270 :param fix_all_faces: if ``true`` all faces of the domain are fixed
271 otherwise only the top surface
272 :type fix_all_faces: ``bool``
273 :param tol: tolerance of underlying PDE
274 :type tol: positive ``float``
275 """
276 super(MagneticModel, self).__init__(domain, chi, B, fix_all_faces, tol)
277 self.__background_field=interpolate(background_field, Function(self.getDomain()))
278 self.getPDE().setValue(A=kronecker(self.getDomain()))
279
280 def getArguments(self, k):
281 """
282 Returns precomputed values shared by getValue() and getGradient().
283
284 :param k: susceptibility
285 :type k: `Scalar`
286 :return: scalar magnetic potential and corresponding magnetic field
287 :rtype: ``Scalar``, ``Vector``
288 """
289 phi = self.getPotential(k)
290 magnetic_field = (1+k) * self.__background_field -grad(phi)
291 return phi, magnetic_field
292
293 def getPotential(self, k):
294 """
295 Calculates the magnetic potential for a given susceptibility.
296
297 :param k: susceptibility
298 :type k: `Scalar`
299 :return: magnetic potential
300 :rtype: ``Scalar``
301 """
302 pde=self.getPDE()
303
304 pde.resetRightHandSideCoefficients()
305 pde.setValue(X = (1+k)* self.__background_field)
306 phi=pde.getSolution()
307
308 return phi
309
310 def getValue(self, k, phi, magnetic_field):
311 """
312 Returns the value of the defect.
313
314 :param k: susceptibility
315 :type k: `Scalar`
316 :param phi: corresponding potential
317 :type phi: `Scalar`
318 :param magnetic_field: magnetic field
319 :type magnetic_field: `Vector`
320 :rtype: ``float``
321 """
322 return self.getDefect(magnetic_field)
323
324 def getGradient(self, k, phi, magnetic_field):
325 """
326 Returns the gradient of the defect with respect to susceptibility.
327
328 :param k: susceptibility
329 :type k: `Scalar`
330 :param phi: corresponding potential
331 :type phi: `Scalar`
332 :param magnetic_field: magnetic field
333 :type magnetic_field: `Vector`
334 :rtype: `Scalar`
335 """
336 pde=self.getPDE()
337 B, chi = self.getSurvey()
338
339 Y=0.
340 for s in xrange(len(chi)):
341 Y = chi[s] * (magnetic_field-B[s]) + Y
342
343 pde.resetRightHandSideCoefficients()
344 pde.setValue(X=Y)
345 YT=pde.getSolution()
346 return inner(Y-grad(YT),self.__background_field)
347

  ViewVC Help
Powered by ViewVC 1.1.26