Moving Problem1D/3D folders to NSEM file

Looking at making Jvec more general.
2026-07-04 01:58:52 +08:00 · 2016-05-10 16:07:45 -07:00
parent 1ab67b5790
commit 6165e619ed
8 changed files with 559 additions and 441 deletions
@@ -4,6 +4,7 @@ import sys
 from numpy.lib import recfunctions as recFunc
 from SimPEG.EM.Utils import omega

+
 ##############
 ### Fields ###
 ##############
@@ -12,8 +13,10 @@ class BaseNSEMFields(Problem.Fields):
    knownFields = {}
    dtype = complex

-
-class Fields1D_e(BaseNSEMFields):
+###########
+# 1D Fields
+###########
+class Fields1D_ePrimSec(BaseNSEMFields):
    """
    Fields storage for the 1D NSEM solution.
    """
@@ -44,6 +47,118 @@ class Fields1D_e(BaseNSEMFields):
    def _e(self, eSolution, srcList):
        return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)

+    def _eDeriv_u(self, src, du_dm_v, adjoint = False):
+
+
+        return Utils.Identity()*du_dm_v
+
+    def _eDeriv_m(self, src, v, adjoint = False):
+        # assuming primary does not depend on the model
+        return Utils.Zero()
+
+    def _bPrimary(self, eSolution, srcList):
+        bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
+        for i, src in enumerate(srcList):
+            bp = src.bPrimary(self.survey.prob)
+            if bp is not None:
+                bPrimary[:,i] += bp[:,-1]
+        return bPrimary
+
+    def _bSecondary(self, eSolution, srcList):
+        C = self.mesh.nodalGrad
+        b = (C * eSolution)
+        for i, src in enumerate(srcList):
+            b[:,i] *= - 1./(1j*omega(src.freq))
+            # There is no magnetic source in the MT problem
+            # S_m, _ = src.eval(self.survey.prob)
+            # if S_m is not None:
+            #     b[:,i] += 1./(1j*omega(src.freq)) * S_m
+        return b
+
+    def _b(self, eSolution, srcList):
+        return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
+
+    def _bSecondaryDeriv_u(self, src, v, adjoint = False):
+        C = self.mesh.nodalGrad
+        if adjoint:
+            return - 1./(1j*omega(src.freq)) * (C.T * v)
+        return - 1./(1j*omega(src.freq)) * (C * v)
+
+    def _bSecondaryDeriv_m(self, src, v, adjoint = False):
+        # Doesn't depend on m
+        # _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
+        # S_eDeriv = S_eDeriv(v)
+        # if S_eDeriv is not None:
+        #     return 1./(1j * omega(src.freq)) * S_eDeriv
+        return None
+
+    def _bDeriv_u(self, src, v, adjoint=False):
+        # Primary does not depend on u
+        return self._bSecondaryDeriv_u(src, v, adjoint)
+
+    def _bDeriv_m(self, src, v, adjoint=False):
+        # Assuming the primary does not depend on the model
+        return self._bSecondaryDeriv_m(src, v, adjoint)
+
+    def _fDeriv_u(self, src, v, adjoint=False):
+        """
+        Derivative of the fields object wrt u.
+
+        :param NSEMsrc src: NSEM source
+        :param numpy.ndarray v: random vector of f_sol.size
+        This function stacks the fields derivatives appropriately
+
+        return a vector of size (nreEle+nrbEle)
+        """
+
+        de_du = v #Utils.spdiag(np.ones((self.nF,)))
+        db_du = self._bDeriv_u(src, v, adjoint)
+        # Return the stack
+        # This doesn't work...
+        return np.vstack((de_du,db_du))
+
+    def _fDeriv_m(self, src, v, adjoint=False):
+        """
+        Derivative of the fields object wrt m.
+
+        This function stacks the fields derivatives appropriately
+        """
+        return None
+
+
+class Fields1D_eTotal(BaseNSEMFields):
+    """
+    Fields storage for the 1D NSEM solution solved with for a total domain formulation.
+
+    Used in conjuction with Problem1D_eTotal.
+    """
+    knownFields = {'e_1dSolution':'F'}
+    aliasFields = {
+                    'e_1d' : ['e_1dSolution','F','_e'],
+                    'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
+                    'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
+                    'b_1d' : ['e_1dSolution','E','_b'],
+                    'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
+                    'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
+                  }
+
+    def __init__(self,mesh,survey,**kwargs):
+        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
+
+    def _ePrimary(self, eSolution, srcList):
+        ePrimary = np.zeros_like(eSolution)
+        for i, src in enumerate(srcList):
+            ep = src.ePrimary(self.survey.prob)
+            if ep is not None:
+                ePrimary[:,i] = ep[:,-1]
+        return ePrimary
+
+    def _eSecondary(self, eSolution, srcList):
+        return eSolution
+
+    def _e(self, eSolution, srcList):
+        return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
+
    def _eDeriv_u(self, src, v, adjoint = False):
        return v

