/[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 4098 - (hide annotations)
Fri Dec 7 03:45:20 2012 UTC (6 years, 11 months ago) by caltinay
File MIME type: text/x-python
File size: 5899 byte(s)
Unbroke my breakage.

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 4098 def __init__(self, domain, z0=None, rho0=None, drho=None, beta=None):
100 caltinay 4097 """
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 caltinay 4098 if rho0 is None: rho0=0.
115 caltinay 4097 if drho is None: drho=2750*U.kg/U.m**3
116 caltinay 4098 if beta is None: beta=2.
117 caltinay 4097
118     self.domain=domain
119     if z0 is not None:
120     DIM=self.domain.getDim()
121     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
122     a = drho * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
123     else:
124     a = drho
125     super(DensityMapping,self).__init__(a=a, p0=rho0)
126    
127 gross 4074 class SusceptibilityMapping(LinearMapping):
128     """
129     Susceptibility mapping with depth weighting
130 caltinay 4097
131     *k = k0 + dk * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
132    
133 gross 4074 """
134 caltinay 4098 def __init__(self, domain, z0=None, k0=None, dk=None, beta=None):
135 caltinay 4097 """
136     set up mapping
137    
138     :param domain: domain of the mapping
139     :type domain: ``Domain``
140     :param z0: depth weighting offset. If not present no depth scaling is applied.
141     :type z0: scalar
142     :param k0: reference density, defaults to 0
143     :type k0: scalar
144     :param dk: susceptibility scale, defaults to 1
145     :type dk: scalar
146     :param beta: depth weighting exponent, defaults to 2
147     :type beta: ``float``
148     """
149 caltinay 4098 if k0 is None: k0=0.
150     if beta is None: beta=2.
151     if dk is None: dk=1.
152 caltinay 4097 self.domain=domain
153 caltinay 4098
154 caltinay 4097 if z0 is not None:
155     DIM=self.domain.getDim()
156     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
157     a = dk * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
158     else:
159     a = dk
160     super(SusceptibilityMapping,self).__init__(a=a, p0=k0)
161    
162     # needs REVISION
163 caltinay 3946 class BoundedRangeMapping(Mapping):
164     """
165     Maps an unbounded parameter to a bounded range. The mapping is smooth and
166     continuous.
167     """
168    
169     def __init__(self, s_min=0, s_max=1):
170     if not s_min < s_max:
171     raise ValueError("value for s_min must be less than the value for s_max.")
172     self.s_min=s_min
173     self.s_max=s_max
174    
175     def getValue(self, m):
176     """
177     returns the value of the mapping for m
178     """
179     return (self.s_max+self.s_min)/2. + (self.s_max-self.s_min)/2. * tanh(m)
180    
181     def getDerivative(self, m):
182     """
183     returns the value for the derivative of the mapping for m
184     """
185     return ((self.s_max-self.s_min)/2.) * (1.-tanh(m)**2.)
186    
187     def getInverse(self, s):
188     """
189     returns the value of the inverse of the mapping for s
190     """
191     if not (inf(s) > self.s_min and sup(s) < self.s_max):
192     raise ValueError("s is out of range [%f,%f]"%(inf(s),sup(s)))
193    
194     return 1./2. * log( (s-self.s_min)/(self.s_max-s) )
195    

  ViewVC Help
Powered by ViewVC 1.1.26