Working on MT problem and survey. Fixed bugs and completed example

2026-07-01 13:11:59 +08:00 · 2015-02-17 00:09:24 -08:00
parent fb717b5f31
commit 9b650041d1
4 changed files with 126 additions and 47 deletions
@@ -1,9 +1,36 @@
 # Test script to use simpegMT platform to forward model synthetic data.

 # Import 
-import simpegMT as simpegMT, SimPEG as simpeg 
-
+import simpegMT as simpegmt, SimPEG as simpeg 
+import numpy as np

+# Make a mesh
+M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
+# Setup the model
+conds = [1e-2,1]
+sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds)
+sig[M.gridCC[:,2]>0] = 1e-8
+sig[M.gridCC[:,2]<-600] = 1e-1
+sigBG = np.zeros(M.nC) + conds[0]
+sigBG[M.gridCC[:,2]>0] = 1e-8

+## Setup the the survey object
+# Receiver locations
+rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50))
+rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1))))
 # Make a receiver list
-for 
+rxList = []
+for loc in rx_loc:
+    for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
+        rxList.append(simpegmt.SurveyMT.RxMT(loc,rxType))
+# Source list
+srcList =[]
+for freq in np.logspace(3,-1,21):
+    srcList.append(simpegmt.SurveyMT.srcMT(freq,rxList))
+# Survey MT 
+survey = simpegmt.SurveyMT.SurveyMT(srcList)
+
+## Setup the problem object
+problem = simpegmt.ProblemMT.MTProblem(M)
+problem.pair(survey)
+
@@ -1,4 +1,4 @@
-from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver
+from SimPEG import Survey, Problem, Utils, Models, np, sp, Solver as SimpegSolver
 from scipy.constants import mu_0
 from SurveyMT import SurveyMT, FieldsMT

@@ -47,15 +47,6 @@ class MTProblem(Problem.BaseProblem):
            self._MeSigma = self.mesh.getEdgeInnerProduct(sigma)
        return self._MeSigma

-    # TODO:
-    # MeSigmaBG
-    @property
-    def MeSigmaBG(self):
-        #TODO: hardcoded to sigma as the model
-        if getattr(self, '_MeSigma', None) is None:
-            sigmaBG = self.curTModelBG
-            self._MeSigmaBG = self.mesh.getEdgeInnerProduct(sigmaBG)
-        return self._MeSigmaBG

    @property
    def MeSigmaI(self):
@@ -70,38 +61,68 @@ class MTProblem(Problem.BaseProblem):
    @property
    def curTModel(self):
        if getattr(self, '_curTModel', None) is None:
-            self._curTModel = self.mapping.transform(self.curModel)
+            self._curTModel = self.mapping*self.curModel
        return self._curTModel

    @property
    def curTModelDeriv(self):
        if getattr(self, '_curTModelDeriv', None) is None:
-            self._curTModelDeriv = self.mapping.transformDeriv(self.curModel)
+            self._curTModelDeriv = self.mapping*self.curModel
        return self._curTModelDeriv

-    def fields(self, m):
+
+    # Background model
+    @property
+    def backModel(self):
+        """
+            Sets the model, and removes dependent mass matrices.
+        """
+        return getattr(self, '_backModel', None)
+
+    @backModel.setter
+    def backModel(self, value):
+        if value is self.backModel:
+            return # it is the same!
+        self._backModel = Models.Model(value, self.mapping)
+        for prop in self.deleteTheseOnModelUpdate:
+            if hasattr(self, prop):
+                delattr(self, prop)
+
+    @property
+    def MeSigmaBG(self):
+        #TODO: hardcoded to sigma as the model
+        if getattr(self, '_MeSigmaBG', None) is None:
+            sigmaBG = self.backModel
+            self._MeSigmaBG = self.mesh.getEdgeInnerProduct(sigmaBG)
+        return self._MeSigmaBG
+
+    def fields(self, m, m_back):
        '''
        Function to calculate all the fields for the model m.

-        :param np.ndarray (nC,) m: f
+        :param np.ndarray (nC,) m: Conductivity model
+        :param np.ndarray (nC,) m_back: Background conductivity model
        '''
        self.curModel = m
