/[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 4079 - (show annotations)
Fri Nov 16 07:59:01 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 11275 byte(s)
some modifications to scaling in downunder. still not perfect.
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
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, weight, data, tol=1e-8):
66 """
67 initialization.
68
69 :param domain: domain of the model
70 :type domain: `esys.escript.Domain`
71 :param weight: data weighting factors
72 :type weight: `Vector` or list of `Vector`
73 :param data: data
74 :type data: `Vector` or list of `Vector`
75 :param tol: tolerance of underlying PDE
76 :type tol: positive ``float``
77
78 """
79 super(ForwardModelWithPotential, self).__init__()
80 self.__domain = domain
81
82 try:
83 n=len(weight)
84 m=len(data)
85 if m != n:
86 raise ValueError("Length of weight and data must be the same.")
87 self.__weight = weight
88 self.__data = data
89 except TypeError:
90 self.__weight = [weight]
91 self.__data = [data]
92
93 A=0
94 for s in xrange(len(self.__weight)):
95 A += integrate(inner(self.__weight[s], self.__data[s]**2))
96 if A > 0:
97 A=vol(domain)/A
98 for s in xrange(len(self.__weight)): self.__weight[s]*=A
99 print "FORWARD SCALE = ",self.__weight[0]
100 else:
101 raise ValueError("No non-zero data contribution found.")
102
103
104 BX = boundingBox(domain)
105 DIM = domain.getDim()
106 x = domain.getX()
107 self.__pde=LinearSinglePDE(domain)
108 self.__pde.getSolverOptions().setTolerance(tol)
109 self.__pde.setSymmetryOn()
110 #self.__pde.setValue(q=whereZero(x[DIM-1]-BX[DIM-1][1]))
111 self.__pde.setValue(q=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0]))
112 def getDomain(self):
113 """
114 Returns the domain of the forward model.
115 """
116 return self.__domain
117
118 def getPDE(self):
119 """
120 Return the underlying PDE.
121
122 :rtype: `LinearPDE`
123 """
124 return self.__pde
125
126 def getDefect(self, result):
127 """
128 Returns the defect value.
129
130 :param result: a result vector
131 :type result: `Vector`
132 :rtype: ``float``
133 """
134 A=0.
135 for s in xrange(len(self.__weight)):
136 A += integrate(inner(self.__weight[s], (result-self.__data[s])**2))
137 return A/2
138
139 def getSurvey(self, index=None):
140 """
141 Returns the pair (data_index, weight_index), where data_i is the data
142 of survey i, weight_i is the weighting factor for survey i.
143 If index is None, all surveys will be returned in a pair of lists.
144 """
145 if index is None:
146 return self.__data, self.__weight
147 if index>=len(self.__data):
148 raise IndexError("Forward model only has %d surveys"%len(self.__data))
149 return self.__data[index], self.__weight[index]
150
151
152
153 class GravityModel(ForwardModelWithPotential):
154 """
155 Forward Model for gravity inversion as described in the inversion
156 cookbook.
157 """
158 def __init__(self, domain, chi, g, gravity_constant=U.Gravitational_Constant, tol=1e-8):
159 """
160 Creates a new gravity model on the given domain with one or more
161 surveys (chi, g).
162
163 :param domain: domain of the model
164 :type domain: `esys.escript.Domain`
165 :param chi: data weighting factors
166 :type chi: `Vector` or list of `Vector`
167 :param g: gravity anomaly data
168 :type g: `Vector` or list of `Vector`
169 :param fix_all_faces: if ``true`` all faces of the domain are fixed
170 otherwise only the top surface
171 :type fix_all_faces: ``bool``
172 :param tol: tolerance of underlying PDE
173 :type tol: positive ``float``
174 """
175 super(GravityModel, self).__init__(domain, chi, g, tol)
176
177 self.__G = gravity_constant
178 self.getPDE().setValue(A=kronecker(self.getDomain()))
179
180 def getArguments(self, rho):
181 """
182 Returns precomputed values shared by getValue() and getGradient().
183
184 :param rho: a suggestion for the density distribution
185 :type rho: `Scalar`
186 :return: gravity potential and corresponding gravity field.
187 :rtype: ``Scalar``, ``Vector``
188 """
189 phi = self.getPotential(rho)
190 gravity_force = -grad(phi)
191 return phi, gravity_force
192
193 def getPotential(self, rho):
194 """
195 Calculates the gravity potential for a given density distribution.
196
197 :param rho: a suggestion for the density distribution
198 :type rho: `Scalar`
199 :return: gravity potential
200 :rtype: ``Scalar``
201 """
202 pde=self.getPDE()
203
204 pde.resetRightHandSideCoefficients()
205 pde.setValue(Y=-4.*PI*self.__G*rho)
206 phi=pde.getSolution()
207
208 return phi
209
210 def getValue(self, rho, phi, gravity_force):
211 """
212 Returns the value of the defect
213
214 :param rho: density distribution
215 :type rho: `Scalar`
216 :param phi: corresponding potential
217 :type phi: `Scalar`
218 :param gravity_force: gravity force
219 :type gravity_force: `Vector`
220 :rtype: ``float``
221 """
222 return self.getDefect(gravity_force)
223
224 def getGradient(self, rho, phi, gravity_force):
225 """
226 Returns the gradient of the defect with respect to density.
227
228 :param rho: density distribution
229 :type rho: `Scalar`
230 :param phi: corresponding potential
231 :type phi: `Scalar`
232 :param gravity_force: gravity force
233 :type gravity_force: `Vector`
234 :rtype: `Scalar`
235 """
236 pde=self.getPDE()
237 g, chi = self.getSurvey()
238
239 Z=0.
240 for s in xrange(len(chi)):
241 Z = chi[s] * (g[s]-gravity_force) + Z
242
243 pde.resetRightHandSideCoefficients()
244 pde.setValue(X=Z)
245 ZT=pde.getSolution()
246 return ZT*(-4*PI*self.__G)
247
248
249 class MagneticModel(ForwardModelWithPotential):
250 """
251 Forward Model for magnetic inversion as described in the inversion
252 cookbook.
253 """
254 def __init__(self, domain, chi, B, background_field, tol=1e-8):
255 """
256 Creates a new magnetic model on the given domain with one or more
257 surveys (chi, B).
258
259 :param domain: domain of the model
260 :type domain: `Domain`
261 :param chi: data weighting factors
262 :type chi: `Vector` or list of `Vector`
263 :param B: magnetic field data
264 :type B: `Vector` or list of `Vector`
265 :param fix_all_faces: if ``true`` all faces of the domain are fixed
266 otherwise only the top surface
267 :type fix_all_faces: ``bool``
268 :param tol: tolerance of underlying PDE
269 :type tol: positive ``float``
270 """
271 super(MagneticModel, self).__init__(domain, chi, B, tol)
272 print "B = ",B[0]
273 print "background_field = ",background_field
274
275 self.__background_field=interpolate(background_field, B[0].getFunctionSpace())
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 print "YT = ",YT
350 return inner(Y-grad(YT),self.__background_field)
351

  ViewVC Help
Powered by ViewVC 1.1.26