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

Contents of /trunk/downunder/py_src/forwardmodels/magnetic.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 6496 - (show annotations)
Tue Feb 14 00:38:08 2017 UTC (2 years, 2 months ago) by gross
File MIME type: text/x-python
File size: 17037 byte(s)
bug to obtain function space fixed.
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2016 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Apache License, version 2.0
8 # http://www.apache.org/licenses/LICENSE-2.0
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15 from __future__ import division, print_function
16
17 """Forward models for magnetic fields"""
18
19 __copyright__="""Copyright (c) 2003-2016 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Apache License, version 2.0
23 http://www.apache.org/licenses/LICENSE-2.0"""
24 __url__="https://launchpad.net/escript-finley"
25
26 __all__ = ['MagneticModel', 'SelfDemagnetizationModel', 'MagneticIntensityModel']
27
28 from .base import ForwardModelWithPotential
29 from esys.escript import Scalar
30 from esys.escript.util import *
31
32
33 class MagneticModel(ForwardModelWithPotential):
34 """
35 Forward Model for magnetic inversion as described in the inversion
36 cookbook.
37 """
38 def __init__(self, domain, w, B, background_magnetic_flux_density,
39 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
40 """
41 Creates a new magnetic model on the given domain with one or more
42 surveys (w, B).
43
44 :param domain: domain of the model
45 :type domain: `Domain`
46 :param w: data weighting factors
47 :type w: ``Vector`` or list of ``Vector``
48 :param B: magnetic field data
49 :type B: ``Vector`` or list of ``Vector``
50 :param tol: tolerance of underlying PDE
51 :type tol: positive ``float``
52 :param background_magnetic_flux_density: background magnetic flux
53 density (in Tesla) with components (B_east, B_north, B_vertical)
54 :type background_magnetic_flux_density: ``Vector`` or list of `float`
55 :param coordinates: defines coordinate system to be used
56 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
57 :param fixPotentialAtBottom: if true potential is fixed to zero at the
58 bottom of the domain in addition to the top
59 :type fixPotentialAtBottom: ``bool``
60 """
61 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
62 background_magnetic_flux_density=interpolate(background_magnetic_flux_density, self.getDataFunctionSpace() )
63 if not self.getCoordinateTransformation().isCartesian():
64 s = self.getCoordinateTransformation().getScalingFactors()
65 v = self.getCoordinateTransformation().getVolumeFactor()
66 self.__B_r = background_magnetic_flux_density * s * v
67 self.__B_b = background_magnetic_flux_density / s
68
69 A = self.getPDE().createCoefficient("A")
70 fw = s**2 * v
71 for i in range(self.getDomain().getDim()):
72 A[i,i]=fw[i]
73 self.getPDE().setValue(A=A)
74 else: # cartesian
75 self.getPDE().setValue(A=kronecker(self.getDomain()))
76 self.__B_r = background_magnetic_flux_density
77 self.__B_b = background_magnetic_flux_density
78
79 def rescaleWeights(self, scale=1., k_scale=1.):
80 """
81 rescales the weights such that
82
83 *sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale*
84
85 :param scale: scale of data weighting factors
86 :type scale: positive ``float``
87 :param k_scale: scale of susceptibility.
88 :type k_scale: ``Scalar``
89 """
90 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
91
92 def getArguments(self, k):
93 """
94 Returns precomputed values shared by `getDefect()` and `getGradient()`.
95
96 :param k: susceptibility
97 :type k: ``Scalar``
98 :return: scalar magnetic potential and corresponding magnetic field
99 :rtype: ``Scalar``, ``Vector``
100 """
101 phi = self.getPotential(k)
102 magnetic_flux_density = k * self.__B_b -grad(phi)
103 return phi, magnetic_flux_density
104
105 def getPotential(self, k):
106 """
107 Calculates the magnetic potential for a given susceptibility.
108
109 :param k: susceptibility
110 :type k: ``Scalar``
111 :return: magnetic potential
112 :rtype: ``Scalar``
113 """
114 pde=self.getPDE()
115 pde.resetRightHandSideCoefficients()
116 pde.setValue(X = k* self.__B_r)
117 phi=pde.getSolution()
118
119 return phi
120
121 def getDefect(self, k, phi, magnetic_flux_density):
122 """
123 Returns the value of the defect.
124
125 :param k: susceptibility
126 :type k: ``Scalar``
127 :param phi: corresponding potential
128 :type phi: ``Scalar``
129 :param magnetic_flux_density: magnetic field
130 :type magnetic_flux_density: ``Vector``
131 :rtype: ``float``
132 """
133 return self._getDefect(magnetic_flux_density)
134
135 def getGradient(self, k, phi, magnetic_flux_density):
136 """
137 Returns the gradient of the defect with respect to susceptibility.
138
139 :param k: susceptibility
140 :type k: ``Scalar``
141 :param phi: corresponding potential
142 :type phi: ``Scalar``
143 :param magnetic_flux_density: magnetic field
144 :type magnetic_flux_density: ``Vector``
145 :rtype: ``Scalar``
146 """
147 Y=self.getDefectGradient(magnetic_flux_density)
148 pde=self.getPDE()
149 pde.resetRightHandSideCoefficients()
150 pde.setValue(X=Y)
151 YT=pde.getSolution()
152 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
153
154 class SelfDemagnetizationModel(ForwardModelWithPotential):
155 """
156 Forward Model for magnetic inversion with self-demagnetization as
157 described in the inversion cookbook.
158 """
159 def __init__(self, domain, w, B, background_magnetic_flux_density,
160 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
161 """
162 Creates a new magnetic model on the given domain with one or more
163 surveys (w, B).
164
165 :param domain: domain of the model
166 :type domain: `Domain`
167 :param w: data weighting factors
168 :type w: ``Vector`` or list of ``Vector``
169 :param B: magnetic field data
170 :type B: ``Vector`` or list of ``Vector``
171 :param background_magnetic_flux_density: background magnetic flux
172 density (in Tesla) with components (B_east, B_north, B_vertical)
173 :type background_magnetic_flux_density: ``Vector`` or list of `float`
174 :param coordinates: defines coordinate system to be used
175 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
176 :param fixPotentialAtBottom: if true potential is fixed to zero at the
177 bottom of the domain in addition to the top
178 :type fixPotentialAtBottom: ``bool``
179 :param tol: tolerance of underlying PDE
180 :type tol: positive ``float``
181 """
182 super(SelfDemagnetizationModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
183 #=========================================================
184 background_magnetic_flux_density = interpolate(background_magnetic_flux_density, self.getDataFunctionSpace())
185 if not self.getCoordinateTransformation().isCartesian():
186 s = self.getCoordinateTransformation().getScalingFactors()
187 v = self.getCoordinateTransformation().getVolumeFactor()
188 self.__B_r = background_magnetic_flux_density * s * v
189 self.__B_b = background_magnetic_flux_density / s
190
191 self.__fw = s**2 * v
192 else: # cartesian
193 self.__fw = 1
194 self.__B_r = background_magnetic_flux_density
195 self.__B_b = background_magnetic_flux_density
196
197 # keep track of k used to build PDE:
198 self.__last_k = None
199 # this is just the initial set_up
200 A=self.getPDE().createCoefficient("A")
201 for i in range(self.getDomain().getDim()):
202 A[i,i]=1.
203 self.getPDE().setValue(A=A)
204
205 def rescaleWeights(self, scale=1., k_scale=1.):
206 """
207 rescales the weights such that
208
209 *sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale*
210
211 :param scale: scale of data weighting factors
212 :type scale: positive ``float``
213 :param k_scale: scale of susceptibility.
214 :type k_scale: ``Scalar``
215 """
216 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
217
218 def getArguments(self, k):
219 """
220 Returns precomputed values shared by `getDefect()` and `getGradient()`.
221
222 :param k: susceptibility
223 :type k: ``Scalar``
224 :return: scalar magnetic potential and corresponding magnetic field
225 :rtype: ``Scalar``, ``Vector``
226 """
227 phi = self.getPotential(k)
228 grad_phi=grad(phi)
229 magnetic_flux_density = k * self.__B_b -(1+k)*grad_phi
230 return phi, grad_phi, magnetic_flux_density
231
232 def __updateA(self, k):
233 """
234 updates PDE coefficient if PDE is used with new k
235 """
236 pde=self.getPDE()
237 if self.__last_k is not k:
238 A=pde.getCoefficient('A')
239 if self.getCoordinateTransformation().isCartesian():
240 for i in range(self.getDomain().getDim()):
241 A[i,i] = 1+k
242 else:
243 for i in range(self.getDomain().getDim()):
244 A[i,i] = (1+k)*self.__fw[i]
245
246 self.__last_k = k
247 pde.setValue(A=A)
248
249 def getPotential(self, k):
250 """
251 Calculates the magnetic potential for a given susceptibility.
252
253 :param k: susceptibility
254 :type k: ``Scalar``
255 :return: magnetic potential
256 :rtype: ``Scalar``
257 """
258 self.__updateA(k)
259 pde=self.getPDE()
260 pde.resetRightHandSideCoefficients()
261 pde.setValue(X = k*self.__B_r)
262 phi=pde.getSolution()
263 return phi
264
265 def getDefect(self, k, phi, grad_phi, magnetic_flux_density):
266 """
267 Returns the value of the defect.
268
269 :param k: susceptibility
270 :type k: ``Scalar``
271 :param phi: corresponding potential
272 :type phi: ``Scalar``
273 :param magnetic_flux_density: magnetic field
274 :type magnetic_flux_density: ``Vector``
275 :rtype: ``float``
276 """
277 return self._getDefect(magnetic_flux_density)
278
279 def getGradient(self, k, phi, grad_phi, magnetic_flux_density):
280 """
281 Returns the gradient of the defect with respect to susceptibility.
282
283 :param k: susceptibility
284 :type k: ``Scalar``
285 :param phi: corresponding potential
286 :type phi: ``Scalar``
287 :param magnetic_flux_density: magnetic field
288 :type magnetic_flux_density: ``Vector``
289 :rtype: ``Scalar``
290 """
291 self.__updateA(k)
292 Y=self.getDefectGradient(magnetic_flux_density)
293 pde=self.getPDE()
294 pde.resetRightHandSideCoefficients()
295 pde.setValue(X=(1+k)*Y)
296 grad_YT=grad(pde.getSolution())
297
298 if self.getCoordinateTransformation().isCartesian(): # then b_r=B_b
299 return inner(grad_YT-Y, self.__B_r-grad_phi)
300 else:
301 return inner(grad_YT,self.__B_r-grad_phi)+inner(Y,grad_phi-self.__B_b)
302
303 class MagneticIntensityModel(ForwardModelWithPotential):
304 """
305 Forward Model for magnetic intensity inversion as described in the inversion
306 cookbook.
307 """
308 def __init__(self, domain, w, b, background_magnetic_flux_density,
309 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
310 """
311 Creates a new magnetic intensity model on the given domain with one or more
312 surveys (w, b).
313
314 :param domain: domain of the model
315 :type domain: `Domain`
316 :param w: data weighting factors
317 :type w: ``Scalar`` or list of ``Scalar``
318 :param b: magnetic intensity field data
319 :type b: ``Scalar`` or list of ``Scalar``
320 :param tol: tolerance of underlying PDE
321 :type tol: positive ``float``
322 :param background_magnetic_flux_density: background magnetic flux
323 density (in Tesla) with components (B_east, B_north, B_vertical)
324 :type background_magnetic_flux_density: ``Vector`` or list of `float`
325 :param coordinates: defines coordinate system to be used
326 :type coordinates: None
327 :param fixPotentialAtBottom: if true potential is fixed to zero at the
328 bottom of the domain in addition to the top
329 :type fixPotentialAtBottom: ``bool``
330 """
331 super(MagneticIntensityModel, self).__init__(domain, w, b, coordinates, fixPotentialAtBottom, tol)
332 background_magnetic_flux_density=interpolate(background_magnetic_flux_density, self.getDataFunctionSpace())
333 if not self.getCoordinateTransformation().isCartesian(): # need to be checked!
334 s = self.getCoordinateTransformation().getScalingFactors()
335 v = self.getCoordinateTransformation().getVolumeFactor()
336 self.__B_r = background_magnetic_flux_density * s * v
337 self.__B_b = background_magnetic_flux_density / s
338
339 A = self.getPDE().createCoefficient("A")
340 fw = s**2 * v
341 for i in range(self.getDomain().getDim()):
342 A[i,i]=fw[i]
343 self.getPDE().setValue(A=A)
344 else: # cartesian
345 self.getPDE().setValue(A=kronecker(self.getDomain()))
346 self.__B_r = background_magnetic_flux_density
347 self.__B_b = background_magnetic_flux_density
348 self.__normalized_B_b=normalize(self.__B_b)
349
350 def rescaleWeights(self, scale=1., k_scale=1.):
351 """
352 rescales the weights such that
353
354 *sum_s integrate( ( w_i[s] *B_i[s]) (w_j[s]*1/L_j) * L**2 * (background_magnetic_flux_density_j[s]*1/L_j) * k_scale )=scale*
355
356 :param scale: scale of data weighting factors
357 :type scale: positive ``float``
358 :param k_scale: scale of susceptibility.
359 :type k_scale: ``Scalar``
360 """
361 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
362
363 def getArguments(self, k):
364 """
365 Returns precomputed values shared by `getDefect()` and `getGradient()`.
366
367 :param k: susceptibility
368 :type k: ``Scalar``
369 :return: scalar magnetic potential and corresponding magnetic field
370 :rtype: ``Scalar``, ``Vector``
371 """
372 phi = self.getPotential(k)
373 magnetic_flux_density = k * self.__B_b -grad(phi)
374 return phi, magnetic_flux_density
375
376 def getPotential(self, k):
377 """
378 Calculates the magnetic potential for a given susceptibility.
379
380 :param k: susceptibility
381 :type k: ``Scalar``
382 :return: magnetic potential
383 :rtype: ``Scalar``
384 """
385 pde=self.getPDE()
386 pde.resetRightHandSideCoefficients()
387 pde.setValue(X = k* self.__B_r)
388 phi=pde.getSolution()
389
390 return phi
391
392 def getDefect(self, k, phi, magnetic_flux_density):
393 """
394 Returns the value of the defect.
395
396 :param k: susceptibility
397 :type k: ``Scalar``
398 :param phi: corresponding potential
399 :type phi: ``Scalar``
400 :param magnetic_flux_density: magnetic field
401 :type magnetic_flux_density: ``Vector``
402 :rtype: ``float``
403 """
404 weights=self.getMisfitWeights()
405 data=self.getData()
406 A=0.
407 for s in xrange(len(weights)):
408 A += integrate( (weights[s]*(inner(self.__normalized_B_b, magnetic_flux_density)-data[s]) )**2 )
409 return A/2
410
411 def getGradient(self, k, phi, magnetic_flux_density):
412 """
413 Returns the gradient of the defect with respect to susceptibility.
414
415 :param k: susceptibility
416 :type k: ``Scalar``
417 :param phi: corresponding potential
418 :type phi: ``Scalar``
419 :param magnetic_flux_density: magnetic field
420 :type magnetic_flux_density: ``Vector``
421 :rtype: ``Scalar``
422 """
423 weights=self.getMisfitWeights()
424 data=self.getData()
425 Y=Scalar(0.,magnetic_flux_density.getFunctionSpace())
426 for s in xrange(len(weights)):
427 Y+=weights[s]**2*(data[s]-inner(self.__normalized_B_b, magnetic_flux_density))
428 pde=self.getPDE()
429 pde.resetRightHandSideCoefficients()
430 pde.setValue(X=Y*self.__normalized_B_b)
431 YT=pde.getSolution()
432 return inner(grad(YT),self.__B_r) -Y*inner(self.__normalized_B_b,self.__B_b)

  ViewVC Help
Powered by ViewVC 1.1.26