/[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 4433 - (show annotations)
Fri May 31 12:09:58 2013 UTC (6 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 16102 byte(s)
some clarifications on geodetic coordinates. 
order of background magnetic flux density component has been corrected: input is now B_east, B_north, B_vertical.


1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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-2013 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 import Data, Vector, Scalar, Function
29 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
30 from .coordinates import makeTranformation
31 from esys.escript.util import *
32 from math import pi as PI
33 import numpy as np
34
35
36 class ForwardModel(object):
37 """
38 An abstract forward model that can be plugged into a cost function.
39 Subclasses need to implement `getDefect()`, `getGradient()`, and possibly
40 `getArguments()` and 'getCoordinateTransformation'.
41 """
42 def __init__(self):
43 pass
44
45 def getArguments(self, x):
46 return ()
47
48 def getDefect(self, x, *args):
49 raise NotImplementedError
50
51 def getGradient(self, x, *args):
52 raise NotImplementedError
53
54 def getCoordinateTransformation(self):
55 return None
56
57
58 class ForwardModelWithPotential(ForwardModel):
59 """
60 Base class for a forward model using a potential such as magnetic or
61 gravity. It defines a cost function::
62
63 defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) ) ** 2 )
64
65 where s runs over the survey, weight_i are weighting factors, data_i are
66 the data, and r_i are the results produced by the forward model.
67 It is assumed that the forward model is produced through postprocessing
68 of the solution of a potential PDE.
69 """
70 def __init__(self, domain, w, data, coordinates=None,
71 fixPotentialAtBottom=False,
72 tol=1e-8):
73 """
74 initializes a new forward model with potential.
75
76 :param domain: domain of the model
77 :type domain: `Domain`
78 :param w: data weighting factors
79 :type w: ``Vector`` or list of ``Vector``
80 :param data: data
81 :type data: ``Vector`` or list of ``Vector``
82 :param coordinates: defines coordinate system to be used
83 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
84 :param tol: tolerance of underlying PDE
85 :type tol: positive ``float``
86 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
87 in addition to the top.
88 :type fixPotentialAtBottom: ``bool``
89 """
90 super(ForwardModelWithPotential, self).__init__()
91 self.__domain = domain
92 self.__trafo=makeTranformation(domain, coordinates)
93
94 try:
95 n=len(w)
96 m=len(data)
97 if not m == n:
98 raise ValueError("Length of weight and data must be the same.")
99 self.__weight = w
100 self.__data = data
101 except TypeError:
102 self.__weight = [w]
103 self.__data = [data]
104
105 BX = boundingBox(domain)
106 DIM = domain.getDim()
107 x = domain.getX()
108 self.__pde=LinearSinglePDE(domain)
109 self.__pde.getSolverOptions().setTolerance(tol)
110 self.__pde.setSymmetryOn()
111 z=x[DIM-1]
112 q0=whereZero(z-BX[DIM-1][1])
113 if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
114 self.__pde.setValue(q=q0)
115
116 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
117 self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
118
119 if not self.__trafo.isCartesian():
120 fd=1./self.__trafo.getScalingFactors()
121 fw=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
122 for s in range(len(self.__weight)):
123 self.__weight[s] = fw * self.__weight[s]
124 self.__data[s] = fd * self.__data[s]
125
126 def _rescaleWeights(self, scale=1., fetch_factor=1.):
127 """
128 rescales the weights such that
129
130 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
131 """
132 if not scale > 0:
133 raise ValueError("Value for scale must be positive.")
134 A=0
135 for s in range(len(self.__weight)):
136 A += integrate( abs(inner(self.__weight[s], self.__data[s])* inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
137 if A > 0:
138 A=sqrt(scale/A)/self.diameter
139 for s in range(len(self.__weight)): self.__weight[s]*=A
140 else:
141 raise ValueError("Rescaling of weights failed.")
142
143 def getDomain(self):
144 """
145 Returns the domain of the forward model.
146
147 :rtype: `Domain`
148 """
149 return self.__domain
150
151 def getCoordinateTransformation(self):
152 """
153 returns the coordinate transformation being used
154
155 :rtype: ``CoordinateTransformation``
156 """
157 return self.__trafo
158
159 def getPDE(self):
160 """
161 Return the underlying PDE.
162
163 :rtype: `LinearPDE`
164 """
165 return self.__pde
166
167 def _getDefect(self, result):
168 """
169 Returns the defect value.
170
171 :param result: a result vector
172 :type result: `Vector`
173 :rtype: ``float``
174 """
175 A=0.
176 for s in range(len(self.__weight)):
177 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
178 return A/2
179
180 def getDefectGradient(self, result):
181 Y=0.
182 for s in range(len(self.__weight)):
183 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
184 return Y
185
186 def getSurvey(self, index=None):
187 """
188 Returns the pair (data_index, weight_index), where data_i is the data
189 of survey i, weight_i is the weighting factor for survey i.
190 If index is None, all surveys will be returned in a pair of lists.
191 """
192 if index is None:
193 return self.__data, self.__weight
194 if index>=len(self.__data):
195 raise IndexError("Forward model only has %d surveys"%len(self.__data))
196 return self.__data[index], self.__weight[index]
197
198
199
200 class GravityModel(ForwardModelWithPotential):
201 """
202 Forward Model for gravity inversion as described in the inversion
203 cookbook.
204 """
205 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
206 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
207 """
208 Creates a new gravity model on the given domain with one or more
209 surveys (w, g).
210
211 :param domain: domain of the model
212 :type domain: `Domain`
213 :param w: data weighting factors
214 :type w: ``Vector`` or list of ``Vector``
215 :param g: gravity anomaly data
216 :type g: ``Vector`` or list of ``Vector``
217 :param coordinates: defines coordinate system to be used
218 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
219 :param tol: tolerance of underlying PDE
220 :type tol: positive ``float``
221 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
222 in addition to the top.
223 :type fixPotentialAtBottom: ``bool``
224
225
226 :note: It is advisable to call rescaleWeights() to rescale weights before start inversion.
227 """
228 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
229
230 if not self.getCoordinateTransformation().isCartesian():
231 self.__G = 4*PI*gravity_constant * self.getCoordinateTransformation().getVolumeFactor()
232
233 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
234 A=self.getPDE().createCoefficient("A")
235 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
236 self.getPDE().setValue(A=A)
237 else: # cartesian
238 self.__G = 4*PI*gravity_constant
239 self.getPDE().setValue(A=kronecker(self.getDomain()))
240
241 def rescaleWeights(self, scale=1., rho_scale=1.):
242 """
243 rescales the weights such that
244
245 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
246
247 :param scale: scale of data weighting factors
248 :type scale: positive ``float``
249 :param rho_scale: scale of density.
250 :type rho_scale: ``Scalar``
251 """
252 self._rescaleWeights(scale, self.__G*rho_scale)
253
254 def getArguments(self, rho):
255 """
256 Returns precomputed values shared by `getDefect()` and `getGradient()`.
257
258 :param rho: a suggestion for the density distribution
259 :type rho: ``Scalar``
260 :return: gravity potential and corresponding gravity field.
261 :rtype: ``Scalar``, ``Vector``
262 """
263 phi = self.getPotential(rho)
264 gravity_force = -grad(phi)
265 return phi, gravity_force
266
267 def getPotential(self, rho):
268 """
269 Calculates the gravity potential for a given density distribution.
270
271 :param rho: a suggestion for the density distribution
272 :type rho: ``Scalar``
273 :return: gravity potential
274 :rtype: ``Scalar``
275 """
276 pde=self.getPDE()
277
278 pde.resetRightHandSideCoefficients()
279 pde.setValue(Y=-self.__G*rho)
280 phi=pde.getSolution()
281
282 return phi
283
284 def getDefect(self, rho, phi, gravity_force):
285 """
286 Returns the value of the defect
287
288 :param rho: density distribution
289 :type rho: ``Scalar``
290 :param phi: corresponding potential
291 :type phi: ``Scalar``
292 :param gravity_force: gravity force
293 :type gravity_force: ``Vector``
294 :rtype: ``float``
295 """
296 return self._getDefect(gravity_force)
297
298 def getGradient(self, rho, phi, gravity_force):
299 """
300 Returns the gradient of the defect with respect to density.
301
302 :param rho: density distribution
303 :type rho: ``Scalar``
304 :param phi: corresponding potential
305 :type phi: ``Scalar``
306 :param gravity_force: gravity force
307 :type gravity_force: ``Vector``
308 :rtype: ``Scalar``
309 """
310 pde=self.getPDE()
311 pde.resetRightHandSideCoefficients()
312 pde.setValue(X=self.getDefectGradient(gravity_force))
313 ZT=pde.getSolution()
314 return ZT*(-self.__G)
315
316
317 class MagneticModel(ForwardModelWithPotential):
318 """
319 Forward Model for magnetic inversion as described in the inversion
320 cookbook.
321 """
322 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
323 """
324 Creates a new magnetic model on the given domain with one or more
325 surveys (w, B).
326
327 :param domain: domain of the model
328 :type domain: `Domain`
329 :param w: data weighting factors
330 :type w: ``Vector`` or list of ``Vector``
331 :param B: magnetic field data
332 :type B: ``Vector`` or list of ``Vector``
333 :param tol: tolerance of underlying PDE
334 :type tol: positive ``float``
335 :param background_magnetic_flux_density: background magnetic flux density (in Teslar) with components (B_east, B_north, B_vertical)
336 :type background_magnetic_flux_density: ``Vector`` or list of `float`
337 :param coordinates: defines coordinate system to be used
338 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
339 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
340 in addition to the top.
341 :type fixPotentialAtBottom: ``bool``
342 """
343 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
344 self.__background_magnetic_flux_density = \
345 interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
346 if not self.getCoordinateTransformation().isCartesian():
347 self.__F = self.getCoordinateTransformation().getVolumeFactor()
348 self.__B_r=self.__background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
349 self.__B_b=self.__background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
350
351 A=self.getPDE().createCoefficient("A")
352 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
353 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
354 self.getPDE().setValue(A=A)
355 else: # cartesian
356 self.getPDE().setValue(A=kronecker(self.getDomain()))
357 self.__B_r=self.__background_magnetic_flux_density
358 self.__B_b=self.__background_magnetic_flux_density
359
360 def rescaleWeights(self, scale=1., k_scale=1.):
361 """
362 rescales the weights such that
363
364 *sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale*
365
366 :param scale: scale of data weighting factors
367 :type scale: positive ``float``
368 :param k_scale: scale of susceptibility.
369 :type k_scale: ``Scalar``
370 """
371 self._rescaleWeights(scale, inner(self.__background_magnetic_flux_density,1/self.edge_lengths ) * k_scale)
372
373 def getArguments(self, k):
374 """
375 Returns precomputed values shared by `getDefect()` and `getGradient()`.
376
377 :param k: susceptibility
378 :type k: ``Scalar``
379 :return: scalar magnetic potential and corresponding magnetic field
380 :rtype: ``Scalar``, ``Vector``
381 """
382 phi = self.getPotential(k)
383 magnetic_flux_density = k * self.__B_b -grad(phi)
384 return phi, magnetic_flux_density
385
386 def getPotential(self, k):
387 """
388 Calculates the magnetic potential for a given susceptibility.
389
390 :param k: susceptibility
391 :type k: ``Scalar``
392 :return: magnetic potential
393 :rtype: ``Scalar``
394 """
395 pde=self.getPDE()
396
397 pde.resetRightHandSideCoefficients()
398 pde.setValue(X = k* self.__B_r)
399 phi=pde.getSolution()
400
401 return phi
402
403 def getDefect(self, k, phi, magnetic_flux_density):
404 """
405 Returns the value of the defect.
406
407 :param k: susceptibility
408 :type k: ``Scalar``
409 :param phi: corresponding potential
410 :type phi: ``Scalar``
411 :param magnetic_flux_density: magnetic field
412 :type magnetic_flux_density: ``Vector``
413 :rtype: ``float``
414 """
415 return self._getDefect(magnetic_flux_density)
416
417 def getGradient(self, k, phi, magnetic_flux_density):
418 """
419 Returns the gradient of the defect with respect to susceptibility.
420
421 :param k: susceptibility
422 :type k: ``Scalar``
423 :param phi: corresponding potential
424 :type phi: ``Scalar``
425 :param magnetic_flux_density: magnetic field
426 :type magnetic_flux_density: ``Vector``
427 :rtype: ``Scalar``
428 """
429 Y=self.getDefectGradient(magnetic_flux_density)
430 pde=self.getPDE()
431 pde.resetRightHandSideCoefficients()
432 pde.setValue(X=Y)
433 YT=pde.getSolution()
434 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
435

  ViewVC Help
Powered by ViewVC 1.1.26