/[escript]/trunk/downunder/py_src/forwardmodels/gravity.py
ViewVC logotype

Contents of /trunk/downunder/py_src/forwardmodels/gravity.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 5707 - (show annotations)
Mon Jun 29 03:59:06 2015 UTC (3 years, 6 months ago) by sshaw
File MIME type: text/x-python
File size: 5447 byte(s)
adding copyright headers to files without copyright info, moved header to top of file in some cases where it wasn't
1 ##############################################################################
2 #
3 # Copyright (c) 2003-2015 by The University of Queensland
4 # http://www.uq.edu.au
5 #
6 # Primary Business: Queensland, Australia
7 # Licensed under the Open Software License version 3.0
8 # http://www.opensource.org/licenses/osl-3.0.php
9 #
10 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
11 # Development 2012-2013 by School of Earth Sciences
12 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
13 #
14 ##############################################################################
15
16 """Forward model for gravity (Bouguer) anomaly"""
17 from __future__ import division, print_function
18
19 __copyright__="""Copyright (c) 2003-2015 by The 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__ = ['GravityModel']
27
28 from .base import ForwardModelWithPotential
29 from esys.escript import unitsSI as U
30 from esys.escript.util import *
31 from math import pi as PI
32
33
34 class GravityModel(ForwardModelWithPotential):
35 """
36 Forward Model for gravity inversion as described in the inversion
37 cookbook.
38 """
39 def __init__(self, domain, w, g, gravity_constant=U.Gravitational_Constant,
40 coordinates=None, fixPotentialAtBottom=False, tol=1e-8):
41 """
42 Creates a new gravity model on the given domain with one or more
43 surveys (w, g).
44
45 :param domain: domain of the model
46 :type domain: `Domain`
47 :param w: data weighting factors
48 :type w: ``Vector`` or list of ``Vector``
49 :param g: gravity anomaly data
50 :type g: ``Vector`` or list of ``Vector``
51 :param coordinates: defines coordinate system to be used
52 :type coordinates: ReferenceSystem` or `SpatialCoordinateTransformation`
53 :param tol: tolerance of underlying PDE
54 :type tol: positive ``float``
55 :param fixPotentialAtBottom: if true potential is fixed to zero at the
56 base of the domain in addition to the top
57 :type fixPotentialAtBottom: ``bool``
58
59 :note: It is advisable to call rescaleWeights() to rescale weights
60 before starting the inversion.
61 """
62 super(GravityModel, self).__init__(domain, w, g, coordinates, fixPotentialAtBottom, tol)
63
64 trafo = self.getCoordinateTransformation()
65 if not trafo.isCartesian():
66 self.__G = 4*PI*gravity_constant * trafo.getVolumeFactor()
67 #* trafo.getReferenceSystem().getHeightUnit()**(-3)
68
69 fw = trafo.getScalingFactors()**2 * trafo.getVolumeFactor()
70 A=self.getPDE().createCoefficient("A")
71 for i in range(self.getDomain().getDim()): A[i,i]=fw[i]
72 self.getPDE().setValue(A=A)
73 else: # cartesian
74 self.__G = 4*PI*gravity_constant
75 self.getPDE().setValue(A=kronecker(self.getDomain()))
76
77 def rescaleWeights(self, scale=1., rho_scale=1.):
78 """
79 rescales the weights such that
80
81 *sum_s integrate( ( w_i[s] *g_i[s]) (w_j[s]*1/L_j) * L**2 * 4*pi*G*rho_scale )=scale*
82
83 :param scale: scale of data weighting factors
84 :type scale: positive ``float``
85 :param rho_scale: scale of density.
86 :type rho_scale: ``Scalar``
87 """
88 self._rescaleWeights(scale, self.__G*rho_scale)
89
90 def getArguments(self, rho):
91 """
92 Returns precomputed values shared by `getDefect()` and `getGradient()`.
93
94 :param rho: a suggestion for the density distribution
95 :type rho: ``Scalar``
96 :return: gravity potential and corresponding gravity field.
97 :rtype: ``Scalar``, ``Vector``
98 """
99 phi = self.getPotential(rho)
100 gravity_force = -grad(phi)
101 return phi, gravity_force
102
103 def getPotential(self, rho):
104 """
105 Calculates the gravity potential for a given density distribution.
106
107 :param rho: a suggestion for the density distribution
108 :type rho: ``Scalar``
109 :return: gravity potential
110 :rtype: ``Scalar``
111 """
112 pde = self.getPDE()
113 pde.resetRightHandSideCoefficients()
114 pde.setValue(Y=-self.__G*rho)
115 phi=pde.getSolution()
116
117 return phi
118
119 def getDefect(self, rho, phi, gravity_force):
120 """
121 Returns the value of the defect
122
123 :param rho: density distribution
124 :type rho: ``Scalar``
125 :param phi: corresponding potential
126 :type phi: ``Scalar``
127 :param gravity_force: gravity force
128 :type gravity_force: ``Vector``
129 :rtype: ``float``
130 """
131 return self._getDefect(gravity_force)
132
133 def getGradient(self, rho, phi, gravity_force):
134 """
135 Returns the gradient of the defect with respect to density.
136
137 :param rho: density distribution
138 :type rho: ``Scalar``
139 :param phi: corresponding potential
140 :type phi: ``Scalar``
141 :param gravity_force: gravity force
142 :type gravity_force: ``Vector``
143 :rtype: ``Scalar``
144 """
145 pde=self.getPDE()
146 pde.resetRightHandSideCoefficients()
147 pde.setValue(X=self.getDefectGradient(gravity_force))
148 ZT=pde.getSolution()
149 return ZT*(-self.__G)
150

  ViewVC Help
Powered by ViewVC 1.1.26