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

  ViewVC Help
Powered by ViewVC 1.1.26