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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4097 - (hide annotations)
Fri Dec 7 01:18:35 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 5731 byte(s)
Minor edits.

1 caltinay 3946
2 jfenwick 3981 ##############################################################################
3 caltinay 3946 #
4     # Copyright (c) 2003-2012 by University of Queensland
5 jfenwick 3981 # http://www.uq.edu.au
6 caltinay 3946 #
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 jfenwick 3981 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
12     # Development since 2012 by School of Earth Sciences
13     #
14     ##############################################################################
15 caltinay 3946
16 caltinay 4007 """Collection of parametrizations that map physical values to model parameters
17     and back"""
18    
19 caltinay 3946 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
20 jfenwick 3981 http://www.uq.edu.au
21 caltinay 3946 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 gross 4074 __all__ = ['Mapping', 'DensityMapping', 'SusceptibilityMapping', 'BoundedRangeMapping', 'LinearMapping']
27 caltinay 3947
28 gross 4074 from esys.escript import inf, sup, log, tanh, boundingBoxEdgeLengths
29     import esys.escript.unitsSI as U
30 caltinay 3946
31     class Mapping(object):
32     """
33 caltinay 4097 An abstract mapping class to map level set functions *m* to physical
34     parameters *p*.
35 caltinay 3946 """
36    
37     def __init__(self, *args):
38     pass
39    
40     def __call__(self, m):
41     """
42     short for getValue(m).
43     """
44     return self.getValue(m)
45    
46     def getValue(self, m):
47     """
48     returns the value of the mapping for m
49     """
50     raise NotImplementedError
51    
52     def getDerivative(self, m):
53     """
54     returns the value for the derivative of the mapping for m
55     """
56     raise NotImplementedError
57    
58     def getInverse(self, s):
59     """
60 gross 4074 returns the value of the inverse of the mapping for physical parameter p
61 caltinay 3946 """
62     raise NotImplementedError
63    
64 gross 4074 class LinearMapping(Mapping):
65     """
66 jfenwick 4080 Maps a parameter by a linear transformation p = a * m + p0
67 gross 4074 """
68 caltinay 3946
69 gross 4074 def __init__(self, a=1, p0=0):
70     self.__a=a
71     self.__p0=p0
72     self.__a_inv = 1/a
73 caltinay 4097
74 gross 4074 def getValue(self, m):
75     """
76     returns the value of the mapping for m
77     """
78     return self.__a*m + self.__p0
79    
80     def getDerivative(self, m):
81     """
82     returns the value for the derivative of the mapping for m
83     """
84     return self.__a
85    
86     def getInverse(self, p):
87     """
88     returns the value of the inverse of the mapping for s
89     """
90     return self.__a_inv * ( p - self.__p0)
91    
92     class DensityMapping(LinearMapping):
93     """
94     Density mapping with depth weighting
95 caltinay 4097
96     *rho = rho0 + drho * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
97    
98 gross 4074 """
99 caltinay 4097 def __init__(self, domain, z0=None, rho0=0., drho=None, beta=2.):
100     """
101     initializes the mapping
102    
103     :param domain: domain of the mapping
104     :type domain: ``Domain``
105     :param z0: depth weighting offset. If not present no depth scaling is applied.
106     :type z0: scalar
107     :param rho0: reference density, defaults to 0
108     :type rho0: scalar
109     :param drho: density scale. By default density of granite = 2750kg/m**3 is used.
110     :type drho: scalar
111     :param beta: depth weighting exponent, defaults to 2
112     :type beta: ``float``
113     """
114     if drho is None: drho=2750*U.kg/U.m**3
115    
116     self.domain=domain
117     if z0 is not None:
118     DIM=self.domain.getDim()
119     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
120     a = drho * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
121     else:
122     a = drho
123     super(DensityMapping,self).__init__(a=a, p0=rho0)
124    
125 gross 4074 class SusceptibilityMapping(LinearMapping):
126     """
127     Susceptibility mapping with depth weighting
128 caltinay 4097
129     *k = k0 + dk * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
130    
131 gross 4074 """
132 caltinay 4097 def __init__(self, domain, z0=None, k0=0., dk=1., beta=2.):
133     """
134     set up mapping
135    
136     :param domain: domain of the mapping
137     :type domain: ``Domain``
138     :param z0: depth weighting offset. If not present no depth scaling is applied.
139     :type z0: scalar
140     :param k0: reference density, defaults to 0
141     :type k0: scalar
142     :param dk: susceptibility scale, defaults to 1
143     :type dk: scalar
144     :param beta: depth weighting exponent, defaults to 2
145     :type beta: ``float``
146     """
147     self.domain=domain
148     if z0 is not None:
149     DIM=self.domain.getDim()
150     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
151     a = dk * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
152     else:
153     a = dk
154     super(SusceptibilityMapping,self).__init__(a=a, p0=k0)
155    
156     # needs REVISION
157 caltinay 3946 class BoundedRangeMapping(Mapping):
158     """
159     Maps an unbounded parameter to a bounded range. The mapping is smooth and
160     continuous.
161     """
162    
163     def __init__(self, s_min=0, s_max=1):
164     if not s_min < s_max:
165     raise ValueError("value for s_min must be less than the value for s_max.")
166     self.s_min=s_min
167     self.s_max=s_max
168    
169     def getValue(self, m):
170     """
171     returns the value of the mapping for m
172     """
173     return (self.s_max+self.s_min)/2. + (self.s_max-self.s_min)/2. * tanh(m)
174    
175     def getDerivative(self, m):
176     """
177     returns the value for the derivative of the mapping for m
178     """
179     return ((self.s_max-self.s_min)/2.) * (1.-tanh(m)**2.)
180    
181     def getInverse(self, s):
182     """
183     returns the value of the inverse of the mapping for s
184     """
185     if not (inf(s) > self.s_min and sup(s) < self.s_max):
186     raise ValueError("s is out of range [%f,%f]"%(inf(s),sup(s)))
187    
188     return 1./2. * log( (s-self.s_min)/(self.s_max-s) )
189    

  ViewVC Help
Powered by ViewVC 1.1.26