@@ -120,7 +235,16 @@ class Fields1D_e(BaseNSEMFields):
        """
        return None

-class Fields3D_e(BaseNSEMFields):
+
+###########
+# 2D Fields
+###########
+
+
+###########
+# 3D Fields
+###########
+class Fields3D_ePrimSec(BaseNSEMFields):
    """
    Fields storage for the 3D NSEM solution. Labels polarizations by px and py.

@@ -1,8 +1,9 @@
+from SimPEG.EM.Utils import omega, mu_0
 from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
 from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
 from SurveyNSEM import Survey, Data
-from FieldsNSEM import BaseNSEMFields
-
+from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
+from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields

 class BaseNSEMProblem(BaseFDEMProblem):
    """
@@ -56,9 +57,9 @@ class BaseNSEMProblem(BaseFDEMProblem):
                # We need fDeriv_m = df/du*du/dm + df/dm
                # Construct du/dm, it requires a solve
                # NOTE: need to account for the 2 polarizations in the derivatives.
-                f_src = f[src,:]
+                u_src = f[src,:] # u should be a vector by definition. Need to fix this...
                # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
-                dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns.
+                dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
                dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
                if dRHS_dm is None:
                    du_dm = dA_duI * ( -dA_dm )
@@ -80,7 +81,7 @@ class BaseNSEMProblem(BaseFDEMProblem):

            :param numpy.ndarray m (nC, 1) - conductive model
            :param numpy.ndarray v (nD, 1) - vector
-            :param NSEMfields object u (optional) - NSEM fields object, if not given it is calculated
+            :param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
            :rtype: NSEMdata object
            :return: Data sensitivities wrt m
        """
@@ -103,7 +104,7 @@ class BaseNSEMProblem(BaseFDEMProblem):

            for src in self.survey.getSrcByFreq(freq):
                ftype = self._fieldType + 'Solution'
-                f_src = f[src, :]
+                f_src = f[src, :] # Need to fix this...

                for rx in src.rxList:
                    # Get the adjoint evalDeriv
@@ -130,3 +131,428 @@ class BaseNSEMProblem(BaseFDEMProblem):
            # Clean the factorization, clear memory.
            ATinv.clean()
        return Jtv