-        RHS, CalcFields = self.getRHS(freq), self.calcFields
+        self.backModel = m_back
+        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields

        F = FieldsMT(self.mesh, self.survey)

        for freq in self.survey.freqs:
            A = self.getA(freq)
-            rhs = RHS(freq)
+            rhs = self.getRHS(freq,m_back)
            Ainv = self.Solver(A, **self.solverOpts)
            e = Ainv * rhs 

            # Store the fields
            Src = self.survey.getSources(freq)
-            # Stroe the fields
-            F[Src, 'e'] = e
-            F[Src, 'b'] = self.mesh.edgeCurl * e 
-
+            # Store the fields
+            F[Src, 'e_px'] = e[:,0]
+            F[Src, 'e_py'] = e[:,1]
+            b = self.mesh.edgeCurl * e     
+            F[Src, 'b_px'] = b[:,0]
+            F[Src, 'b_py'] = b[:,1]
        return F


@@ -157,8 +178,8 @@ class MTProblem(Problem.BaseProblem):
        # assert len()
        # Get the background electric fields
        from simpegMT.Sources import homo1DModelSource
-        eBG_bp = home1DModelSource(self.mesh,freq,backSigma)
-        Abg = getAbg(freq)
+        eBG_bp = homo1DModelSource(self.mesh,freq,backSigma)
+        Abg = self.getAbg(freq)

        return Abg*eBG_bp

@@ -1,36 +1,40 @@
-def homo1DModelSource(mesh,freq,bgMod):
+import SimPEG as simpeg, numpy as np
+from simpegMT.Utils import get1DEfields
+
+def homo1DModelSource(mesh,freq,m_back):
    '''
        Function that calculates and return background fields

        :param Simpeg mesh object mesh: Holds information on the discretization
        :param float freq: The frequency to solve at
-        :param np.array bgMod: Background model to base the calculations on.
+        :param np.array m_back: Background model of conductivity to base the calculations on.
        :rtype: numpy.ndarray (mesh.nE,2)
        :return: eBG_bp, E fields for the background model at both polarizations.

    '''

    # import
+    import SimPEG as simpeg
    from simpegMT.Utils import get1DEfields
    # Get a 1d solution for a halfspace background
    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
-    e0_1d = get1DEfields(mesh1d,mesh.r(sigBG,'CC','CC','M')[0,0,:],freq)
+    e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
    # Setup x (east) polarization (_x)
    ex_px = np.zeros(mesh.vnEx,dtype=complex)
    ey_px = np.zeros((mesh.nEy,1),dtype=complex)
    ez_px = np.zeros((mesh.nEz,1),dtype=complex)
    # Assign the source to ex_x
-    for i in arange(mesh.vnEx[0]):
-        for j in arange(mesh.vnEx[1]):
+    for i in np.arange(mesh.vnEx[0]):
+        for j in np.arange(mesh.vnEx[1]):
            ex_px[i,j,:] = e0_1d
    eBG_px = np.vstack((simpeg.Utils.mkvc(mesh.r(ex_px,'Ex','Ex','V'),2),ey_px,ez_px))
    # Setup y (north) polarization (_py)
-    ex_py = np.zeros(mesh.nEx, dtype='complex128')
-    ey_py = np.zeros((mesh.vnEy), dtype='complex128')
-    ez_py = np.zeros(mesh.nEz, dtype='complex128')
+    ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
+    ey_py = np.zeros(mesh.vnEy, dtype='complex128')
+    ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
    # Assign the source to ey_py
-    for i in arange(mesh.vnEy[0]):
-        for j in arange(mesh.vnEy[1]):
+    for i in np.arange(mesh.vnEy[0]):
+        for j in np.arange(mesh.vnEy[1]):
            ey_py[i,j,:] = e0_1d 
    eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))

@@ -1,4 +1,5 @@
 from SimPEG import Survey, Utils, Problem, np, sp
+from scipy.constants import mu_0

 class RxMT(Survey.BaseRx):

