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

  ViewVC Help
Powered by ViewVC 1.1.26