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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5088 - (show annotations)
Tue Jul 1 00:32:22 2014 UTC (5 years, 1 month ago) by jduplessis
File MIME type: text/x-python
File size: 43826 byte(s)
option for direct solver with mtInversion
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2014 by University of Queensland
5 # http://www.uq.edu.au
6 #
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12 # Development 2012-2013 by School of Earth Sciences
13 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
14 #
15 ##############################################################################
16
17 """Collection of forward models that define the inversion problem"""
18
19 __copyright__="""Copyright (c) 2003-2014 by 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__ = ['ForwardModel','ForwardModelWithPotential','GravityModel','MagneticModel', 'SelfDemagnetizationModel', 'AcousticWaveForm', 'MT2DModelTEMode']
27
28 from esys.escript import unitsSI as U
29 from esys.escript import Data, Vector, Scalar, Function, DiracDeltaFunctions, FunctionOnBoundary, Solution, length, exp
30 from esys.escript.linearPDEs import LinearSinglePDE, LinearPDE, SolverOptions
31 from .coordinates import makeTranformation
32 from esys.escript.util import *
33 from math import pi as PI
34 from esys.weipa import saveSilo
35 import numpy as np
36
37
38 class ForwardModel(object):
39 """
40 An abstract forward model that can be plugged into a cost function.
41 Subclasses need to implement `getDefect()`, `getGradient()`, and possibly
42 `getArguments()` and 'getCoordinateTransformation'.
43 """
44 def __init__(self):
45 pass
46
47 def getArguments(self, x):
48 return ()
49
50 def getDefect(self, x, *args):
51 raise NotImplementedError
52
53 def getGradient(self, x, *args):
54 raise NotImplementedError
55
56 def getCoordinateTransformation(self):
57 return None
58
59
60 class ForwardModelWithPotential(ForwardModel):
61 """
62 Base class for a forward model using a potential such as magnetic or
63 gravity. It defines a cost function:
64
65 defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) ) ** 2 )
66
67 where s runs over the survey, weight_i are weighting factors, data_i are
68 the data, and r_i are the results produced by the forward model.
69 It is assumed that the forward model is produced through postprocessing
70 of the solution of a potential PDE.
71 """
72 def __init__(self, domain, w, data, coordinates=None,
73 fixPotentialAtBottom=False,
74 tol=1e-8):
75 """
76 initializes a new forward model with potential.
77
78 :param domain: domain of the model
79 :type domain: `Domain`
80 :param w: data weighting factors
81 :type w: ``Vector`` or list of ``Vector``
82 :param data: data
83 :type data: ``Vector`` or list of ``Vector``
84 :param coordinates: defines coordinate system to be used
85 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
86 :param tol: tolerance of underlying PDE
87 :type tol: positive ``float``
88 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
89 in addition to the top.
90 :type fixPotentialAtBottom: ``bool``
91 """
92 super(ForwardModelWithPotential, self).__init__()
93 self.__domain = domain
94 self.__trafo=makeTranformation(domain, coordinates)
95
96 try:
97 n=len(w)
98 m=len(data)
99 if not m == n:
100 raise ValueError("Length of weight and data must be the same.")
101 self.__weight = w
102 self.__data = data
103 except TypeError:
104 self.__weight = [w]
105 self.__data = [data]
106
107 BX = boundingBox(domain)
108 DIM = domain.getDim()
109 x = domain.getX()
110 self.__pde=LinearSinglePDE(domain)
111 self.__pde.getSolverOptions().setTolerance(tol)
112 self.__pde.setSymmetryOn()
113 z=x[DIM-1]
114 q0=whereZero(z-BX[DIM-1][1])
115 if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
116 self.__pde.setValue(q=q0)
117
118 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
119 self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
120
121 self.__origweight=[]
122 for s in range(len(self.__weight)):
123 # save a copy of the original weights in case of rescaling
124 self.__origweight.append(1.*self.__weight[s])
125
126 if not self.__trafo.isCartesian():
127 fd=1./self.__trafo.getScalingFactors()
128 fw=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
129 for s in range(len(self.__weight)):
130 self.__weight[s] = fw * self.__weight[s]
131 self.__data[s] = fd * self.__data[s]
132
133 def _rescaleWeights(self, scale=1., fetch_factor=1.):
134 """
135 rescales the weights such that
136
137 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
138 """
139 if not scale > 0:
140 raise ValueError("Value for scale must be positive.")
141 A=0
142 # copy back original weights before rescaling
143 self.__weight=[1.*ow for ow in self.__origweight]
144
145 for s in range(len(self.__weight)):
146 A += integrate(abs(inner(self.__weight[s], self.__data[s]) * inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
147 if A > 0:
148 A=sqrt(scale/A)/self.diameter
149 if not self.__trafo.isCartesian():
150 A*=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
151 for s in range(len(self.__weight)):
152 self.__weight[s]*=A
153 else:
154 raise ValueError("Rescaling of weights failed.")
155
156 def getDomain(self):
157 """
158 Returns the domain of the forward model.
159
160 :rtype: `Domain`
161 """
162 return self.__domain
163
164 def getCoordinateTransformation(self):
165 """
166 returns the coordinate transformation being used
167
168 :rtype: ``CoordinateTransformation``
169 """
170 return self.__trafo
171
172 def getPDE(self):
173 """
174 Return the underlying PDE.
175
176 :rtype: `LinearPDE`
177 """
178 return self.__pde
179
180 def _getDefect(self, result):
181 """
182 Returns the defect value.
183
184 :param result: a result vector
185 :type result: `Vector`
186 :rtype: ``float``
187 """
188 A=0.
189 for s in range(len(self.__weight)):
190 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
191 return A/2
192
193 def getDefectGradient(self, result):
194 Y=0.
195 for s in range(len(self.__weight)):
196 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
197 return Y
198
199 def getSurvey(self, index=None):
200 """
201 Returns the pair (data_index, weight_index), where data_i is the data
202 of survey i, weight_i is the weighting factor for survey i.
203 If index is None, all surveys will be returned in a pair of lists.
204 """
205 if index is None:
206 return self.__data, self.__weight
207 if index>=len(self.__data):
208 raise IndexError("Forward model only has %d surveys"%len(self.__data))
209 return self.__data[index], self.__weight[index]
210
211
212
213 class GravityModel(ForwardModelWithPotential):
214 """
215 Forward Model for gravity inversion as described in the inversion
216 cookbook.
217 """
218 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
219 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
220 """
221 Creates a new gravity model on the given domain with one or more
222 surveys (w, g).
223
224 :param domain: domain of the model
225 :type domain: `Domain`
226 :param w: data weighting factors
227 :type w: ``Vector`` or list of ``Vector``
228 :param g: gravity anomaly data
229 :type g: ``Vector`` or list of ``Vector``
230 :param coordinates: defines coordinate system to be used
231 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
232 :param tol: tolerance of underlying PDE
233 :type tol: positive ``float``
234 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
235 in addition to the top.
236 :type fixPotentialAtBottom: ``bool``
237
238
239 :note: It is advisable to call rescaleWeights() to rescale weights before start inversion.
240 """
241 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
242
243 trafo = self.getCoordinateTransformation()
244 if not trafo.isCartesian():
245 self.__G = 4*PI*gravity_constant * trafo.getVolumeFactor()
246 #* trafo.getReferenceSystem().getHeightUnit()**(-3)
247
248 fw = trafo.getScalingFactors()**2 * trafo.getVolumeFactor()
249 A=self.getPDE().createCoefficient("A")
250 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
251 self.getPDE().setValue(A=A)
252 else: # cartesian
253 self.__G = 4*PI*gravity_constant
254 self.getPDE().setValue(A=kronecker(self.getDomain()))
255
256 def rescaleWeights(self, scale=1., rho_scale=1.):
257 """
258 rescales the weights such that
259
260 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
261
262 :param scale: scale of data weighting factors
263 :type scale: positive ``float``
264 :param rho_scale: scale of density.
265 :type rho_scale: ``Scalar``
266 """
267 self._rescaleWeights(scale, self.__G*rho_scale)
268
269 def getArguments(self, rho):
270 """
271 Returns precomputed values shared by `getDefect()` and `getGradient()`.
272
273 :param rho: a suggestion for the density distribution
274 :type rho: ``Scalar``
275 :return: gravity potential and corresponding gravity field.
276 :rtype: ``Scalar``, ``Vector``
277 """
278 phi = self.getPotential(rho)
279 gravity_force = -grad(phi)
280 return phi, gravity_force
281
282 def getPotential(self, rho):
283 """
284 Calculates the gravity potential for a given density distribution.
285
286 :param rho: a suggestion for the density distribution
287 :type rho: ``Scalar``
288 :return: gravity potential
289 :rtype: ``Scalar``
290 """
291 pde=self.getPDE()
292
293 pde.resetRightHandSideCoefficients()
294 pde.setValue(Y=-self.__G*rho)
295 phi=pde.getSolution()
296
297 return phi
298
299 def getDefect(self, rho, phi, gravity_force):
300 """
301 Returns the value of the defect
302
303 :param rho: density distribution
304 :type rho: ``Scalar``
305 :param phi: corresponding potential
306 :type phi: ``Scalar``
307 :param gravity_force: gravity force
308 :type gravity_force: ``Vector``
309 :rtype: ``float``
310 """
311 return self._getDefect(gravity_force)
312
313 def getGradient(self, rho, phi, gravity_force):
314 """
315 Returns the gradient of the defect with respect to density.
316
317 :param rho: density distribution
318 :type rho: ``Scalar``
319 :param phi: corresponding potential
320 :type phi: ``Scalar``
321 :param gravity_force: gravity force
322 :type gravity_force: ``Vector``
323 :rtype: ``Scalar``
324 """
325 pde=self.getPDE()
326 pde.resetRightHandSideCoefficients()
327 pde.setValue(X=self.getDefectGradient(gravity_force))
328 ZT=pde.getSolution()
329 return ZT*(-self.__G)
330
331
332 class MagneticModel(ForwardModelWithPotential):
333 """
334 Forward Model for magnetic inversion as described in the inversion
335 cookbook.
336 """
337 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
338 """
339 Creates a new magnetic model on the given domain with one or more
340 surveys (w, B).
341
342 :param domain: domain of the model
343 :type domain: `Domain`
344 :param w: data weighting factors
345 :type w: ``Vector`` or list of ``Vector``
346 :param B: magnetic field data
347 :type B: ``Vector`` or list of ``Vector``
348 :param tol: tolerance of underlying PDE
349 :type tol: positive ``float``
350 :param background_magnetic_flux_density: background magnetic flux
351 density (in Tesla) with components (B_east, B_north, B_vertical)
352 :type background_magnetic_flux_density: ``Vector`` or list of `float`
353 :param coordinates: defines coordinate system to be used
354 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
355 :param fixPotentialAtBottom: if true potential is fixed to zero at the
356 bottom of the domain in addition to the top
357 :type fixPotentialAtBottom: ``bool``
358 """
359 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
360 background_magnetic_flux_density =interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
361 if not self.getCoordinateTransformation().isCartesian():
362 #self.__F = self.getCoordinateTransformation().getVolumeFactor()
363 self.__B_r=background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
364 self.__B_b=background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
365
366 A=self.getPDE().createCoefficient("A")
367 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
368 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
369 self.getPDE().setValue(A=A)
370 else: # cartesian
371 self.getPDE().setValue(A=kronecker(self.getDomain()))
372 self.__B_r=background_magnetic_flux_density
373 self.__B_b=background_magnetic_flux_density
374
375 def rescaleWeights(self, scale=1., k_scale=1.):
376 """
377 rescales the weights such that
378
379 *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*
380
381 :param scale: scale of data weighting factors
382 :type scale: positive ``float``
383 :param k_scale: scale of susceptibility.
384 :type k_scale: ``Scalar``
385 """
386 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
387
388 def getArguments(self, k):
389 """
390 Returns precomputed values shared by `getDefect()` and `getGradient()`.
391
392 :param k: susceptibility
393 :type k: ``Scalar``
394 :return: scalar magnetic potential and corresponding magnetic field
395 :rtype: ``Scalar``, ``Vector``
396 """
397 phi = self.getPotential(k)
398 magnetic_flux_density = k * self.__B_b -grad(phi)
399 return phi, magnetic_flux_density
400
401 def getPotential(self, k):
402 """
403 Calculates the magnetic potential for a given susceptibility.
404
405 :param k: susceptibility
406 :type k: ``Scalar``
407 :return: magnetic potential
408 :rtype: ``Scalar``
409 """
410 pde=self.getPDE()
411
412 pde.resetRightHandSideCoefficients()
413 pde.setValue(X = k* self.__B_r)
414 phi=pde.getSolution()
415
416 return phi
417
418 def getDefect(self, k, phi, magnetic_flux_density):
419 """
420 Returns the value of the defect.
421
422 :param k: susceptibility
423 :type k: ``Scalar``
424 :param phi: corresponding potential
425 :type phi: ``Scalar``
426 :param magnetic_flux_density: magnetic field
427 :type magnetic_flux_density: ``Vector``
428 :rtype: ``float``
429 """
430 return self._getDefect(magnetic_flux_density)
431
432 def getGradient(self, k, phi, magnetic_flux_density):
433 """
434 Returns the gradient of the defect with respect to susceptibility.
435
436 :param k: susceptibility
437 :type k: ``Scalar``
438 :param phi: corresponding potential
439 :type phi: ``Scalar``
440 :param magnetic_flux_density: magnetic field
441 :type magnetic_flux_density: ``Vector``
442 :rtype: ``Scalar``
443 """
444 Y=self.getDefectGradient(magnetic_flux_density)
445 pde=self.getPDE()
446 pde.resetRightHandSideCoefficients()
447 pde.setValue(X=Y)
448 YT=pde.getSolution()
449 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
450
451 class SelfDemagnetizationModel(ForwardModelWithPotential):
452 """
453 Forward Model for magnetic inversion with self-demagnetization as
454 described in the inversion cookbook.
455 """
456 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
457 """
458 Creates a new magnetic model on the given domain with one or more
459 surveys (w, B).
460
461 :param domain: domain of the model
462 :type domain: `Domain`
463 :param w: data weighting factors
464 :type w: ``Vector`` or list of ``Vector``
465 :param B: magnetic field data
466 :type B: ``Vector`` or list of ``Vector``
467 :param background_magnetic_flux_density: background magnetic flux
468 density (in Tesla) with components (B_east, B_north, B_vertical)
469 :type background_magnetic_flux_density: ``Vector`` or list of `float`
470 :param coordinates: defines coordinate system to be used
471 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
472 :param fixPotentialAtBottom: if true potential is fixed to zero at the
473 bottom of the domain in addition to the top
474 :type fixPotentialAtBottom: ``bool``
475 :param tol: tolerance of underlying PDE
476 :type tol: positive ``float``
477 """
478 super(SelfDemagnetizationModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
479 #=========================================================
480 background_magnetic_flux_density = interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
481 if not self.getCoordinateTransformation().isCartesian():
482 #self.__F = self.getCoordinateTransformation().getVolumeFactor()
483 self.__B_r=background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
484 self.__B_b=background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
485
486 self.__fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
487 else: # cartesian
488 self.__fw=1
489 self.__B_r=background_magnetic_flux_density
490 self.__B_b=background_magnetic_flux_density
491
492 # keep track of k used to build PDE:
493 self.__last_k=None
494 # this is just the initial set_up
495 A=self.getPDE().createCoefficient("A")
496 for i in range(self.getDomain().getDim()): A[i,i]=1.
497 self.getPDE().setValue(A=A)
498
499 def rescaleWeights(self, scale=1., k_scale=1.):
500 """
501 rescales the weights such that
502
503 *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*
504
505 :param scale: scale of data weighting factors
506 :type scale: positive ``float``
507 :param k_scale: scale of susceptibility.
508 :type k_scale: ``Scalar``
509 """
510 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
511
512 def getArguments(self, k):
513 """
514 Returns precomputed values shared by `getDefect()` and `getGradient()`.
515
516 :param k: susceptibility
517 :type k: ``Scalar``
518 :return: scalar magnetic potential and corresponding magnetic field
519 :rtype: ``Scalar``, ``Vector``
520 """
521 phi = self.getPotential(k)
522 grad_phi=grad(phi)
523 magnetic_flux_density = k * self.__B_b -(1+k)*grad_phi
524 return phi, grad_phi, magnetic_flux_density
525
526 def __updateA(self, k):
527 """
528 updates PDE coefficient if PDE is used with new k
529 """
530 pde=self.getPDE()
531 if self.__last_k is not k:
532 A=pde.getCoefficient('A')
533 if self.getCoordinateTransformation().isCartesian():
534 for i in range(self.getDomain().getDim()): A[i,i]=(1+k)
535 else:
536 for i in range(self.getDomain().getDim()): A[i,i]=(1+k)*self.__fw[i]
537 self.__last_k = k
538 pde.setValue(A=A)
539
540 def getPotential(self, k):
541 """
542 Calculates the magnetic potential for a given susceptibility.
543
544 :param k: susceptibility
545 :type k: ``Scalar``
546 :return: magnetic potential
547 :rtype: ``Scalar``
548 """
549 self.__updateA(k)
550 pde=self.getPDE()
551 pde.resetRightHandSideCoefficients()
552 pde.setValue(X = k* self.__B_r)
553 phi=pde.getSolution()
554 return phi
555
556 def getDefect(self, k, phi, grad_phi, magnetic_flux_density):
557 """
558 Returns the value of the defect.
559
560 :param k: susceptibility
561 :type k: ``Scalar``
562 :param phi: corresponding potential
563 :type phi: ``Scalar``
564 :param magnetic_flux_density: magnetic field
565 :type magnetic_flux_density: ``Vector``
566 :rtype: ``float``
567 """
568 return self._getDefect(magnetic_flux_density)
569
570 def getGradient(self, k, phi, grad_phi, magnetic_flux_density):
571 """
572 Returns the gradient of the defect with respect to susceptibility.
573
574 :param k: susceptibility
575 :type k: ``Scalar``
576 :param phi: corresponding potential
577 :type phi: ``Scalar``
578 :param magnetic_flux_density: magnetic field
579 :type magnetic_flux_density: ``Vector``
580 :rtype: ``Scalar``
581 """
582 self.__updateA(k)
583 Y=self.getDefectGradient(magnetic_flux_density)
584 pde=self.getPDE()
585 pde.resetRightHandSideCoefficients()
586 pde.setValue(X=(1+k)*Y)
587 grad_YT=grad(pde.getSolution())
588
589 if self.getCoordinateTransformation().isCartesian(): # then b_r=B_b
590 return inner(grad_YT-Y, self.__B_r-grad_phi)
591 else:
592 return inner(grad_YT,self.__B_r-grad_phi)+inner(Y,grad_phi-self.__B_b)
593
594 class AcousticWaveForm(ForwardModel):
595 """
596 Forward Model for acoustic waveform inversion in the frequency domain.
597 It defines a cost function:
598
599 defect = 1/2 integrate( ( w * ( a * u - data ) ) ** 2 )
600
601 where w are weighting factors, data are the measured data (as a 2-comp
602 vector of real and imaginary part) for real frequency omega, and u is
603 the corresponding result produced by the forward model.
604 u (as a 2-comp vector) is the solution of the complex Helmholtz equation
605 for frequency omega, source F and complex, inverse, squared p-velocity
606 sigma:
607
608 * -u_{ii} - omega**2 * sigma * u = F
609
610 It is assumed that the exact scale of source F is unknown and the scaling
611 factor a of F is calculated by minimizing the defect.
612 """
613 def __init__(self, domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True):
614 """
615 initializes a new forward model with acoustic wave form inversion.
616
617 :param domain: domain of the model
618 :type domain: `Domain`
619 :param w: weighting factors
620 :type w: ``Scalar``
621 :param data: real and imaginary part of data
622 :type data: ``Data`` of shape (2,)
623 :param F: real and imaginary part of source given at Dirac points,
624 on surface or at volume.
625 :type F: ``Data`` of shape (2,)
626 :param coordinates: defines coordinate system to be used (not supported yet)
627 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
628 :param tol: tolerance of underlying PDE
629 :type tol: positive ``float``
630 :param saveMemory: if true stiffness matrix is deleted after solution
631 of PDE to minimize memory requests. This will
632 require more compute time as the matrix needs to be
633 reallocated.
634 :type saveMemory: ``bool``
635 :param scaleF: if true source F is scaled to minimize defect.
636 :type scaleF: ``bool``
637 :param fixAtBottom: if true pressure is fixed to zero at the bottom of
638 the domain
639 :type fixAtBottom: ``bool``
640 """
641 super(AcousticWaveForm, self).__init__()
642 self.__domain = domain
643 self.__omega=omega
644 self.scaleF=scaleF
645 self.__trafo=makeTranformation(domain, coordinates)
646 if not self.getCoordinateTransformation().isCartesian():
647 raise ValueError("Non-Cartesian Coordinates are not supported yet.")
648 if not isinstance(data, Data):
649 raise ValueError("data must be an escript.Data object.")
650 if not data.getFunctionSpace() == FunctionOnBoundary(domain):
651 raise ValueError("data must be defined on boundary")
652 if not data.getShape() == (2,):
653 raise ValueError("data must have shape 2 (real and imaginary part).")
654 if w is None:
655 w = 1.
656 if not isinstance(w, Data):
657 w=Data(w, FunctionOnBoundary(domain))
658 else:
659 if not w.getFunctionSpace() == FunctionOnBoundary(domain):
660 raise ValueError("Weights must be defined on boundary.")
661 if not w.getShape() == ():
662 raise ValueError("Weights must be scalar.")
663 self.__weight = w
664 self.__data = data
665 if scaleF:
666 A = integrate(self.__weight*length(self.__data)**2)
667 if A>0:
668 self.__data*=1./sqrt(A)
669
670 self.__BX = boundingBox(domain)
671 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
672
673 self.__F=Data()
674 self.__f=Data()
675 self.__f_dirac=Data()
676
677 if not isinstance(F, Data):
678 F=interpolate(F, DiracDeltaFunctions(domain))
679 if not F.getShape() == (2,):
680 raise ValueError("Source must have shape (2,) (real and imaginary part).")
681
682 if F.getFunctionSpace() == DiracDeltaFunctions(domain):
683 self.__f_dirac=F
684 elif F.getFunctionSpace() == FunctionOnBoundary(domain):
685 self.__f=F
686 else:
687 self.__F=F
688 self.__tol=tol
689 self.__fixAtBottom=fixAtBottom
690 self.__pde=None
691 if not saveMemory:
692 self.__pde=self.setUpPDE()
693
694 def getSurvey(self, index=None):
695 """
696 Returns the pair (data, weight)
697
698 If argument index is ignored.
699 """
700 return self.__data, self.__weight
701
702 def rescaleWeights(self, scale=1., sigma_scale=1.):
703 """
704 rescales the weights such that
705
706 *integrate( ( w omega**2 * sigma_scale * data * ((1/L_j)**2)**-1) +1 )/(data*omega**2 * ((1/L_j)**2)**-1) * sigma_scale )=scale*
707
708 :param scale: scale of data weighting factors
709 :type scale: positive ``float``
710 :param sigma_scale: scale of 1/vp**2 velocity.
711 :type sigma_scale: ``Scalar``
712 """
713 raise Warning("rescaleWeights is not tested yet.")
714 if not scale > 0:
715 raise ValueError("Value for scale must be positive.")
716 if not sigma_scale*omega**2*d > 0:
717 raise ValueError("Rescaling of weights failed due to zero denominator.")
718 # copy back original weights before rescaling
719 #self.__weight=[1.*ow for ow in self.__origweight]
720 L2=1/length(1/self.edge_length)**2
721 d=Lsup(length(data))
722 A=integrate(self.__weight*(sigma_scale*omega**2*d+1)/(sigma_scale*omega**2*d) )
723 if A > 0:
724 self.__weight*=1./A
725 if self.scaleF: self.__data*=sqrt(A)
726 else:
727 raise ValueError("Rescaling of weights failed.")
728
729 def getDomain(self):
730 """
731 Returns the domain of the forward model.
732
733 :rtype: `Domain`
734 """
735 return self.__domain
736
737 def getCoordinateTransformation(self):
738 """
739 returns the coordinate transformation being used
740
741 :rtype: ``CoordinateTransformation``
742 """
743 return self.__trafo
744
745 def setUpPDE(self):
746 """
747 Creates and returns the underlying PDE.
748
749 :rtype: `LinearPDE`
750 """
751 from esys.escript import getEscriptParamInt
752 if self.__pde is None:
753 pde=LinearPDE(self.__domain, numEquations=2)
754 D=pde.createCoefficient('D')
755 A=pde.createCoefficient('A')
756 A[0,:,0,:]=kronecker(self.__domain.getDim())
757 A[1,:,1,:]=kronecker(self.__domain.getDim())
758 pde.setValue(A=A, D=D)
759 if self.__fixAtBottom:
760 DIM=self.__domain.getDim()
761 z = self.__domain.getX()[DIM-1]
762 pde.setValue(q=whereZero(z-self.__BX[DIM-1][0])*[1,1])
763
764 if getEscriptParamInt("PASO_DIRECT")==0:
765 raise ValueError("Either this build of escript or the current MPI configuration does not support direct solvers.")
766 pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
767 pde.getSolverOptions().setTolerance(self.__tol)
768 pde.setSymmetryOff()
769 else:
770 pde=self.__pde
771 pde.resetRightHandSideCoefficients()
772 return pde
773
774 def getSourceScaling(self, u):
775 """
776 returns the scaling factor s required to rescale source F to minimize defect |s * u- data|^2
777 :param u: value of pressure solution (real and imaginary part)
778 :type u: ``Data`` of shape (2,)
779 :rtype: `complex`
780 """
781 uTu=integrate(self.__weight * length(u)**2)
782 uTar=integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
783 uTai=integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
784 if uTu > 0:
785 return complex(uTar/uTu, uTai/uTu)
786 else:
787 return complex(1.,0)
788
789 def getArguments(self, sigma):
790 """
791 Returns precomputed values shared by `getDefect()` and `getGradient()`.
792
793 :param sigma: a suggestion for complex 1/V**2
794 :type sigma: ``Data`` of shape (2,)
795 :return: solution, uTar, uTai, uTu
796 :rtype: ``Data`` of shape (2,), 3 x `float`
797 """
798 pde=self.setUpPDE()
799 D=pde.getCoefficient('D')
800 D[0,0]=-self.__omega**2 * sigma[0]
801 D[0,1]= self.__omega**2 * sigma[1]
802 D[1,0]=-self.__omega**2 * sigma[1]
803 D[1,1]=-self.__omega**2 * sigma[0]
804 pde.setValue(D=D, Y=self.__F, y=self.__f, y_dirac=self.__f_dirac)
805 u=pde.getSolution()
806
807 uTar=integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
808 uTai=integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
809 uTu = integrate( self.__weight * length(u)**2 )
810 return u, uTar, uTai, uTu
811
812 def getDefect(self, sigma, u, uTar, uTai, uTu):
813 """
814 Returns the defect value.
815
816 :param sigma: a suggestion for complex 1/V**2
817 :type sigma: ``Data`` of shape (2,)
818 :param u: a u vector
819 :type u: ``Data`` of shape (2,)
820 :param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1]))
821 :type uTar: float
822 :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1]))
823 :type uTa: float
824 :param uTu : = integrate( w * (u,u))
825 :type uTu: float
826
827 :rtype: ``float``
828 """
829 # assuming integrate(w * length(data)**2) =1
830 if self.scaleF and abs(uTu) >0:
831 A = 1.-(uTar**2 + uTai**2)/uTu
832 else:
833 A = integrate(self.__weight*length(self.__data)**2)- 2 * uTar + uTu
834 return A/2
835
836 def getGradient(self, sigma, u, uTar, uTai, uTu):
837 """
838 Returns the gradient of the defect with respect to density.
839
840 :param sigma: a suggestion for complex 1/V**2
841 :type sigma: ``Data`` of shape (2,)
842 :param u: a u vector
843 :type u: ``Data`` of shape (2,)
844 :param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1]))
845 :type uTar: float
846 :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1]))
847 :type uTa: float
848 :param uTu : = integrate( w * (u,u))
849 :type uTu: float
850 """
851 pde=self.setUpPDE()
852
853 if self.scaleF and abs(uTu) >0:
854 Z=((uTar**2+uTai**2)/uTu**2) *interpolate(u, self.__data.getFunctionSpace())
855 Z[0]+= (-uTar/uTu) * self.__data[0]+ (-uTai/uTu) * self.__data[1]
856 Z[1]+= (-uTar/uTu) * self.__data[1]+ uTai/uTu * self.__data[0]
857
858 else:
859 Z= u - self.__data
860 if Z.getFunctionSpace() == DiracDeltaFunctions(self.getDomain()):
861 pde.setValue(y_dirac=self.__weight * Z)
862 else:
863 pde.setValue(y=self.__weight * Z)
864 D=pde.getCoefficient('D')
865 D[0,0]=-self.__omega**2 * sigma[0]
866 D[0,1]=-self.__omega**2 * sigma[1]
867 D[1,0]= self.__omega**2 * sigma[1]
868 D[1,1]=-self.__omega**2 * sigma[0]
869 pde.setValue(D=D)
870 ZTo2=pde.getSolution()*self.__omega**2
871 return inner(ZTo2,u)*[1,0]+(ZTo2[1]*u[0]-ZTo2[0]*u[1])*[0,1]
872
873 class MT2DModelTEMode(ForwardModel):
874 """
875 Forward Model for two dimensional MT model in the TE mode for a given
876 frequency omega.
877 It defines a cost function:
878
879 * defect = 1/2 integrate( sum_s w^s * ( E_x - Z_XY^s * H_y ) ) ** 2 *
880
881 where E_x is the horizontal electric field perpendicular to the YZ-domain,
882 horizontal magnetic field H_y=1/(i*omega*mu) * E_{x,z} with complex unit
883 i and permeability mu. The weighting factor w^s is set to
884
885 * w^s(X) = w_0^s/(2pi*eta**2) * exp(-length(X-X^s)**2/(2*eta**2)) *
886
887 where X^s is the location of impedance measurement Z_XY^s, w_0 is the level
888 of confidence (eg. 1/measurement error) and eta is level of spatial
889 confidence.
890
891 E_x is given as solution of the PDE
892
893 * -E_{x,ii} - i omega * mu * sigma * E_x = 0
894
895 where E_x is set to E_0 on the top of the domain. Homogeneous Neuman
896 conditions are assumed elsewhere.
897 """
898 def __init__(self, domain, omega, x, Z_XY, eta, w0=1., mu=4*PI*1e-7, E_x0=1,
899 coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, directSolver=False):
900 """
901 initializes a new forward model.
902
903 :param data: real and imaginary part of data
904 :type data: ``Data`` of shape (2,)
905 :param F: real and imaginary part of source given at Dirac points, on surface or at volume.
906 :type F: ``Data`` of shape (2,)
907
908 :param domain: domain of the model
909 :type domain: `Domain`
910 :param omega: frequency
911 :type omega: positive ``float``
912 :param x: coordinates of measurements
913 :type x: ``list`` of ``tuple`` with ``float``
914 :param Z_XY: measured impedance
915 :type Z_XY: ``list`` of ``complex``
916 :param eta: spatial confidence radius
917 :type eta: positive ``float`` or ``list`` of positive ``float``
918 :param w0: confidence factors for meassurements. If not set one is assumed.
919 :type w0: ``None`` or a list of positive ``float``
920 :param E_x0: value of E_x at top the domain and if `fixAtBottom` is set at the bottom of the domain.
921 :type E_x0: ``float``, ``complex`` or ``Data`` of shape (2,)
922 :param mu: permeability
923 :type mu: ``float``
924 :param coordinates: defines coordinate system to be used (not supported yet)
925 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
926 :param tol: tolerance of underlying PDE
927 :type tol: positive ``float``
928 :param fixAtBottom: if true E_x is set E_x0 at the bottom of the domain
929 :type fixAtBottom: ``bool``
930 :param saveMemory: if true stiffness matrix is deleted after solution
931 of the PDE to minimize memory requests. This will
932 require more compute time as the matrix needs to be
933 reallocated.
934 :type saveMemory: ``bool``
935 """
936 super(MT2DModelTEMode, self).__init__()
937 self.__domain = domain
938 DIM=domain.getDim()
939 self.__trafo=makeTranformation(domain, coordinates)
940 if not self.getCoordinateTransformation().isCartesian():
941 raise ValueError("Non-Cartesian coordinates are not supported yet.")
942 self.__omega=omega
943 self.__mu=mu
944
945 #====================================
946 if not len(x) == len(Z_XY):
947 raise ValueError("Number of data points and number of impedance values don't match.")
948 self.__x=x
949 f=-1./(complex(0,1)*omega*mu)
950 self.__scaledZxy=[ z*f for z in Z_XY ]
951 #====================================
952 if isinstance(eta, float) or isinstance(eta, int) :
953 self.__eta =[ float(eta) for z in Z_XY]
954 else:
955 self.__eta =eta
956 if not len(self.__eta) == len(Z_XY):
957 raise ValueError("Number of confidence radii and number of impedance values don't match.")
958 #====================================
959 if isinstance(w0, float) or isinstance(w0, int):
960 w0 =[ float(w0) for z in Z_XY]
961 if not len(w0) == len(Z_XY):
962 raise ValueError("Number of confidence factors and number of impedance values don't match.")
963
964 self.__wx0=[0.] * len(Z_XY)
965 x=Function(domain).getX()
966 totalS=0
967 s = 0
968 while s < len(self.__scaledZxy):
969 totalS+=w0[s]
970 f=integrate(self.getWeightingFactor(x, 1., self.__x[s], self.__eta[s]))
971 if f < 2*PI*self.__eta[s]**2 * 0.1 :
972 raise ValueError("Zero weight (almost) for data point %s. Change eta or refine mesh."%(s,))
973 self.__wx0[s]=w0[s]/f
974 s += 1
975 if not totalS >0 :
976 raise ValueError("Scaling of weight factors failed as sum is zero.")
977
978 self.__wx0=[ w/totalS for w in self.__wx0 ]
979 #====================================
980 if isinstance(E_x0, float) or isinstance(E_x0, int) :
981 self.__E_x0 =Data((E_x0,0), Solution(domain))
982 elif isinstance(E_x0, tuple):
983 self.__E_x0 =Data((E_x0[0],E_x0[1]), Solution(domain))
984 elif isinstance(E_x0, complex):
985 self.__E_x0 =Data((E_x0.real,E_x0.imag), Solution(domain))
986 else:
987 if not E_x0.getShape() == (2,):
988 raise ValueError("Expected shape of E_x0 is (2,)")
989 self.__E_x0= E_x0
990 #=======================
991 self.__tol=tol
992 self.__fixAtBottom=fixAtBottom
993 self.__directSolver=directSolver
994 self.__pde=None
995 if not saveMemory:
996 self.__pde=self.setUpPDE()
997
998 def getDomain(self):
999 """
1000 Returns the domain of the forward model.
1001
1002 :rtype: `Domain`
1003 """
1004 return self.__domain
1005
1006 def getCoordinateTransformation(self):
1007 """
1008 returns the coordinate transformation being used
1009
1010 :rtype: ``CoordinateTransformation``
1011 """
1012 return self.__trafo
1013
1014 def setUpPDE(self):
1015 """
1016 Return the underlying PDE.
1017
1018 :rtype: `LinearPDE`
1019 """
1020 if self.__pde is None:
1021 pde=LinearPDE(self.__domain, numEquations=2)
1022 if(directSolver == True):
1023 pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
1024 D=pde.createCoefficient('D')
1025 A=pde.createCoefficient('A')
1026 Y=pde.createCoefficient('Y')
1027 X=pde.createCoefficient('X')
1028 A[0,:,0,:]=kronecker(self.__domain.getDim())
1029 A[1,:,1,:]=kronecker(self.__domain.getDim())
1030 pde.setValue(A=A, D=D, X=X, Y=Y)
1031 DIM=self.__domain.getDim()
1032 z = self.__domain.getX()[DIM-1]
1033 q=whereZero(z-sup(z))
1034 if self.__fixAtBottom:
1035 q+=whereZero(z-inf(z))
1036 pde.setValue(q=q*[1,1])
1037 pde.getSolverOptions().setTolerance(self.__tol)
1038 pde.setSymmetryOff()
1039 else:
1040 pde=self.__pde
1041 pde.resetRightHandSideCoefficients()
1042 return pde
1043
1044 def getArguments(self, sigma):
1045 """
1046 Returns precomputed values shared by `getDefect()` and `getGradient()`.
1047
1048 :param sigma: conductivity
1049 :type sigma: ``Data`` of shape (2,)
1050 :return: E_x, E_x,z
1051 :rtype: ``Data`` of shape (2,)
1052 """
1053 pde=self.setUpPDE()
1054 D=pde.getCoefficient('D')
1055 f=self.__omega * self.__mu * sigma
1056 D[0,1]= - f
1057 D[1,0]= f
1058 pde.setValue(D=D, r=self.__E_x0)
1059 u=pde.getSolution()
1060 return u, grad(u)[:,1]
1061
1062 def getWeightingFactor(self, x, wx0, x0, eta):
1063 """
1064 returns the weighting factor
1065 """
1066 return wx0 * exp(length(x-x0)**2*(-0.5/eta**2))
1067
1068 def getDefect(self, sigma, E_x, dE_xdz):
1069 """
1070 Returns the defect value.
1071
1072 :param sigma: a suggestion for complex 1/V**2
1073 :type sigma: ``Data`` of shape (2,)
1074 :param E_x: electric field
1075 :type E_x: ``Data`` of shape (2,)
1076 :param dE_xdz: vertical derivative of electric field
1077 :type dE_xdz: ``Data`` of shape (2,)
1078
1079 :rtype: ``float``
1080 """
1081 x=Function(self.getDomain()).getX()
1082 A=0
1083 u0=E_x[0]
1084 u1=E_x[1]
1085 u01=dE_xdz[0]
1086 u11=dE_xdz[1]
1087
1088 # this cane be done faster!
1089 s = 0
1090 while s < len(self.__scaledZxy):
1091 ws=self.getWeightingFactor(x, self.__wx0[s], self.__x[s], self.__eta[s])
1092 A+=integrate( ws * ( (u0-u01*self.__scaledZxy[s].real+u11*self.__scaledZxy[s].imag)**2 + (u1-u01*self.__scaledZxy[s].imag-u11*self.__scaledZxy[s].real)**2 ) )
1093 s += 1
1094 return A/2
1095
1096 def getGradient(self, sigma, E_x, dE_xdz):
1097 """
1098 Returns the gradient of the defect with respect to density.
1099
1100 :param sigma: a suggestion for complex 1/V**2
1101 :type sigma: ``Data`` of shape (2,)
1102 :param E_x: electric field
1103 :type E_x: ``Data`` of shape (2,)
1104 :param dE_xdz: vertical derivative of electric field
1105 :type dE_xdz: ``Data`` of shape (2,)
1106 """
1107 x=Function(self.getDomain()).getX()
1108 pde=self.setUpPDE()
1109
1110 u0=E_x[0]
1111 u1=E_x[1]
1112 u01=dE_xdz[0]
1113 u11=dE_xdz[1]
1114
1115 D=pde.getCoefficient('D')
1116 Y=pde.getCoefficient('Y')
1117 X=pde.getCoefficient('X')
1118
1119 f=self.__omega * self.__mu * sigma
1120 D[0,1]= f
1121 D[1,0]= -f
1122
1123 s = 0
1124 while s < len(self.__scaledZxy):
1125 ws=self.getWeightingFactor(x, self.__wx0[s], self.__x[s], self.__eta[s])
1126 Y[0]+=ws * (u0 - u01* self.__scaledZxy[s].real + u11* self.__scaledZxy[s].imag)
1127 Y[1]+=ws * (u1 - u01* self.__scaledZxy[s].imag - u11* self.__scaledZxy[s].real)
1128 X[0,1]+=ws * (u01* abs(self.__scaledZxy[s])**2 - u0* self.__scaledZxy[s].real - u1* self.__scaledZxy[s].imag )
1129 X[1,1]+=ws * (u11* abs(self.__scaledZxy[s])**2 + u0* self.__scaledZxy[s].imag - u1* self.__scaledZxy[s].real )
1130 s += 1
1131 pde.setValue(D=D, X=X, Y=Y)
1132 Zstar=pde.getSolution()
1133 return self.__omega * self.__mu * (Zstar[0]*u1-Zstar[1]*u0)
1134

  ViewVC Help
Powered by ViewVC 1.1.26