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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5575 - (show annotations)
Wed Apr 1 03:17:22 2015 UTC (4 years ago) by sshaw
File MIME type: text/x-python
File size: 12419 byte(s)
sphinx now builds python doco including module-level docstrings and class __init__ docstrings, also cleared up some doco errors and added multires documentation
1 from __future__ import print_function
2 from __future__ import division
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2015 by University of Queensland
6 # http://www.uq.edu.au
7 #
8 # Primary Business: Queensland, Australia
9 # Licensed under the Open Software License version 3.0
10 # http://www.opensource.org/licenses/osl-3.0.php
11 #
12 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 """Forward model for acoustic wave forms"""
19
20 __copyright__="""Copyright (c) 2003-2015 by 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__ = ['AcousticWaveForm']
28
29 from .base import ForwardModel
30 from esys.downunder.coordinates import makeTranformation
31 from esys.escript import Data, DiracDeltaFunctions, FunctionOnBoundary
32 from esys.escript.linearPDEs import LinearPDE, SolverOptions
33 from esys.escript.util import *
34 import numpy as np
35
36
37 class AcousticWaveForm(ForwardModel):
38 """
39 Forward Model for acoustic waveform inversion in the frequency domain.
40 It defines a cost function:
41
42 :math: `defect = 1/2 integrate( ( w * ( a * u - data ) ) ** 2 )`
43
44 where w are weighting factors, data are the measured data (as a 2-comp
45 vector of real and imaginary part) for real frequency omega, and u is
46 the corresponding result produced by the forward model.
47 u (as a 2-comp vector) is the solution of the complex Helmholtz equation
48 for frequency omega, source F and complex, inverse, squared p-velocity
49 sigma:
50
51 :math: `-u_{ii} - omega**2 * sigma * u = F`
52
53 It is assumed that the exact scale of source F is unknown and the scaling
54 factor a of F is calculated by minimizing the defect.
55 """
56 def __init__(self, domain, omega, w, data, F, coordinates=None,
57 fixAtBottom=False, tol=1e-8, saveMemory=True, scaleF=True):
58 """
59 initializes a new forward model with acoustic wave form inversion.
60
61 :param domain: domain of the model
62 :type domain: `Domain`
63 :param w: weighting factors
64 :type w: ``Scalar``
65 :param data: real and imaginary part of data
66 :type data: ``Data`` of shape (2,)
67 :param F: real and imaginary part of source given at Dirac points,
68 on surface or at volume.
69 :type F: ``Data`` of shape (2,)
70 :param coordinates: defines coordinate system to be used (not supported yet)
71 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
72 :param tol: tolerance of underlying PDE
73 :type tol: positive ``float``
74 :param saveMemory: if true stiffness matrix is deleted after solution
75 of PDE to minimize memory requests. This will
76 require more compute time as the matrix needs to be
77 reallocated.
78 :type saveMemory: ``bool``
79 :param scaleF: if true source F is scaled to minimize defect.
80 :type scaleF: ``bool``
81 :param fixAtBottom: if true pressure is fixed to zero at the bottom of
82 the domain
83 :type fixAtBottom: ``bool``
84 """
85 super(AcousticWaveForm, self).__init__()
86 self.__trafo = makeTranformation(domain, coordinates)
87 if not self.getCoordinateTransformation().isCartesian():
88 raise ValueError("Non-Cartesian Coordinates are not supported yet.")
89 if not isinstance(data, Data):
90 raise ValueError("data must be an escript.Data object.")
91 if not data.getFunctionSpace() == FunctionOnBoundary(domain):
92 raise ValueError("data must be defined on boundary")
93 if not data.getShape() == (2,):
94 raise ValueError("data must have shape (2,) (real and imaginary part).")
95 if w is None:
96 w = 1.
97 if not isinstance(w, Data):
98 w = Data(w, FunctionOnBoundary(domain))
99 else:
100 if not w.getFunctionSpace() == FunctionOnBoundary(domain):
101 raise ValueError("Weights must be defined on boundary.")
102 if not w.getShape() == ():
103 raise ValueError("Weights must be scalar.")
104
105 self.__domain = domain
106 self.__omega = omega
107 self.__weight = w
108 self.__data = data
109 self.scaleF = scaleF
110 if scaleF:
111 A = integrate(self.__weight*length(self.__data)**2)
112 if A > 0:
113 self.__data*=1./sqrt(A)
114
115 self.__BX = boundingBox(domain)
116 self.edge_lengths = np.asarray(boundingBoxEdgeLengths(domain))
117
118 if not isinstance(F, Data):
119 F=interpolate(F, DiracDeltaFunctions(domain))
120 if not F.getShape() == (2,):
121 raise ValueError("Source must have shape (2,) (real and imaginary part).")
122
123 self.__F=Data()
124 self.__f=Data()
125 self.__f_dirac=Data()
126
127 if F.getFunctionSpace() == DiracDeltaFunctions(domain):
128 self.__f_dirac=F
129 elif F.getFunctionSpace() == FunctionOnBoundary(domain):
130 self.__f=F
131 else:
132 self.__F=F
133 self.__tol=tol
134 self.__fixAtBottom=fixAtBottom
135 self.__pde=None
136 if not saveMemory:
137 self.__pde=self.setUpPDE()
138
139 def getSurvey(self, index=None):
140 """
141 Returns the pair (data, weight)
142
143 If argument index is ignored.
144 """
145 return self.__data, self.__weight
146
147 def rescaleWeights(self, scale=1., sigma_scale=1.):
148 """
149 rescales the weights such that
150
151 :math: integrate( ( w omega**2 * sigma_scale * data * ((1/L_j)**2)**-1) +1 )/(data*omega**2 * ((1/L_j)**2)**-1) * sigma_scale )=scale
152
153 :param scale: scale of data weighting factors
154 :type scale: positive ``float``
155 :param sigma_scale: scale of 1/vp**2 velocity.
156 :type sigma_scale: ``Scalar``
157 """
158 raise Warning("rescaleWeights is not tested yet.")
159 if not scale > 0:
160 raise ValueError("Value for scale must be positive.")
161 if not sigma_scale*omega**2*d > 0:
162 raise ValueError("Rescaling of weights failed due to zero denominator.")
163 # copy back original weights before rescaling
164 #self.__weight=[1.*ow for ow in self.__origweight]
165 L2=1/length(1/self.edge_length)**2
166 d=Lsup(length(data))
167 A=integrate(self.__weight*(sigma_scale*omega**2*d+1)/(sigma_scale*omega**2*d) )
168 if A > 0:
169 self.__weight*=1./A
170 if self.scaleF:
171 self.__data*=sqrt(A)
172 else:
173 raise ValueError("Rescaling of weights failed.")
174
175 def getDomain(self):
176 """
177 Returns the domain of the forward model.
178
179 :rtype: `Domain`
180 """
181 return self.__domain
182
183 def getCoordinateTransformation(self):
184 """
185 returns the coordinate transformation being used
186
187 :rtype: ``CoordinateTransformation``
188 """
189 return self.__trafo
190
191 def setUpPDE(self):
192 """
193 Creates and returns the underlying PDE.
194
195 :rtype: `LinearPDE`
196 """
197 from esys.escript import getEscriptParamInt
198 if self.__pde is None:
199 if getEscriptParamInt("PASO_DIRECT")==0:
200 raise ValueError("Either this build of escript or the current MPI configuration does not support direct solvers.")
201 pde=LinearPDE(self.__domain, numEquations=2)
202 D=pde.createCoefficient('D')
203 A=pde.createCoefficient('A')
204 A[0,:,0,:]=kronecker(self.__domain.getDim())
205 A[1,:,1,:]=kronecker(self.__domain.getDim())
206 pde.setValue(A=A, D=D)
207 if self.__fixAtBottom:
208 DIM=self.__domain.getDim()
209 z = self.__domain.getX()[DIM-1]
210 pde.setValue(q=whereZero(z-self.__BX[DIM-1][0])*[1,1])
211
212 pde.getSolverOptions().setSolverMethod(SolverOptions.DIRECT)
213 pde.getSolverOptions().setTolerance(self.__tol)
214 pde.setSymmetryOff()
215 else:
216 pde=self.__pde
217 pde.resetRightHandSideCoefficients()
218 return pde
219
220 def getSourceScaling(self, u):
221 """
222 returns the scaling factor s required to rescale source F to minimize defect ``|s * u- data|^2``
223
224 :param u: value of pressure solution (real and imaginary part)
225 :type u: ``Data`` of shape (2,)
226 :rtype: `complex`
227 """
228 uTu = integrate(self.__weight * length(u)**2)
229 uTar = integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
230 uTai = integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
231 if uTu > 0:
232 return complex(uTar/uTu, uTai/uTu)
233 else:
234 return complex(1.,0)
235
236 def getArguments(self, sigma):
237 """
238 Returns precomputed values shared by `getDefect()` and `getGradient()`.
239
240 :param sigma: a suggestion for complex 1/V**2
241 :type sigma: ``Data`` of shape (2,)
242 :return: solution, uTar, uTai, uTu
243 :rtype: ``Data`` of shape (2,), 3 x `float`
244 """
245 pde=self.setUpPDE()
246 D=pde.getCoefficient('D')
247 D[0,0]=-self.__omega**2 * sigma[0]
248 D[0,1]= self.__omega**2 * sigma[1]
249 D[1,0]=-self.__omega**2 * sigma[1]
250 D[1,1]=-self.__omega**2 * sigma[0]
251 pde.setValue(D=D, Y=self.__F, y=self.__f, y_dirac=self.__f_dirac)
252 u=pde.getSolution()
253
254 uTar=integrate(self.__weight * ( u[0]*self.__data[0]+u[1]*self.__data[1]) )
255 uTai=integrate(self.__weight * ( u[0]*self.__data[1]-u[1]*self.__data[0]) )
256 uTu = integrate( self.__weight * length(u)**2 )
257 return u, uTar, uTai, uTu
258
259 def getDefect(self, sigma, u, uTar, uTai, uTu):
260 """
261 Returns the defect value.
262
263 :param sigma: a suggestion for complex 1/V**2
264 :type sigma: ``Data`` of shape (2,)
265 :param u: a u vector
266 :type u: ``Data`` of shape (2,)
267 :param uTar: equals `integrate( w * (data[0]*u[0]+data[1]*u[1]))`
268 :type uTar: `float`
269 :param uTai: equals `integrate( w * (data[1]*u[0]-data[0]*u[1]))`
270 :type uTa: `float`
271 :param uTu: equals `integrate( w * (u,u))`
272 :type uTu: `float`
273
274 :rtype: ``float``
275 """
276 # assuming integrate(w * length(data)**2) =1
277 if self.scaleF and abs(uTu) >0:
278 A = 1.-(uTar**2 + uTai**2)/uTu
279 else:
280 A = integrate(self.__weight*length(self.__data)**2)- 2 * uTar + uTu
281 return A/2
282
283 def getGradient(self, sigma, u, uTar, uTai, uTu):
284 """
285 Returns the gradient of the defect with respect to density.
286
287 :param sigma: a suggestion for complex 1/V**2
288 :type sigma: ``Data`` of shape (2,)
289 :param u: a u vector
290 :type u: ``Data`` of shape (2,)
291 :param uTar: equals `integrate( w * (data[0]*u[0]+data[1]*u[1]))`
292 :type uTar: `float`
293 :param uTai: equals `integrate( w * (data[1]*u[0]-data[0]*u[1]))`
294 :type uTa: `float`
295 :param uTu: equals `integrate( w * (u,u))`
296 :type uTu: `float`
297 """
298 pde=self.setUpPDE()
299
300 if self.scaleF and abs(uTu) >0:
301 Z=((uTar**2+uTai**2)/uTu**2) *interpolate(u, self.__data.getFunctionSpace())
302 Z[0]+= (-uTar/uTu) * self.__data[0]+ (-uTai/uTu) * self.__data[1]
303 Z[1]+= (-uTar/uTu) * self.__data[1]+ uTai/uTu * self.__data[0]
304
305 else:
306 Z = u - self.__data
307 if Z.getFunctionSpace() == DiracDeltaFunctions(self.getDomain()):
308 pde.setValue(y_dirac=self.__weight * Z)
309 else:
310 pde.setValue(y=self.__weight * Z)
311 D=pde.getCoefficient('D')
312 D[0,0]=-self.__omega**2 * sigma[0]
313 D[0,1]=-self.__omega**2 * sigma[1]
314 D[1,0]= self.__omega**2 * sigma[1]
315 D[1,1]=-self.__omega**2 * sigma[0]
316 pde.setValue(D=D)
317 ZTo2=pde.getSolution()*self.__omega**2
318 return inner(ZTo2,u)*[1,0]+(ZTo2[1]*u[0]-ZTo2[0]*u[1])*[0,1]
319

  ViewVC Help
Powered by ViewVC 1.1.26