+
+###################################
+## 1D problems
+###################################
+
+class Problem1D_ePrimSec(BaseNSEMProblem):
+    """
+    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly half space ).
+
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = 'e_1dSolution'
+    _formulation  = 'EF'
+    fieldsPair = Fields1D_ePrimSec
+
+    # Initiate properties
+    _sigmaPrimary = None
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+        self.fieldsPair = Fields1D_e
+        # self._sigmaPrimary = sigmaPrimary
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        return A
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        The derivative of A wrt sigma
+        """
+
+        dsig_dm = self.curModel.sigmaDeriv
+        MeMui = self.MeMui
+        #
+        u_src = u['e_1dSolution']
+        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
+        if adjoint:
+            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
+        # Note: output has to be nN/nF, not nC/nE.
+        # v should be nC
+        return 1j * omega(freq) * ( dMfSigma_dm * v )
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
+            :return: RHS for 1 polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS wrt sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+        # Make the fields object
+        F = self.fieldsPair(self.mesh, self.survey)
+        # Loop over the frequencies
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
+            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
+
+            # Note curl e = -iwb so b = -curl e /iw
+            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
+            # F[Src, 'b_1d'] = b[:,1]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F
+
+# Note this is not fully functional.
+# Missing:
+# Fields class corresponding to the fields
+# Update Jvec and Jtvec to include all the derivatives components
+# Other things ...
+class Problem1D_eTotal(BaseNSEMProblem):
+    """
+    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
+
+    Solves the equation:
+
+    Math:
+        Have to do this...
+        Not implement correctly.......
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = 'e_1dSolution'
+    _formulation  = 'EF'
+    fieldsPair = Fields1D_eTotal
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+    @property
+    def MeMui(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MeMui', None) is None:
+            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
+        return self._MeMui
+
+    @property
+    def MfSigma(self):
+        """
+            Edge inner product matrix
+        """
+        if getattr(self, '_MfSigma', None) is None:
+            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
+        return self._MfSigma
+
+    def getA(self, freq, full=False):
+        """
+            Function to get the A matrix.
+
+            :param float freq: Frequency
+            :param logic full: Return full A or the inner part
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+
+        MeMui = self.MeMui
+        MfSigma = self.MfSigma
+        # Note: need to use the code above since in the 1D problem I want
+        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
+        # Possible that _fieldType and _eqLocs can fix this
+        # MeMui = self.MfMui
+        # MfSigma = self.MfSigma
+        C = self.mesh.nodalGrad
+        # Make A
+        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
+        # Either return full or only the inner part of A
+        if full:
+            return A
+        else:
+            return A[1:-1,1:-1]
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        raise NotImplementedError('getADeriv is not implemented')
+
+    def getRHS(self, freq):
+        """
+            Function to return the right hand side for the system.
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+        # Get sources for the frequency
+        # NOTE: Need to use the source information, doesn't really apply in 1D
+        src = self.survey.getSrcByFreq(freq)
+        # Get the full A
+        A = self.getA(freq,full=True)
+        # Define the outer part of the solution matrix
+        Aio = A[1:-1,[0,-1]]
+        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
+        Etot = (Ed + Eu)
+        sourceAmp = 1.0
+        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
+        ## Note: The analytic solution is derived with e^iwt
+        eBC = np.r_[Etot[0],Etot[-1]]
+        # The right hand side
+
+        return -Aio*eBC, eBC
+
+    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
+        raise NotImplementedError('getRHSDeriv not implemented yet')
+        return None
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        :param np.ndarray (nC,) m_back: Background conductivity model
+        '''
+        self.curModel = m
+        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
+
+        F = Fields1D_eTotal(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs, e_o = self.getRHS(freq)
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_i = Ainv * rhs
+            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)
+            # NOTE: only store e fields
+            F[Src, 'e_1dSolution'] = e[:,0]
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+        return F
+
+
+###################################
+## 3D problems
+###################################
+class Problem3D_ePrimSec(BaseNSEMProblem):
+    """
+    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
+
+    By eliminating the magnetic flux density using
+
+        .. math ::
+
+            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
+
+
+    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
+
+    .. math ::
+        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
+    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
+
+    The primary field is estimated from a background model (commonly as a 1D model).
+
+    """
+
+    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
+    _solutionType = [ 'e_pxSolution', 'e_pySolution']  # Forces order on the object
+    _formulation  = 'EB'
+    fieldsPair = Fields3D_ePrimSec
+
+    # Initiate properties
+    _sigmaPrimary = None
+
+    def __init__(self, mesh, **kwargs):
+        BaseNSEMProblem.__init__(self, mesh, **kwargs)
+
+    @property
+    def sigmaPrimary(self):
+        """
+        A background model, use for the calculation of the primary fields.
+
+        """
+        return self._sigmaPrimary
+    @sigmaPrimary.setter
+    def sigmaPrimary(self, val):
+        # Note: TODO add logic for val, make sure it is the correct size.
+        self._sigmaPrimary = val
+
+    def getA(self, freq):
+        """
+            Function to get the A system.
+
+            :param float freq: Frequency
+            :rtype: scipy.sparse.csr_matrix
+            :return: A
+        """
+        Mmui = self.MfMui
+        Msig = self.MeSigma
+        C = self.mesh.edgeCurl
+
+        return C.T*Mmui*C + 1j*omega(freq)*Msig
+
+    def getADeriv_m(self, freq, u, v, adjoint=False):
+        """
+        Calculate the derivative of A wrt m.
+
+        """
+        # Fix u to be a matrix nE,2
+        # This considers both polarizations and returns a nE,2 matrix for each polarization
+        if adjoint:
+            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
+        else:
+            # Need a nE,2 matrix to be returned
+            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
+        return 1j * omega(freq) * dMe_dsigV
+
+
+    def getRHS(self, freq):
+        """
+        Function to return the right hand side for the system.
+
+            :param float freq: Frequency
+            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
+            :return: RHS for both polarizations, primary fields
+        """
+
+        # Get sources for the frequncy(polarizations)
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_e = Src.S_e(self)
+        return -1j * omega(freq) * S_e
+
+    def getRHSDeriv_m(self, freq, v, adjoint=False):
+        """
+        The derivative of the RHS with respect to sigma
+        """
+
+        Src = self.survey.getSrcByFreq(freq)[0]
+        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
+        return -1j * omega(freq) * S_eDeriv
+
+    def fields(self, m):
+        '''
+        Function to calculate all the fields for the model m.
+
+        :param np.ndarray (nC,) m: Conductivity model
+        '''
+        # Set the current model
+        self.curModel = m
+
+        F = self.fieldsPair(self.mesh, self.survey)
+        for freq in self.survey.freqs:
+            if self.verbose:
+                startTime = time.time()
+                print 'Starting work for {:.3e}'.format(freq)
+                sys.stdout.flush()
+            A = self.getA(freq)
+            rhs  = self.getRHS(freq)
+            # Solve the system
+            Ainv = self.Solver(A, **self.solverOpts)
+            e_s = Ainv * rhs
+
+            # Store the fields
+            Src = self.survey.getSrcByFreq(freq)[0]
+            # Store the fields
+            # Use self._solutionType
+            F[Src, 'e_pxSolution'] = e_s[:,0]
+            F[Src, 'e_pySolution'] = e_s[:,1]
+            # Note curl e = -iwb so b = -curl/iw
+
+            if self.verbose:
+                print 'Ran for {:f} seconds'.format(time.time()-startTime)
+                sys.stdout.flush()
+            Ainv.clean()
+        return F
@@ -1,291 +0,0 @@
-from SimPEG.EM.Utils import omega
-from SimPEG import mkvc
-from scipy.constants import mu_0
-from SimPEG.NSEM.NSEM import BaseNSEMProblem
-from SimPEG.NSEM.SurveyNSEM import Survey, Data
-from SimPEG.NSEM.FieldsNSEM import Fields1D_e
-from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
-import numpy as np
-import multiprocessing, sys, time
-
-
-class eForm_psField(BaseNSEMProblem):
-    """
-    A NSEM problem soving a e formulation and primary/secondary fields decomposion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly half space ).
-
-
-    """
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e_1d'
-    _eqLocs    = 'EF'
-    _sigmaPrimary = None
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-        self.fieldsPair = Fields1D_e
-        # self._sigmaPrimary = sigmaPrimary
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        return A
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        The derivative of A wrt sigma
-        """
-
-        dsig_dm = self.curModel.sigmaDeriv
-        MeMui = self.MeMui
-        #
-        u_src = u['e_1dSolution']
-        dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
-        if adjoint:
-            return 1j * omega(freq) * (  dMfSigma_dm.T * v )
-        # Note: output has to be nN/nF, not nC/nE.
-        # v should be nC
-        return 1j * omega(freq) * ( dMfSigma_dm * v )
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
-            :return: RHS for 1 polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS wrt sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # NOTE: only store the e_solution(secondary), all other components calculated in the fields object
-            F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
-
-            # Note curl e = -iwb so b = -curl e /iw
-            # b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
-            # F[Src, 'b_1d'] = b[:,1]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
-
-# Note this is not fully functional.
-# Missing:
-# Fields class corresponding to the fields
-# Update Jvec and Jtvec to include all the derivatives components
-# Other things ...
-class eForm_TotalField(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a Total bondary domain decompostion.
-
-    Solves the equation:
-
-    Math:
-
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'EF'
-
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-    @property
-    def MeMui(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MeMui', None) is None:
-            self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
-        return self._MeMui
-
-    @property
-    def MfSigma(self):
-        """
-            Edge inner product matrix
-        """
-        if getattr(self, '_MfSigma', None) is None:
-            self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
-        return self._MfSigma
-
-    def getA(self, freq, full=False):
-        """
-            Function to get the A matrix.
-
-            :param float freq: Frequency
-            :param logic full: Return full A or the inner part
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-
-        MeMui = self.MeMui
-        MfSigma = self.MfSigma
-        # Note: need to use the code above since in the 1D problem I want
-        # e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
-        # Possible that _fieldType and _eqLocs can fix this
-        # MeMui = self.MfMui
-        # MfSigma = self.MfSigma
-        C = self.mesh.nodalGrad
-        # Make A
-        A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
-        # Either return full or only the inner part of A
-        if full:
-            return A
-        else:
-            return A[1:-1,1:-1]
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        raise NotImplementedError('getADeriv is not implemented')
-
-    def getRHS(self, freq):
-        """
-            Function to return the right hand side for the system.
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-        # Get sources for the frequency
-        # NOTE: Need to use the source information, doesn't really apply in 1D
-        src = self.survey.getSrcByFreq(freq)
-        # Get the full A
-        A = self.getA(freq,full=True)
-        # Define the outer part of the solution matrix
-        Aio = A[1:-1,[0,-1]]
-        Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
-        Etot = (Ed + Eu)
-        sourceAmp = 1.0
-        Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
-        ## Note: The analytic solution is derived with e^iwt
-        eBC = np.r_[Etot[0],Etot[-1]]
-        # The right hand side
-
-        return -Aio*eBC, eBC
-
-    def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
-        raise NotImplementedError('getRHSDeriv not implemented yet')
-        return None
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        :param np.ndarray (nC,) m_back: Background conductivity model
-        '''
-        self.curModel = m
-        # RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
-
-        F = Fields1D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs, e_o = self.getRHS(freq)
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_i = Ainv * rhs
-            e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)
-            # NOTE: only store e fields
-            F[Src, 'e_1dSolution'] = e[:,0]
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-        return F
@@ -1 +0,0 @@
-from Probs import eForm_TotalField, eForm_psField
@@ -1 +0,0 @@
-pass
@@ -1,138 +0,0 @@
-from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
-from SimPEG.EM.Utils import omega
-from scipy.constants import mu_0
-from SimPEG.NSEM.NSEM import BaseNSEMProblem
-from SimPEG.NSEM.SurveyNSEM import Survey, Data
-from SimPEG.NSEM.FieldsNSEM import Fields3D_e
-import multiprocessing, sys, time
-
-
-
-class eForm_ps(BaseNSEMProblem):
-    """
-    A NSEM problem solving a e formulation and a primary/secondary fields decompostion.
-
-    By eliminating the magnetic flux density using
-
-        .. math ::
-
-            \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
-
-
-    we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
-
-    .. math ::
-        \\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
-    which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
-
-    The primary field is estimated from a background model (commonly as a 1D model).
-
-    """
-
-    # From FDEMproblem: Used to project the fields. Currently not used for NSEMproblem.
-    _fieldType = 'e'
-    _eqLocs    = 'FE'
-    fieldsPair = Fields3D_e
-    _sigmaPrimary = None
-
-    def __init__(self, mesh, **kwargs):
-        BaseNSEMProblem.__init__(self, mesh, **kwargs)
-
-    @property
-    def sigmaPrimary(self):
-        """
-        A background model, use for the calculation of the primary fields.
-
-        """
-        return self._sigmaPrimary
-    @sigmaPrimary.setter
-    def sigmaPrimary(self, val):
-        # Note: TODO add logic for val, make sure it is the correct size.
-        self._sigmaPrimary = val
-
-    def getA(self, freq):
-        """
-            Function to get the A system.
-
-            :param float freq: Frequency
-            :rtype: scipy.sparse.csr_matrix
-            :return: A
-        """
-        Mmui = self.MfMui
-        Msig = self.MeSigma
-        C = self.mesh.edgeCurl
-
-        return C.T*Mmui*C + 1j*omega(freq)*Msig
-
-    def getADeriv_m(self, freq, u, v, adjoint=False):
-        """
-        Calculate the derivative of A wrt m.
-
-        """
-
-        # This considers both polarizations and returns a nE,2 matrix for each polarization
-        if adjoint:
-            dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
-        else:
-            # Need a nE,2 matrix to be returned
-            dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
-        return 1j * omega(freq) * dMe_dsigV
-
-
-    def getRHS(self, freq):
-        """
-        Function to return the right hand side for the system.
-
-            :param float freq: Frequency
-            :rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
-            :return: RHS for both polarizations, primary fields
-        """
-
-        # Get sources for the frequncy(polarizations)
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_e = Src.S_e(self)
-        return -1j * omega(freq) * S_e
-
-    def getRHSDeriv_m(self, freq, v, adjoint=False):
-        """
-        The derivative of the RHS with respect to sigma
-        """
-
-        Src = self.survey.getSrcByFreq(freq)[0]
-        S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
-        return -1j * omega(freq) * S_eDeriv
-
-    def fields(self, m):
-        '''
-        Function to calculate all the fields for the model m.
-
-        :param np.ndarray (nC,) m: Conductivity model
-        '''
-        # Set the current model
-        self.curModel = m
-
-        F = Fields3D_e(self.mesh, self.survey)
-        for freq in self.survey.freqs:
-            if self.verbose:
-                startTime = time.time()
-                print 'Starting work for {:.3e}'.format(freq)
-                sys.stdout.flush()
-            A = self.getA(freq)
-            rhs  = self.getRHS(freq)
-            # Solve the system
-            Ainv = self.Solver(A, **self.solverOpts)
-            e_s = Ainv * rhs
-
-            # Store the fields
-            Src = self.survey.getSrcByFreq(freq)[0]
-            # Store the fieldss
-            F[Src, 'e_pxSolution'] = e_s[:,0]
-            F[Src, 'e_pySolution'] = e_s[:,1]
-            # Note curl e = -iwb so b = -curl/iw
-
-            if self.verbose:
-                print 'Ran for {:f} seconds'.format(time.time()-startTime)
-                sys.stdout.flush()
-            Ainv.clean()
-        return F
-
@@ -1 +0,0 @@
-from Probs import eForm_ps
				`@@ -1 +0,0 @@`
				`from Probs import eForm_TotalField, eForm_psField`