/[escript]/trunk/downunder/py_src/forwardmodels.py
ViewVC logotype

Annotation of /trunk/downunder/py_src/forwardmodels.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4080 - (hide annotations)
Mon Nov 19 01:45:38 2012 UTC (6 years, 11 months ago) by jfenwick
File MIME type: text/x-python
File size: 10925 byte(s)
Epydoc fixes

1 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4     # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16 caltinay 4007 """Collection of forward models that define the inversion problem"""
17    
18 caltinay 3946 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 jfenwick 3981 http://www.uq.edu.au
20 caltinay 3946 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 caltinay 4058 __all__ = ['ForwardModel','ForwardModelWithPotential','GravityModel','MagneticModel']
26 caltinay 3947
27 caltinay 3946 from esys.escript import unitsSI as U
28 gross 4037 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE
29 caltinay 4007 from esys.escript.util import *
30 gross 4037 from esys.escript import Data, Vector, Scalar, Function
31     from math import pi as PI
32 caltinay 3946
33 gross 4074
34 caltinay 3946 class ForwardModel(object):
35     """
36     An abstract forward model that can be plugged into a cost function.
37 caltinay 4007 Subclasses need to implement `getValue()`, `getGradient()`, and possibly
38     `getArguments()`.
39 caltinay 3946 """
40     def __init__(self):
41     pass
42    
43     def getArguments(self, x):
44 caltinay 4044 return ()
45 caltinay 3946
46     def getValue(self, x, *args):
47 caltinay 4044 raise NotImplementedError
48 caltinay 3946
49     def getGradient(self, x, *args):
50 caltinay 4044 raise NotImplementedError
51 gross 4076
52 caltinay 4077
53 caltinay 4044 class ForwardModelWithPotential(ForwardModel):
54 caltinay 3946 """
55 caltinay 4044 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 caltinay 3946 """
65 gross 4074 def __init__(self, domain, weight, data, tol=1e-8):
66 caltinay 3946 """
67 caltinay 4044 initialization.
68    
69 gross 4037 :param domain: domain of the model
70 caltinay 4044 :type domain: `esys.escript.Domain`
71     :param weight: data weighting factors
72     :type weight: `Vector` or list of `Vector`
73 gross 4037 :param data: data
74     :type data: `Vector` or list of `Vector`
75     :param tol: tolerance of underlying PDE
76     :type tol: positive ``float``
77 caltinay 4044
78 caltinay 3946 """
79 gross 4037 super(ForwardModelWithPotential, self).__init__()
80 caltinay 3946 self.__domain = domain
81 caltinay 4044
82 caltinay 3946 try:
83 gross 4037 n=len(weight)
84     m=len(data)
85 caltinay 3946 if m != n:
86 caltinay 4077 raise ValueError("Length of weight and data must be the same.")
87 gross 4037 self.__weight = weight
88     self.__data = data
89 caltinay 3946 except TypeError:
90 gross 4037 self.__weight = [weight]
91     self.__data = [data]
92 caltinay 3946
93     A=0
94 gross 4037 for s in xrange(len(self.__weight)):
95 caltinay 4077 A += integrate(inner(self.__weight[s], self.__data[s]**2))
96 gross 4074 if A > 0:
97 gross 4079 A=vol(domain)/A
98     for s in xrange(len(self.__weight)): self.__weight[s]*=A
99     print "FORWARD SCALE = ",self.__weight[0]
100 caltinay 4077 else:
101 gross 4074 raise ValueError("No non-zero data contribution found.")
102 caltinay 4044
103 caltinay 4077
104 caltinay 3946 BX = boundingBox(domain)
105     DIM = domain.getDim()
106     x = domain.getX()
107     self.__pde=LinearSinglePDE(domain)
108     self.__pde.getSolverOptions().setTolerance(tol)
109 caltinay 4058 self.__pde.setSymmetryOn()
110 gross 4079 #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 gross 4037 def getDomain(self):
113     """
114 caltinay 4044 Returns the domain of the forward model.
115 gross 4037 """
116     return self.__domain
117 caltinay 4044
118 gross 4037 def getPDE(self):
119     """
120 caltinay 4044 Return the underlying PDE.
121    
122 gross 4037 :rtype: `LinearPDE`
123     """
124     return self.__pde
125 caltinay 4044
126 gross 4037 def getDefect(self, result):
127     """
128 caltinay 4044 Returns the defect value.
129    
130 gross 4037 :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 gross 4079 A += integrate(inner(self.__weight[s], (result-self.__data[s])**2))
137     return A/2
138 caltinay 4044
139 caltinay 3960 def getSurvey(self, index=None):
140     """
141 caltinay 4063 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 caltinay 3960 If index is None, all surveys will be returned in a pair of lists.
144     """
145     if index is None:
146 gross 4037 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 caltinay 4044
151    
152 gross 4074
153 gross 4037 class GravityModel(ForwardModelWithPotential):
154     """
155 caltinay 4044 Forward Model for gravity inversion as described in the inversion
156     cookbook.
157 gross 4037 """
158 gross 4074 def __init__(self, domain, chi, g, gravity_constant=U.Gravitational_Constant, tol=1e-8):
159 gross 4037 """
160     Creates a new gravity model on the given domain with one or more
161 caltinay 4044 surveys (chi, g).
162    
163 gross 4037 :param domain: domain of the model
164 caltinay 4044 :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 gross 4037 :type g: `Vector` or list of `Vector`
169     :param tol: tolerance of underlying PDE
170     :type tol: positive ``float``
171     """
172 gross 4074 super(GravityModel, self).__init__(domain, chi, g, tol)
173 caltinay 4044
174 gross 4037 self.__G = gravity_constant
175     self.getPDE().setValue(A=kronecker(self.getDomain()))
176 caltinay 3960
177 caltinay 3946 def getArguments(self, rho):
178     """
179     Returns precomputed values shared by getValue() and getGradient().
180 caltinay 4044
181 gross 4037 :param rho: a suggestion for the density distribution
182     :type rho: `Scalar`
183     :return: gravity potential and corresponding gravity field.
184     :rtype: ``Scalar``, ``Vector``
185 caltinay 3946 """
186     phi = self.getPotential(rho)
187 gross 4037 gravity_force = -grad(phi)
188     return phi, gravity_force
189 caltinay 3946
190     def getPotential(self, rho):
191 gross 4037 """
192 caltinay 4044 Calculates the gravity potential for a given density distribution.
193    
194 gross 4037 :param rho: a suggestion for the density distribution
195 caltinay 4044 :type rho: `Scalar`
196 gross 4037 :return: gravity potential
197     :rtype: ``Scalar``
198     """
199     pde=self.getPDE()
200 caltinay 4044
201 gross 4037 pde.resetRightHandSideCoefficients()
202     pde.setValue(Y=-4.*PI*self.__G*rho)
203     phi=pde.getSolution()
204 caltinay 4044
205 caltinay 3946 return phi
206 caltinay 4044
207 caltinay 3946 def getValue(self, rho, phi, gravity_force):
208 gross 4037 """
209 caltinay 4044 Returns the value of the defect
210    
211 gross 4037 :param rho: density distribution
212 caltinay 4044 :type rho: `Scalar`
213 gross 4037 :param phi: corresponding potential
214 caltinay 4044 :type phi: `Scalar`
215 gross 4037 :param gravity_force: gravity force
216 caltinay 4044 :type gravity_force: `Vector`
217 gross 4037 :rtype: ``float``
218     """
219     return self.getDefect(gravity_force)
220 caltinay 3946
221     def getGradient(self, rho, phi, gravity_force):
222 gross 4037 """
223 caltinay 4044 Returns the gradient of the defect with respect to density.
224    
225 gross 4037 :param rho: density distribution
226 caltinay 4044 :type rho: `Scalar`
227 gross 4037 :param phi: corresponding potential
228 caltinay 4044 :type phi: `Scalar`
229 gross 4037 :param gravity_force: gravity force
230 caltinay 4044 :type gravity_force: `Vector`
231 gross 4037 :rtype: `Scalar`
232 caltinay 4044 """
233 gross 4037 pde=self.getPDE()
234     g, chi = self.getSurvey()
235 caltinay 4044
236 caltinay 3946 Z=0.
237 gross 4037 for s in xrange(len(chi)):
238     Z = chi[s] * (g[s]-gravity_force) + Z
239 caltinay 3946
240 gross 4037 pde.resetRightHandSideCoefficients()
241     pde.setValue(X=Z)
242     ZT=pde.getSolution()
243 caltinay 3946 return ZT*(-4*PI*self.__G)
244    
245 caltinay 4044
246 gross 4037 class MagneticModel(ForwardModelWithPotential):
247     """
248 caltinay 4044 Forward Model for magnetic inversion as described in the inversion
249     cookbook.
250 gross 4037 """
251 gross 4074 def __init__(self, domain, chi, B, background_field, tol=1e-8):
252 gross 4037 """
253 caltinay 4044 Creates a new magnetic model on the given domain with one or more
254     surveys (chi, B).
255    
256 gross 4037 :param domain: domain of the model
257     :type domain: `Domain`
258 caltinay 4044 :param chi: data weighting factors
259     :type chi: `Vector` or list of `Vector`
260 gross 4037 :param B: magnetic field data
261     :type B: `Vector` or list of `Vector`
262     :param tol: tolerance of underlying PDE
263     :type tol: positive ``float``
264     """
265 gross 4074 super(MagneticModel, self).__init__(domain, chi, B, tol)
266 gross 4079 print "B = ",B[0]
267     print "background_field = ",background_field
268    
269 caltinay 4077 self.__background_field=interpolate(background_field, B[0].getFunctionSpace())
270 gross 4037 self.getPDE().setValue(A=kronecker(self.getDomain()))
271    
272     def getArguments(self, k):
273     """
274     Returns precomputed values shared by getValue() and getGradient().
275 caltinay 4044
276     :param k: susceptibility
277 gross 4037 :type k: `Scalar`
278 caltinay 4044 :return: scalar magnetic potential and corresponding magnetic field
279 gross 4037 :rtype: ``Scalar``, ``Vector``
280     """
281     phi = self.getPotential(k)
282 gross 4074 magnetic_field = k * self.__background_field -grad(phi)
283 gross 4037 return phi, magnetic_field
284    
285     def getPotential(self, k):
286     """
287 caltinay 4044 Calculates the magnetic potential for a given susceptibility.
288    
289 gross 4037 :param k: susceptibility
290 caltinay 4044 :type k: `Scalar`
291     :return: magnetic potential
292 gross 4037 :rtype: ``Scalar``
293     """
294     pde=self.getPDE()
295 caltinay 4044
296 gross 4037 pde.resetRightHandSideCoefficients()
297 gross 4074 pde.setValue(X = k* self.__background_field)
298 gross 4037 phi=pde.getSolution()
299 caltinay 4044
300 gross 4037 return phi
301 caltinay 4044
302 gross 4037 def getValue(self, k, phi, magnetic_field):
303     """
304 caltinay 4044 Returns the value of the defect.
305    
306 gross 4037 :param k: susceptibility
307 caltinay 4044 :type k: `Scalar`
308 gross 4037 :param phi: corresponding potential
309 caltinay 4044 :type phi: `Scalar`
310 gross 4037 :param magnetic_field: magnetic field
311 caltinay 4044 :type magnetic_field: `Vector`
312 gross 4037 :rtype: ``float``
313     """
314     return self.getDefect(magnetic_field)
315    
316     def getGradient(self, k, phi, magnetic_field):
317     """
318 caltinay 4044 Returns the gradient of the defect with respect to susceptibility.
319    
320 gross 4037 :param k: susceptibility
321 caltinay 4044 :type k: `Scalar`
322 gross 4037 :param phi: corresponding potential
323 caltinay 4044 :type phi: `Scalar`
324 gross 4037 :param magnetic_field: magnetic field
325 caltinay 4044 :type magnetic_field: `Vector`
326 gross 4037 :rtype: `Scalar`
327 caltinay 4044 """
328 gross 4037 pde=self.getPDE()
329     B, chi = self.getSurvey()
330 caltinay 4044
331 gross 4037 Y=0.
332     for s in xrange(len(chi)):
333     Y = chi[s] * (magnetic_field-B[s]) + Y
334 caltinay 4044
335 gross 4037 pde.resetRightHandSideCoefficients()
336     pde.setValue(X=Y)
337 caltinay 4077
338     #print "test VTK file generated"
339     #from esys.weipa import saveVTK
340     #saveVTK("Y.vtu",RHS=pde.getRightHandSide(), Y=Y)
341    
342 gross 4037 YT=pde.getSolution()
343 gross 4079 print "YT = ",YT
344 gross 4039 return inner(Y-grad(YT),self.__background_field)
345 gross 4037

  ViewVC Help
Powered by ViewVC 1.1.26