/[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 3947 - (show annotations)
Wed Aug 22 23:19:10 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/x-python
File size: 3783 byte(s)
Compiling and installing downunder module now. Adjusted import statements
accordingly. Added a gravity test run.

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

  ViewVC Help
Powered by ViewVC 1.1.26