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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4012 - (show annotations)
Wed Oct 3 02:25:10 2012 UTC (7 years, 1 month ago) by caltinay
File MIME type: text/x-python
File size: 4516 byte(s)
Initialize X_reduced in forward model to avoid issues later on.

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 forward models that define the inversion problem"""
17
18 __copyright__="""Copyright (c) 2003-2012 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 __all__ = ['ForwardModel','GravityModel']
26
27 from esys.escript import unitsSI as U
28 from esys.escript.linearPDEs import LinearSinglePDE
29 from esys.escript.util import *
30 from esys.escript import Data
31
32 PI = 3.14159265358979323846
33 G = 6.6742e-11*U.m**3/(U.kg*U.sec**2)
34
35 class ForwardModel(object):
36 """
37 An abstract forward model that can be plugged into a cost function.
38 Subclasses need to implement `getValue()`, `getGradient()`, and possibly
39 `getArguments()`.
40 """
41 def __init__(self):
42 pass
43
44 def getArguments(self, x):
45 return None
46
47 def getValue(self, x, *args):
48 pass
49
50 def getGradient(self, x, *args):
51 pass
52
53 class GravityModel(ForwardModel):
54 """
55 Forward Model for gravity inversion as described in the inversion cookbook.
56 """
57 def __init__(self, domain, chi, g, constrain_top=True, gravity_constant=G, tol=1e-8):
58 """
59 Creates a new gravity model on the given domain with one or more
60 surveys (chi, g).
61 """
62 self.__domain = domain
63 try:
64 n=len(chi)
65 m=len(g)
66 if m != n:
67 raise ValueError("Length of chi and g must be the same.")
68 self.__chi = chi
69 self.__g = g
70 except TypeError:
71 self.__chi = [chi]
72 self.__g = [g]
73
74 A=0
75 for s in xrange(len(self.__chi)):
76 A2 = integrate(inner(self.__chi[s], self.__g[s]**2))
77 if A2 < 0:
78 raise ValueError("Negative weighting factor for survey %s"%s)
79 A=max(A2, A)
80 if not A > 0:
81 raise ValueError("No reference data set.")
82 self.__G = gravity_constant
83 BX = boundingBox(domain)
84 DIM = domain.getDim()
85 x = domain.getX()
86 self.__pde=LinearSinglePDE(domain)
87 self.__pde.getSolverOptions().setTolerance(tol)
88 if constrain_top:
89 constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])
90 else:
91 constraint=whereZero(x[DIM-1]-BX[DIM-1][0])
92 for i in xrange(DIM-1):
93 constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])
94 self.__pde.setValue(A=kronecker(DIM), q=constraint)
95
96 def getSurvey(self, index=None):
97 """
98 Returns the pair (g_index, chi_index), where g_i is the gravity
99 anomaly of survey i, chi_i is the weighting factor for survey i.
100 If index is None, all surveys will be returned in a pair of lists.
101 """
102 if index is None:
103 return self.__g, self.__chi
104 if index>=len(self.__g):
105 raise IndexError("Forward model only has %d surveys"%len(self.__g))
106 return self.__g[index], self.__chi[index]
107
108 def getArguments(self, rho):
109 """
110 Returns precomputed values shared by getValue() and getGradient().
111 """
112 phi = self.getPotential(rho)
113 return phi, -grad(phi)
114
115 def getPotential(self, rho):
116 self.__pde.setValue(Y=(-4*PI*self.__G) * rho, X=Data(), X_reduced=Data())
117 phi=self.__pde.getSolution()
118 return phi
119
120 def getValue(self, rho, phi, gravity_force):
121 A=0.
122 for s in xrange(len(self.__chi)):
123 A = inner(self.__chi[s], (gravity_force-self.__g[s])**2) + A
124 return 0.5*integrate(A)
125
126 def getGradient(self, rho, phi, gravity_force):
127 Z=0.
128 for s in xrange(len(self.__chi)):
129 Z = self.__chi[s] * (-gravity_force+self.__g[s]) + Z
130
131 self.__pde.setValue(Y=Data(), X=Z)
132 ZT=self.__pde.getSolution()
133 return ZT*(-4*PI*self.__G)
134

  ViewVC Help
Powered by ViewVC 1.1.26