/[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 4373 - (show annotations)
Mon Apr 22 03:16:24 2013 UTC (6 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 15991 byte(s)
more on geodetic coordinates
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 `getValue()`, `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 getValue(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
123 for s in range(len(self.__weight)):
124 self.__weight[s] = fw * self.__weight[s]
125 self.__data[s] = fd * self.__data[s]
126
127
128 def _rescaleWeights(self, scale=1., fetch_factor=1.):
129 """
130 rescales the weights such that
131
132 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
133 """
134 if not scale > 0:
135 raise ValueError("Value for scale must be positive.")
136 A=0
137 for s in range(len(self.__weight)):
138 A += integrate( abs(inner(self.__weight[s], self.__data[s])* inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
139 if A > 0:
140 A=sqrt(scale/A)/self.diameter
141 for s in range(len(self.__weight)): self.__weight[s]*=A
142 else:
143 raise ValueError("Rescaling of weights failed.")
144
145 def getDomain(self):
146 """
147 Returns the domain of the forward model.
148
149 :rtype: `Domain`
150 """
151 return self.__domain
152
153 def getCoordinateTransformation(self):
154 """
155 returns the coordinate transformation being used
156
157 :rtype: ``CoordinateTransformation``
158 """
159 return self.__trafo
160
161 def getPDE(self):
162 """
163 Return the underlying PDE.
164
165 :rtype: `LinearPDE`
166 """
167 return self.__pde
168
169 def getDefect(self, result):
170 """
171 Returns the defect value.
172
173 :param result: a result vector
174 :type result: `Vector`
175 :rtype: ``float``
176 """
177 A=0.
178 for s in range(len(self.__weight)):
179 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
180 return A/2
181
182 def getDefectGradient(self, result):
183 Y=0.
184 for s in range(len(self.__weight)):
185 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
186 return Y
187
188 def getSurvey(self, index=None):
189 """
190 Returns the pair (data_index, weight_index), where data_i is the data
191 of survey i, weight_i is the weighting factor for survey i.
192 If index is None, all surveys will be returned in a pair of lists.
193 """
194 if index is None:
195 return self.__data, self.__weight
196 if index>=len(self.__data):
197 raise IndexError("Forward model only has %d surveys"%len(self.__data))
198 return self.__data[index], self.__weight[index]
199
200
201
202 class GravityModel(ForwardModelWithPotential):
203 """
204 Forward Model for gravity inversion as described in the inversion
205 cookbook.
206 """
207 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
208 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
209 """
210 Creates a new gravity model on the given domain with one or more
211 surveys (w, g).
212
213 :param domain: domain of the model
214 :type domain: `Domain`
215 :param w: data weighting factors
216 :type w: ``Vector`` or list of ``Vector``
217 :param g: gravity anomaly data
218 :type g: ``Vector`` or list of ``Vector``
219 :param coordinates: defines coordinate system to be used
220 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
221 :param tol: tolerance of underlying PDE
222 :type tol: positive ``float``
223 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
224 in addition to the top.
225 :type fixPotentialAtBottom: ``bool``
226
227
228 :note: It is advisable to call rescaleWeights() to rescale weights before start inversion.
229 """
230 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
231
232 if not self.getCoordinateTransformation().isCartesian():
233 self.__G = 4*PI*gravity_constant * self.getCoordinateTransformation().getVolumeFactor()
234
235 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
236 A=self.getPDE().createCoeffiecient("A")
237 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
238 self.getPDE().setValue(A=A)
239
240 else:
241 self.__G = 4*PI*gravity_constant
242 self.getPDE().setValue(A=kronecker(self.getDomain()))
243
244 def rescaleWeights(self, scale=1., rho_scale=1.):
245 """
246 rescales the weights such that
247
248 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
249
250 :param scale: scale of data weighting factors
251 :type scale: positive ``float``
252 :param rho_scale: scale of density.
253 :type rho_scale: ``Scalar``
254 """
255 self._rescaleWeights(scale, self.__G*rho_scale)
256
257 def getArguments(self, rho):
258 """
259 Returns precomputed values shared by `getValue()` and `getGradient()`.
260
261 :param rho: a suggestion for the density distribution
262 :type rho: ``Scalar``
263 :return: gravity potential and corresponding gravity field.
264 :rtype: ``Scalar``, ``Vector``
265 """
266 phi = self.getPotential(rho)
267 gravity_force = -grad(phi)
268 return phi, gravity_force
269
270 def getPotential(self, rho):
271 """
272 Calculates the gravity potential for a given density distribution.
273
274 :param rho: a suggestion for the density distribution
275 :type rho: ``Scalar``
276 :return: gravity potential
277 :rtype: ``Scalar``
278 """
279 pde=self.getPDE()
280
281 pde.resetRightHandSideCoefficients()
282 pde.setValue(Y=-self.__G*rho)
283 phi=pde.getSolution()
284
285 return phi
286
287 def getValue(self, rho, phi, gravity_force):
288 """
289 Returns the value of the defect
290
291 :param rho: density distribution
292 :type rho: ``Scalar``
293 :param phi: corresponding potential
294 :type phi: ``Scalar``
295 :param gravity_force: gravity force
296 :type gravity_force: ``Vector``
297 :rtype: ``float``
298 """
299 return self.getDefect(gravity_force)
300
301 def getGradient(self, rho, phi, gravity_force):
302 """
303 Returns the gradient of the defect with respect to density.
304
305 :param rho: density distribution
306 :type rho: ``Scalar``
307 :param phi: corresponding potential
308 :type phi: ``Scalar``
309 :param gravity_force: gravity force
310 :type gravity_force: ``Vector``
311 :rtype: ``Scalar``
312 """
313 pde=self.getPDE()
314 pde.resetRightHandSideCoefficients()
315 pde.setValue(X=self.getDefectGradient(gravity_force))
316 ZT=pde.getSolution()
317 return ZT*(-self.__G)
318
319
320 class MagneticModel(ForwardModelWithPotential):
321 """
322 Forward Model for magnetic inversion as described in the inversion
323 cookbook.
324 """
325 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
326 """
327 Creates a new magnetic model on the given domain with one or more
328 surveys (w, B).
329
330 :param domain: domain of the model
331 :type domain: `Domain`
332 :param w: data weighting factors
333 :type w: ``Vector`` or list of ``Vector``
334 :param B: magnetic field data
335 :type B: ``Vector`` or list of ``Vector``
336 :param tol: tolerance of underlying PDE
337 :type tol: positive ``float``
338 :param coordinates: defines coordinate system to be used
339 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
340 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
341 in addition to the top.
342 :type fixPotentialAtBottom: ``bool``
343 """
344 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
345 self.__background_magnetic_flux_density=interpolate(background_magnetic_flux_density,
346 B[0].getFunctionSpace())
347 if not self.getCoordinateTransformation().isCartesian():
348 self.__F = self.getCoordinateTransformation().getVolumeFactor()
349 self.__B_r=self.__background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
350 self.__B_b=self.__background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
351
352 A=self.getPDE().createCoeffiecient("A")
353 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
354 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
355 self.getPDE().setValue(A=A)
356
357 else:
358 self.getPDE().setValue(A=kronecker(self.getDomain()))
359 self.__B_r=self.__background_magnetic_flux_density
360 self.__B_b=self.__background_magnetic_flux_density
361
362 def rescaleWeights(self, scale=1., k_scale=1.):
363 """
364 rescales the weights such that
365
366 *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*
367
368 :param scale: scale of data weighting factors
369 :type scale: positive ``float``
370 :param k_scale: scale of susceptibility.
371 :type k_scale: ``Scalar``
372 """
373 self._rescaleWeights(scale, inner(self.__background_magnetic_flux_density,1/self.edge_lengths ) * k_scale)
374
375 def getArguments(self, k):
376 """
377 Returns precomputed values shared by `getValue()` and `getGradient()`.
378
379 :param k: susceptibility
380 :type k: ``Scalar``
381 :return: scalar magnetic potential and corresponding magnetic field
382 :rtype: ``Scalar``, ``Vector``
383 """
384 phi = self.getPotential(k)
385 magnetic_flux_density = k * self.__B_b -grad(phi)
386 return phi, magnetic_flux_density
387
388 def getPotential(self, k):
389 """
390 Calculates the magnetic potential for a given susceptibility.
391
392 :param k: susceptibility
393 :type k: ``Scalar``
394 :return: magnetic potential
395 :rtype: ``Scalar``
396 """
397 pde=self.getPDE()
398
399 pde.resetRightHandSideCoefficients()
400 pde.setValue(X = k* self.__B_r)
401 phi=pde.getSolution()
402
403 return phi
404
405 def getValue(self, k, phi, magnetic_flux_density):
406 """
407 Returns the value of the defect.
408
409 :param k: susceptibility
410 :type k: ``Scalar``
411 :param phi: corresponding potential
412 :type phi: ``Scalar``
413 :param magnetic_flux_density: magnetic field
414 :type magnetic_flux_density: ``Vector``
415 :rtype: ``float``
416 """
417 return self.getDefect(magnetic_flux_density)
418
419 def getGradient(self, k, phi, magnetic_flux_density):
420 """
421 Returns the gradient of the defect with respect to susceptibility.
422
423 :param k: susceptibility
424 :type k: ``Scalar``
425 :param phi: corresponding potential
426 :type phi: ``Scalar``
427 :param magnetic_flux_density: magnetic field
428 :type magnetic_flux_density: ``Vector``
429 :rtype: ``Scalar``
430 """
431 Y=self.getDefectGradient(magnetic_flux_density)
432 pde=self.getPDE()
433 pde.resetRightHandSideCoefficients()
434 pde.setValue(X=Y)
435 YT=pde.getSolution()
436 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
437

  ViewVC Help
Powered by ViewVC 1.1.26