/[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 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 months ago) by sshaw
File MIME type: text/x-python
File size: 11695 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
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-2015 by The University of Queensland
20 http://www.uq.edu.au
21 Primary Business: Queensland, Australia"""
22 __license__="""Licensed under the Open Software License version 3.0
23 http://www.opensource.org/licenses/osl-3.0.php"""
24 __url__="https://launchpad.net/escript-finley"
25
26 __all__ = ['MagneticModel', 'SelfDemagnetizationModel']
27
28 from .base import ForwardModelWithPotential
29 from esys.escript.util import *
30
31
32 class MagneticModel(ForwardModelWithPotential):
33 """
34 Forward Model for magnetic inversion as described in the inversion
35 cookbook.
36 """
37 def __init__(self, domain, w, B, background_magnetic_flux_density,
38 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
39 """
40 Creates a new magnetic model on the given domain with one or more
41 surveys (w, B).
42
43 :param domain: domain of the model
44 :type domain: `Domain`
45 :param w: data weighting factors
46 :type w: ``Vector`` or list of ``Vector``
47 :param B: magnetic field data
48 :type B: ``Vector`` or list of ``Vector``
49 :param tol: tolerance of underlying PDE
50 :type tol: positive ``float``
51 :param background_magnetic_flux_density: background magnetic flux
52 density (in Tesla) with components (B_east, B_north, B_vertical)
53 :type background_magnetic_flux_density: ``Vector`` or list of `float`
54 :param coordinates: defines coordinate system to be used
55 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
56 :param fixPotentialAtBottom: if true potential is fixed to zero at the
57 bottom of the domain in addition to the top
58 :type fixPotentialAtBottom: ``bool``
59 """
60 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
61 background_magnetic_flux_density=interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
62 if not self.getCoordinateTransformation().isCartesian():
63 s = self.getCoordinateTransformation().getScalingFactors()
64 v = self.getCoordinateTransformation().getVolumeFactor()
65 self.__B_r = background_magnetic_flux_density * s * v
66 self.__B_b = background_magnetic_flux_density / s
67
68 A = self.getPDE().createCoefficient("A")
69 fw = s**2 * v
70 for i in range(self.getDomain().getDim()):
71 A[i,i]=fw[i]
72 self.getPDE().setValue(A=A)
73 else: # cartesian
74 self.getPDE().setValue(A=kronecker(self.getDomain()))
75 self.__B_r = background_magnetic_flux_density
76 self.__B_b = background_magnetic_flux_density
77
78 def rescaleWeights(self, scale=1., k_scale=1.):
79 """
80 rescales the weights such that
81
82 *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*
83
84 :param scale: scale of data weighting factors
85 :type scale: positive ``float``
86 :param k_scale: scale of susceptibility.
87 :type k_scale: ``Scalar``
88 """
89 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
90
91 def getArguments(self, k):
92 """
93 Returns precomputed values shared by `getDefect()` and `getGradient()`.
94
95 :param k: susceptibility
96 :type k: ``Scalar``
97 :return: scalar magnetic potential and corresponding magnetic field
98 :rtype: ``Scalar``, ``Vector``
99 """
100 phi = self.getPotential(k)
101 magnetic_flux_density = k * self.__B_b -grad(phi)
102 return phi, magnetic_flux_density
103
104 def getPotential(self, k):
105 """
106 Calculates the magnetic potential for a given susceptibility.
107
108 :param k: susceptibility
109 :type k: ``Scalar``
110 :return: magnetic potential
111 :rtype: ``Scalar``
112 """
113 pde=self.getPDE()
114 pde.resetRightHandSideCoefficients()
115 pde.setValue(X = k* self.__B_r)
116 phi=pde.getSolution()
117
118 return phi
119
120 def getDefect(self, k, phi, magnetic_flux_density):
121 """
122 Returns the value of the defect.
123
124 :param k: susceptibility
125 :type k: ``Scalar``
126 :param phi: corresponding potential
127 :type phi: ``Scalar``
128 :param magnetic_flux_density: magnetic field
129 :type magnetic_flux_density: ``Vector``
130 :rtype: ``float``
131 """
132 return self._getDefect(magnetic_flux_density)
133
134 def getGradient(self, k, phi, magnetic_flux_density):
135 """
136 Returns the gradient of the defect with respect to susceptibility.
137
138 :param k: susceptibility
139 :type k: ``Scalar``
140 :param phi: corresponding potential
141 :type phi: ``Scalar``
142 :param magnetic_flux_density: magnetic field
143 :type magnetic_flux_density: ``Vector``
144 :rtype: ``Scalar``
145 """
146 Y=self.getDefectGradient(magnetic_flux_density)
147 pde=self.getPDE()
148 pde.resetRightHandSideCoefficients()
149 pde.setValue(X=Y)
150 YT=pde.getSolution()
151 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
152
153 class SelfDemagnetizationModel(ForwardModelWithPotential):
154 """
155 Forward Model for magnetic inversion with self-demagnetization as
156 described in the inversion cookbook.
157 """
158 def __init__(self, domain, w, B, background_magnetic_flux_density,
159 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
160 """
161 Creates a new magnetic model on the given domain with one or more
162 surveys (w, B).
163
164 :param domain: domain of the model
165 :type domain: `Domain`
166 :param w: data weighting factors
167 :type w: ``Vector`` or list of ``Vector``
168 :param B: magnetic field data
169 :type B: ``Vector`` or list of ``Vector``
170 :param background_magnetic_flux_density: background magnetic flux
171 density (in Tesla) with components (B_east, B_north, B_vertical)
172 :type background_magnetic_flux_density: ``Vector`` or list of `float`
173 :param coordinates: defines coordinate system to be used
174 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
175 :param fixPotentialAtBottom: if true potential is fixed to zero at the
176 bottom of the domain in addition to the top
177 :type fixPotentialAtBottom: ``bool``
178 :param tol: tolerance of underlying PDE
179 :type tol: positive ``float``
180 """
181 super(SelfDemagnetizationModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
182 #=========================================================
183 background_magnetic_flux_density = interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
184 if not self.getCoordinateTransformation().isCartesian():
185 s = self.getCoordinateTransformation().getScalingFactors()
186 v = self.getCoordinateTransformation().getVolumeFactor()
187 self.__B_r = background_magnetic_flux_density * s * v
188 self.__B_b = background_magnetic_flux_density / s
189
190 self.__fw = s**2 * v
191 else: # cartesian
192 self.__fw = 1
193 self.__B_r = background_magnetic_flux_density
194 self.__B_b = background_magnetic_flux_density
195
196 # keep track of k used to build PDE:
197 self.__last_k = None
198 # this is just the initial set_up
199 A=self.getPDE().createCoefficient("A")
200 for i in range(self.getDomain().getDim()):
201 A[i,i]=1.
202 self.getPDE().setValue(A=A)
203
204 def rescaleWeights(self, scale=1., k_scale=1.):
205 """
206 rescales the weights such that
207
208 *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*
209
210 :param scale: scale of data weighting factors
211 :type scale: positive ``float``
212 :param k_scale: scale of susceptibility.
213 :type k_scale: ``Scalar``
214 """
215 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
216
217 def getArguments(self, k):
218 """
219 Returns precomputed values shared by `getDefect()` and `getGradient()`.
220
221 :param k: susceptibility
222 :type k: ``Scalar``
223 :return: scalar magnetic potential and corresponding magnetic field
224 :rtype: ``Scalar``, ``Vector``
225 """
226 phi = self.getPotential(k)
227 grad_phi=grad(phi)
228 magnetic_flux_density = k * self.__B_b -(1+k)*grad_phi
229 return phi, grad_phi, magnetic_flux_density
230
231 def __updateA(self, k):
232 """
233 updates PDE coefficient if PDE is used with new k
234 """
235 pde=self.getPDE()
236 if self.__last_k is not k:
237 A=pde.getCoefficient('A')
238 if self.getCoordinateTransformation().isCartesian():
239 for i in range(self.getDomain().getDim()):
240 A[i,i] = 1+k
241 else:
242 for i in range(self.getDomain().getDim()):
243 A[i,i] = (1+k)*self.__fw[i]
244
245 self.__last_k = k
246 pde.setValue(A=A)
247
248 def getPotential(self, k):
249 """
250 Calculates the magnetic potential for a given susceptibility.
251
252 :param k: susceptibility
253 :type k: ``Scalar``
254 :return: magnetic potential
255 :rtype: ``Scalar``
256 """
257 self.__updateA(k)
258 pde=self.getPDE()
259 pde.resetRightHandSideCoefficients()
260 pde.setValue(X = k*self.__B_r)
261 phi=pde.getSolution()
262 return phi
263
264 def getDefect(self, k, phi, grad_phi, magnetic_flux_density):
265 """
266 Returns the value of the defect.
267
268 :param k: susceptibility
269 :type k: ``Scalar``
270 :param phi: corresponding potential
271 :type phi: ``Scalar``
272 :param magnetic_flux_density: magnetic field
273 :type magnetic_flux_density: ``Vector``
274 :rtype: ``float``
275 """
276 return self._getDefect(magnetic_flux_density)
277
278 def getGradient(self, k, phi, grad_phi, magnetic_flux_density):
279 """
280 Returns the gradient of the defect with respect to susceptibility.
281
282 :param k: susceptibility
283 :type k: ``Scalar``
284 :param phi: corresponding potential
285 :type phi: ``Scalar``
286 :param magnetic_flux_density: magnetic field
287 :type magnetic_flux_density: ``Vector``
288 :rtype: ``Scalar``
289 """
290 self.__updateA(k)
291 Y=self.getDefectGradient(magnetic_flux_density)
292 pde=self.getPDE()
293 pde.resetRightHandSideCoefficients()
294 pde.setValue(X=(1+k)*Y)
295 grad_YT=grad(pde.getSolution())
296
297 if self.getCoordinateTransformation().isCartesian(): # then b_r=B_b
298 return inner(grad_YT-Y, self.__B_r-grad_phi)
299 else:
300 return inner(grad_YT,self.__B_r-grad_phi)+inner(Y,grad_phi-self.__B_b)
301

  ViewVC Help
Powered by ViewVC 1.1.26