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

Contents of /trunk/downunder/py_src/forwardmodels/subsidence.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: 4548 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 Subsidence modelling"""
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__ = ['Subsidence']
27
28 from .base import ForwardModel
29 from esys.escript import Data, FunctionOnBoundary
30 from esys.escript.linearPDEs import LinearPDESystem
31 from esys.escript.util import *
32
33
34 class Subsidence(ForwardModel):
35 """
36 Forward Model for subsidence inversion minimizing
37 integrate( (inner(w,u)-d)**2)
38 where u is the surface displacement due to a pressure change P
39 """
40 def __init__(self, domain, w, d, lam, mu, coordinates=None, tol=1e-8):
41 """
42 Creates a new subsidence on the given domain
43
44 :param domain: domain of the model
45 :type domain: `Domain`
46 :param w: data weighting factors and direction
47 :type w: ``Vector`` with ``FunctionOnBoundary``
48 :param d: displacement measured at surface
49 :type d: ``Scalar`` with ``FunctionOnBoundary``
50 :param lam: 1st Lame coefficient
51 :type lam: ``Scalar`` with ``Function``
52 :param lam: 2st Lame coefficient/Shear modulus
53 :type lam: ``Scalar`` with ``Function``
54 :param coordinates: defines coordinate system to be used (not supported yet))
55 :type coordinates: `ReferenceSystem` or `SpatialCoordinateTransformation`
56 :param tol: tolerance of underlying PDE
57 :type tol: positive ``float``
58 """
59 super(Subsidence, self).__init__()
60 DIM=domain.getDim()
61
62 self.__pde=LinearPDESystem(domain)
63 self.__pde.setSymmetryOn()
64 self.__pde.getSolverOptions().setTolerance(tol)
65 #... set coefficients ...
66 C=self.__pde.createCoefficient('A')
67 for i in range(DIM):
68 for j in range(DIM):
69 C[i,i,j,j]+=lam
70 C[i,j,i,j]+=mu
71 C[i,j,j,i]+=mu
72 x=domain.getX()
73 msk=whereZero(x[DIM-1])*kronecker(DIM)[DIM-1]
74 for i in range(DIM-1):
75 xi=x[i]
76 msk+=(whereZero(xi-inf(xi))+whereZero(xi-sup(xi))) *kronecker(DIM)[i]
77 self.__pde.setValue(A=C,q=msk)
78
79 self.__w=interpolate(w, FunctionOnBoundary(domain))
80 self.__d=interpolate(d, FunctionOnBoundary(domain))
81
82 def rescaleWeights(self, scale=1., P_scale=1.):
83 """
84 rescales the weights
85
86 :param scale: scale of data weighting factors
87 :type scale: positive ``float``
88 :param P_scale: scale of pressure increment
89 :type P_scale: ``Scalar``
90 """
91 pass
92
93 def getArguments(self, P):
94 """
95 Returns precomputed values shared by `getDefect()` and `getGradient()`.
96
97 :param P: pressure
98 :type P: ``Scalar``
99 :return: displacement u
100 :rtype: ``Vector``
101 """
102 DIM=self.__pde.getDim()
103 self.__pde.setValue(y=Data(),X=P*kronecker(DIM))
104 u= self.__pde.getSolution()
105 return u,
106
107 def getDefect(self, P,u):
108 """
109 Returns the value of the defect.
110
111 :param P: pressure
112 :type P: ``Scalar``
113 :param u: corresponding displacement
114 :type u: ``Vector``
115 :rtype: ``float``
116 """
117 return 0.5*integrate((inner(u,self.__w)-self.__d)**2)
118
119 def getGradient(self, P, u):
120 """
121 Returns the gradient of the defect with respect to susceptibility.
122
123 :param P: pressure
124 :type P: ``Scalar``
125 :param u: corresponding displacement
126 :type u: ``Vector``
127 :rtype: ``Scalar``
128 """
129 d=inner(u,self.__w)-self.__d
130 self.__pde.setValue(y=d*self.__w,X=Data())
131 ustar=self.__pde.getSolution()
132
133 return div(ustar)
134

  ViewVC Help
Powered by ViewVC 1.1.26