/[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 4821 - (show annotations)
Tue Apr 1 04:58:33 2014 UTC (5 years, 4 months ago) by sshaw
File MIME type: text/x-python
File size: 33659 byte(s)
moved SolverOptions to c++, split into SolverOptions for the options and SolverBuddy as the state as a precursor to per-pde solving... does break some use cases (e.g. pde.getSolverOptions().DIRECT will now fail, new value access is with SolverOptions.DIRECT), examples and documentation updated to match
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']
27
28 from esys.escript import unitsSI as U
29 from esys.escript import Data, Vector, Scalar, Function, DiracDeltaFunctions, FunctionOnBoundary
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 import numpy as np
35
36
37 class ForwardModel(object):
38 """
39 An abstract forward model that can be plugged into a cost function.
40 Subclasses need to implement `getDefect()`, `getGradient()`, and possibly
41 `getArguments()` and 'getCoordinateTransformation'.
42 """
43 def __init__(self):
44 pass
45
46 def getArguments(self, x):
47 return ()
48
49 def getDefect(self, x, *args):
50 raise NotImplementedError
51
52 def getGradient(self, x, *args):
53 raise NotImplementedError
54
55 def getCoordinateTransformation(self):
56 return None
57
58
59 class ForwardModelWithPotential(ForwardModel):
60 """
61 Base class for a forward model using a potential such as magnetic or
62 gravity. It defines a cost function:
63
64 defect = 1/2 sum_s integrate( ( weight_i[s] * ( r_i - data_i[s] ) ) ** 2 )
65
66 where s runs over the survey, weight_i are weighting factors, data_i are
67 the data, and r_i are the results produced by the forward model.
68 It is assumed that the forward model is produced through postprocessing
69 of the solution of a potential PDE.
70 """
71 def __init__(self, domain, w, data, coordinates=None,
72 fixPotentialAtBottom=False,
73 tol=1e-8):
74 """
75 initializes a new forward model with potential.
76
77 :param domain: domain of the model
78 :type domain: `Domain`
79 :param w: data weighting factors
80 :type w: ``Vector`` or list of ``Vector``
81 :param data: data
82 :type data: ``Vector`` or list of ``Vector``
83 :param coordinates: defines coordinate system to be used
84 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
85 :param tol: tolerance of underlying PDE
86 :type tol: positive ``float``
87 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
88 in addition to the top.
89 :type fixPotentialAtBottom: ``bool``
90 """
91 super(ForwardModelWithPotential, self).__init__()
92 self.__domain = domain
93 self.__trafo=makeTranformation(domain, coordinates)
94
95 try:
96 n=len(w)
97 m=len(data)
98 if not m == n:
99 raise ValueError("Length of weight and data must be the same.")
100 self.__weight = w
101 self.__data = data
102 except TypeError:
103 self.__weight = [w]
104 self.__data = [data]
105
106 BX = boundingBox(domain)
107 DIM = domain.getDim()
108 x = domain.getX()
109 self.__pde=LinearSinglePDE(domain)
110 self.__pde.getSolverOptions().setTolerance(tol)
111 self.__pde.setSymmetryOn()
112 z=x[DIM-1]
113 q0=whereZero(z-BX[DIM-1][1])
114 if fixPotentialAtBottom: q0+=whereZero(z-BX[DIM-1][0])
115 self.__pde.setValue(q=q0)
116
117 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
118 self.diameter=1./sqrt(sum(1./self.edge_lengths**2))
119
120 self.__origweight=[]
121 for s in range(len(self.__weight)):
122 # save a copy of the original weights in case of rescaling
123 self.__origweight.append(1.*self.__weight[s])
124
125 if not self.__trafo.isCartesian():
126 fd=1./self.__trafo.getScalingFactors()
127 fw=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
128 for s in range(len(self.__weight)):
129 self.__weight[s] = fw * self.__weight[s]
130 self.__data[s] = fd * self.__data[s]
131
132 def _rescaleWeights(self, scale=1., fetch_factor=1.):
133 """
134 rescales the weights such that
135
136 *sum_s integrate( ( weight_i[s] *data_i[s]) (weight_j[s]*1/L_j) * L**2 * fetch_factor )=scale*
137 """
138 if not scale > 0:
139 raise ValueError("Value for scale must be positive.")
140 A=0
141 # copy back original weights before rescaling
142 self.__weight=[1.*ow for ow in self.__origweight]
143
144 for s in range(len(self.__weight)):
145 A += integrate(abs(inner(self.__weight[s], self.__data[s]) * inner(self.__weight[s], 1/self.edge_lengths) * fetch_factor))
146 if A > 0:
147 A=sqrt(scale/A)/self.diameter
148 if not self.__trafo.isCartesian():
149 A*=self.__trafo.getScalingFactors()*sqrt(self.__trafo.getVolumeFactor())
150 for s in range(len(self.__weight)):
151 self.__weight[s]*=A
152 else:
153 raise ValueError("Rescaling of weights failed.")
154
155 def getDomain(self):
156 """
157 Returns the domain of the forward model.
158
159 :rtype: `Domain`
160 """
161 return self.__domain
162
163 def getCoordinateTransformation(self):
164 """
165 returns the coordinate transformation being used
166
167 :rtype: ``CoordinateTransformation``
168 """
169 return self.__trafo
170
171 def getPDE(self):
172 """
173 Return the underlying PDE.
174
175 :rtype: `LinearPDE`
176 """
177 return self.__pde
178
179 def _getDefect(self, result):
180 """
181 Returns the defect value.
182
183 :param result: a result vector
184 :type result: `Vector`
185 :rtype: ``float``
186 """
187 A=0.
188 for s in range(len(self.__weight)):
189 A += integrate( inner(self.__weight[s], self.__data[s]-result)**2 )
190 return A/2
191
192 def getDefectGradient(self, result):
193 Y=0.
194 for s in range(len(self.__weight)):
195 Y = inner(self.__weight[s], self.__data[s]-result) * self.__weight[s] + Y
196 return Y
197
198 def getSurvey(self, index=None):
199 """
200 Returns the pair (data_index, weight_index), where data_i is the data
201 of survey i, weight_i is the weighting factor for survey i.
202 If index is None, all surveys will be returned in a pair of lists.
203 """
204 if index is None:
205 return self.__data, self.__weight
206 if index>=len(self.__data):
207 raise IndexError("Forward model only has %d surveys"%len(self.__data))
208 return self.__data[index], self.__weight[index]
209
210
211
212 class GravityModel(ForwardModelWithPotential):
213 """
214 Forward Model for gravity inversion as described in the inversion
215 cookbook.
216 """
217 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
218 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
219 """
220 Creates a new gravity model on the given domain with one or more
221 surveys (w, g).
222
223 :param domain: domain of the model
224 :type domain: `Domain`
225 :param w: data weighting factors
226 :type w: ``Vector`` or list of ``Vector``
227 :param g: gravity anomaly data
228 :type g: ``Vector`` or list of ``Vector``
229 :param coordinates: defines coordinate system to be used
230 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
231 :param tol: tolerance of underlying PDE
232 :type tol: positive ``float``
233 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
234 in addition to the top.
235 :type fixPotentialAtBottom: ``bool``
236
237
238 :note: It is advisable to call rescaleWeights() to rescale weights before start inversion.
239 """
240 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
241
242 if not self.getCoordinateTransformation().isCartesian():
243 self.__G = 4*PI*gravity_constant * self.getCoordinateTransformation().getVolumeFactor()
244
245 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
246 A=self.getPDE().createCoefficient("A")
247 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
248 self.getPDE().setValue(A=A)
249 else: # cartesian
250 self.__G = 4*PI*gravity_constant
251 self.getPDE().setValue(A=kronecker(self.getDomain()))
252
253 def rescaleWeights(self, scale=1., rho_scale=1.):
254 """
255 rescales the weights such that
256
257 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
258
259 :param scale: scale of data weighting factors
260 :type scale: positive ``float``
261 :param rho_scale: scale of density.
262 :type rho_scale: ``Scalar``
263 """
264 self._rescaleWeights(scale, self.__G*rho_scale)
265
266 def getArguments(self, rho):
267 """
268 Returns precomputed values shared by `getDefect()` and `getGradient()`.
269
270 :param rho: a suggestion for the density distribution
271 :type rho: ``Scalar``
272 :return: gravity potential and corresponding gravity field.
273 :rtype: ``Scalar``, ``Vector``
274 """
275 phi = self.getPotential(rho)
276 gravity_force = -grad(phi)
277 return phi, gravity_force
278
279 def getPotential(self, rho):
280 """
281 Calculates the gravity potential for a given density distribution.
282
283 :param rho: a suggestion for the density distribution
284 :type rho: ``Scalar``
285 :return: gravity potential
286 :rtype: ``Scalar``
287 """
288 pde=self.getPDE()
289
290 pde.resetRightHandSideCoefficients()
291 pde.setValue(Y=-self.__G*rho)
292 phi=pde.getSolution()
293
294 return phi
295
296 def getDefect(self, rho, phi, gravity_force):
297 """
298 Returns the value of the defect
299
300 :param rho: density distribution
301 :type rho: ``Scalar``
302 :param phi: corresponding potential
303 :type phi: ``Scalar``
304 :param gravity_force: gravity force
305 :type gravity_force: ``Vector``
306 :rtype: ``float``
307 """
308 return self._getDefect(gravity_force)
309
310 def getGradient(self, rho, phi, gravity_force):
311 """
312 Returns the gradient of the defect with respect to density.
313
314 :param rho: density distribution
315 :type rho: ``Scalar``
316 :param phi: corresponding potential
317 :type phi: ``Scalar``
318 :param gravity_force: gravity force
319 :type gravity_force: ``Vector``
320 :rtype: ``Scalar``
321 """
322 pde=self.getPDE()
323 pde.resetRightHandSideCoefficients()
324 pde.setValue(X=self.getDefectGradient(gravity_force))
325 ZT=pde.getSolution()
326 return ZT*(-self.__G)
327
328
329 class MagneticModel(ForwardModelWithPotential):
330 """
331 Forward Model for magnetic inversion as described in the inversion
332 cookbook.
333 """
334 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
335 """
336 Creates a new magnetic model on the given domain with one or more
337 surveys (w, B).
338
339 :param domain: domain of the model
340 :type domain: `Domain`
341 :param w: data weighting factors
342 :type w: ``Vector`` or list of ``Vector``
343 :param B: magnetic field data
344 :type B: ``Vector`` or list of ``Vector``
345 :param tol: tolerance of underlying PDE
346 :type tol: positive ``float``
347 :param background_magnetic_flux_density: background magnetic flux density (in Teslar) with components (B_east, B_north, B_vertical)
348 :type background_magnetic_flux_density: ``Vector`` or list of `float`
349 :param coordinates: defines coordinate system to be used
350 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
351 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
352 in addition to the top.
353 :type fixPotentialAtBottom: ``bool``
354 """
355 super(MagneticModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
356 background_magnetic_flux_density =interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
357 if not self.getCoordinateTransformation().isCartesian():
358 #self.__F = self.getCoordinateTransformation().getVolumeFactor()
359 self.__B_r=background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
360 self.__B_b=background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
361
362 A=self.getPDE().createCoefficient("A")
363 fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
364 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
365 self.getPDE().setValue(A=A)
366 else: # cartesian
367 self.getPDE().setValue(A=kronecker(self.getDomain()))
368 self.__B_r=background_magnetic_flux_density
369 self.__B_b=background_magnetic_flux_density
370
371 def rescaleWeights(self, scale=1., k_scale=1.):
372 """
373 rescales the weights such that
374
375 *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*
376
377 :param scale: scale of data weighting factors
378 :type scale: positive ``float``
379 :param k_scale: scale of susceptibility.
380 :type k_scale: ``Scalar``
381 """
382 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
383
384 def getArguments(self, k):
385 """
386 Returns precomputed values shared by `getDefect()` and `getGradient()`.
387
388 :param k: susceptibility
389 :type k: ``Scalar``
390 :return: scalar magnetic potential and corresponding magnetic field
391 :rtype: ``Scalar``, ``Vector``
392 """
393 phi = self.getPotential(k)
394 magnetic_flux_density = k * self.__B_b -grad(phi)
395 return phi, magnetic_flux_density
396
397 def getPotential(self, k):
398 """
399 Calculates the magnetic potential for a given susceptibility.
400
401 :param k: susceptibility
402 :type k: ``Scalar``
403 :return: magnetic potential
404 :rtype: ``Scalar``
405 """
406 pde=self.getPDE()
407
408 pde.resetRightHandSideCoefficients()
409 pde.setValue(X = k* self.__B_r)
410 phi=pde.getSolution()
411
412 return phi
413
414 def getDefect(self, k, phi, magnetic_flux_density):
415 """
416 Returns the value of the defect.
417
418 :param k: susceptibility
419 :type k: ``Scalar``
420 :param phi: corresponding potential
421 :type phi: ``Scalar``
422 :param magnetic_flux_density: magnetic field
423 :type magnetic_flux_density: ``Vector``
424 :rtype: ``float``
425 """
426 return self._getDefect(magnetic_flux_density)
427
428 def getGradient(self, k, phi, magnetic_flux_density):
429 """
430 Returns the gradient of the defect with respect to susceptibility.
431
432 :param k: susceptibility
433 :type k: ``Scalar``
434 :param phi: corresponding potential
435 :type phi: ``Scalar``
436 :param magnetic_flux_density: magnetic field
437 :type magnetic_flux_density: ``Vector``
438 :rtype: ``Scalar``
439 """
440 Y=self.getDefectGradient(magnetic_flux_density)
441 pde=self.getPDE()
442 pde.resetRightHandSideCoefficients()
443 pde.setValue(X=Y)
444 YT=pde.getSolution()
445 return inner(grad(YT),self.__B_r) -inner(Y,self.__B_b)
446
447 class SelfDemagnetizationModel(ForwardModelWithPotential):
448 """
449 Forward Model for magnetic inversion with self-demagnetization as described in the inversion
450 cookbook.
451 """
452 def __init__(self, domain, w, B, background_magnetic_flux_density, coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
453 """
454 Creates a new magnetic model on the given domain with one or more
455 surveys (w, B).
456
457 :param domain: domain of the model
458 :type domain: `Domain`
459 :param w: data weighting factors
460 :type w: ``Vector`` or list of ``Vector``
461 :param B: magnetic field data
462 :type B: ``Vector`` or list of ``Vector``
463 :param tol: tolerance of underlying PDE
464 :type tol: positive ``float``
465 :param background_magnetic_flux_density: background magnetic flux density (in Teslar) with components (B_east, B_north, B_vertical)
466 :type background_magnetic_flux_density: ``Vector`` or list of `float`
467 :param coordinates: defines coordinate system to be used
468 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
469 :param fixPotentialAtBottom: if true potential is fixed to zero at the bottom of the domain
470 in addition to the top.
471 :type fixPotentialAtBottom: ``bool``
472 """
473 super(SelfDemagnetizationModel, self).__init__(domain, w, B, coordinates, fixPotentialAtBottom, tol)
474 #=========================================================
475 background_magnetic_flux_density = interpolate(background_magnetic_flux_density, B[0].getFunctionSpace())
476 if not self.getCoordinateTransformation().isCartesian():
477 #self.__F = self.getCoordinateTransformation().getVolumeFactor()
478 self.__B_r=background_magnetic_flux_density*self.getCoordinateTransformation().getScalingFactors()*self.getCoordinateTransformation().getVolumeFactor()
479 self.__B_b=background_magnetic_flux_density/self.getCoordinateTransformation().getScalingFactors()
480
481 self.__fw=self.getCoordinateTransformation().getScalingFactors()**2*self.getCoordinateTransformation().getVolumeFactor()
482 else: # cartesian
483 self.__fw=1
484 self.__B_r=background_magnetic_flux_density
485 self.__B_b=background_magnetic_flux_density
486
487 # keep track on k used to build PDE:
488 self.__last_k=None
489 # this is just the initial set_up
490 A=self.getPDE().createCoefficient("A")
491 for i in range(self.getDomain().getDim()): A[i,i]=1
492 self.getPDE().setValue(A=A)
493
494
495 def rescaleWeights(self, scale=1., k_scale=1.):
496 """
497 rescales the weights such that
498
499 *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*
500
501 :param scale: scale of data weighting factors
502 :type scale: positive ``float``
503 :param k_scale: scale of susceptibility.
504 :type k_scale: ``Scalar``
505 """
506 self._rescaleWeights(scale, inner(self.__B_r,1/self.edge_lengths ) * k_scale)
507
508 def getArguments(self, k):
509 """
510 Returns precomputed values shared by `getDefect()` and `getGradient()`.
511
512 :param k: susceptibility
513 :type k: ``Scalar``
514 :return: scalar magnetic potential and corresponding magnetic field
515 :rtype: ``Scalar``, ``Vector``
516 """
517 phi = self.getPotential(k)
518 grad_phi=grad(phi)
519 magnetic_flux_density = k * self.__B_b -(1+k)*grad_phi
520 return phi, grad_phi, magnetic_flux_density
521
522 def __updateA(self, k):
523 """
524 updates PDE coefficient if PDE is used with new k
525 """
526 pde=self.getPDE()
527 if not self.__last_k == k:
528 A=pde.getCoefficient('A')
529 if self.getCoordinateTransformation().isCartesian():
530 for i in range(self.getDomain().getDim()): A[i,i]=(1+k)
531 else:
532 for i in range(self.getDomain().getDim()): A[i,i]=(1+k)*self.__fw[i]
533 self.__last_k = k
534 pde.setValue(A=A)
535
536 def getPotential(self, k):
537 """
538 Calculates the magnetic potential for a given susceptibility.
539
540 :param k: susceptibility
541 :type k: ``Scalar``
542 :return: magnetic potential
543 :rtype: ``Scalar``
544 """
545 self.__updateA(k)
546 pde=self.getPDE()
547 pde.resetRightHandSideCoefficients()
548 pde.setValue(X = k* self.__B_r)
549 phi=pde.getSolution()
550 return phi
551
552 def getDefect(self, k, phi, grad_phi, magnetic_flux_density):
553 """
554 Returns the value of the defect.
555
556 :param k: susceptibility
557 :type k: ``Scalar``
558 :param phi: corresponding potential
559 :type phi: ``Scalar``
560 :param magnetic_flux_density: magnetic field
561 :type magnetic_flux_density: ``Vector``
562 :rtype: ``float``
563 """
564 return self._getDefect(magnetic_flux_density)
565
566 def getGradient(self, k, phi, grad_phi, magnetic_flux_density):
567 """
568 Returns the gradient of the defect with respect to susceptibility.
569
570 :param k: susceptibility
571 :type k: ``Scalar``
572 :param phi: corresponding potential
573 :type phi: ``Scalar``
574 :param magnetic_flux_density: magnetic field
575 :type magnetic_flux_density: ``Vector``
576 :rtype: ``Scalar``
577 """
578 self.__updateA(k)
579 Y=self.getDefectGradient(magnetic_flux_density)
580 pde=self.getPDE()
581 pde.resetRightHandSideCoefficients()
582 pde.setValue(X=(1+k)*Y)
583 grad_YT=grad(pde.getSolution())
584
585 if self.getCoordinateTransformation().isCartesian(): # then b_r=B_b
586 return inner(grad_YT-Y, self.__B_r-grad_phi)
587 else:
588 return inner(grad_YT,self.__B_r-grad_phi)+inner(Y,grad_phi-self.__B_b)
589
590 class AcousticWaveForm(ForwardModel):
591 """
592 Forward Model for acoustic waveform inversion in the frequence domain
593 It defines a cost function:
594
595 defect = 1/2 integrate( ( w * ( a * u - data ) ) ** 2 )
596
597 where w are weighting factors, data are the measured data (as a 2-comp vector of real and imaginary part) for real frequency omega,
598 and u is the coresponding result produced by the forward model. u (as a 2-comp vector) is the solution of the
599 complex Helmholtz equation for frequency omega, source F and complex, inverse, squared p-velocity sigma:
600 * -u_{ii} - omega**2 * sigma * u = F
601 It is assumed that the exact scale of source F is unknown and the scaling factor a of F is calculated by minimizing the
602 defect
603 """
604 def __init__(self, domain, omega, w, data, F, coordinates=None, fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True):
605 """
606 initializes a new forward model with acoustic wave form inversion.
607
608 :param domain: domain of the model
609 :type domain: `Domain`
610 :param w: weighting factors
611 :type w: ``Scalar``
612 :param data: real and imaginary part of data
613 :type data: ``Data`` of shape (2,)
614 :param F: real and imaginary part of source given at Dirac points, on surface or at volume.
615 :type F: ``Data`` of shape (2,)
616 :param coordinates: defines coordinate system to be used (not supported yet)
617 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
618 :param tol: tolerance of underlying PDE
619 :type tol: positive ``float``
620 :param saveMemory: if true stiffness matrix is deleted after solution of PDE to
621 minimize memory requests. This will require more compute time as
622 the matrix needs to be reallocated.
623 :type saveMemory: ``bool``
624 :param scaleF: if true source F is scaled to minimize defect.
625 :type scaleF: ``bool``
626 :param fixAtBottom: if true pressure is fixed to zero at the bottom of the domain
627 :type fixAtBottom: ``bool``
628 """
629 super(AcousticWaveForm, self).__init__()
630 self.__domain = domain
631 self.__omega=omega
632 self.scaleF=scaleF
633 self.__trafo=makeTranformation(domain, coordinates)
634 if not self.getCoordinateTransformation().isCartesian():
635 raise ValueError("Non-Cartesian Coordinates are not supported yet.")
636 if not isinstance(data, Data):
637 raise ValueError("data must be escript.Data object.")
638 if not data.getFunctionSpace() == FunctionOnBoundary(domain):
639 raise ValueError("data must be defined on boundary")
640 if not data.getShape() == (2,):
641 raise ValueError("data must have shape 2 (real and imaginary part).")
642 if w == None:
643 w = 1.
644 if not isinstance(w, Data):
645 w=Data(w, FunctionOnBoundary(domain))
646 else:
647 if not w.getFunctionSpace() == FunctionOnBoundary(domain):
648 raise ValueError("Weights must be defined on boundary.")
649 if not w.getShape() == ():
650 raise ValueError("Weights must be scalar.")
651 self.__weight = w
652 self.__data = data
653 if scaleF:
654 A = integrate(self.__weight*length(self.__data)**2)
655 if A>0:
656 self.__data*=1./sqrt(A)
657
658 self.__BX = boundingBox(domain)
659 self.edge_lengths=np.asarray(boundingBoxEdgeLengths(domain))
660
661 self.__F=Data()
662 self.__f=Data()
663 self.__f_dirac=Data()
664
665 if not isinstance(F, Data):
666 F=interpolate(F, DiracDeltaFunctions(domain))
667 if not F.getShape() == (2,):
668 raise ValueError("Sourcemust have shape 2 (real and imaginary part).")
669
670 if F.getFunctionSpace() == DiracDeltaFunctions(domain):
671 self.__f_dirac=F
672 elif F.getFunctionSpace() == FunctionOnBoundary(domain):
673 self.__f=F
674 else:
675 self.__F=F
676 self.__tol=tol
677 self.__fixAtBottom=fixAtBottom
678 self.__pde=None
679 if not saveMemory:
680 self.__pde=self.setUpPDE()
681 def getSurvey(self, index=None):
682 """
683 Returns the pair (data, weight)
684
685 If argument index is ignored.
686 """
687 return self.__data, self.__weight
688
689
690 def rescaleWeights(self, scale=1., sigma_scale=1.):
691 """
692 rescales the weights such that
693
694 *integrate( ( w omega**2 * sigma_scale * data * ((1/L_j)**2)**-1) +1 )/(data*omega**2 * ((1/L_j)**2)**-1) * sigma_scale )=scale*
695
696 :param scale: scale of data weighting factors
697 :type scale: positive ``float``
698 :param sigma_scale: scale of 1/vp**2 velocity.
699 :type sigma_scale: ``Scalar``
700 """
701 raise Warning("rescaleWeights is not tested yet.")
702 if not scale > 0:
703 raise ValueError("Value for scale must be positive.")
704 if not sigma_scale*omega**2*d > 0:
705 raise ValueError("Rescaling of weights failed due to zero denominator.")
706 # copy back original weights before rescaling
707 #self.__weight=[1.*ow for ow in self.__origweight]
708 L2=1/length(1/self.edge_length)**2
709 d=Lsup(length(data))
710 A=integrate(self.__weight*(sigma_scale*omega**2*d+1)/(sigma_scale*omega**2*d) )
711 if A > 0:
712 self.__weight*=1./A
713 if self.scaleF: self.__data*=sqrt(A)
714 else:
715 raise ValueError("Rescaling of weights failed.")
716
717 def getDomain(self):
718 """
719 Returns the domain of the forward model.
720
721 :rtype: `Domain`
722 """
723 return self.__domain
724
725 def getCoordinateTransformation(self):
726 """
727 returns the coordinate transformation being used
728
729 :rtype: ``CoordinateTransformation``
730 """
731 return self.__trafo
732
733 def setUpPDE(self):
734 """
735 Return the underlying PDE.
736
737 :rtype: `LinearPDE`
738 """
739 from esys.escript import getEscriptParamInt
740 if self.__pde is None:
741 pde=LinearPDE(self.__domain, numEquations=2)
742 D=pde.createCoefficient('D')
743 A=pde.createCoefficient('A')
744 A[0,:,0,:]=kronecker(self.__domain.getDim())
745 A[1,:,1,:]=kronecker(self.__domain.getDim())
746 pde.setValue(A=A, D=D)
747 if self.__fixAtBottom:
748 DIM=self.__domain.getDim()
749 z = self.__domain.getX()[DIM-1]
750 pde.setValue(q=whereZero(z-self.__BX[DIM-1][0])*[1,1])
751
752 if getEscriptParamInt("PASO_DIRECT")==0:
753 raise ValueError("Either this build of escript or the current MPI configuration does not support direct solvers.")
754 pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
755 pde.getSolverOptions().setTolerance(self.__tol)
756 pde.setSymmetryOff()
757 else:
758 pde=self.__pde
759 pde.resetRightHandSideCoefficients()
760 return pde
761
762 def getSourceScaling(self, u):
763 """
764 returns the scaling factor s required to rescale source F to minimize defect |s * u- data|^2
765 :param u: value of pressure solution (real and imaginary part)
766 :type u: ``Data`` of shape (2,)
767 :rtype: `complex`
768 """
769 uTu=integrate(self.__weight * length(u)**2)
770 uTar=integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
771 uTai=integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
772 if uTu > 0:
773 return complex(uTar/uTu, uTai/uTu)
774 else:
775 return complex(1.,0)
776
777 def getArguments(self, sigma):
778 """
779 Returns precomputed values shared by `getDefect()` and `getGradient()`.
780
781 :param sigma: a suggestion for complex 1/V**2
782 :type sigma: ``Data`` of shape (2,)
783 :return: solution, uTar, uTai, uTu
784 :rtype: ``Data`` of shape (2,), 3 x `float`
785 """
786 pde=self.setUpPDE()
787 D=pde.getCoefficient('D')
788 D[0,0]=-self.__omega**2 * sigma[0]
789 D[0,1]= self.__omega**2 * sigma[1]
790 D[1,0]=-self.__omega**2 * sigma[1]
791 D[1,1]=-self.__omega**2 * sigma[0]
792 pde.setValue(D=D, Y=self.__F, y=self.__f, y_dirac=self.__f_dirac)
793 u=pde.getSolution()
794
795 uTar=integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
796 uTai=integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
797 uTu = integrate( self.__weight * length(u)**2 )
798 return u, uTar, uTai, uTu
799
800 def getDefect(self, sigma, u, uTar, uTai, uTu):
801 """
802 Returns the defect value.
803
804 :param sigma: a suggestion for complex 1/V**2
805 :type sigma: ``Data`` of shape (2,)
806 :param u: a u vector
807 :type u: ``Data`` of shape (2,)
808 :param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1]))
809 :type uTar: float
810 :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1]))
811 :type uTa: float
812 :param uTu : = integrate( w * (u,u))
813 :type uTu: float
814
815 :rtype: ``float``
816 """
817 # assummuing integrate(w * length(data)**2) =1
818 if self.scaleF and abs(uTu) >0:
819 A=1.-(uTar**2 + uTai**2)/uTu
820 else:
821 A = integrate(self.__weight*length(self.__data)**2)- 2 * uTar + uTu
822 return A/2
823
824 def getGradient(self,sigma, u,uTar, uTai, uTu ):
825 """
826 Returns the gradient of the defect with respect to density.
827
828 :param sigma: a suggestion for complex 1/V**2
829 :type sigma: ``Data`` of shape (2,)
830 :param u: a u vector
831 :type u: ``Data`` of shape (2,)
832 :param uTar : = integrate( w * (data[0]*u[0]+data[1]*u[1]))
833 :type uTar: float
834 :param uTai : = integrate( w * (data[1]*u[0]-data[0]*u[1]))
835 :type uTa: float
836 :param uTu : = integrate( w * (u,u))
837 :type uTu: float
838 """
839 pde=self.setUpPDE()
840
841 if self.scaleF and abs(uTu) >0:
842 Z=((uTar**2+uTai**2)/uTu**2) *interpolate(u, self.__data.getFunctionSpace())
843 Z[0]+= (-uTar/uTu) * self.__data[0]+ (-uTai/uTu) * self.__data[1]
844 Z[1]+= (-uTar/uTu) * self.__data[1]+ uTai/uTu * self.__data[0]
845
846 else:
847 Z= u - self.__data
848 if Z.getFunctionSpace() == DiracDeltaFunctions(self.getDomain()):
849 pde.setValue(y_dirac=self.__weight * Z)
850 else:
851 pde.setValue(y=self.__weight * Z)
852 D=pde.getCoefficient('D')
853 D[0,0]=-self.__omega**2 * sigma[0]
854 D[0,1]=-self.__omega**2 * sigma[1]
855 D[1,0]= self.__omega**2 * sigma[1]
856 D[1,1]=-self.__omega**2 * sigma[0]
857 pde.setValue(D=D)
858 ZTo2=pde.getSolution()*self.__omega**2
859 return inner(ZTo2,u)*[1,0]+(ZTo2[1]*u[0]-ZTo2[0]*u[1])*[0,1]
860
861
862
863
864
865
866

  ViewVC Help
Powered by ViewVC 1.1.26