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

Contents of /trunk/downunder/py_src/forwardmodels/magnetotelluric2d.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: 20766 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
16 """Forward models for 2D MT (TE and TM mode)"""
17
18 from __future__ import division, print_function
19
20 __copyright__="""Copyright (c) 2003-2015 by The University of Queensland
21 http://www.uq.edu.au
22 Primary Business: Queensland, Australia"""
23 __license__="""Licensed under the Open Software License version 3.0
24 http://www.opensource.org/licenses/osl-3.0.php"""
25 __url__="https://launchpad.net/escript-finley"
26
27 __all__ = ['MT2DModelTEMode', 'MT2DModelTMMode']
28
29 from .base import ForwardModel
30 from esys.downunder.coordinates import makeTranformation
31 from esys.escript import Data, Scalar, Vector, Function, FunctionOnBoundary, Solution
32 from esys.escript.linearPDEs import LinearPDE, SolverOptions
33 from esys.escript.util import *
34 from math import pi as PI
35
36
37 class MT2DBase(ForwardModel):
38 """
39 Base class for 2D MT forward models. See `MT2DModelTEMode` and
40 `MT2DModelTMMode` for actual implementations.
41 """
42 def __init__(self, domain, omega, x, Z, eta=None, w0=1., mu=4*PI*1e-7,
43 Ftop=1., fixAtTop=False, fixAboveLevel=None, Fbottom=0.,
44 fixAtBottom=False, coordinates=None, tol=1e-8,
45 saveMemory=False, directSolver=True):
46 """
47 initializes a new forward model.
48
49 :param domain: domain of the model
50 :type domain: `Domain`
51 :param omega: frequency
52 :type omega: positive ``float``
53 :param x: coordinates of measurements
54 :type x: ``list`` of ``tuple`` with ``float``
55 :param Z: measured impedance (possibly scaled)
56 :type Z: ``list`` of ``complex``
57 :param eta: spatial confidence radius
58 :type eta: positive ``float`` or ``list`` of positive ``float``
59 :param w0: confidence factors for meassurements.
60 :type w0: ``None`` or a list of positive ``float``
61 :param mu: permeability
62 :type mu: ``float``
63 :param Ftop: value of field at top of the domain, see `fixAtTop` and
64 `fixAboveLevel`
65 :type Ftop: ``float``, ``complex`` or ``Data`` of shape (2,)
66 :param fixAtTop: if true F is set to Ftop at the top of the domain.
67 If both `fixAtTop` and `fixAboveLevel` are set, then
68 `fixAboveLevel` takes precedence.
69 :type fixAtTop: ``bool``
70 :param fixAboveLevel: level above which F is set to Ftop (typically
71 the level of the air layer).
72 Use `fixAtTop` *or* `fixAboveLevel`, not both.
73 :type fixAboveLevel : ``float`` or ``None``
74 :param Fbottom: value of field at base of the domain
75 :type Fbottom: ``float``, ``complex`` or ``Data`` of shape (2,)
76 :param fixAtBottom: if true F is set to Fbottom at the bottom of the domain
77 :type fixAtBottom: ``bool``
78 :param coordinates: defines coordinate system to be used (not supported yet)
79 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
80 :param tol: tolerance of underlying PDE
81 :type tol: positive ``float``
82 :param saveMemory: if true stiffness matrix is deleted after solution
83 of the PDE to minimize memory use. This will
84 require more compute time as the matrix needs to be
85 reallocated at each iteration.
86 :type saveMemory: ``bool``
87 :param directSolver: if true a direct solver (rather than an iterative
88 solver) will be used to solve the PDE
89 :type directSolver: ``bool``
90 """
91 super(MT2DBase, self).__init__()
92 self.__trafo = makeTranformation(domain, coordinates)
93 if not self.getCoordinateTransformation().isCartesian():
94 raise ValueError("Non-Cartesian coordinates are not supported yet.")
95 if len(x) != len(Z):
96 raise ValueError("Number of data points and number of impedance values don't match.")
97
98 if eta is None:
99 eta = sup(domain.getSize())*0.45
100
101 if isinstance(eta, float) or isinstance(eta, int):
102 eta = [float(eta)]*len(Z)
103 elif not len(eta) == len(Z):
104 raise ValueError("Number of confidence radii and number of impedance values don't match.")
105
106 if isinstance(w0, float) or isinstance(w0, int):
107 w0 =[float(w0)]*len(Z)
108 elif not len(w0) == len(Z):
109 raise ValueError("Number of confidence factors and number of impedance values don't match.")
110
111 self.__domain = domain
112 self._omega_mu = omega * mu
113
114 xx=Function(domain).getX()
115 totalS=0
116 self._Z = [ Scalar(0., Function(domain)), Scalar(0., Function(domain)) ]
117 self._weight = Scalar(0., Function(domain))
118
119 for s in range(len(Z)):
120 chi = self.getWeightingFactor(xx, 1., x[s], eta[s])
121 f = integrate(chi)
122 if f < eta[s]**2 * 0.01 :
123 raise ValueError("Zero weight (almost) for data point %s. Change eta or refine mesh."%(s,))
124 w02 = w0[s]/f
125 totalS += w02
126 self._Z[0] += chi*Z[s].real
127 self._Z[1] += chi*Z[s].imag
128 self._weight += chi*w02/(abs(Z[s])**2)
129
130 if not totalS > 0:
131 raise ValueError("Scaling of weight factors failed as sum is zero.")
132
133 DIM = domain.getDim()
134 z = domain.getX()[DIM-1]
135 self._ztop = sup(z)
136 self._zbottom = inf(z)
137 self._q=Vector(0.,Solution(domain))
138 self._r=Vector(0.,Solution(domain))
139 #====================================
140 if fixAtTop or fixAboveLevel is not None:
141 if fixAboveLevel is not None:
142 m=whereNonNegative(z-fixAboveLevel)
143 else:
144 m=whereZero(z-self._ztop)
145 if isinstance(Ftop, float) or isinstance(Ftop, int):
146 d = Data((Ftop,0), Solution(domain))
147 elif isinstance(Ftop, tuple):
148 d = Data((Ftop[0],Ftop[1]), Solution(domain))
149 elif isinstance(Ftop, complex):
150 d = Data((Ftop.real,Ftop.imag), Solution(domain))
151 else:
152 if not Ftop.getShape() == (2,):
153 raise ValueError("Expected shape of top value is (2,)")
154 d = Ftop
155 self._r+=m*d
156 self._q+=m*[1.,1.]
157 if fixAtBottom:
158 m=whereZero(z-self._zbottom)
159 if isinstance(Fbottom, float) or isinstance(Fbottom, int) :
160 d = Data((Fbottom,0), Solution(domain))
161 elif isinstance(Fbottom, tuple):
162 d = Data((Fbottom[0],Fbottom[1]), Solution(domain))
163 elif isinstance(Fbottom, complex):
164 d = Data((Fbottom.real,Fbottom.imag), Solution(domain))
165 else:
166 if not Fbottom.getShape() == (2,):
167 raise ValueError("Expected shape of top value is (2,)")
168 d = Fbottom
169 self._r+=m*d
170 self._q+=m*[1.,1.]
171 #====================================
172 self.__tol = tol
173 self._directSolver = directSolver
174 self._saveMemory = saveMemory
175 self.__pde = None
176 if not saveMemory:
177 self.__pde = self.setUpPDE()
178
179 def getDomain(self):
180 """
181 Returns the domain of the forward model.
182
183 :rtype: `Domain`
184 """
185 return self.__domain
186
187 def getCoordinateTransformation(self):
188 """
189 returns the coordinate transformation being used
190
191 :rtype: ``CoordinateTransformation``
192 """
193 return self.__trafo
194
195 def setUpPDE(self):
196 """
197 Return the underlying PDE.
198
199 :rtype: `LinearPDE`
200 """
201 if self.__pde is None:
202 DIM=self.__domain.getDim()
203 pde=LinearPDE(self.__domain, numEquations=2)
204 if self._directSolver == True:
205 pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
206 D=pde.createCoefficient('D')
207 A=pde.createCoefficient('A')
208 pde.setValue(A=A, D=D, q=self._q)
209 pde.getSolverOptions().setTolerance(self.__tol)
210 pde.setSymmetryOff()
211 else:
212 pde=self.__pde
213 pde.resetRightHandSideCoefficients()
214 pde.setValue(X=pde.createCoefficient('X'), Y=pde.createCoefficient('Y'))
215 return pde
216
217 def getWeightingFactor(self, x, wx0, x0, eta):
218 """
219 returns the weighting factor
220 """
221 try:
222 origin, spacing, NE = x.getDomain().getGridParameters()
223 cell = [int((x0[i]-origin[i])/spacing[i]) for i in range(2)]
224 midpoint = [origin[i]+cell[i]*spacing[i]+spacing[i]/2. for i in range(2)]
225 return wx0 * whereNegative(length(x-midpoint)-eta)
226 except:
227 return wx0 * whereNegative(length(x-x0)-eta)
228
229 def getArguments(self, x):
230 """
231 Returns precomputed values shared by `getDefect()` and `getGradient()`.
232 Needs to be implemented in subclasses.
233 """
234 raise NotImplementedError
235
236 def getDefect(self, x, Ex, dExdz):
237 """
238 Returns the defect value. Needs to be implemented in subclasses.
239 """
240 raise NotImplementedError
241
242 def getGradient(self, x, Ex, dExdz):
243 """
244 Returns the gradient. Needs to be implemented in subclasses.
245 """
246 raise NotImplementedError
247
248
249 class MT2DModelTEMode(MT2DBase):
250 """
251 Forward Model for two dimensional MT model in the TE mode for a given
252 frequency omega.
253 It defines a cost function:
254
255 * defect = 1/2 integrate( sum_s w^s * ( E_x/H_y - Z_XY^s ) ) ** 2 *
256
257 where E_x is the horizontal electric field perpendicular to the YZ-domain,
258 horizontal magnetic field H_y=1/(i*omega*mu) * E_{x,z} with complex unit
259 i and permeability mu. The weighting factor w^s is set to
260 * w^s(X) = w_0^s *
261 if length(X-X^s) <= eta and zero otherwise. X^s is the location of
262 impedance measurement Z_XY^s, w_0^s is the level
263 of confidence (eg. 1/measurement error) and eta is level of spatial
264 confidence.
265
266 E_x is given as solution of the PDE
267
268 * -E_{x,ii} - i omega * mu * sigma * E_x = 0
269
270 where by default the normal derivative of E_x at the top of the domain
271 is set to Ex_n=1 and E_x is set to zero at the bottom. Homogeneous Neuman
272 conditions are assumed elsewhere. If fixAtTop is set E_x is set to Ex_top.
273 """
274 def __init__(self, domain, omega, x, Z_XY, eta=None, w0=1., mu=4*PI*1e-7,
275 Ex_n=1, coordinates=None, Ex_top=1, fixAtTop=False, tol=1e-8,
276 saveMemory=False, directSolver=True):
277 """
278 initializes a new forward model. See base class for a description of
279 the arguments.
280 """
281
282 f = -1./(complex(0,1)*omega*mu)
283 scaledZXY = [ z*f for z in Z_XY ]
284 self.__Ex_n=complex(Ex_n)
285 super(MT2DModelTEMode, self).__init__(domain=domain, omega=omega, x=x,
286 Z=scaledZXY, eta=eta, w0=w0, mu=mu,
287 Ftop=Ex_top, fixAtTop=fixAtTop, fixAboveLevel= None, Fbottom=0., fixAtBottom=True,
288 coordinates=coordinates, tol=tol, saveMemory=saveMemory, directSolver=directSolver)
289
290
291 def getArguments(self, sigma):
292 """
293 Returns precomputed values shared by `getDefect()` and `getGradient()`.
294
295 :param sigma: conductivity
296 :type sigma: ``Data`` of shape (2,)
297 :return: Ex_, Ex_,z
298 :rtype: ``Data`` of shape (2,)
299 """
300 DIM = self.getDomain().getDim()
301 pde = self.setUpPDE()
302 D = pde.getCoefficient('D')
303 f = self._omega_mu * sigma
304 D[0,1] = -f
305 D[1,0] = f
306 A= pde.getCoefficient('A')
307 A[0,:,0,:]=kronecker(DIM)
308 A[1,:,1,:]=kronecker(DIM)
309
310 Z = FunctionOnBoundary(self.getDomain()).getX()[DIM-1]
311 pde.setValue(A=A, D=D, y=whereZero(Z-self._ztop)*[self.__Ex_n.real,self.__Ex_n.imag], r=self._r)
312 u = pde.getSolution()
313 return u, grad(u)[:,1]
314
315 def getDefect(self, sigma, Ex, dExdz):
316 """
317 Returns the defect value.
318
319 :param sigma: a suggestion for conductivity
320 :type sigma: ``Data`` of shape ()
321 :param Ex: electric field
322 :type Ex: ``Data`` of shape (2,)
323 :param dExdz: vertical derivative of electric field
324 :type dExdz: ``Data`` of shape (2,)
325
326 :rtype: ``float``
327 """
328 x=dExdz.getFunctionSpace().getX()
329 Ex=interpolate(Ex, x.getFunctionSpace())
330 u0=Ex[0]
331 u1=Ex[1]
332 u01=dExdz[0]
333 u11=dExdz[1]
334 scale = self._weight / ( u01**2 + u11**2 )
335
336 Z = self._Z
337 A = integrate(scale * ( (Z[0]**2+Z[1]**2)*(u01**2+u11**2)
338 + 2*Z[1]*(u0*u11-u01*u1)
339 - 2*Z[0]*(u0*u01+u11*u1)
340 + u0**2 + u1**2 ))
341
342 return A/2
343
344 def getGradient(self, sigma, Ex, dExdz):
345 """
346 Returns the gradient of the defect with respect to density.
347
348 :param sigma: a suggestion for conductivity
349 :type sigma: ``Data`` of shape ()
350 :param Ex: electric field
351 :type Ex: ``Data`` of shape (2,)
352 :param dExdz: vertical derivative of electric field
353 :type dExdz: ``Data`` of shape (2,)
354 """
355 pde=self.setUpPDE()
356 DIM = self.getDomain().getDim()
357
358 x=dExdz.getFunctionSpace().getX()
359 Ex=interpolate(Ex, x.getFunctionSpace())
360 u0 = Ex[0]
361 u1 = Ex[1]
362 u01 = dExdz[0]
363 u11 = dExdz[1]
364
365 D=pde.getCoefficient('D')
366 Y=pde.getCoefficient('Y')
367 X=pde.getCoefficient('X')
368 A= pde.getCoefficient('A')
369
370 A[0,:,0,:]=kronecker(DIM)
371 A[1,:,1,:]=kronecker(DIM)
372
373 f = self._omega_mu * sigma
374 D[0,1] = f
375 D[1,0] = -f
376
377 Z = self._Z
378 scale = 1./( u01**2 + u11**2 )
379 scale2 = scale**2
380 scale *= self._weight
381 scale2 *= self._weight
382
383 Y[0] = scale * (u0 - u01*Z[0] + u11*Z[1])
384 Y[1] = scale * (u1 - u01*Z[1] - u11*Z[0])
385 X[0,1] = scale2 * (2*u01*u11*(Z[0]*u1-Z[1]*u0) \
386 + (Z[0]*u0+Z[1]*u1)*(u01**2-u11**2)
387 - u01*(u0**2 + u1**2))
388 X[1,1] = scale2 * (2*u01*u11*(Z[1]*u1+Z[0]*u0) \
389 + (Z[1]*u0-Z[0]*u1)*(u01**2-u11**2)
390 - u11*(u0**2 + u1**2))
391
392 pde.setValue(A=A, D=D, X=X, Y=Y)
393 Zstar=pde.getSolution()
394 return (-self._omega_mu)* (Zstar[1]*u0-Zstar[0]*u1)
395
396
397 class MT2DModelTMMode(MT2DBase):
398 """
399 Forward Model for two-dimensional MT model in the TM mode for a given
400 frequency omega.
401 It defines a cost function:
402
403 * defect = 1/2 integrate( sum_s w^s * ( rho*H_x/Hy - Z_YX^s ) ) ** 2 *
404
405 where H_x is the horizontal magnetic field perpendicular to the YZ-domain,
406 horizontal magnetic field H_y=1/(i*omega*mu) * E_{x,z} with complex unit
407 i and permeability mu. The weighting factor w^s is set to
408 * w^s(X) = w_0^s *
409 if length(X-X^s) <= eta and zero otherwise. X^s is the location of
410 impedance measurement Z_XY^s, w_0^s is the level
411 of confidence (eg. 1/measurement error) and eta is level of spatial
412 confidence.
413
414 H_x is given as solution of the PDE
415
416 * -(rho*H_{x,i})_{,i} + i omega * mu * H_x = 0
417
418 where H_x is set to Hx_top on the top of the domain. Homogeneous Neuman
419 conditions are assumed elsewhere. If fixAtBottom is set H_x is set to Hx_bottom
420 at the bottom of the domain overwrtining the Neuman condition.
421 If fixAboveLevel is set then H_x is set to Hx_top for any location above and including fixAboveLevel
422 typically including teh top boundary.
423 """
424 def __init__(self, domain, omega, x, Z_YX, eta=None, w0=1., mu=4*PI*1e-7,
425 fixAboveLevel=None, Hx_top=1, coordinates=None, Hx_bottom=1.,
426 fixAtBottom=False, tol=1e-8, saveMemory=False,
427 directSolver=True):
428 """
429 initializes a new forward model. See base class for a description of
430 the arguments.
431 """
432 super(MT2DModelTMMode, self).__init__(domain=domain, omega=omega, x=x,
433 Z=Z_YX, eta=eta, w0=w0, mu=mu, Ftop=Hx_top, fixAtTop=True,
434 fixAboveLevel=fixAboveLevel, Fbottom=Hx_bottom,
435 fixAtBottom=fixAtBottom, coordinates=coordinates, tol=tol,
436 saveMemory=saveMemory, directSolver=directSolver)
437
438 def getArguments(self, rho):
439 """
440 Returns precomputed values shared by `getDefect()` and `getGradient()`.
441
442 :param rho: resistivity
443 :type rho: ``Data`` of shape (2,)
444 :return: Hx, grad(Hx)
445 :rtype: ``tuple`` of ``Data``
446 """
447 DIM = self.getDomain().getDim()
448 pde = self.setUpPDE()
449
450 D = pde.getCoefficient('D')
451 f = self._omega_mu
452 D[0,1] = -f
453 D[1,0] = f
454
455 A= pde.getCoefficient('A')
456 for i in range(DIM):
457 A[0,i,0,i]=rho
458 A[1,i,1,i]=rho
459
460 pde.setValue(A=A, D=D, r=self._r)
461 u = pde.getSolution()
462 return u, grad(u)
463
464 def getDefect(self, rho, Hx, g_Hx):
465 """
466 Returns the defect value.
467
468 :param rho: a suggestion for resistivity
469 :type rho: ``Data`` of shape ()
470 :param Hx: magnetic field
471 :type Hx: ``Data`` of shape (2,)
472 :param g_Hx: gradient of magnetic field
473 :type g_Hx: ``Data`` of shape (2,2)
474
475 :rtype: ``float``
476 """
477 x = g_Hx.getFunctionSpace().getX()
478 Hx = interpolate(Hx, x.getFunctionSpace())
479 u0 = Hx[0]
480 u1 = Hx[1]
481 u01 = g_Hx[0,1]
482 u11 = g_Hx[1,1]
483 scale = rho / ( u0**2 + u1**2 )
484
485 Z = self._Z
486 A = integrate(self._weight * ( Z[0]**2 + Z[1]**2
487 + scale*(-2*Z[0]*(u0*u01 + u1*u11)
488 +2*Z[1]*(u1*u01 - u0*u11)
489 +rho*(u01**2 + u11**2)) ))
490 return A/2
491
492 def getGradient(self, rho, Hx, g_Hx):
493 """
494 Returns the gradient of the defect with respect to resistivity.
495
496 :param rho: a suggestion for resistivity
497 :type rho: ``Data`` of shape ()
498 :param Hx: magnetic field
499 :type Hx: ``Data`` of shape (2,)
500 :param g_Hx: gradient of magnetic field
501 :type g_Hx: ``Data`` of shape (2,2)
502 """
503 pde=self.setUpPDE()
504 DIM = self.getDomain().getDim()
505
506 x=g_Hx.getFunctionSpace().getX()
507 Hx=interpolate(Hx, x.getFunctionSpace())
508 u0 = Hx[0]
509 u1 = Hx[1]
510 u00 = g_Hx[0,0]
511 u10 = g_Hx[1,0]
512 u01 = g_Hx[0,1]
513 u11 = g_Hx[1,1]
514
515 A=pde.getCoefficient('A')
516 D=pde.getCoefficient('D')
517 Y=pde.getCoefficient('Y')
518 X=pde.getCoefficient('X')
519
520 for i in range(DIM):
521 A[0,i,0,i]=rho
522 A[1,i,1,i]=rho
523
524 f = self._omega_mu
525 D[0,1] = f
526 D[1,0] = -f
527
528 Z = self._Z
529 scale = 1./( u0**2 + u1**2 )
530 scale2 = scale**2
531 scale *= self._weight
532 scale2 *= rho*self._weight
533 rho_scale = rho*scale
534
535 gscale = u01**2 + u11**2
536
537 Y[0] = scale2 * ( (Z[0]*u01+Z[1]*u11)*(u0**2-u1**2)
538 + 2*u0*u1*(Z[0]*u11-Z[1]*u01)
539 - rho*u0*gscale )
540 Y[1] = scale2 * ( (Z[0]*u11-Z[1]*u01)*(u1**2-u0**2)
541 + 2*u0*u1*(Z[0]*u01+Z[1]*u11)
542 - rho*u1*gscale )
543 X[0,1] = rho_scale * (-Z[0]*u0 + Z[1]*u1 + rho*u01)
544 X[1,1] = rho_scale * (-Z[0]*u1 - Z[1]*u0 + rho*u11)
545
546 pde.setValue(A=A, D=D, X=X, Y=Y)
547 g=grad(pde.getSolution())
548
549 Hstarr_x = g[0,0]
550 Hstari_x = g[1,0]
551 Hstarr_z = g[0,1]
552 Hstari_z = g[1,1]
553 return -scale*(u0*(Z[0]*u01+Z[1]*u11)+u1*(Z[0]*u11-Z[1]*u01)-rho*gscale)\
554 - Hstarr_x*u00 - Hstarr_z*u01 - Hstari_x*u10 - Hstari_z*u11
555

  ViewVC Help
Powered by ViewVC 1.1.26