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

  ViewVC Help
Powered by ViewVC 1.1.26