/[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 4058 - (hide annotations)
Fri Nov 2 06:47:06 2012 UTC (7 years ago) by caltinay
File MIME type: text/x-python
File size: 11447 byte(s)
-Ripley now checks for input fs when calling gradient
-interactive mode works again when filename is supplied (typo)
-forward model PDE sets symmetry on now

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     class ForwardModel(object):
34     """
35     An abstract forward model that can be plugged into a cost function.
36 caltinay 4007 Subclasses need to implement `getValue()`, `getGradient()`, and possibly
37     `getArguments()`.
38 caltinay 3946 """
39     def __init__(self):
40     pass
41    
42     def getArguments(self, x):
43 caltinay 4044 return ()
44 caltinay 3946
45     def getValue(self, x, *args):
46 caltinay 4044 raise NotImplementedError
47 caltinay 3946
48     def getGradient(self, x, *args):
49 caltinay 4044 raise NotImplementedError
50    
51     class ForwardModelWithPotential(ForwardModel):
52 caltinay 3946 """
53 caltinay 4044 Base class for a forward model using a potential such as magnetic or
54     gravity. It defines a cost function
55    
56     defect = 1/2 sum_s integrate( weight_i[s] * ( r_i - data_i[s] )**2 )
57    
58     where s runs over the survey, weight_i are weighting factors, data_i are
59     the data, and r_i are the results produced by the forward model.
60     It is assumed that the forward model is produced through postprocessing
61     of the solution of a potential PDE.
62 caltinay 3946 """
63 caltinay 4044 def __init__(self, domain, weight, data, fix_all_faces=False, tol=1e-8):
64 caltinay 3946 """
65 caltinay 4044 initialization.
66    
67 gross 4037 :param domain: domain of the model
68 caltinay 4044 :type domain: `esys.escript.Domain`
69     :param weight: data weighting factors
70     :type weight: `Vector` or list of `Vector`
71 gross 4037 :param data: data
72     :type data: `Vector` or list of `Vector`
73 caltinay 4044 :param fix_all_faces: if ``true`` all faces of the domain are fixed
74     otherwise only the top surface
75 gross 4037 :type fix_all_faces: ``bool``
76     :param tol: tolerance of underlying PDE
77     :type tol: positive ``float``
78 caltinay 4044
79 caltinay 3946 """
80 gross 4037 super(ForwardModelWithPotential, self).__init__()
81 caltinay 3946 self.__domain = domain
82 caltinay 4044
83 caltinay 3946 try:
84 gross 4037 n=len(weight)
85     m=len(data)
86 caltinay 3946 if m != n:
87 gross 4037 raise ValueError("Length of weight and g must be the same.")
88     self.__weight = weight
89     self.__data = data
90 caltinay 3946 except TypeError:
91 gross 4037 self.__weight = [weight]
92     self.__data = [data]
93 caltinay 3946
94     A=0
95 gross 4037 for s in xrange(len(self.__weight)):
96     A2 = integrate(inner(self.__weight[s], self.__data[s]**2))
97 caltinay 3946 if A2 < 0:
98     raise ValueError("Negative weighting factor for survey %s"%s)
99     A=max(A2, A)
100     if not A > 0:
101     raise ValueError("No reference data set.")
102 caltinay 4044
103 caltinay 3946 BX = boundingBox(domain)
104     DIM = domain.getDim()
105     x = domain.getX()
106     self.__pde=LinearSinglePDE(domain)
107     self.__pde.getSolverOptions().setTolerance(tol)
108 caltinay 4058 self.__pde.setSymmetryOn()
109 caltinay 4044
110 gross 4037 if fix_all_faces:
111 caltinay 4044 constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])
112     for i in xrange(DIM-1):
113     constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])
114     else:
115     constraint=whereZero(x[DIM-1]-BX[DIM-1][1])
116     self.__pde.setValue(q=constraint)
117    
118 gross 4037 def getDomain(self):
119     """
120 caltinay 4044 Returns the domain of the forward model.
121 gross 4037 """
122     return self.__domain
123 caltinay 4044
124 gross 4037 def getPDE(self):
125     """
126 caltinay 4044 Return the underlying PDE.
127    
128 gross 4037 :rtype: `LinearPDE`
129     """
130     return self.__pde
131 caltinay 4044
132 gross 4037 def getDefect(self, result):
133     """
134 caltinay 4044 Returns the defect value.
135    
136 gross 4037 :param result: a result vector
137     :type result: `Vector`
138     :rtype: ``float``
139     """
140     A=0.
141     for s in xrange(len(self.__weight)):
142     A = inner(self.__weight[s], (result-self.__data[s])**2) + A
143     return 0.5*integrate(A)
144 caltinay 4044
145 caltinay 3960 def getSurvey(self, index=None):
146     """
147 gross 4037 Returns the pair (g_index, weight_index), where g_i is the gravity
148     anomaly of survey i, weight_i is the weighting factor for survey i.
149 caltinay 3960 If index is None, all surveys will be returned in a pair of lists.
150     """
151     if index is None:
152 gross 4037 return self.__data, self.__weight
153     if index>=len(self.__data):
154     raise IndexError("Forward model only has %d surveys"%len(self.__data))
155     return self.__data[index], self.__weight[index]
156 caltinay 4044
157    
158 gross 4037 class GravityModel(ForwardModelWithPotential):
159     """
160 caltinay 4044 Forward Model for gravity inversion as described in the inversion
161     cookbook.
162 gross 4037 """
163 caltinay 4053 def __init__(self, domain, chi, g, fix_all_faces=False, gravity_constant=U.Gravitational_Constant, tol=1e-8):
164 gross 4037 """
165     Creates a new gravity model on the given domain with one or more
166 caltinay 4044 surveys (chi, g).
167    
168 gross 4037 :param domain: domain of the model
169 caltinay 4044 :type domain: `esys.escript.Domain`
170     :param chi: data weighting factors
171     :type chi: `Vector` or list of `Vector`
172     :param g: gravity anomaly data
173 gross 4037 :type g: `Vector` or list of `Vector`
174 caltinay 4044 :param fix_all_faces: if ``true`` all faces of the domain are fixed
175     otherwise only the top surface
176 gross 4037 :type fix_all_faces: ``bool``
177     :param tol: tolerance of underlying PDE
178     :type tol: positive ``float``
179     """
180     super(GravityModel, self).__init__(domain, chi, g, fix_all_faces, tol)
181 caltinay 4044
182 gross 4037 self.__G = gravity_constant
183     self.getPDE().setValue(A=kronecker(self.getDomain()))
184 caltinay 3960
185 caltinay 3946 def getArguments(self, rho):
186     """
187     Returns precomputed values shared by getValue() and getGradient().
188 caltinay 4044
189 gross 4037 :param rho: a suggestion for the density distribution
190     :type rho: `Scalar`
191     :return: gravity potential and corresponding gravity field.
192     :rtype: ``Scalar``, ``Vector``
193 caltinay 3946 """
194     phi = self.getPotential(rho)
195 gross 4037 gravity_force = -grad(phi)
196     return phi, gravity_force
197 caltinay 3946
198     def getPotential(self, rho):
199 gross 4037 """
200 caltinay 4044 Calculates the gravity potential for a given density distribution.
201    
202 gross 4037 :param rho: a suggestion for the density distribution
203 caltinay 4044 :type rho: `Scalar`
204 gross 4037 :return: gravity potential
205     :rtype: ``Scalar``
206     """
207     pde=self.getPDE()
208 caltinay 4044
209 gross 4037 pde.resetRightHandSideCoefficients()
210     pde.setValue(Y=-4.*PI*self.__G*rho)
211     phi=pde.getSolution()
212 caltinay 4044
213 caltinay 3946 return phi
214 caltinay 4044
215 caltinay 3946 def getValue(self, rho, phi, gravity_force):
216 gross 4037 """
217 caltinay 4044 Returns the value of the defect
218    
219 gross 4037 :param rho: density distribution
220 caltinay 4044 :type rho: `Scalar`
221 gross 4037 :param phi: corresponding potential
222 caltinay 4044 :type phi: `Scalar`
223 gross 4037 :param gravity_force: gravity force
224 caltinay 4044 :type gravity_force: `Vector`
225 gross 4037 :rtype: ``float``
226     """
227     return self.getDefect(gravity_force)
228 caltinay 3946
229     def getGradient(self, rho, phi, gravity_force):
230 gross 4037 """
231 caltinay 4044 Returns the gradient of the defect with respect to density.
232    
233 gross 4037 :param rho: density distribution
234 caltinay 4044 :type rho: `Scalar`
235 gross 4037 :param phi: corresponding potential
236 caltinay 4044 :type phi: `Scalar`
237 gross 4037 :param gravity_force: gravity force
238 caltinay 4044 :type gravity_force: `Vector`
239 gross 4037 :rtype: `Scalar`
240 caltinay 4044 """
241 gross 4037 pde=self.getPDE()
242     g, chi = self.getSurvey()
243 caltinay 4044
244 caltinay 3946 Z=0.
245 gross 4037 for s in xrange(len(chi)):
246     Z = chi[s] * (g[s]-gravity_force) + Z
247 caltinay 3946
248 gross 4037 pde.resetRightHandSideCoefficients()
249     pde.setValue(X=Z)
250     ZT=pde.getSolution()
251 caltinay 3946 return ZT*(-4*PI*self.__G)
252    
253 caltinay 4044
254 gross 4037 class MagneticModel(ForwardModelWithPotential):
255     """
256 caltinay 4044 Forward Model for magnetic inversion as described in the inversion
257     cookbook.
258 gross 4037 """
259 caltinay 4053 def __init__(self, domain, chi, B, background_field, fix_all_faces=False, tol=1e-8):
260 gross 4037 """
261 caltinay 4044 Creates a new magnetic model on the given domain with one or more
262     surveys (chi, B).
263    
264 gross 4037 :param domain: domain of the model
265     :type domain: `Domain`
266 caltinay 4044 :param chi: data weighting factors
267     :type chi: `Vector` or list of `Vector`
268 gross 4037 :param B: magnetic field data
269     :type B: `Vector` or list of `Vector`
270 caltinay 4044 :param fix_all_faces: if ``true`` all faces of the domain are fixed
271     otherwise only the top surface
272 gross 4037 :type fix_all_faces: ``bool``
273     :param tol: tolerance of underlying PDE
274     :type tol: positive ``float``
275     """
276 gross 4039 super(MagneticModel, self).__init__(domain, chi, B, fix_all_faces, tol)
277 gross 4037 self.__background_field=interpolate(background_field, Function(self.getDomain()))
278     self.getPDE().setValue(A=kronecker(self.getDomain()))
279    
280     def getArguments(self, k):
281     """
282     Returns precomputed values shared by getValue() and getGradient().
283 caltinay 4044
284     :param k: susceptibility
285 gross 4037 :type k: `Scalar`
286 caltinay 4044 :return: scalar magnetic potential and corresponding magnetic field
287 gross 4037 :rtype: ``Scalar``, ``Vector``
288     """
289     phi = self.getPotential(k)
290     magnetic_field = (1+k) * self.__background_field -grad(phi)
291     return phi, magnetic_field
292    
293     def getPotential(self, k):
294     """
295 caltinay 4044 Calculates the magnetic potential for a given susceptibility.
296    
297 gross 4037 :param k: susceptibility
298 caltinay 4044 :type k: `Scalar`
299     :return: magnetic potential
300 gross 4037 :rtype: ``Scalar``
301     """
302     pde=self.getPDE()
303 caltinay 4044
304 gross 4037 pde.resetRightHandSideCoefficients()
305     pde.setValue(X = (1+k)* self.__background_field)
306     phi=pde.getSolution()
307 caltinay 4044
308 gross 4037 return phi
309 caltinay 4044
310 gross 4037 def getValue(self, k, phi, magnetic_field):
311     """
312 caltinay 4044 Returns the value of the defect.
313    
314 gross 4037 :param k: susceptibility
315 caltinay 4044 :type k: `Scalar`
316 gross 4037 :param phi: corresponding potential
317 caltinay 4044 :type phi: `Scalar`
318 gross 4037 :param magnetic_field: magnetic field
319 caltinay 4044 :type magnetic_field: `Vector`
320 gross 4037 :rtype: ``float``
321     """
322     return self.getDefect(magnetic_field)
323    
324     def getGradient(self, k, phi, magnetic_field):
325     """
326 caltinay 4044 Returns the gradient of the defect with respect to susceptibility.
327    
328 gross 4037 :param k: susceptibility
329 caltinay 4044 :type k: `Scalar`
330 gross 4037 :param phi: corresponding potential
331 caltinay 4044 :type phi: `Scalar`
332 gross 4037 :param magnetic_field: magnetic field
333 caltinay 4044 :type magnetic_field: `Vector`
334 gross 4037 :rtype: `Scalar`
335 caltinay 4044 """
336 gross 4037 pde=self.getPDE()
337     B, chi = self.getSurvey()
338 caltinay 4044
339 gross 4037 Y=0.
340     for s in xrange(len(chi)):
341     Y = chi[s] * (magnetic_field-B[s]) + Y
342 caltinay 4044
343 gross 4037 pde.resetRightHandSideCoefficients()
344     pde.setValue(X=Y)
345     YT=pde.getSolution()
346 gross 4039 return inner(Y-grad(YT),self.__background_field)
347 gross 4037

  ViewVC Help
Powered by ViewVC 1.1.26