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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4063 by caltinay, Fri Nov 9 06:15:00 2012 UTC revision 4074 by gross, Thu Nov 15 03:30:59 2012 UTC
# Line 30  from esys.escript.util import * Line 30  from esys.escript.util import *
30  from esys.escript import Data, Vector, Scalar, Function  from esys.escript import Data, Vector, Scalar, Function
31  from math import pi as PI  from math import pi as PI
32    
33    
34  class ForwardModel(object):  class ForwardModel(object):
35      """      """
36      An abstract forward model that can be plugged into a cost function.      An abstract forward model that can be plugged into a cost function.
# Line 47  class ForwardModel(object): Line 48  class ForwardModel(object):
48    
49      def getGradient(self, x, *args):      def getGradient(self, x, *args):
50          raise NotImplementedError          raise NotImplementedError
51        def getDirectionalDerivative(self, x, d, *args):
52            """
53            return directional derivatiave
54            """
55            Y=self.getGradient(x, *args)
56            return integrate(p * Y)
57            
58  class ForwardModelWithPotential(ForwardModel):  class ForwardModelWithPotential(ForwardModel):
59      """      """
60      Base class for a forward model using a potential such as magnetic or      Base class for a forward model using a potential such as magnetic or
# Line 60  class ForwardModelWithPotential(ForwardM Line 67  class ForwardModelWithPotential(ForwardM
67      It is assumed that the forward model is produced through postprocessing      It is assumed that the forward model is produced through postprocessing
68      of the solution of a potential PDE.      of the solution of a potential PDE.
69      """      """
70      def __init__(self, domain, weight, data, fix_all_faces=False, tol=1e-8):      def __init__(self, domain, weight, data,  tol=1e-8):
71          """          """
72          initialization.          initialization.
73    
# Line 70  class ForwardModelWithPotential(ForwardM Line 77  class ForwardModelWithPotential(ForwardM
77          :type weight: `Vector` or list of `Vector`          :type weight: `Vector` or list of `Vector`
78          :param data: data          :param data: data
79          :type data: `Vector` or list of `Vector`          :type data: `Vector` or list of `Vector`
         :param fix_all_faces: if ``true`` all faces of the domain are fixed  
                               otherwise only the top surface  
         :type fix_all_faces: ``bool``  
80          :param tol: tolerance of underlying PDE          :param tol: tolerance of underlying PDE
81          :type tol: positive ``float``          :type tol: positive ``float``
82    
# Line 93  class ForwardModelWithPotential(ForwardM Line 97  class ForwardModelWithPotential(ForwardM
97    
98          A=0          A=0
99          for s in xrange(len(self.__weight)):          for s in xrange(len(self.__weight)):
100              A2 = integrate(inner(self.__weight[s], self.__data[s]**2))              A += integrate(inner(self.__weight[s], self.__data[s]**2))      
101              if A2 < 0:          if A > 0:
102                  raise ValueError("Negative weighting factor for survey %s"%s)          for s in xrange(len(self.__weight)):  self.__weight[s]*=1./A
103              A=max(A2, A)      else:
104          if not A > 0:              raise ValueError("No non-zero data contribution found.")
105              raise ValueError("No reference data set.")              
106    
107          BX = boundingBox(domain)          BX = boundingBox(domain)
108          DIM = domain.getDim()          DIM = domain.getDim()
# Line 106  class ForwardModelWithPotential(ForwardM Line 110  class ForwardModelWithPotential(ForwardM
110          self.__pde=LinearSinglePDE(domain)          self.__pde=LinearSinglePDE(domain)
111          self.__pde.getSolverOptions().setTolerance(tol)          self.__pde.getSolverOptions().setTolerance(tol)
112          self.__pde.setSymmetryOn()          self.__pde.setSymmetryOn()
113            self.__pde.setValue(q=whereZero(x[DIM-1]-BX[DIM-1][1]))
         if fix_all_faces:  
             constraint=whereZero(x[DIM-1]-BX[DIM-1][1])+whereZero(x[DIM-1]-BX[DIM-1][0])  
             for i in xrange(DIM-1):  
                 constraint=constraint+whereZero(x[i]-BX[i][1])+whereZero(x[i]-BX[i][0])  
         else:  
             constraint=whereZero(x[DIM-1]-BX[DIM-1][1])  
         self.__pde.setValue(q=constraint)  
