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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3947 - (hide 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 caltinay 3946
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 caltinay 3947 __all__ = ['ForwardModel','GravityModel']
23    
24 caltinay 3946 from esys.escript import unitsSI as U
25     from esys.escript import *
26     from esys.escript.linearPDEs import LinearSinglePDE
27    
28 caltinay 3947 PI = 3.14159265358979323846
29 caltinay 3946 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