/[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 4074 - (hide annotations)
Thu Nov 15 03:30:59 2012 UTC (7 years ago) by gross
File MIME type: text/x-python
File size: 5832 byte(s)
this will break the downundermodule: major revision of the classes.
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 gross 4074 An abstract mapping class to map level set functions *m* to physical parameters *p*
34 caltinay 3946 """
35    
36     def __init__(self, *args):
37     pass
38    
39     def __call__(self, m):
40     """
41     short for getValue(m).
42     """
43     return self.getValue(m)
44    
45     def getValue(self, m):
46     """
47     returns the value of the mapping for m
48     """
49     raise NotImplementedError
50    
51     def getDerivative(self, m):
52     """
53     returns the value for the derivative of the mapping for m
54     """
55     raise NotImplementedError
56    
57     def getInverse(self, s):
58     """
59 gross 4074 returns the value of the inverse of the mapping for physical parameter p
60 caltinay 3946 """
61     raise NotImplementedError
62    
63 gross 4074 class LinearMapping(Mapping):
64     """
65     Maps a parameter by a linear transformation p = a *m + p0
66     """
67 caltinay 3946
68 gross 4074 def __init__(self, a=1, p0=0):
69     self.__a=a
70     self.__p0=p0
71     self.__a_inv = 1/a
72    
73     def getValue(self, m):
74     """
75     returns the value of the mapping for m
76     """
77     return self.__a*m + self.__p0
78    
79     def getDerivative(self, m):
80     """
81     returns the value for the derivative of the mapping for m
82     """
83     return self.__a
84    
85     def getInverse(self, p):
86     """
87     returns the value of the inverse of the mapping for s
88     """
89     return self.__a_inv * ( p - self.__p0)
90    
91     class DensityMapping(LinearMapping):
92     """
93     Density mapping with depth weighting
94    
95     *rho = rho0 + drho * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
96    
97     """
98     def __init__(self, domain, z0=None, rho0=None, drho=None, beta=None):
99     """
100     set up mapping
101    
102     :param domain: domain of the mapping
103     :param z0: depth weighting offset. If not present no depth scaling is applied.
104     :param drho: density scale. By default density of granite = 2750kg/m**3 is used.
105     :param rho0: reference density.
106     :param beta: depth weighting exponent.
107     :type domain: ``Domain``
108     :type z0: scalar
109     :type drho: scalar
110     :type rho0: scalar
111     :type beta: ``float``
112     """
113     if rho0 == None: rho0=0.
114     if drho == None: drho=2750*U.kg/U.m**3
115     if beta == None: beta=2.
116    
117     self.domain=domain
118     if z0:
119     DIM=self.domain.getDim()
120     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
121     a = drho * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
122     else:
123     a = drho
124     super(DensityMapping,self).__init__(a=a, p0=rho0)
125    
126     class SusceptibilityMapping(LinearMapping):
127     """
128     Susceptibility mapping with depth weighting
129    
130     *k = k0 + dk * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
131    
132     """
133     def __init__(self, domain, z0=None, k0=None, dk=None, beta=None):
134     """
135     set up mapping
136    
137     :param domain: domain of the mapping
138     :param z0: depth weighting offset. If not present no depth scaling is applied.
139     :param dk: susceptibility scale. Default 1..
140     :param k0: reference density. Default 0.
141     :param beta: depth weighting exponent.
142     :type domain: ``Domain``
143     :type z0: scalar
144     :type dk: scalar
145     :type k0: scalar
146     :type beta: ``float``
147     """
148     if k0==None: k0=0
149     if beta==None: beta=2.
150     if dk==None: dk=1
151     self.domain=domain
152     if z0:
153     DIM=self.domain.getDim()
154     l_z=boundingBoxEdgeLengths(domain)[DIM-1]
155     a = dk * ( (domain.getX()[DIM-1]-z0)/l_z )**(beta/2)
156     else:
157     a = dk
158     super(SusceptibilityMapping,self).__init__(a=a, p0=k0)
159    
160     # needs REVISION
161 caltinay 3946 class BoundedRangeMapping(Mapping):
162     """
163     Maps an unbounded parameter to a bounded range. The mapping is smooth and
164     continuous.
165     """
166    
167     def __init__(self, s_min=0, s_max=1):
168     if not s_min < s_max:
169     raise ValueError("value for s_min must be less than the value for s_max.")
170     self.s_min=s_min
171     self.s_max=s_max
172    
173     def getValue(self, m):
174     """
175     returns the value of the mapping for m
176     """
177     return (self.s_max+self.s_min)/2. + (self.s_max-self.s_min)/2. * tanh(m)
178    
179     def getDerivative(self, m):
180     """
181     returns the value for the derivative of the mapping for m
182     """
183     return ((self.s_max-self.s_min)/2.) * (1.-tanh(m)**2.)
184    
185     def getInverse(self, s):
186     """
187     returns the value of the inverse of the mapping for s
188     """
189     if not (inf(s) > self.s_min and sup(s) < self.s_max):
190     raise ValueError("s is out of range [%f,%f]"%(inf(s),sup(s)))
191    
192     return 1./2. * log( (s-self.s_min)/(self.s_max-s) )
193    
194    
195    
196    

  ViewVC Help
Powered by ViewVC 1.1.26