/[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 4074 - (show annotations)
Thu Nov 15 03:30:59 2012 UTC (6 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 11215 byte(s)
this will break the downundermodule: major revision of the classes.
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.linearPDEs import LinearSinglePDE, LinearPDE
29 from esys.escript.util import *
30 from esys.escript import Data, Vector, Scalar, Function
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 def getDirectionalDerivative(self, x, d, *args):
52 """
53 return directional derivatiave
54 """
55 Y=self.getGradient(x, *args)
56 return integrate(p * Y)
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, weight, data, tol=1e-8):
71 """
72 initialization.
73
74 :param domain: domain of the model
75 :type domain: `esys.escript.Domain`
76 :param weight: data weighting factors
77 :type weight: `Vector` or list of `Vector`
78 :param data: data
79 :type data: `Vector` or list of `Vector`
80 :param tol: tolerance of underlying PDE
81 :type tol: positive ``float``
82
83 """
84 super(ForwardModelWithPotential, self).__init__()
85 self.__domain = domain
86
87 try:
88 n=len(weight)
89 m=len(data)
90 if m != n:
91 raise ValueError("Length of weight and g must be the same.")
92 self.__weight = weight
93 self.__data = data
94 except TypeError:
95 self.__weight = [weight]
96 self.__data = [data]
97
98 A=0
99 for s in xrange(len(self.__weight)):
100 A += integrate(inner(self.__weight[s], self.__data[s]**2))
101 if A > 0:
102 for s in xrange(len(self.__weight)): self.__weight[s]*=1./A
103 else:
104 raise ValueError("No non-zero data contribution found.")
105
106
107 BX = boundingBox(domain)
108 DIM = domain.getDim()
109 x = domain.getX()
110 self.__pde=LinearSinglePDE(domain)
111 self.__pde.getSolverOptions().setTolerance(tol)
112 self.__pde.setSymmetryOn()
113 self.__pde.setValue(q=whereZero(x[DIM-1]-BX[DIM-1][1]))
114
115 def getDomain(self):
116 """
117 Returns the domain of the forward model.
118 """
119 return self.__domain
120
121 def getPDE(self):
122 """
123 Return the underlying PDE.
124
125 :rtype: `LinearPDE`
126 """
127 return self.__pde
128
129 def getDefect(self, result):
130 """
131 Returns the defect value.
132
133 :param result: a result vector
134 :type result: `Vector`
135 :rtype: ``float``
136 """
137 A=0.
138 for s in xrange(len(self.__weight)):
139 A = inner(self.__weight[s], (result-self.__data[s])**2) + A
140 return 0.5*integrate(A)
141
142 def getSurvey(self, index=None):
143 """
144 Returns the pair (data_index, weight_index), where data_i is the data
145 of survey i, weight_i is the weighting factor for survey i.
146 If index is None, all surveys will be returned in a pair of lists.
147 """
148 if index is None:
149 return self.__data, self.__weight
150 if index>=len(self.__data):
151 raise IndexError("Forward model only has %d surveys"%len(self.__data))
152 return self.__data[index], self.__weight[index]
153
154
155
156 class GravityModel(ForwardModelWithPotential):
157 """
158 Forward Model for gravity inversion as described in the inversion
159 cookbook.
160 """
161 def __init__(self, domain, chi, g, gravity_constant=U.Gravitational_Constant, tol=1e-8):
162 """
163 Creates a new gravity model on the given domain with one or more
164 surveys (chi, g).
165
166 :param domain: domain of the model
167 :type domain: `esys.escript.Domain`
168 :param chi: data weighting factors
169 :type chi: `Vector` or list of `Vector`
170 :param g: gravity anomaly data
171 :type g: `Vector` or list of `Vector`
172 :param fix_all_faces: if ``true`` all faces of the domain are fixed
173 otherwise only the top surface
174 :type fix_all_faces: ``bool``
175 :param tol: tolerance of underlying PDE
176 :type tol: positive ``float``
177 """
178 super(GravityModel, self).__init__(domain, chi, g, tol)
179
180 self.__G = gravity_constant
181 self.getPDE().setValue(A=kronecker(self.getDomain()))
182
183 def getArguments(self, rho):
184 """
185 Returns precomputed values shared by getValue() and getGradient().
186
187 :param rho: a suggestion for the density distribution
188 :type rho: `Scalar`
189 :return: gravity potential and corresponding gravity field.
190 :rtype: ``Scalar``, ``Vector``
191 """
192 phi = self.getPotential(rho)
193 gravity_force = -grad(phi)
194 return phi, gravity_force
195
196 def getPotential(self, rho):
197 """
198 Calculates the gravity potential for a given density distribution.
199
200 :param rho: a suggestion for the density distribution
201 :type rho: `Scalar`
202 :return: gravity potential
203 :rtype: ``Scalar``
204 """
205 pde=self.getPDE()
206
207 pde.resetRightHandSideCoefficients()
208 pde.setValue(Y=-4.*PI*self.__G*rho)
209 phi=pde.getSolution()
210
211 return phi
212
213 def getValue(self, rho, phi, gravity_force):
214 """
215 Returns the value of the defect
216
217 :param rho: density distribution
218 :type rho: `Scalar`
219 :param phi: corresponding potential
220 :type phi: `Scalar`
221 :param gravity_force: gravity force
222 :type gravity_force: `Vector`
223 :rtype: ``float``
224 """
225 return self.getDefect(gravity_force)
226
227 def getGradient(self, rho, phi, gravity_force):
228 """
229 Returns the gradient of the defect with respect to density.
230
231 :param rho: density distribution
232 :type rho: `Scalar`
233 :param phi: corresponding potential
234 :type phi: `Scalar`
235 :param gravity_force: gravity force
236 :type gravity_force: `Vector`
237 :rtype: `Scalar`
238 """
239 pde=self.getPDE()
240 g, chi = self.getSurvey()
241
242 Z=0.
243 for s in xrange(len(chi)):
244 Z = chi[s] * (g[s]-gravity_force) + Z
245
246 pde.resetRightHandSideCoefficients()
247 pde.setValue(X=Z)
248 ZT=pde.getSolution()
249 return ZT*(-4*PI*self.__G)
250
251
252 class MagneticModel(ForwardModelWithPotential):
253 """
254 Forward Model for magnetic inversion as described in the inversion
255 cookbook.
256 """
257 def __init__(self, domain, chi, B, background_field, tol=1e-8):
258 """
259 Creates a new magnetic model on the given domain with one or more
260 surveys (chi, B).
261
262 :param domain: domain of the model
263 :type domain: `Domain`
264 :param chi: data weighting factors
265 :type chi: `Vector` or list of `Vector`
266 :param B: magnetic field data
267 :type B: `Vector` or list of `Vector`
268 :param fix_all_faces: if ``true`` all faces of the domain are fixed
269 otherwise only the top surface
270 :type fix_all_faces: ``bool``
271 :param tol: tolerance of underlying PDE
272 :type tol: positive ``float``
273 """
274 super(MagneticModel, self).__init__(domain, chi, B, tol)
275 self.__background_field=interpolate(background_field, Function(self.getDomain()))
276 self.getPDE().setValue(A=kronecker(self.getDomain()))
277
278 def getArguments(self, k):
279 """
280 Returns precomputed values shared by getValue() and getGradient().
281
282 :param k: susceptibility
283 :type k: `Scalar`
284 :return: scalar magnetic potential and corresponding magnetic field
285 :rtype: ``Scalar``, ``Vector``
286 """
287 phi = self.getPotential(k)
288 magnetic_field = k * self.__background_field -grad(phi)
289 return phi, magnetic_field
290
291 def getPotential(self, k):
292 """
293 Calculates the magnetic potential for a given susceptibility.
294
295 :param k: susceptibility
296 :type k: `Scalar`
297 :return: magnetic potential
298 :rtype: ``Scalar``
299 """
300 pde=self.getPDE()
301
302 pde.resetRightHandSideCoefficients()
303 pde.setValue(X = k* self.__background_field)
304 phi=pde.getSolution()
305
306 return phi
307
308 def getValue(self, k, phi, magnetic_field):
309 """
310 Returns the value of the defect.
311
312 :param k: susceptibility
313 :type k: `Scalar`
314 :param phi: corresponding potential
315 :type phi: `Scalar`
316 :param magnetic_field: magnetic field
317 :type magnetic_field: `Vector`
318 :rtype: ``float``
319 """
320 return self.getDefect(magnetic_field)
321
322 def getGradient(self, k, phi, magnetic_field):
323 """
324 Returns the gradient of the defect with respect to susceptibility.
325
326 :param k: susceptibility
327 :type k: `Scalar`
328 :param phi: corresponding potential
329 :type phi: `Scalar`
330 :param magnetic_field: magnetic field
331 :type magnetic_field: `Vector`
332 :rtype: `Scalar`
333 """
334 pde=self.getPDE()
335 B, chi = self.getSurvey()
336
337 Y=0.
338 for s in xrange(len(chi)):
339 Y = chi[s] * (magnetic_field-B[s]) + Y
340
341 pde.resetRightHandSideCoefficients()
342 pde.setValue(X=Y)
343
344 print "test VTK file generated"
345 from esys.weipa import saveVTK
346 saveVTK("Y.vtu",RHS=pde.getRightHandSide(), Y=Y)
347
348 YT=pde.getSolution()
349 return inner(Y-grad(YT),self.__background_field)
350
351

  ViewVC Help
Powered by ViewVC 1.1.26