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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4344 - (show annotations)
Wed Mar 27 07:58:34 2013 UTC (6 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 6217 byte(s)
fix potential at botttom added.
1
2 ##############################################################################
3 #
4 # Copyright (c) 2003-2013 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 """Collection of parametrizations that map physical values to model parameters
17 and back"""
18
19 __copyright__="""Copyright (c) 2003-2013 by University of Queensland
20 http://www.uq.edu.au
21 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 __all__ = ['Mapping', 'DensityMapping', 'SusceptibilityMapping', 'BoundedRangeMapping', 'LinearMapping']
27
28 from esys.escript import inf, sup, log, tanh, boundingBoxEdgeLengths, clip
29 import esys.escript.unitsSI as U
30
31 class Mapping(object):
32 """
33 An abstract mapping class to map level set functions *m* to physical
34 parameters *p*.
35 """
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 returns the value of the inverse of the mapping for physical parameter p
61 """
62 raise NotImplementedError
63
64 def getTypicalDerivative(self):
65 """
66 returns a typical value for the derivative
67 """
68 raise NotImplementedError
69
70
71 class LinearMapping(Mapping):
72 """
73 Maps a parameter by a linear transformation p = a * m + p0
74 """
75
76 def __init__(self, a=1, p0=0):
77 self.__a=a
78 self.__p0=p0
79 self.__a_inv = 1/a
80
81 def getValue(self, m):
82 """
83 returns the value of the mapping for m
84 """
85 return self.__a*m + self.__p0
86
87 def getDerivative(self, m):
88 """
89 returns the value for the derivative of the mapping for m
90 """
91 return self.__a
92
93 def getInverse(self, p):
94 """
95 returns the value of the inverse of the mapping for s
96 """
97 return self.__a_inv * ( p - self.__p0)
98
99 def getTypicalDerivative(self):
100 """
101 returns a typical value for the derivative
102 """
103 return self.__a
104
105 class DensityMapping(LinearMapping):
106 """
107 Density mapping with depth weighting
108
109 *rho = rho0 + drho * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
110
111 """
112 def __init__(self, domain, z0=None, rho0=None, drho=None, beta=None):
113 """
114 initializes the mapping
115
116 :param domain: domain of the mapping
117 :type domain: ``Domain``
118 :param z0: depth weighting offset. If not present no depth scaling is applied.
119 :type z0: scalar
120 :param rho0: reference density, defaults to 0
121 :type rho0: scalar
122 :param drho: density scale. By default density of granite = 2750kg/m**3 is used.
123 :type drho: scalar
124 :param beta: depth weighting exponent, defaults to 2
125 :type beta: ``float``
126 """
127 if rho0 is None: rho0=0.
128 if drho is None: drho=2750*U.kg/U.m**3
129 if beta is None: beta=2.
130
131 self.domain=domain
132 if z0 is not None:
133 DIM=self.domain.getDim()
134 l_z=boundingBoxEdgeLengths(domain)[DIM-1]
135 a = drho * ( clip(z0-domain.getX()[DIM-1], minval=0)/l_z )**(beta/2)
136 else:
137 a = drho
138 super(DensityMapping,self).__init__(a=a, p0=rho0)
139
140 class SusceptibilityMapping(LinearMapping):
141 """
142 Susceptibility mapping with depth weighting
143
144 *k = k0 + dk * ( (x_2 - z0)/l_z)^(beta/2) ) * m*
145
146 """
147 def __init__(self, domain, z0=None, k0=None, dk=None, beta=None):
148 """
149 set up mapping
150
151 :param domain: domain of the mapping
152 :type domain: ``Domain``
153 :param z0: depth weighting offset. If not present no depth scaling is applied.
154 :type z0: scalar
155 :param k0: reference density, defaults to 0
156 :type k0: scalar
157 :param dk: susceptibility scale, defaults to 1
158 :type dk: scalar
159 :param beta: depth weighting exponent, defaults to 2
160 :type beta: ``float``
161 """
162 if k0 is None: k0=0.
163 if beta is None: beta=2.
164 if dk is None: dk=1.
165 self.domain=domain
166
167 if z0 is not None:
168 DIM=self.domain.getDim()
169 l_z=boundingBoxEdgeLengths(domain)[DIM-1]
170 a = dk * ( clip(z0-domain.getX()[DIM-1] , minval=0)/l_z )**(beta/2)
171 else:
172 a = dk
173 super(SusceptibilityMapping,self).__init__(a=a, p0=k0)
174
175 # needs REVISION
176 class BoundedRangeMapping(Mapping):
177 """
178 Maps an unbounded parameter to a bounded range. The mapping is smooth and
179 continuous.
180 """
181
182 def __init__(self, s_min=0, s_max=1):
183 if not s_min < s_max:
184 raise ValueError("value for s_min must be less than the value for s_max.")
185 self.s_min=s_min
186 self.s_max=s_max
187
188 def getValue(self, m):
189 """
190 returns the value of the mapping for m
191 """
192 return (self.s_max+self.s_min)/2. + (self.s_max-self.s_min)/2. * tanh(m)
193
194 def getDerivative(self, m):
195 """
196 returns the value for the derivative of the mapping for m
197 """
198 return ((self.s_max-self.s_min)/2.) * (1.-tanh(m)**2.)
199
200 def getInverse(self, s):
201 """
202 returns the value of the inverse of the mapping for s
203 """
204 if not (inf(s) > self.s_min and sup(s) < self.s_max):
205 raise ValueError("s is out of range [%f,%f]"%(inf(s),sup(s)))
206
207 return 1./2. * log( (s-self.s_min)/(self.s_max-s) )
208

  ViewVC Help
Powered by ViewVC 1.1.26