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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4120 - (show annotations)
Tue Dec 18 04:49:37 2012 UTC (6 years, 10 months ago) by gross
File MIME type: text/x-python
File size: 12637 byte(s)
some improvements to the robutness of the minimizer:
  a break down in the orthogonalization is hnadles via a restart.
  a restart can now be set manually
  the iteration throughs now an exception in case of failure.

It is also possible to set the density anomaly to zero for regions below a certain depth. 




1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2012 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 since 2012 by School of Earth Sciences
13 #
14 ##############################################################################
15
16 """Domain construction from survey data for inversions"""
17
18 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = ['DomainBuilder']
26
27 import logging
28 import numpy as np
29 from esys.escript.util import *
30 from esys.escript import unitsSI as U
31 from esys.ripley import Brick, Rectangle
32 from .datasources import DataSource
33
34 class DomainBuilder(object):
35 """
36 This class is responsible for constructing an escript Domain object with
37 suitable extents and resolution for survey data (`DataSource` objects)
38 that is added to it.
39 """
40
41 def __init__(self, dim=3):
42 """
43 Constructor.
44
45 :param dim: Dimensionality (2 or 3) of the target domain.
46 This has implications for the survey data than can be
47 added. By default a 3D domain is created.
48 :type dim: ``int``
49 """
50 self.logger = logging.getLogger('inv.%s'%self.__class__.__name__)
51 if dim not in (2,3):
52 raise ValueError("Number of dimensions must be 2 or 3")
53 self._domain=None
54 self._dim=dim
55 self._gravity_surveys=[]
56 self._magnetic_surveys=[]
57 self._sources=[]
58 self.setPadding()
59 self.setVerticalExtents()
60 self.__background_magnetic_field = None
61 self.fixDensityBelow()
62 self.fixSusceptibilityBelow()
63
64 def addSource(self, source):
65 """
66 Adds a survey data provider to the domain builder.
67
68 :param source: The data source to be added.
69 :type source: `DataSource`
70 """
71 if not isinstance(source, DataSource):
72 raise TypeError("source is not a DataSource")
73 self._sources.append(source)
74
75 def setFractionalPadding(self, pad_x=None, pad_y=None):
76 """
77 Sets the amount of padding around the dataset as a fraction of the
78 dataset side lengths.
79
80 For example, calling ``setFractionalPadding(0.2, 0.1)`` to a data
81 source with size 10x20 will result in the padded data set size
82 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
83
84 :param pad_x: Padding per side in x direction (default: no padding)
85 :type pad_x: ``float``
86 :param pad_y: Padding per side in y direction (default: no padding)
87 :type pad_y: ``float``
88 :note: `pad_y` is ignored for 2-dimensional datasets.
89 """
90 if pad_x is not None:
91 if pad_x < 0:
92 raise ValueError("setFractionalPadding: Arguments must be non-negative")
93 if pad_x > 10:
94 raise ValueError("setFractionalPadding: Argument too large")
95 if pad_y is not None:
96 if pad_y < 0:
97 raise ValueError("setFractionalPadding: Arguments must be non-negative")
98 if pad_y > 10:
99 raise ValueError("setFractionalPadding: Argument too large")
100 self._padding = [pad_x,pad_y], 'f'
101
102 def setPadding(self, pad_x=None, pad_y=None):
103 """
104 Sets the amount of padding around the dataset in absolute length units.
105
106 The final domain size will be the length in x (in y) of the dataset
107 plus twice the value of `pad_x` (`pad_y`). The arguments must be
108 non-negative.
109
110 :param pad_x: Padding per side in x direction (default: no padding)
111 :type pad_x: ``float``
112 :param pad_y: Padding per side in y direction (default: no padding)
113 :type pad_y: ``float``
114 :note: `pad_y` is ignored for 2-dimensional datasets.
115 """
116 if pad_x is not None:
117 if pad_x < 0:
118 raise ValueError("setPadding: Arguments must be non-negative")
119 if pad_y is not None:
120 if pad_y < 0:
121 raise ValueError("setPadding: Arguments must be non-negative")
122 self._padding = [pad_x,pad_y], 'l'
123
124 def setElementPadding(self, pad_x=None, pad_y=None):
125 """
126 Sets the amount of padding around the dataset in number of elements.
127
128 When the domain is constructed `pad_x` (`pad_y`) elements are added
129 on each side of the x- (y-) dimension. The arguments must be
130 non-negative.
131
132 :param pad_x: Padding per side in x direction (default: no padding)
133 :type pad_x: ``int``
134 :param pad_y: Padding per side in y direction (default: no padding)
135 :type pad_y: ``int``
136 :note: `pad_y` is ignored for 2-dimensional datasets.
137 """
138 if pad_x is not None:
139 if type(pad_x) is not int:
140 raise TypeError("setElementPadding expects integer arguments")
141 if pad_x < 0:
142 raise ValueError("setElementPadding: Arguments must be non-negative")
143 if pad_y is not None:
144 if type(pad_y) is not int:
145 raise TypeError("setElementPadding expects integer arguments")
146 if pad_y < 0:
147 raise ValueError("setElementPadding: Arguments must be non-negative")
148 self._padding = [pad_x,pad_y], 'e'
149
150 def getGravitySurveys(self):
151 """
152 Returns a list of gravity surveys, see `` getSurveys`` for details
153 """
154 return self.getSurveys(DataSource.GRAVITY)
155
156 def getMagneticSurveys(self):
157 """
158 Returns a list of magnetic surveys, see `` getSurveys`` for details
159 """
160 return self.getSurveys(DataSource.MAGNETIC)
161 def fixDensityBelow(self,depth=None):
162 """
163 defines the depth below which density anomaly is set to zero.
164 """
165 self.__fix_density_below=depth
166
167 def fixSusceptibilityBelow(self,depth=None):
168 """
169 defines the depth below which density anomaly is set to zero.
170 """
171 self.__fix_susceptibility_below=depth
172
173
174
175 def getSurveys(self, datatype):
176 """
177 Returns a list of `Data` objects for all gravity surveys available to
178 this domain builder.
179
180 :return: List of gravity surveys which are tuples (anomaly,error).
181 :rtype: ``list``
182 """
183 if len(self._gravity_surveys)==0:
184 for src in self._sources:
185 if src.getDataType()==datatype:
186 survey=src.getSurveyData(self.getDomain(), self._dom_origin, self._dom_NE, self._spacing)
187 self._gravity_surveys.append(survey)
188 return self._gravity_surveys
189
190 def setBackgroundMagneticFluxDensity(self, B):
191 """
192 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
193 """
194 self.__background_magnetic_field=B
195
196 def getBackgroundMagneticFluxDensity(self):
197 """
198 returns the back ground magnetic flux density.
199 """
200 B = self.__background_magnetic_field
201 # this is for Cartesian (FIXME ?)
202 if self._dim<3:
203 return np.array([-B[2], -B[0]])
204 else:
205 return np.array([-B[1], -B[2], -B[0]])
206
207 def getSetDensityMask(self):
208 z=self.getDomain().getX()[self._dim-1]
209 m = whereNonNegative(z)
210 if self.__fix_density_below:
211 m += whereNonPositive(z+self.__fix_density_below)
212
213 return m
214
215 def getSetSusceptibilityMask(self):
216 z=self.getDomain().getX()[self._dim-1]
217 m = whereNonNegative(z)
218 if self.__fix_susceptibility_below:
219 m += whereNonPositive(z+self.__fix_susceptibility_below)
220 return m
221
222
223 def getDomain(self):
224 """
225 Returns a domain that spans the data area plus padding.
226 The domain is created the first time this method is called, subsequent
227 calls return the same domain so anything that affects the domain
228 (such as padding) needs to be set beforehand.
229
230 :return: The escript domain for this data source.
231 :rtype: `esys.escript.Domain`
232 """
233 if self._domain is None:
234 self._domain=self.__createDomain()
235 return self._domain
236
237 def setVerticalExtents(self, depth=40000., air_layer=10000., num_cells=25):
238 """
239 This method sets the target domain parameters for the vertical
240 dimension.
241
242 :param depth: Depth in meters of the domain.
243 :type depth: ``float``
244 :param air_layer: Depth of the layer above sea level in meters
245 :type air_layer: ``float``
246 :param num_cells: Number of domain elements for the entire vertical
247 dimension
248 :type num_cells: ``int``
249 """
250 self._v_depth=depth
251 self._v_air_layer=air_layer
252 self._v_num_cells=num_cells
253
254 def __getTotalExtentsWithPadding(self):
255 """
256 Helper method that computes origin and number of elements
257 after adding padding to the bounding box of all available survey data.
258 """
259 X0, NX, DX = self.__getTotalExtents()
260 DIM=len(X0)
261 frac=[]
262 # padding is applied to each side so multiply by 2 to get the total
263 # amount of padding per dimension
264 pad, pt = self._padding
265 for i in range(DIM):
266 if pad[i] is None: continue
267 if pt == 'f': # fraction of side length
268 frac.append(2.*pad[i])
269 elif pt == 'e': # number of elements
270 frac.append(2.*pad[i]/float(NX[i]))
271 else: # absolute length
272 f=pad[i]/DX[i]
273 frac.append(2.*f/float(NX[i]))
274
275 # calculate new number of elements
276 NX_padded=[int(round(NX[i]*(1+frac[i]))) for i in range(DIM)]
277 NXdiff=[NX_padded[i]-NX[i] for i in range(DIM)]
278 X0_padded=[X0[i]-NXdiff[i]/2.*DX[i] for i in range(DIM)]
279 return X0_padded, NX_padded, DX
280
281 def __getTotalExtents(self):
282 """
283 Helper method that computes the origin, number of elements and
284 minimal element spacing taking into account all available survey data.
285 """
286 if len(self._sources)==0:
287 raise ValueError("No data")
288 X0, NE, DX = self._sources[0].getDataExtents()
289 XN=[X0[i]+NE[i]*DX[i] for i in range(len(NE))]
290
291 for src in self._sources[1:]:
292 d_x0, d_ne, d_dx = src.getDataExtents()
293 for i in range(len(d_x0)):
294 X0[i]=min(X0[i], d_x0[i])
295 for i in range(len(d_dx)):
296 DX[i]=min(DX[i], d_dx[i])
297 for i in range(len(d_ne)):
298 XN[i]=max(XN[i], d_x0[i]+d_ne[i]*d_dx[i])
299 NE=[int((XN[i]-X0[i])/DX[i]) for i in range(len(XN))]
300 return X0, NE, DX
301
302 def __createDomain(self):
303 """
304 Creates and returns an escript domain that spans the entire area of
305 available data plus a padding zone. This method is called only once
306 the first time `getDomain()` is invoked.
307
308 :return: The escript domain
309 :rtype: `esys.escript.Domain`
310 """
311 X0, NX, DX = self.__getTotalExtentsWithPadding()
312
313 # number of domain elements
314 NE = NX + [self._v_num_cells]
315
316 # origin of domain
317 origin = X0 + [-self._v_depth]
318 self._dom_origin = [np.floor(oi) for oi in origin]
319
320 # cell size / point spacing
321 spacing = DX + [(self._v_depth+self._v_air_layer)/self._v_num_cells]
322 self._spacing = [float(np.floor(si)) for si in spacing]
323
324 lo=[(self._dom_origin[i], self._dom_origin[i]+NE[i]*self._spacing[i]) for i in range(self._dim)]
325 if self._dim==3:
326 dom=Brick(*NE, l0=lo[0], l1=lo[1], l2=lo[2])
327 else:
328 dom=Rectangle(*NE, l0=lo[0], l1=lo[1])
329
330 # ripley may internally adjust NE and length, so recompute
331 self._dom_len=[sup(dom.getX()[i])-inf(dom.getX()[i]) for i in range(self._dim)]
332 self._dom_NE=[int(self._dom_len[i]/self._spacing[i]) for i in range(self._dim)]
333
334 self.logger.debug("Domain size: "+str(self._dom_NE))
335 self.logger.debug(" length: "+str(self._dom_len))
336 self.logger.debug(" origin: "+str(self._dom_origin))
337 return dom
338

  ViewVC Help
Powered by ViewVC 1.1.26