114    
115      def getDomain(self):      def getDomain(self):
116          """          """
# Line 155  class ForwardModelWithPotential(ForwardM Line 152  class ForwardModelWithPotential(ForwardM
152          return self.__data[index], self.__weight[index]          return self.__data[index], self.__weight[index]
153    
154    
155    
156  class GravityModel(ForwardModelWithPotential):  class GravityModel(ForwardModelWithPotential):
157      """      """
158      Forward Model for gravity inversion as described in the inversion      Forward Model for gravity inversion as described in the inversion
159      cookbook.      cookbook.
160      """      """
161      def __init__(self, domain, chi, g, fix_all_faces=False, gravity_constant=U.Gravitational_Constant, tol=1e-8):      def __init__(self, domain, chi, g,  gravity_constant=U.Gravitational_Constant, tol=1e-8):
162          """          """
163          Creates a new gravity model on the given domain with one or more          Creates a new gravity model on the given domain with one or more
164          surveys (chi, g).          surveys (chi, g).
# Line 177  class GravityModel(ForwardModelWithPoten Line 175  class GravityModel(ForwardModelWithPoten
175          :param tol: tolerance of underlying PDE          :param tol: tolerance of underlying PDE
176          :type tol: positive ``float``          :type tol: positive ``float``
177          """          """
178          super(GravityModel, self).__init__(domain, chi, g, fix_all_faces, tol)          super(GravityModel, self).__init__(domain, chi, g, tol)
179    
180          self.__G = gravity_constant          self.__G = gravity_constant
181          self.getPDE().setValue(A=kronecker(self.getDomain()))          self.getPDE().setValue(A=kronecker(self.getDomain()))
# Line 256  class MagneticModel(ForwardModelWithPote Line 254  class MagneticModel(ForwardModelWithPote
254      Forward Model for magnetic inversion as described in the inversion      Forward Model for magnetic inversion as described in the inversion
255      cookbook.      cookbook.
256      """      """
257      def __init__(self, domain, chi, B, background_field, fix_all_faces=False, tol=1e-8):      def __init__(self, domain, chi, B, background_field,  tol=1e-8):
258          """          """
259          Creates a new magnetic model on the given domain with one or more          Creates a new magnetic model on the given domain with one or more
260          surveys (chi, B).          surveys (chi, B).
# Line 273  class MagneticModel(ForwardModelWithPote Line 271  class MagneticModel(ForwardModelWithPote
271          :param tol: tolerance of underlying PDE          :param tol: tolerance of underlying PDE
272          :type tol: positive ``float``          :type tol: positive ``float``
273          """          """
274          super(MagneticModel, self).__init__(domain, chi, B, fix_all_faces, tol)          super(MagneticModel, self).__init__(domain, chi, B, tol)
275          self.__background_field=interpolate(background_field, Function(self.getDomain()))          self.__background_field=interpolate(background_field, Function(self.getDomain()))
276          self.getPDE().setValue(A=kronecker(self.getDomain()))          self.getPDE().setValue(A=kronecker(self.getDomain()))
277    
# Line 287  class MagneticModel(ForwardModelWithPote Line 285  class MagneticModel(ForwardModelWithPote
285          :rtype: ``Scalar``, ``Vector``          :rtype: ``Scalar``, ``Vector``
286          """          """
287          phi = self.getPotential(k)          phi = self.getPotential(k)
288          magnetic_field = (1+k) * self.__background_field -grad(phi)          magnetic_field = k * self.__background_field -grad(phi)
289          return phi, magnetic_field          return phi, magnetic_field
290    
291      def getPotential(self, k):      def getPotential(self, k):
# Line 302  class MagneticModel(ForwardModelWithPote Line 300  class MagneticModel(ForwardModelWithPote
300          pde=self.getPDE()          pde=self.getPDE()
301    
302          pde.resetRightHandSideCoefficients()          pde.resetRightHandSideCoefficients()
303          pde.setValue(X = (1+k)* self.__background_field)          pde.setValue(X = k* self.__background_field)
304          phi=pde.getSolution()          phi=pde.getSolution()
305    
306          return phi          return phi
# Line 342  class MagneticModel(ForwardModelWithPote Line 340  class MagneticModel(ForwardModelWithPote
340    
341          pde.resetRightHandSideCoefficients()          pde.resetRightHandSideCoefficients()
342          pde.setValue(X=Y)          pde.setValue(X=Y)
343            
344            print "test VTK file generated"
345            from esys.weipa import saveVTK
346            saveVTK("Y.vtu",RHS=pde.getRightHandSide(), Y=Y)
347            
348          YT=pde.getSolution()          YT=pde.getSolution()
349          return inner(Y-grad(YT),self.__background_field)          return inner(Y-grad(YT),self.__background_field)
350    
351            

Legend:
Removed from v.4063  
changed lines
  Added in v.4074

  ViewVC Help
Powered by ViewVC 1.1.26