/[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 4114 - (show annotations)
Fri Dec 14 04:24:46 2012 UTC (6 years, 10 months ago) by caltinay
File MIME type: text/x-python
File size: 13199 byte(s)
Time to remove deprecated saveVTK/DX methods from Data and Domain.

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

  ViewVC Help
Powered by ViewVC 1.1.26