# Contents of /trunk/downunder/py_src/forwardmodels/magnetotelluric2d.py

Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 8 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 5 # 6 # Primary Business: Queensland, Australia 7 # Licensed under the Open Software License version 3.0 8 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 25 __url__= 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