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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4120 - (hide 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 caltinay 4060
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 caltinay 4097 from esys.escript import unitsSI as U
31 caltinay 4060 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 gross 4106 self.__background_magnetic_field = None
61 gross 4120 self.fixDensityBelow()
62     self.fixSusceptibilityBelow()
63 caltinay 4060
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 caltinay 4108 def setFractionalPadding(self, pad_x=None, pad_y=None):
76 caltinay 4060 """
77 caltinay 4108 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 caltinay 4060 source with size 10x20 will result in the padded data set size
82 caltinay 4108 14x24 (10*(1+2*0.2), 20*(1+2*0.1))
83 caltinay 4060
84 caltinay 4108 :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 caltinay 4060 """
90 caltinay 4108 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 caltinay 4060
106 caltinay 4108 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 gross 4115
150 caltinay 4060 def getGravitySurveys(self):
151 gross 4115 """
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 gross 4120 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 gross 4115
167 gross 4120 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 gross 4115 def getSurveys(self, datatype):
176 caltinay 4060 """
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 gross 4115 if src.getDataType()==datatype:
186 caltinay 4060 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 gross 4115 def setBackgroundMagneticFluxDensity(self, B):
191 caltinay 4060 """
192 gross 4115 sets the back ground magnetic flux density B=(B_r,B_theta, B_phi)
193 caltinay 4060 """
194 gross 4106 self.__background_magnetic_field=B
195 caltinay 4108
196 gross 4115 def getBackgroundMagneticFluxDensity(self):
197 gross 4106 """
198 gross 4115 returns the back ground magnetic flux density.
199 caltinay 4108 """
200 gross 4106 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 caltinay 4060 else:
205 gross 4106 return np.array([-B[1], -B[2], -B[0]])
206 caltinay 4060
207     def getSetDensityMask(self):
208 gross 4120 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 caltinay 4060 def getSetSusceptibilityMask(self):
216 gross 4120 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 caltinay 4060
222 gross 4120
223 caltinay 4060 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 caltinay 4108 pad, pt = self._padding
265 caltinay 4060 for i in range(DIM):
266 caltinay 4108 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 caltinay 4060
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 caltinay 4108
330 caltinay 4060 # 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