@@ -65,14 +66,37 @@ class RxMT(Survey.BaseRx):
        '''
        Project the fields and return the 
        '''
-        # Get the numerator information
-        P_num = self.getP(mesh,self.projGLoc('numerator'))
-        u_num_complex = u[src, self.projField('numerator')]
-        # Get the denominator information
-        P_den = self.getP(mesh,self.projGLoc('denominator'))
-        u_den_complex = u[src, self.projField('denominator')]
-        # Calculate the fraction
-        f_part_complex = (P_num*u_num_complex)/(P_den*u_den_complex)
+
+        # Get the projection
+        Pex = self.getP(mesh,'Ex')
+        Pey = self.getP(mesh,'Ey')
+        Pbx = self.getP(mesh,'Fx')
+        Pby = self.getP(mesh,'Fy')
+        # Get the fields at location
+        ex_px = Pex*u[src,'e_px']
+        ey_px = Pey*u[src,'e_px']
+        ex_py = Pex*u[src,'e_py']
+        ey_py = Pey*u[src,'e_py']
+        hx_px = Pbx*u[src,'b_px']/mu_0
+        hy_px = Pby*u[src,'b_px']/mu_0
+        hx_py = Pbx*u[src,'b_py']/mu_0
+        hy_py = Pby*u[src,'b_py']/mu_0
+        if 'zxx' in self.rxType:
+            f_part_complex = (ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
+        elif 'zxy' in self.rxType:
+            f_part_complex  = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
+        elif 'zyx' in self.rxType:
+            f_part_complex  = (ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
+        elif 'zyy' in self.rxType:
+            f_part_complex  = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
+
+        # P_num = self.getP(mesh,self.projGLoc('numerator'))
+        # u_num_complex = u[src, self.projField('numerator')]
+        # # Get the denominator information
+        # P_den = self.getP(mesh,self.projGLoc('denominator'))
+        # u_den_complex = u[src, self.projField('denominator')]
+        # # Calculate the fraction
+        # f_part_complex = (P_num*u_num_complex)/(P_den*u_den_complex)
        # get the real or imag component
        real_or_imag = self.projComp
        f_part = getattr(f_part_complex, real_or_imag)
@@ -100,6 +124,7 @@ class RxMT(Survey.BaseRx):


 # Call this Source or polarization or something...?
+# Note: Might need to add tests to make sure that both polarization have the same rxList. 
 class srcMT(Survey.BaseTx):
    '''
    Sources for the MT problem. 
@@ -107,24 +132,25 @@ class srcMT(Survey.BaseTx):

    :param float freq: The frequency of the source
    :param list rxList: A list of receivers associated with the source
+    :param str srcPol: The polarization of the source
    '''

    freq = None #: Frequency (float)

    rxPair = RxMT

-    knownTxTypes = ['ORTPOL'] # ORThogonal POLarization
+    knownTxTypes = ['pol_xy','pol_x','pol_y'] # ORThogonal POLarization

-    def __init__(self, freq, rxList): # remove txType? hardcode to one thing. always polarizations
+    def __init__(self, freq, rxList, srcPol = 'pol_xy'): # remove txType? hardcode to one thing. always polarizations
        self.freq = float(freq)
-        Survey.BaseTx.__init__(self, None, 'ORTPOL', rxList)
+        Survey.BaseTx.__init__(self, None, srcPol, rxList)
        # Survey.BaseTx.__init__(self, loc, 'polarization', rxList)



 class FieldsMT(Problem.Fields):
    """Fancy Field Storage for a MT survey."""
-    knownFields = {'b': 'F', 'e': 'E'}
+    knownFields = {'b_px': 'F','b_py': 'F', 'e_px': 'E','e_py': 'E'}
    dtype = complex


@@ -141,6 +167,7 @@ class SurveyMT(Survey.BaseSurvey):
    def __init__(self, srcList, **kwargs):
        # Sort these by frequency
        self.srcList = srcList
+        self.txList = self.srcList # Hack - make txList index srcList, since it is used in the backend. 
        Survey.BaseSurvey.__init__(self, **kwargs)

        _freqDict = {}