/[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 4074 - (show 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
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 """Collection of parametrizations that map physical values to model parameters
17 and back"""
18
19 __copyright__="""Copyright (c) 2003-2012 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
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 parameters *p*
34 """
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 returns the value of the inverse of the mapping for physical parameter p
60 """
61 raise NotImplementedError
62
63 class LinearMapping(Mapping):
64 """
65 Maps a parameter by a linear transformation p = a *m + p0
66 """
67
68 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 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