/[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 4095 - (show annotations)
Wed Dec 5 05:32:22 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 12728 byte(s)
A bit of doco cleanup.

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
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 = scale/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, scale=1., tol=1e-8):
66 """
67 initializes a new forward model with potential.
68
69 :param domain: domain of the model
70 :type domain: `esys.escript.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 :param scale: constant scaling factor
80 :type scale: positive ``float``
81
82 :note: The weighting factors are rescaled such that
83 *sum_s integrate( ( weight_i[s] *data_i[s]) **2 )=1*
84 """
85 super(ForwardModelWithPotential, self).__init__()
86 self.__domain = domain
87 if scale > 0:
88 self.__scale = scale
89 else:
90 raise ValueError("Scaling factor must be positive.")
91
92
93 if useSphericalCoordinates:
94 raise ValueError("Spherical coordinates are not supported yet.")
95 else:
96 self.__useSphericalCoordinates=useSphericalCoordinates
97 try:
98 n=len(w)
99 m=len(data)
100 if m != n:
101 raise ValueError("Length of weight and data must be the same.")
102 self.__weight = w
103 self.__data = data
104 except TypeError:
105 self.__weight = [w]
106 self.__data = [data]
107
108 A=0
109 for s in xrange(len(self.__weight)):
110 A += integrate( inner(self.__weight[s], self.__data[s]) **2 )
111 if A > 0:
112 A=vol(domain)/A
113 for s in xrange(len(self.__weight)): self.__weight[s]*=A
114 else:
115 raise ValueError("No non-zero data contribution found.")
116
117
118 BX = boundingBox(domain)
119 DIM = domain.getDim()
120 x = domain.getX()
121 self.__pde=LinearSinglePDE(domain)
122 self.__pde.getSolverOptions().setTolerance(tol)
123 self.__pde.setSymmetryOn()
124 self.__pde.setValue(q=whereZero(x[DIM-1]-BX[DIM-1][1]))
125
126 def useSphericalCoordinates(self):
127 """
128 Returns ``True`` if spherical coordinates are used.
129 """
130 return self.__useSphericalCoordinates
131
132 def getScalingFactor(self):
133 """
134 Returns the scaling factor.
135
136 :rtype: ``float``
137 """
138 return self.__scale
139
140 def getDomain(self):
141 """
142 Returns the domain of the forward model.
143
144 :rtype: `esys.escript.Domain`
145 """
146 return self.__domain
147
148 def getPDE(self):
149 """
150 Return the underlying PDE.
151
152 :rtype: `LinearPDE`
153 """
154 return self.__pde
155
156 def getDefect(self, result):
157 """
158 Returns the defect value.
159
160 :param result: a result vector
161 :type result: `Vector`
162 :rtype: ``float``
163 """
164 A=0.
165 for s in xrange(len(self.__weight)):
166 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
167 return self.__scale*A/2
168
169 def getDefectGradient(self, result):
170 Y=0.
171 for s in range(len(self.__weight)):
172 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
173 print "self.__scale*Y =",self.__scale*Y
174 return self.__scale*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 """
196 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
197 useSphericalCoordinates=False, scale=1., tol=1e-8):
198 """
199 Creates a new gravity model on the given domain with one or more
200 surveys (w, g).
201
202 :param domain: domain of the model
203 :type domain: `esys.escript.Domain`
204 :param w: data weighting factors
205 :type w: `Vector` or list of `Vector`
206 :param g: gravity anomaly data
207 :type g: `Vector` or list of `Vector`
208 :param useSphericalCoordinates: if set spherical coordinates are used.
209 :type useSphericalCoordinates: ``bool``
210 :param tol: tolerance of underlying PDE
211 :type tol: positive ``float``
212 :param scale: constant scaling factor
213 :type scale: positive ``float``
214
215 :note: The weighting factors are rescaled such that
216 *sum_s integrate( ( w_i[s] *g_i[s]) **2 )=1*
217 """
218 super(GravityModel, self).__init__(domain, w, g, useSphericalCoordinates, scale, tol)
219
220 self.__G = gravity_constant
221 self.getPDE().setValue(A=kronecker(self.getDomain()))
222
223 def getArguments(self, rho):
224 """
225 Returns precomputed values shared by getValue() and getGradient().
226
227 :param rho: a suggestion for the density distribution
228 :type rho: `Scalar`
229 :return: gravity potential and corresponding gravity field.
230 :rtype: ``Scalar``, ``Vector``
231 """
232 phi = self.getPotential(rho)
233 gravity_force = -grad(phi)
234 return phi, gravity_force
235
236 def getPotential(self, rho):
237 """
238 Calculates the gravity potential for a given density distribution.
239
240 :param rho: a suggestion for the density distribution
241 :type rho: `Scalar`
242 :return: gravity potential
243 :rtype: ``Scalar``
244 """
245 pde=self.getPDE()
246
247 pde.resetRightHandSideCoefficients()
248 pde.setValue(Y=-4.*PI*self.__G*rho)
249 phi=pde.getSolution()
250
251 return phi
252
253 def getValue(self, rho, phi, gravity_force):
254 """
255 Returns the value of the defect
256
257 :param rho: density distribution
258 :type rho: `Scalar`
259 :param phi: corresponding potential
260 :type phi: `Scalar`
261 :param gravity_force: gravity force
262 :type gravity_force: `Vector`
263 :rtype: ``float``
264 """
265 return self.getDefect(gravity_force)
266
267 def getGradient(self, rho, phi, gravity_force):
268 """
269 Returns the gradient of the defect with respect to density.
270
271 :param rho: density distribution
272 :type rho: `Scalar`
273 :param phi: corresponding potential
274 :type phi: `Scalar`
275 :param gravity_force: gravity force
276 :type gravity_force: `Vector`
277 :rtype: `Scalar`
278 """
279 pde=self.getPDE()
280 pde.resetRightHandSideCoefficients()
281 pde.setValue(X=self.getDefectGradient(gravity_force))
282 ZT=pde.getSolution()
283 return ZT*(-4*PI*self.__G)
284
285
286 class MagneticModel(ForwardModelWithPotential):
287 """
288 Forward Model for magnetic inversion as described in the inversion
289 cookbook.
290 """
291 def __init__(self, domain, w, B, background_field, useSphericalCoordinates=False, scale=1., tol=1e-8):
292 """
293 Creates a new magnetic model on the given domain with one or more
294 surveys (w, B).
295
296 :param domain: domain of the model
297 :type domain: `Domain`
298 :param w: data weighting factors
299 :type w: `Vector` or list of `Vector`
300 :param B: magnetic field data
301 :type B: `Vector` or list of `Vector`
302 :param tol: tolerance of underlying PDE
303 :type tol: positive ``float``
304 :param useSphericalCoordinates: if set spherical coordinates are used.
305 :type useSphericalCoordinates: ``bool``
306 :param scale: constant scaling factor
307 :type scale: positive ``float``
308
309 :note: The weighting factors are rescaled such that
310 *sum_s integrate( ( w_i[s] *B_i[s]) **2 )=1*
311 """
312 super(MagneticModel, self).__init__(domain, w, B, useSphericalCoordinates, scale, tol)
313 self.__background_field=interpolate(background_field, B[0].getFunctionSpace())
314 self.getPDE().setValue(A=kronecker(self.getDomain()))
315
316 def getArguments(self, k):
317 """
318 Returns precomputed values shared by getValue() and getGradient().
319
320 :param k: susceptibility
321 :type k: `Scalar`
322 :return: scalar magnetic potential and corresponding magnetic field
323 :rtype: ``Scalar``, ``Vector``
324 """
325 phi = self.getPotential(k)
326 magnetic_field = k * self.__background_field -grad(phi)
327 return phi, magnetic_field
328
329 def getPotential(self, k):
330 """
331 Calculates the magnetic potential for a given susceptibility.
332
333 :param k: susceptibility
334 :type k: `Scalar`
335 :return: magnetic potential
336 :rtype: ``Scalar``
337 """
338 pde=self.getPDE()
339
340 pde.resetRightHandSideCoefficients()
341 pde.setValue(X = k* self.__background_field)
342 phi=pde.getSolution()
343
344 return phi
345
346 def getValue(self, k, phi, magnetic_field):
347 """
348 Returns the value of the defect.
349
350 :param k: susceptibility
351 :type k: `Scalar`
352 :param phi: corresponding potential
353 :type phi: `Scalar`
354 :param magnetic_field: magnetic field
355 :type magnetic_field: `Vector`
356 :rtype: ``float``
357 """
358 return self.getDefect(magnetic_field)
359
360 def getGradient(self, k, phi, magnetic_field):
361 """
362 Returns the gradient of the defect with respect to susceptibility.
363
364 :param k: susceptibility
365 :type k: `Scalar`
366 :param phi: corresponding potential
367 :type phi: `Scalar`
368 :param magnetic_field: magnetic field
369 :type magnetic_field: `Vector`
370 :rtype: `Scalar`
371 """
372 Y=self.getDefectGradient(magnetic_field)
373 pde=self.getPDE()
374 pde.resetRightHandSideCoefficients()
375 pde.setValue(X=Y)
376
377 #print "test VTK file generated"
378 #from esys.weipa import saveVTK
379 #saveVTK("Y.vtu",RHS=pde.getRightHandSide(), Y=Y)
380
381 YT=pde.getSolution()
382 return inner(grad(YT)-Y,self.__background_field)
383

  ViewVC Help
Powered by ViewVC 1.1.26