minor fix

- make PTv out from the loop
- do not use Data class in the loop for Jv

.travis.yml

@@ -35,7 +35,7 @@ before_install:
 
 # Install packages
 install:
-  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython ipywidgets nose vtk
+  - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk
   - pip install nose-cov python-coveralls
 
   - git clone https://github.com/rowanc1/pymatsolver.git

SimPEG/Directives.py

@@ -213,33 +213,42 @@ class SaveOutputEveryIteration(_SaveEveryIteration):
         f.close()
 
 class SaveOutputDictEveryIteration(_SaveEveryIteration):
-    """
-    Saves inversion parameters at every iteraion.
-
-
-    """
+    """SaveOutputDictEveryIteration"""
 
     def initialize(self):
         print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName
 
     def endIter(self):
-
-        # Initialize the output dict
-        outDict = {}
         # Save the data.
-        outDict['iter'] = self.opt.iter
-        outDict['beta'] = self.invProb.beta
-        outDict['phi_d'] = self.invProb.phi_d
-        outDict['phi_ms'] = self.reg._evalSmall(self.invProb.curModel)
-        outDict['phi_mx'] = self.reg._evalSmoothx(self.invProb.curModel)
-        outDict['phi_my'] = self.reg._evalSmoothy(self.invProb.curModel) if self.prob.mesh.dim >= 2 else 'NaN'
-        outDict['phi_mz'] = self.reg._evalSmoothz(self.invProb.curModel) if self.prob.mesh.dim==3 else 'NaN'
-        outDict['f'] = self.opt.f
-        outDict['m'] = self.invProb.curModel
-        outDict['dpred'] = self.invProb.dpred
+        ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
+        phi_ms = 0.5*ms.dot(ms)
+        if self.reg.mrefInSmooth == True:
+            mref = self.reg.mref
+        else:
+            mref = 0
+        mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
+        phi_mx = 0.5 * mx.dot(mx)
+        if self.prob.mesh.dim >= 2:
+            my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_my = 0.5 * my.dot(my)
+        else:
+            phi_my = 'NaN'
+        if self.prob.mesh.dim==3:
+            mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
+            phi_mz = 0.5 * mz.dot(mz)
+        else:
+            phi_mz = 'NaN'
+
 
         # Save the file as a npz
-        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), outDict)
+        np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
+
+#         mref = getattr(self, 'm_prev', None)
+#         if mref is None:
+#             if self.debug: print 'UpdateReferenceModel is using mref0'
+#             mref = self.mref0
+#         self.m_prev = self.invProb.m_current
+#         return mref
 
 class Update_IRLS(InversionDirective):
 

SimPEG/EM/FDEM/FieldsFDEM.py

@@ -60,6 +60,20 @@ class Fields(SimPEG.Problem.Fields):
 
         return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList)
 
+    def _bSecondary(self, solution, srcList):
+        """
+        Total magnetic flux density is sum of primary and secondary
+
+        :param numpy.ndarray solution: field we solved for
+        :param list srcList: list of sources
+        :rtype: numpy.ndarray
+        :return: total magnetic flux density
+        """
+        if getattr(self, '_bSecondary', None) is None:
+            raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0])
+
+        return  self._bSecondary(solution, srcList)
+
     def _h(self, solution, srcList):
         """
         Total magnetic field is sum of primary and secondary
@@ -124,6 +138,21 @@ class Fields(SimPEG.Problem.Fields):
             return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint)
         return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex)
 
+    def _bSecondaryDeriv(self, src, du_dm_v, v, adjoint = False):
+        """
+        Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
+
+        :param Src src: sorce
+        :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint)
+        :param numpy.ndarray v: vector to take sensitivity product with
+        :param bool adjoint: adjoint?
+        :rtype: numpy.ndarray
+        :return: derivative times a vector (or tuple for adjoint)
+        """
+        # TODO: modify when primary field is dependent on m
+
+        return self._bDeriv(src, du_dm_v, v, adjoint = adjoint)
+
     def _hDeriv(self, src, du_dm_v, v, adjoint = False):
         """
         Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint
@@ -471,6 +500,8 @@ class Fields3D_b(Fields):
             return 'E'
         elif fieldType == 'b':
             return 'F'
+        elif fieldType == 'bSecondary':
+            return 'F'
         elif (fieldType == 'h') or (fieldType == 'j'):
             return'CCV'
         else:

SimPEG/EM/FDEM/RxFDEM.py

@@ -97,6 +97,19 @@ class Point_b(BaseRx):
         self.projField = 'b'
         super(Point_b, self).__init__(locs, orientation, component)
 
+class Point_bSecondary(BaseRx):
+    """
+    Magnetic flux FDEM receiver
+
+    :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`)
+    :param string orientation: receiver orientation 'x', 'y' or 'z'
+    :param string component: real or imaginary component 'real' or 'imag'
+    """
+
+    def __init__(self, locs, orientation=None, component=None):
+        self.projField = 'bSecondary'
+        super(Point_bSecondary, self).__init__(locs, orientation, component)
+
 
 class Point_h(BaseRx):
     """

SimPEG/EM/FDEM/SrcFDEM.py

@@ -555,7 +555,7 @@ class CircularLoop(BaseSrc):
             a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu)
 
         else:
-            srcfct = MagneticDipoleVectorPotential
+            srcfct = MagneticLoopVectorPotential
             ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu)
             ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu)
             az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu)

SimPEG/EM/TDEM/FieldsTDEM.py

@@ -0,0 +1,142 @@
+import numpy as np
+import scipy.sparse as sp
+import SimPEG
+from SimPEG import Utils
+from SimPEG.EM.Utils import omega
+from SimPEG.Utils import Zero, Identity
+
+class Fields(SimPEG.Problem.TimeFields):
+    """
+
+    Fancy Field Storage for a TDEM survey. Only one field type is stored for
+    each problem, the rest are computed. The fields obejct acts like an array and is indexed by
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[srcList,'e']
+        b = f[srcList,'b']
+
+    If accessing all sources for a given field, use the :code:`:`
+
+    .. code-block:: python
+
+        f = problem.fields(m)
+        e = f[:,'e']
+        b = f[:,'b']
+
+    The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies)
+    """
+
+    knownFields = {}
+    dtype = float
+
+    def _eDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._eDeriv_u(tInd, src, v, adjoint), self._eDeriv_m(tInd, src, v, adjoint)
+        return self._eDeriv_u(tInd, src, dun_dm_v) + self._eDeriv_m(tInd, src, v)
+
+    def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+        if adjoint is True:
+            return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+        return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+class Fields_Derivs(Fields):
+    knownFields = {
+                    'bDeriv': 'F',
+                    'eDeriv': 'E',
+                    'hDeriv': 'E',
+                    'jDeriv': 'F'
+                  }
+
+
+class Fields_b(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'bSolution': 'F'}
+    aliasFields = {
+                    'b': ['bSolution', 'F', '_b'],
+                    'e': ['bSolution', 'E', '_e'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+    def _b(self, bSolution, srcList, tInd):
+        return bSolution
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        return Identity()*dun_dm_v
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)
+
+    def _e(self, bSolution, srcList, tInd):
+        e = self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * bSolution ) )
+        for i, src in enumerate(srcList):
+            _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+            e[:,i] = e[:,i] - self.MeSigmaI * S_e
+        return e
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        if adjoint is True:
+            return self.MfMui.T * ( self.edgeCurl * ( self.MeSigmaI.T * dun_dm_v ) )
+        return self.MeSigmaI * ( self.edgeCurl.T * ( self.MfMui * dun_dm_v ) )
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        _, S_e = src.eval(self.survey.prob, self.survey.prob.times[tInd])
+        bSolution = self[[src],'bSolution',tInd]
+
+        _, S_eDeriv = src.evalDeriv(self.survey.prob.times[tInd], self, adjoint=adjoint)
+
+        if adjoint is True:
+            return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution ) ).T * v - S_eDeriv(self.MeSigmaI.T * v)
+
+        return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v)
+
+
+
+class Fields_e(Fields):
+    """Fancy Field Storage for a TDEM survey."""
+    knownFields = {'eSolution': 'E'}
+    aliasFields = {
+                    'e': ['eSolution', 'E', '_e'],
+                    'b': ['eSolution', 'F', '_b'],
+                  }
+
+    def startup(self):
+        self.MeSigmaI      = self.survey.prob.MeSigmaI
+        self.MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv
+        self.edgeCurl      = self.survey.prob.mesh.edgeCurl
+        self.MfMui         = self.survey.prob.MfMui
+
+
+    def _e(self, eSolution, srcList, tInd):
+        return eSolution
+
+    def _eDeriv_u(self, tInd, src, dun_dm_v, adjoint = False):
+        return dun_dm_v
+
+    def _eDeriv_m(self, tInd, src, v, adjoint = False):
+        return Zero()
+
+    def _b(self, eSolution, srcList, tInd):
+        raise NotImplementedError
+
+    def _bDeriv_u(self, tInd, src, dun_dm_v, adjoint=False):
+        raise NotImplementedError
+
+    def _bDeriv_m(self, tInd, src, v, adjoint=False):
+        raise NotImplementedError
+
+    # def _bDeriv(self, tInd, src, dun_dm_v, v, adjoint=False):
+    #     if adjoint is True:
+    #         return self._bDeriv_u(tInd, src, v, adjoint), self._bDeriv_m(tInd, src, v, adjoint)
+    #     return self._bDeriv_u(tInd, src, dun_dm_v) + self._bDeriv_m(tInd, src, v)

SimPEG/EM/TDEM/SrcTDEM.py

@@ -0,0 +1,246 @@
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
+from SimPEG.EM.Utils import *
+
+####################################################
+# Sources
+####################################################
+
+class BaseWaveform(object):
+
+    def __init__(self, offTime=0., hasInitialFields=False):
+        self.offTime = offTime
+        self.hasInitialFields = hasInitialFields
+
+    def _assertMatchesPair(self, pair):
+        assert (isinstance(self, pair)
+            ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__)
+
+    def eval(self, time):
+        raise NotImplementedError
+
+    def evalDeriv(self, time):
+        raise NotImplementedError # needed for E-formulation
+
+
+class StepOffWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        return 0.
+
+
+class RawWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('RawWaveform has not been implemented, you should write it!')
+
+
+class TriangularWaveform(BaseWaveform):
+
+    def __init__(self, offTime=0.):
+        BaseWaveform.__init__(self, offTime, hasInitialFields=True)
+
+    def eval(self, time):
+        raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!')
+
+
+
+class BaseSrc(SimPEG.Survey.BaseSrc):
+
+    # rxPair = Rx
+    integrate = True
+    waveformPair = BaseWaveform
+
+    @property
+    def waveform(self):
+        "A waveform instance is not None"
+        return getattr(self, '_waveform', None)
+    @waveform.setter
+    def waveform(self, val):
+        if self.waveform is None:
+            val._assertMatchesPair(self.waveformPair)
+            self._mapping = val
+        else:
+            self._mapping = self.PropMap(val)
+
+
+    def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs):
+        self.waveform = waveform
+        SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
+
+
+    def bInitial(self, prob):
+        return Zero()
+
+    def bInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eInitial(self, prob):
+        return Zero()
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+        return Zero()
+
+    def eval(self, prob, time):
+        S_m = self.S_m(prob, time)
+        S_e = self.S_e(prob, time)
+        return S_m, S_e
+
+    def evalDeriv(self, prob, time, v=None, adjoint=False):
+        if v is not None:
+            return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint)
+        else:
+            return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint)
+
+    def S_m(self, prob, time):
+        return Zero()
+
+    def S_e(self, prob, time):
+        return Zero()
+
+    def S_mDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+    def S_eDeriv(self, prob, time, v=None, adjoint=False):
+        return Zero()
+
+
+class MagDipole(BaseSrc):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    moment = 1.
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+
+        else:
+            srcfct = MagneticDipoleVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+
+
+    def bInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        return self._bfromVectorPotential(prob)
+
+    def eInitial(self, prob):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaI = prob.MeSigmaI
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+
+        return MeSigmaI * (C.T * (MfMui * b))
+
+    def eInitialDeriv(self, prob, v=None, adjoint=False):
+
+        if self.waveform.hasInitialFields is False:
+            return Zero()
+
+        b = self.bInitial(prob)
+        MeSigmaIDeriv = prob.MeSigmaIDeriv
+        MfMui = prob.MfMui
+        C = prob.mesh.edgeCurl
+        S_e = self.S_e(prob, prob.t0)
+
+        # S_e doesn't depend on the model
+
+        if adjoint:
+            return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v
+
+        return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
+
+
+    def S_m(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+    def S_e(self, prob, time):
+        if self.waveform.hasInitialFields is False:
+            raise NotImplementedError
+        return Zero()
+
+class CircularLoop(MagDipole):
+
+    waveform = None
+    loc = None
+    orientation = 'Z'
+    radius = None
+    mu = mu_0
+
+    def __init__(self, rxList, **kwargs):
+        assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
+        self.integrate = False
+        BaseSrc.__init__(self, rxList, **kwargs)
+
+    def _bfromVectorPotential(self, prob):
+        if prob._eqLocs is 'FE':
+            gridX = prob.mesh.gridEx
+            gridY = prob.mesh.gridEy
+            gridZ = prob.mesh.gridEz
+            C = prob.mesh.edgeCurl
+
+        elif prob._eqLocs is 'EF':
+            gridX = prob.mesh.gridFx
+            gridY = prob.mesh.gridFy
+            gridZ = prob.mesh.gridFz
+            C = prob.mesh.edgeCurl.T
+
+
+        if prob.mesh._meshType is 'CYL':
+            if not prob.mesh.isSymmetric:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            a = MagneticLoopVectorPotential(self.loc, gridY, 'y', radius=self.radius, mu=self.mu)
+
+        else:
+            srcfct = MagneticLoopVectorPotential
+            ax = srcfct(self.loc, gridX, 'x', mu=self.mu, radius=self.radius)
+            ay = srcfct(self.loc, gridY, 'y', mu=self.mu, radius=self.radius)
+            az = srcfct(self.loc, gridZ, 'z', mu=self.mu, radius=self.radius)
+            a = np.concatenate((ax, ay, az))
+
+        return C*a
+

SimPEG/EM/TDEM/SurveyTDEM.py

@@ -1,10 +1,16 @@
-from SimPEG import Utils, Survey, np
-from SimPEG.Survey import BaseSurvey
+import SimPEG
+from SimPEG import np, Utils
+from SimPEG.Utils import Zero, Identity
+from scipy.constants import mu_0
 from SimPEG.EM.Utils import *
-from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src
 
 
-class RxTDEM(Survey.BaseTimeRx):
+####################################################
+# Receivers
+####################################################
+
+class Rx(SimPEG.Survey.BaseTimeRx):
 
     knownRxTypes = {
                     'ex':['e', 'Ex', 'N'],
@@ -21,7 +27,7 @@ class RxTDEM(Survey.BaseTimeRx):
                    }
 
     def __init__(self, locs, times, rxType):
-        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+        SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType)
 
     @property
     def projField(self):
@@ -56,144 +62,60 @@ class RxTDEM(Survey.BaseTimeRx):
         u_part = Utils.mkvc(u[src, self.projField, :])
         return P*u_part
 
-    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+    def evalDeriv(self, src, mesh, timeMesh, v, adjoint=False):
         P = self.getP(mesh, timeMesh)
-
         if not adjoint:
-            return P * Utils.mkvc(v[src, self.projField, :])
+            return P * v #Utils.mkvc(v[src, self.projField+'Deriv', :])
         elif adjoint:
-            return P.T * v[src, self]
+            # dP_dF_T = P.T * v #[src, self]
+            # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN )
+            return P.T * v #np.reshape(dP_dF_T, newshape, order='F')
 
 
-class SrcTDEM(Survey.BaseSrc):
-    rxPair = RxTDEM
-    radius = None
+####################################################
+# Survey
+####################################################
 
-    def getInitialFields(self, mesh):
-        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
-        return F0
-
-    def getJs(self, mesh, time):
-        return None
-
-
-class SrcTDEM_VMD_MVP(SrcTDEM):
-
-    def __init__(self,rxList,loc,waveformType="STEPOFF"):
-        self.loc = loc
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Vertical magnetic dipole, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-            else:
-                raise Exception('Unknown mesh for VMD')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
-        else:
-            raise Exception('Unknown mesh for VMD')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SrcTDEM_CircularLoop_MVP(SrcTDEM):
-    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
-        self.loc = loc
-        self.radius = radius
-        self.waveformType = waveformType
-        SrcTDEM.__init__(self,rxList)
-
-    def getInitialFields(self, mesh):
-        """Circular Loop, magnetic vector potential"""
-        if self.waveformType == "STEPOFF":
-            print ">> Step waveform: Non-zero initial condition"
-            if mesh._meshType is 'CYL':
-                if mesh.isSymmetric:
-                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-                else:
-                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-            elif mesh._meshType is 'TENSOR':
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-            else:
-                raise Exception('Unknown mesh for CircularLoop')
-            return {"b": mesh.edgeCurl*MVP}
-        elif self.waveformType == "GENERAL":
-            print ">> General waveform: Zero initial condition"
-            return {"b": np.zeros(mesh.nF)}
-        else:
-            raise NotImplementedError("Only use STEPOFF or GENERAL")
-
-    def getMeS(self, mesh, MfMui):
-        if mesh._meshType is 'CYL':
-            if mesh.isSymmetric:
-                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
-            else:
-                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
-        elif mesh._meshType is 'TENSOR':
-            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
-        else:
-            raise Exception('Unknown mesh for CircularLoop')
-        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
-
-
-class SurveyTDEM(Survey.BaseSurvey):
+class Survey(SimPEG.Survey.BaseSurvey):
     """
-        docstring for SurveyTDEM
+    Time domain electromagnetic survey
     """
-    srcPair = SrcTDEM
+
+    srcPair = Src.BaseSrc
+    rxPair = Rx
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
         self.srcList = srcList
-        Survey.BaseSurvey.__init__(self, **kwargs)
+        SimPEG.Survey.BaseSurvey.__init__(self, **kwargs)
 
     def eval(self, u):
-        data = Survey.Data(self)
+        data = SimPEG.Survey.Data(self)
         for src in self.srcList:
             for rx in src.rxList:
                 data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
         return data
 
     def evalDeriv(self, u, v=None, adjoint=False):
-        assert v is not None, 'v to multiply must be provided.'
+        raise Exception('Use Receivers to project fields deriv.')
+        # assert v is not None, 'v to multiply must be provided.'
 
-        if not adjoint:
-            data = Survey.Data(self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
-            return data
-        else:
-            f = FieldsTDEM(self.mesh, self)
-            for src in self.srcList:
-                for rx in src.rxList:
-                    Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
-                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
-                    if rx.projField not in f: # first time we are projecting
-                        f[src, rx.projField, :] = Ptv
-                    else: # there are already fields, so let's add to them!
-                        f[src, rx.projField, :] += Ptv
-            return f
+        # if not adjoint:
+        #     data = SimPEG.Survey.Data(self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+        #     return data
+        # else:
+        #     f = FieldsTDEM(self.mesh, self)
+        #     for src in self.srcList:
+        #         for rx in src.rxList:
+        #             Ptv = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+        #             Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+        #             if rx.projField not in f: # first time we are projecting
+        #                 f[src, rx.projField, :] = Ptv
+        #             else: # there are already fields, so let's add to them!
+        #                 f[src, rx.projField, :] += Ptv
+        #     return f
 
 

SimPEG/EM/TDEM/TDEM.py

@@ -0,0 +1,553 @@
+from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver
+from SimPEG.EM.Base import BaseEMProblem
+from SimPEG.EM.TDEM.SurveyTDEM import Survey as SurveyTDEM
+from SimPEG.EM.TDEM.FieldsTDEM import *
+from scipy.constants import mu_0
+import time
+
+class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
+    """
+    We start with the first order form of Maxwell's equations
+    """
+    surveyPair = SurveyTDEM
+    fieldsPair = Fields
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        Problem.BaseTimeProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def fields(self, m):
+        """
+        Solve the forward problem for the fields.
+
+        :param numpy.array m: inversion model (nP,)
+        :rtype numpy.array:
+        :return F: fields
+        """
+
+        tic = time.time()
+        self.curModel = m
+
+        F = self.fieldsPair(self.mesh, self.survey)
+
+        # set initial fields
+        F[:,self._fieldType+'Solution',0] = self.getInitialFields()
+
+        # timestep to solve forward
+        if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv = None
+        for tInd, dt in enumerate(self.timeSteps):
+            if Ainv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Ainv.clean()
+                Ainv = None
+
+            if Ainv is None:
+                A = self.getAdiag(tInd)
+                if self.verbose: print 'Factoring...   (dt = %e)'%dt
+                Ainv = self.Solver(A, **self.solverOpts)
+                if self.verbose: print 'Done'
+
+            rhs = self.getRHS(tInd+1) # this is on the nodes of the time mesh
+            Asubdiag = self.getAsubdiag(tInd)
+
+            if self.verbose: print ('    Solving...   (tInd = %i)')% (tInd+1)
+            sol = Ainv * (rhs - Asubdiag * F[:,self._fieldType+'Solution',tInd]) # taking a step
+
+            if self.verbose: print '    Done...'
+
+            if sol.ndim == 1:
+                sol.shape = (sol.size,1)
+            F[:,self._fieldType+'Solution',tInd+1] = sol
+        if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
+        Ainv.clean()
+        return F
+
+
+    def Jvec(self, m, v, f=None):
+        """
+        Jvec computes the sensitivity times a vector
+
+        .. math::
+            \mathbf{J} \mathbf{v} = \\frac{d\mathbf{P}}{d\mathbf{F}} \left( \\frac{d\mathbf{F}}{d\mathbf{u}} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} \\right) \mathbf{v}
+
+        where
+
+        .. math::
+            \mathbf{A} \\frac{d\mathbf{u}}{d\mathbf{m}} + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} = \\frac{d \mathbf{RHS}}{d \mathbf{m}}
+        """
+
+        if f is None:
+           f = self.fields(m)
+
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+        self.curModel = m
+
+        # mat to store previous time-step's solution deriv times a vector for each source
+        # size: nu x nSrc
+
+        # this is a bit silly
+
+        # if self._fieldType is 'b' or self._fieldType is 'j':
+        #     ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        # elif self._fieldType is 'e' or self._fieldType is 'h':
+        #     ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        # for i, src in enumerate(self.survey.srcList):
+        dun_dm_v = np.hstack([Utils.mkvc(self.getInitialFieldsDeriv(src,v),2) for src in self.survey.srcList]) # can over-write this at each timestep
+        #
+        df_dm_v = Fields_Derivs(self.mesh, self.survey) # store the field derivs we need to project to calc full deriv
+
+        Adiaginv = None
+
+        for tInd, dt in zip(range(self.nT), self.timeSteps):
+            if Adiaginv is not None and (tInd > 0 and dt != self.timeSteps[tInd - 1]):# keep factors if dt is the same as previous step b/c A will be the same
+                Adiaginv.clean()
+                Adiaginv = None
+
+            if Adiaginv is None:
+                A = self.getAdiag(tInd)
+                Adiaginv = self.Solver(A, **self.solverOpts)
+
+            Asubdiag = self.getAsubdiag(tInd)
+
+            for i, src in enumerate(self.survey.srcList):
+
+                # here, we are lagging by a timestep, so filling in as we go
+                for projField in set([rx.projField for rx in src.rxList]):
+                    # Seogi: df_duFun?
+                    df_dmFun = getattr(f, '_%sDeriv'%projField, None)
+                    # df_dm_v is dense, but we only need the times at (rx.P.T * ones > 0)
+                    # This should be called rx.footprint
+                    df_dm_v[src, '%sDeriv'%projField , tInd] = df_dmFun(tInd, src, dun_dm_v[:,i], v)
+
+                un_src = f[src,ftype,tInd+1]
+
+                dA_dm_v   = self.getAdiagDeriv(tInd, un_src, v) # cell centered on time mesh
+                dRHS_dm_v = self.getRHSDeriv(tInd+1, src, v) # on nodes of time mesh
+
+                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, f[src,ftype,tInd], v)
+
+                JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
+
+                # step in time and overwrite
+                if tInd != len(self.timeSteps+1):
+                    dun_dm_v[:,i] = Adiaginv * (JRHS - Asubdiag * dun_dm_v[:,i])
+
+        # Seogi: suspcious spot
+        # Jv = self.dataPair(self.survey)
+        Jv = []
+        for src in self.survey.srcList:
+            for rx in src.rxList:
+                # Looping over data class append memory as well!!
+                # Jv[src,rx] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:]))
+                Jv.append(rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(df_dm_v[src,'%sDeriv'%rx.projField,:])))
+        Adiaginv.clean()
+        # del df_dm_v, dun_dm_v, Asubdiag
+        # return Utils.mkvc(Jv)
+        return np.hstack(Jv)
+
+    def Jtvec(self, m, v, f=None):
+
+        """
+        Jvec computes the adjoint of the sensitivity times a vector
+
+        .. math::
+            \mathbf{J}^\\top \mathbf{v} =  \left( \\frac{d\mathbf{u}}{d\mathbf{m}} ^ \\top \\frac{d\mathbf{F}}{d\mathbf{u}} ^ \\top  + \\frac{\partial\mathbf{F}}{\partial\mathbf{m}} ^ \\top \\right) \\frac{d\mathbf{P}}{d\mathbf{F}} ^ \\top \mathbf{v}
+
+        where
+
+        .. math::
+            \\frac{d\mathbf{u}}{d\mathbf{m}} ^\\top \mathbf{A}^\\top  + \\frac{d\mathbf{A}(\mathbf{u})}{d\mathbf{m}} ^ \\top = \\frac{d \mathbf{RHS}}{d \mathbf{m}} ^ \\top
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+        ftype = self._fieldType + 'Solution' # the thing we solved for
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        df_duT_v = Fields_Derivs(self.mesh, self.survey)
+        ATinv_df_duT_v = np.zeros((len(self.survey.srcList), len(f[self.survey.srcList[0],ftype,0])), dtype=float) # same size as fields at a single timestep
+
+        JTv = np.zeros(m.shape, dtype=float)
+
+        # Loop over sources and receivers to create a fields object: PT_v, df_duT_v, df_dmT_v
+        PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+        for src in self.survey.srcList:
+            # Looping over initializing field class is appending memory!
+            # PT_v = Fields_Derivs(self.mesh, self.survey) # initialize storage for PT_v (don't need to preserve over sources)
+            # initialize size
+            df_duT_v[src, '%sDeriv'%self._fieldType, :] = np.zeros_like(f[src, self._fieldType, :])
+
+            for rx in src.rxList:
+                print ('_%sDeriv')%(rx.projField)
+                PT_v[src,'%sDeriv'%rx.projField,:] = rx.evalDeriv(src, self.mesh, self.timeMesh, Utils.mkvc(v[src,rx]), adjoint=True) # this is +=
+
+                # PT_v = np.reshape(curPT_v,(len(curPT_v)/self.timeMesh.nN, self.timeMesh.nN), order='F')
+                df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
+
+                for tInd in range(self.nT+1):
+                    cur = df_duTFun(tInd, src, None, Utils.mkvc(PT_v[src,'%sDeriv'%rx.projField,tInd]), adjoint=True)
+                    df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
+                    JTv = cur[1] + JTv
+
+        del PT_v # no longer need this
+
+        AdiagTinv = None
+
+        # Do the back-solve through time
+        for tIndP in reversed(range(self.nT + 1)):
+            tInd = tIndP - 1
+            if AdiagTinv is not None and (tInd <= self.nT and self.timeSteps[tInd] != self.timeSteps[tInd+1]): # if the previous timestep is the same --> no need to refactor the matrix
+                AdiagTinv.clean()
+                AdiagTinv = None
+
+            # refactor if we need to
+            if AdiagTinv is None and tInd > -1:
+                Adiag = self.getAdiag(tInd)
+                AdiagTinv = self.Solver(Adiag.T, **self.solverOpts)
+
+            dAsubdiag_dm_v = Zero()
+
+            if tInd < self.nT - 1:
+                Asubdiag = self.getAsubdiag(tInd+1)
+
+
+            for isrc, src in enumerate(self.survey.srcList):
+                # solve against df_duT_v
+                if tInd >= self.nT-1:
+                    # last timestep (first to be solved)
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]
+                elif tInd > -1:
+                # else:
+                    ATinv_df_duT_v[isrc,:] = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                else:
+                    # AdiagTinv = I
+                    ATinv_df_duT_v[isrc,:] = Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])
+                    # - Utils.mkvc(Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+                    # (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+
+                if tInd < self.nT - 1:
+                    dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], ATinv_df_duT_v[isrc,:], adjoint = True)
+
+                if tInd > -1:
+                    un_src = f[src,ftype,tInd+1]
+                    dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    dRHST_dm_v = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+
+                    JTv = JTv + Utils.mkvc(- dAT_dm_v - dAsubdiag_dm_v + dRHST_dm_v)
+                else:
+                    # dA_dm_v = self.getInitialFieldsDeriv(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True)
+                    # print np.linalg.norm(self.getInitialFieldsDeriv(src, df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1], adjoint=True))
+                    # print np.linalg.norm(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # vec = - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]) + Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1])
+                    # dAsubdiagT_dm_v = self.getAsubdiagDeriv(tInd+1, f[src,ftype,tInd+1], Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint = True)
+                    dRHST_dm_v = Utils.mkvc(self.getInitialFieldsDeriv(src, Utils.mkvc(ATinv_df_duT_v[isrc,:]) , adjoint=True))
+
+                    JTv = JTv + Utils.mkvc( -dAsubdiagT_dm_v + dRHST_dm_v) #
+
+
+
+                    # # dAT_dm_v = self.getAdiagDeriv(tInd, un_src, ATinv_df_duT_v[isrc,:], adjoint=True) # cell centered on time mesh
+                    # dRHST_dm_v0 = self.getRHSDeriv(tInd+1, src, ATinv_df_duT_v[isrc,:], adjoint=True) # on nodes of time mesh
+                    # dRHST_dm_v1 = self.getInitialFieldsDeriv( Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    # JTv = JTv + Utils.mkvc(dRHST_dm_v0 + dRHST_dm_v1)
+
+                    # print 'here'
+                    # inFields = self.getInitialFieldsDeriv(f[src,ftype,tInd+1], Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,tInd+1]), adjoint=True)
+                    #  # - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]), adjoint=True)
+                    # print inFields.shape
+                    # JTv = JTv + inFields
+                # dAsubdiag_dm_v = 0
+
+
+
+        # Missing the 0 step
+
+        # adding du_dm^T * dF_du^T * P^T vfor time 0 (no dRHS_dm_v at time 0)
+        # Asubdiag = self.getAsubdiag(0)
+        # for src in self.survey.srcList:
+        #     for projField in set(rx.projField):
+        #         v = AdiagTinv * (Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0]) - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:]))
+        #         JTv = JTv - Utils.mkvc(self.getAdiagDeriv(0, f[src, ftype, tInd], v, adjoint = True))
+        #         # JTv = JTv + self.getInitialFieldsDeriv(Utils.mkvc(df_duT_v[src,'%sDeriv'%self._fieldType,0] - Asubdiag.T * Utils.mkvc(ATinv_df_duT_v[isrc,:])), adjoint=True)
+
+        # del df_duT_v, ATinv_df_duT_v, A, Asubdiag
+        if AdiagTinv is not None:
+            AdiagTinv.clean()
+
+        return Utils.mkvc(JTv).astype(float)
+
+
+
+    def getSourceTerm(self, tInd):
+
+        Srcs = self.survey.srcList
+
+        if self._eqLocs is 'FE':
+            S_m = np.zeros((self.mesh.nF,len(Srcs)))
+            S_e = np.zeros((self.mesh.nE,len(Srcs)))
+        elif self._eqLocs is 'EF':
+            S_m = np.zeros((self.mesh.nE,len(Srcs)))
+            S_e = np.zeros((self.mesh.nF,len(Srcs)))
+
+        for i, src in enumerate(Srcs):
+            smi, sei = src.eval(self, self.times[tInd])
+            S_m[:,i] = S_m[:,i] + smi
+            S_e[:,i] = S_e[:,i] + sei
+
+        return S_m, S_e
+
+    def getInitialFields(self):
+
+        Srcs = self.survey.srcList
+
+        if self._fieldType is 'b' or self._fieldType is 'j':
+            ifields = np.zeros((self.mesh.nF, len(Srcs)))
+        elif self._fieldType is 'e' or self._fieldType is 'h':
+            ifields = np.zeros((self.mesh.nE, len(Srcs)))
+
+        for i,src in enumerate(Srcs):
+            ifields[:,i] = ifields[:,i] + getattr(src, '%sInitial'%self._fieldType, None)(self)
+
+        return ifields
+
+    def getInitialFieldsDeriv(self, src, v, adjoint=False):
+
+        if adjoint is False:
+            if self._fieldType is 'b' or self._fieldType is 'j':
+                ifieldsDeriv = np.zeros(self.mesh.nF)
+            elif self._fieldType is 'e' or self._fieldType is 'h':
+                ifieldsDeriv = np.zeros(self.mesh.nE)
+
+        elif adjoint is True:
+            ifieldsDeriv = np.zeros(self.mapping.nP)
+
+        ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+
+            # ifieldsDeriv = Utils.mkvc(getattr(src, '%sInitialDeriv'%self._fieldType, None)(self,v,adjoint)) + ifieldsDeriv
+            # ifieldsDeriv = self.getAdiagDeriv(None, u, v, adjoint)
+            # ifieldsDeriv = ifieldsDeriv.sum()
+
+        return ifieldsDeriv
+
+
+##########################################################################################
+################################ E-B Formulation #########################################
+##########################################################################################
+
+# ------------------------------- Problem_b -------------------------------------------- #
+
+class Problem_b(BaseTDEMProblem):
+    """
+    Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
+
+    .. math::
+
+        \mathbf{C} \mathbf{e} + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{s_m} \\\\
+        \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}
+
+    where :math:`\mathbf{s_e}` is an integrated quantity, we eliminate :math:`\mathbf{e}` using
+
+    .. math::
+        \mathbf{e} = \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}
+
+    to obtain a second order semi-discretized system in :math:`\mathbf{b}`
+
+    .. math::
+        \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}  + \\frac{\partial \mathbf{b}}{\partial t} = \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    and moving everything except the time derivative to the rhs gives
+
+    .. math::
+        \\frac{\partial \mathbf{b}}{\partial t} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e} + \mathbf{s_m}
+
+    For the time discretization, we use backward euler. To solve for the :math:`n+1`th time step, we have
+
+    .. math::
+        \\frac{\mathbf{b}^{n+1} - \mathbf{b}^{n}}{\mathbf{dt}} = -\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b}^{n+1} + \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1}
+
+    re-arranging to put :math:`\mathbf{b}^{n+1}` on the left hand side gives
+
+    .. math::
+        (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}) \mathbf{b}^{n+1} = \mathbf{b}^{n} + \mathbf{dt}(\mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{s_e}^{n+1} + \mathbf{s_m}^{n+1})
+
+    :param Mesh mesh: mesh
+    :param Mapping mapping: mapping
+    """
+
+    _fieldType = 'b'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_b
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        .. math::
+            (\mathbf{I} + \mathbf{dt} \mathbf{C} \mathbf{M_{\sigma}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f})
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+        I = Utils.speye(self.mesh.nF)
+
+        A = 1./dt * I + ( C * ( MeSigmaI * (C.T * MfMui ) ) )
+
+        if self._makeASymmetric is True:
+            return MfMui.T * A
+        return A
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        C = self.mesh.edgeCurl
+        MeSigmaIDeriv = lambda x: self.MeSigmaIDeriv(x)
+        MfMui = self.MfMui
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = MfMui * v
+            return  MeSigmaIDeriv(C.T * ( MfMui * u )).T * ( C.T * v )
+
+        ADeriv = ( C * ( MeSigmaIDeriv(C.T * ( MfMui * u )) * v ) )
+        if self._makeASymmetric is True:
+            return MfMui.T * ADeriv
+        return ADeriv
+
+
+    def getAsubdiag(self, tInd):
+
+        dt = self.timeSteps[tInd]
+        MfMui = self.MfMui
+        Asubdiag = - 1./dt * sp.eye(self.mesh.nF)
+
+        if self._makeASymmetric is True:
+            return MfMui.T * Asubdiag
+
+        return Asubdiag
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        return Zero() * v
+
+
+
+    def getRHS(self, tInd):
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MfMui = self.MfMui
+
+        S_m, S_e = self.getSourceTerm(tInd)
+
+        rhs =  (C * (MeSigmaI * S_e) + S_m)
+        if self._makeASymmetric is True:
+            return MfMui.T * rhs
+        return rhs
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+
+        C = self.mesh.edgeCurl
+        MeSigmaI = self.MeSigmaI
+        MeSigmaIDeriv = lambda u: self.MeSigmaIDeriv(u)
+        MfMui = self.MfMui
+
+        _, S_e = src.eval(tInd, self)
+        S_mDeriv, S_eDeriv = src.evalDeriv(self.times[tInd], self, adjoint=adjoint)
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = self.MfMui * v
+            if isinstance(S_e, Utils.Zero):
+                MeSigmaIDerivT_v = Utils.Zero()
+            else:
+                MeSigmaIDerivT_v = MeSigmaIDeriv(S_e).T * v
+            RHSDeriv = MeSigmaIDerivT_v + S_eDeriv( MeSigmaI.T *  ( C.T * v ) ) + S_mDeriv(v)
+            return RHSDeriv
+
+        if isinstance(S_e, Utils.Zero):
+            MeSigmaIDeriv_v = Utils.Zero()
+        else:
+            MeSigmaIDeriv_v = MeSigmaIDeriv(S_e) * v
+
+        RHSDeriv = (C * (MeSigmaIDeriv_v + MeSigmaI * S_eDeriv(v) + S_mDeriv(v)))
+
+        if self._makeASymmetric is True:
+            return self.MfMui.T * RHSDeriv
+        return RHSDeriv
+
+
+# ------------------------------- Problem_e -------------------------------------------- #
+
+class Problem_e(BaseTDEMProblem):
+
+    _fieldType = 'e'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields_e
+    surveyPair = SurveyTDEM
+
+    def __init__(self, mesh, mapping=None, **kwargs):
+        BaseTDEMProblem.__init__(self, mesh, mapping=mapping, **kwargs)
+
+    def getAdiag(self, tInd):
+        """
+        System matrix at a given time index
+
+        """
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigma = self.MeSigma
+
+        return C.T * ( MfMui * C ) + 1./dt * MeSigma
+
+
+    def getAdiagDeriv(self, tInd, u, v, adjoint=False):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigmaDeriv = self.MeSigmaDeriv(u)
+
+        if adjoint:
+            return 1./dt * MeSigmaDeriv.T * v
+
+        return 1./dt * MeSigmaDeriv * v
+
+
+    def getAsubdiag(self, tInd):
+        assert tInd >= 0 and tInd < self.nT
+
+        dt = self.timeSteps[tInd]
+
+        return - 1./dt * self.MeSigma
+
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        dt = self.timeSteps[tInd]
+
+        if adjoint:
+            return - 1./dt * self.MeSigmaDeriv(u).T * v
+
+        return - 1./dt * self.MeSigmaDeriv(u) * v
+
+    def getRHS(self, tInd):
+        return Zero()
+
+    def getRHSDeriv(self, tInd, src, v, adjoint=False):
+        return Zero()
+
+
+
+
+

SimPEG/EM/TDEM/__init__.py

@@ -1,3 +1,3 @@
-from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
-from BaseTDEM import BaseTDEMProblem, FieldsTDEM
-from TDEM_b import ProblemTDEM_b
+from TDEM import BaseTDEMProblem, Problem_b, Problem_e
+from FieldsTDEM import Fields, Fields_b
+from SurveyTDEM import Survey, Src, Rx

SimPEG/EM/TDEM_old/SurveyTDEM.py

@@ -0,0 +1,199 @@
+from SimPEG import Utils, Survey, np
+from SimPEG.Survey import BaseSurvey
+from SimPEG.EM.Utils import *
+from BaseTDEM import FieldsTDEM
+import SrcTDEM as Src
+
+class RxTDEM(Survey.BaseTimeRx):
+
+    knownRxTypes = {
+                    'ex':['e', 'Ex', 'N'],
+                    'ey':['e', 'Ey', 'N'],
+                    'ez':['e', 'Ez', 'N'],
+
+                    'bx':['b', 'Fx', 'N'],
+                    'by':['b', 'Fy', 'N'],
+                    'bz':['b', 'Fz', 'N'],
+
+                    'dbxdt':['b', 'Fx', 'CC'],
+                    'dbydt':['b', 'Fy', 'CC'],
+                    'dbzdt':['b', 'Fz', 'CC'],
+                   }
+
+    def __init__(self, locs, times, rxType):
+        Survey.BaseTimeRx.__init__(self, locs, times, rxType)
+
+    @property
+    def projField(self):
+        """Field Type projection (e.g. e b ...)"""
+        return self.knownRxTypes[self.rxType][0]
+
+    @property
+    def projGLoc(self):
+        """Grid Location projection (e.g. Ex Fy ...)"""
+        return self.knownRxTypes[self.rxType][1]
+
+    @property
+    def projTLoc(self):
+        """Time Location projection (e.g. CC N)"""
+        return self.knownRxTypes[self.rxType][2]
+
+    def getTimeP(self, timeMesh):
+        """
+            Returns the time projection matrix.
+
+            .. note::
+
+                This is not stored in memory, but is created on demand.
+        """
+        if self.rxType in ['dbxdt','dbydt','dbzdt']:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)*timeMesh.faceDiv
+        else:
+            return timeMesh.getInterpolationMat(self.times, self.projTLoc)
+
+    def eval(self, src, mesh, timeMesh, u):
+        P = self.getP(mesh, timeMesh)
+        u_part = Utils.mkvc(u[src, self.projField, :])
+        return P*u_part
+
+    def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
+        P = self.getP(mesh, timeMesh)
+
+        if not adjoint:
+            return P * Utils.mkvc(v[src, self.projField, :])
+        elif adjoint:
+            return P.T * v[src, self]
+
+
+class SrcTDEM(Survey.BaseSrc):
+    rxPair = RxTDEM
+    radius = None
+
+    def getInitialFields(self, mesh):
+        F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
+        return F0
+
+    def getJs(self, mesh, time):
+        return None
+
+
+class SrcTDEM_VMD_MVP(SrcTDEM):
+
+    def __init__(self,rxList,loc,waveformType="STEPOFF"):
+        self.loc = loc
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Vertical magnetic dipole, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+            else:
+                raise Exception('Unknown mesh for VMD')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
+        else:
+            raise Exception('Unknown mesh for VMD')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SrcTDEM_CircularLoop_MVP(SrcTDEM):
+    def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
+        self.loc = loc
+        self.radius = radius
+        self.waveformType = waveformType
+        SrcTDEM.__init__(self,rxList)
+
+    def getInitialFields(self, mesh):
+        """Circular Loop, magnetic vector potential"""
+        if self.waveformType == "STEPOFF":
+            print ">> Step waveform: Non-zero initial condition"
+            if mesh._meshType is 'CYL':
+                if mesh.isSymmetric:
+                    MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+                else:
+                    raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+            elif mesh._meshType is 'TENSOR':
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+            else:
+                raise Exception('Unknown mesh for CircularLoop')
+            return {"b": mesh.edgeCurl*MVP}
+        elif self.waveformType == "GENERAL":
+            print ">> General waveform: Zero initial condition"
+            return {"b": np.zeros(mesh.nF)}
+        else:
+            raise NotImplementedError("Only use STEPOFF or GENERAL")
+
+    def getMeS(self, mesh, MfMui):
+        if mesh._meshType is 'CYL':
+            if mesh.isSymmetric:
+                MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius)
+            else:
+                raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
+        elif mesh._meshType is 'TENSOR':
+            MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius)
+        else:
+            raise Exception('Unknown mesh for CircularLoop')
+        return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
+
+
+class SurveyTDEM(Survey.BaseSurvey):
+    """
+        docstring for SurveyTDEM
+    """
+    srcPair = SrcTDEM
+
+    def __init__(self, srcList, **kwargs):
+        # Sort these by frequency
+        self.srcList = srcList
+        Survey.BaseSurvey.__init__(self, **kwargs)
+
+    def projectFields(self, u):
+        data = Survey.Data(self)
+        for src in self.srcList:
+            for rx in src.rxList:
+                data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
+        return data
+
+    def projectFieldsDeriv(self, u, v=None, adjoint=False):
+        assert v is not None, 'v to multiply must be provided.'
+
+        if not adjoint:
+            data = Survey.Data(self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
+            return data
+        else:
+            f = FieldsTDEM(self.mesh, self)
+            for src in self.srcList:
+                for rx in src.rxList:
+                    Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
+                    Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
+                    if rx.projField not in f: # first time we are projecting
+                        f[src, rx.projField, :] = Ptv
+                    else: # there are already fields, so let's add to them!
+                        f[src, rx.projField, :] += Ptv
+            return f
+
+

SimPEG/EM/TDEM_old/__init__.py

@@ -0,0 +1,3 @@
+from SurveyTDEM import * #SurveyTDEM, RxTDEM, SrcTDEM
+from BaseTDEM import BaseTDEMProblem, FieldsTDEM
+from TDEM_b import ProblemTDEM_b

SimPEG/Examples/EM_TDEM_1D_Inversion.py

@@ -42,10 +42,10 @@ def run(plotIt=True):
 
 
     rxOffset=1e-3
-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], np.array([0., 0., 80]))
-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80]))
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
 
     prb.Solver = SolverLU
     prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
@@ -53,9 +53,9 @@ def run(plotIt=True):
 
     # create observed data
     std = 0.05
-    
+
     survey.dobs = survey.makeSyntheticData(mtrue,std)
-    survey.std = std 
+    survey.std = std
     survey.eps = 1e-5*np.linalg.norm(survey.dobs)
 
     if plotIt:

SimPEG/Examples/MT_1D_ForwardAndInversion.py

@@ -1,11 +1,9 @@
 import SimPEG as simpeg
 import numpy as np
-from SimPEG import NSEM
+import SimPEG.MT as MT
 from scipy.constants import mu_0
 import matplotlib.pyplot as plt
 
-np.random.seed(1983)
-
 def run(plotIt=True):
     """
         MT: 1D: Inversion
@@ -19,13 +17,13 @@ def run(plotIt=True):
     ## Setup the forward modeling
     # Setting up 1D mesh and conductivity models to forward model data.
     # Frequency
-    nFreq = 26
-    freqs = np.logspace(2,-3,nFreq)
+    nFreq = 31
+    freqs = np.logspace(3,-3,nFreq)
     # Set mesh parameters
-    ct = 10
-    air = simpeg.Utils.meshTensor([(ct,25,1.4)])
+    ct = 20
+    air = simpeg.Utils.meshTensor([(ct,16,1.4)])
     core = np.concatenate( (  np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
-    bot = simpeg.Utils.meshTensor([(core[0],25,-1.4)])
+    bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
     x0 = -np.array([np.sum(np.concatenate((core,bot)))])
     # Make the model
     m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -35,7 +33,7 @@ def run(plotIt=True):
     layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
     layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
     # Set the conductivity values
-    sig_half = 1e-2
+    sig_half = 2e-3
     sig_air = 1e-8
     sig_layer1 = .2
     sig_layer2 = .2
@@ -52,38 +50,38 @@ def run(plotIt=True):
     m_0 = np.log(sigma_0[active])
 
     # Set the mapping
-    actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
+    actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
     mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
 
     ## Setup the layout of the survey, set the sources and the connected receivers
     # Receivers
     rxList = []
     for rxType in ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([-0.5]),2).T,rxType))
+        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
     # Source list
     srcList =[]
     for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
     # Make the survey
-    survey = NSEM.Survey(srcList)
+    survey = MT.Survey(srcList)
     survey.mtrue = m_true
 
     ## Set the problem
-    problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
+    problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
     problem.pair(survey)
 
     ## Forward model data
     # Project the data
     survey.dtrue = survey.dpred(m_true)
-    survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
+    survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
 
     if plotIt:
-        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
         fig.suptitle('Target - smooth true')
 
 
     # Assign uncertainties
-    std = 0.025 # 5% std
+    std = 0.05 # 5% std
     survey.std = np.abs(survey.dobs*std)
     # Assign the data weight
     Wd = 1./survey.std
@@ -92,33 +90,30 @@ def run(plotIt=True):
     # Define a counter
     C =  simpeg.Utils.Counter()
     # Set the optimization
-    opt = simpeg.Optimization.ProjectedGNCG(maxIter = 25)
+    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
     opt.counter = C
-    opt.lower = np.log(1e-4)
-    opt.upper = np.log(5)
-    opt.LSshorten = 0.1
+    opt.LSshorten = 0.5
     opt.remember('xc')
     # Data misfit
     dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
     dmis.Wd = Wd
     # Regularization - with a regularization mesh
-    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[active]],m1d.x0)
+    regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
     reg = simpeg.Regularization.Tikhonov(regMesh)
     reg.mrefInSmooth = True
-    reg.alpha_s = 1e-1
+    reg.alpha_s = 1e-7
     reg.alpha_x = 1.
-
     # Inversion problem
     invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
     invProb.counter = C
     # Beta cooling
     beta = simpeg.Directives.BetaSchedule()
-    beta.coolingRate = 4.
-    beta.coolingFactor = 4.
-    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
-    betaest.beta0 = 1.
+    beta.coolingRate = 4
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
     targmis = simpeg.Directives.TargetMisfit()
     targmis.target = survey.nD
+    saveModel = simpeg.Directives.SaveModelEveryIteration()
+    saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
     # Create an inversion object
     inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
 
@@ -126,9 +121,8 @@ def run(plotIt=True):
     mopt = inv.run(m_0)
 
     if plotIt:
-        fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
         fig.suptitle('Target - smooth true')
-        fig.axes[0].set_ylim([-10000,500])
         plt.show()
 
 if __name__ == '__main__':

SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py

@@ -1,428 +0,0 @@
-from scipy.constants import epsilon_0, mu_0
-import matplotlib.pyplot as plt
-import numpy as np
-from ipywidgets import *
-from SimPEG.EM.Utils import k, omega
-
-"""
-MT1D: n layered earth problem
-*****************************
-
-Author: Thibaut Astic
-Contact: thast@eos.ubc.ca
-Date: January 2016
-
-This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
-
-We start by looking at Maxwell's equations in the electric
-field \\\(\\\mathbf{E}\\) and the magnetic flux
-\\\(\\\mathbf{H}\\) to write the wave equations
-\\(\\ \nabla ^2  \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
-\\(\\ \nabla ^2  \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
-
-Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
-\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
-\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
-
-With U and D the Up and Down components of the E-field.
-
-The iteration from one layer to another is ensure by:
-
-\\(\\  \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
-    P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
-
-And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
-and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
-
-The layer 0 is assumed to be the air layer.
-
-"""
-
-#Define a frquency range for a survey
-frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
-
-#Functions to create random physical Perties for a n-layered earth
-thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
-                                                     np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
-                                                            ,decimals =1))
-
-sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
-                                               np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
-                                                      ,decimals=3))
-
-mu  = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
-                                             np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
-                                                    ,decimals=1))
-
-eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
-                                               np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
-                                                                ,decimals=1))
-
-#Evaluate Impedance Z of a layer
-ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
-
-#Complex Cole-Cole Conductivity - EM utils
-PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
-
-#Converted thickness array into top of layer array
-top = lambda thick: np.cumsum(thick)
-
-#Propagation Matrix and theirs inverses
-
-#matrix T for transition of Up and Down components accross a layer
-T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
-
-Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
-
-#transition of Up and Down components accross a layer
-UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
-
-
-#matrix P relating Up and Down components with E and H fields
-P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
-
-Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
-
-
-#Time Variation of E and H
-E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
-H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
-
-#Plot the configuration of the problem
-def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
-
-    topn = top(thick)
-    widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
-
-    ax.set_ylim([z.min(),z.max()])
-    ax.set_xlim([-widthg,widthg])
-
-    ax.set_ylabel("Depth (m)", fontsize=16.)
-    ax.yaxis.tick_right()
-    ax.yaxis.set_label_position("right")
-
-    #define filling for the different layers
-    hatches=['/'  , '+', 'x', '|' , '\\', '-'  , 'o' , 'O' , '.' , '*' ]
-
-    #Write the physical properties of air
-    ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
-            xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
-            xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
-            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
-            xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\mu_r$= %1i")%(mu[0]),
-            xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
-            xytext=(0, -np.abs(z.max())/3.), textcoords='data',
-            fontsize=14.)
-
-    #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
-    for i in range(1,len(topn)-1,1):
-        if topn[i] == topn[i+1]:
-            pass
-        else:
-            ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
-                xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
-                xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
-                xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
-                xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
-                fontsize=14.)
-
-            ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
-            ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
-
-    #Write the physical properties of the n-th layer and fill it with pattern
-    ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
-    ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
-
-    ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
-            xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
-            xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
-            xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
-            xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
-            fontsize=14.)
-
-    #plot Trees!
-    ax.annotate("",
-            xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(1.2*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
-            )
-
-    ax.annotate("",
-            xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
-            xytext=(1.5*widthg/2., 0.), textcoords='data',
-            arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
-            )
-
-
-    ax.invert_yaxis()
-
-    return ax
-
-#Propagate Up and Down component for a certain frequency & evaluate E and H field
-
-def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
-
-    sigcm = np.zeros_like(sig,dtype='complex_')
-
-    for j in range(1,len(sig)):
-        sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
-
-    K = k(f, sigcm, mu, eps)
-    Z = ImpZ(f,mu,K)
-
-    EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
-    UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
-
-    UD[1,-1] = 1.
-
-    for i in range(-2,-(n+2),-1):
-
-        UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
-        UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
-
-    for j in range(0,n+1):
-        EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
-
-    return UD, EH, Z ,K
-
-
-#Evaluate the apparent resistivity and phase for a frequency range
-def appres(F,H,sig,chg,taux,c,mu,eps,n):
-
-    Res = np.zeros_like(F)
-    Phase = np.zeros_like(F)
-    App_ImpZ= np.zeros_like(F,dtype='complex_')
-
-    for i in range(0,len(F)):
-
-        UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
-
-        App_ImpZ[i] = EH[0,1]/EH[1,1]
-
-        Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
-        Phase[i] = np.angle(App_ImpZ[i], deg = True)
-
-    return Res,Phase
-
-#Evaluate Up, Down components, E and H field, for a frequency range,
-#a discretized depth range and a time range (use to calculate envelope)
-def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
-
-    topc = top(H)
-
-    layer = np.zeros(len(zsample),dtype=np.int)-1
-
-    Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
-    UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
-
-    for i in range(0,n+1,1):
-        layer = layer+(zsample>=topc[i])*1
-
-    for j in range(0,len(F)):
-
-        UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
-
-        for p in range(0,len(zsample)):
-
-            UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
-
-            for q in range(0,len(tsample)):
-
-                Exzt[p,q]  = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
-                Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
-                Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
-                Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
-
-    return  Exzt,Hyzt,Uz,Dz,UDaux,layer
-
-
-#Function to Plot Apparent Resistivity and Phase
-def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
-
-    Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
-
-    fig,ax = plt.subplots(1,2,figsize=(16,10))
-
-    ax[0].scatter(Res,F,color='black')
-    ax[0].set_xscale('Log')
-    ax[0].set_yscale('Log')
-    ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
-    ax[0].set_ylim([F.min(),F.max()])
-    ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
-    ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
-    ax[0].grid(which='major')
-
-    ax0 = ax[0].twiny()
-
-    ax0.set_xlim([0.,90.])
-    ax0.set_ylim([F.min(),F.max()])
-    ax0.scatter(Phase,F,color='purple')
-    ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
-
-    zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
-
-    ax[0].tick_params(labelsize=16)
-    ax[1].tick_params(labelsize=16)
-    ax0.tick_params(labelsize=16)
-
-    if PlotEnvelope:
-
-        widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
-        fenvelope1n=np.ones(100)*fenvelope
-        ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
-
-        tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
-        Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
-
-        ax1=ax[1].twiny()
-
-        ax[1].tick_params(labelsize=16)
-        ax1.tick_params(labelsize=16)
-
-        ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
-
-        ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
-
-        ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
-                      np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
-                      color='blue', alpha=0.1)
-
-        ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
-                      np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
-                      color='red', alpha=0.1)
-
-        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
-        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
-        ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
-    else:
-        print 'No envelop (if True, might be slow)'
-        ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
-        ax[1].get_xaxis().set_ticks([])
-
-    plt.show()
-
-#Interactive MT for Notebook
-def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
-
-    frangn=frange(-5,5,100.)
-    sig3= np.array([0.,0.001,0.1, 0.001])
-    thick3 = np.array([120000.,50.,50.])
-    eps3=np.array([1.,1.,1.,1])
-    mu3=np.array([1.,1.,1.,1])
-    chg3=np.array([0.,0.1,0.,0.2])
-    chg3_0=np.array([0.,0.1,0.,0.])
-    taux3=np.array([0.,0.1,0.,0.1])
-    c3=np.array([1.,1.,1.,1.])
-
-    sig3[1]=sigl1
-    sig3[1]=10.**sig3[1]
-    sig3[2]=sigl2
-    sig3[2]=10.**sig3[2]
-    sig3[3]=sigl3
-    sig3[3]=10.**sig3[3]
-    mu3[1]=mul1
-    mu3[2]=mul2
-    mu3[3]=mul3
-    eps3[1]=epsl1
-    eps3[2]=epsl2
-    eps3[3]=epsl3
-    thick3[1]=h1
-    thick3[2]=h2
-
-    PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
-
-
-def run(n=3,plotIt=True):
-    # something to make a plot
-
-    F = frange(-5.,5.,20)
-    H = thick(50.,100.,n)
-    sign = sig(-5.,0.,n)
-    mun = mu(1.,2.,n)
-    epsn = eps(1.,9.,n)
-    chg = np.zeros_like(sign)
-    taux = np.zeros_like(sign)
-    c = np.zeros_like(sign)
-
-    Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
-
-    if plotIt:
-
-        PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
-
-    return Res, Phase
-
-if __name__ == '__main__':
-    run()
-
-
-
-
-

SimPEG/Examples/MT_3D_Foward.py

@@ -2,7 +2,7 @@
 
 # Import
 import SimPEG as simpeg
-from SimPEG import NSEM
+from SimPEG import MT
 import numpy as np
 try:
     from pymatsolver import MumpsSolver as Solver
@@ -37,16 +37,16 @@ def run(plotIt=True, nFreq=1):
     for loc in rx_loc:
         # NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
         for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
     # Source list
     srcList =[]
     for freq in np.logspace(3,-3,nFreq):
-        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
+        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
     # Survey MT
-    survey = NSEM.Survey(srcList)
+    survey = MT.Survey(srcList)
 
     ## Setup the problem object
-    problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
+    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
     problem.pair(survey)
     problem.Solver = Solver
 
@@ -55,7 +55,7 @@ def run(plotIt=True, nFreq=1):
     dataVec = survey.eval(fields)
 
     # Make the data
-    mtData = NSEM.Data(survey,dataVec)
+    mtData = MT.Data(survey,dataVec)
     # Add plots
     if plotIt:
         pass

SimPEG/Examples/__init__.py

@@ -1,29 +1,28 @@
 # Run this file to add imports.
 
 ##### AUTOIMPORTS #####
-import EM_FDEM_1D_Inversion
-import Mesh_QuadTree_Creation
-import EM_TDEM_1D_Inversion
-import Mesh_QuadTree_FaceDiv
-import Mesh_Tensor_Creation
-import FLOW_Richards_1D_Celia1990
+import DC_Analytic_Dipole
 import DC_Forward_PseudoSection
-import Mesh_Operators_CahnHilliard
-import Mesh_Basic_Types
+import EM_FDEM_1D_Inversion
+import EM_FDEM_Analytic_MagDipoleWholespace
+import EM_Schenkel_Morrison_Casing
+import EM_TDEM_1D_Inversion
+import FLOW_Richards_1D_Celia1990
 import Inversion_IRLS
 import Inversion_Linear
-import EM_Schenkel_Morrison_Casing
-import MT_3D_Foward
 import Mesh_Basic_ForwardDC
-import MT_1D_ForwardAndInversion
-import Utils_surface2ind_topo
-import MT_1D_analytic_nlayer_Earth
-import EM_FDEM_Analytic_MagDipoleWholespace
 import Mesh_Basic_PlotImage
-import DC_Analytic_Dipole
+import Mesh_Basic_Types
+import Mesh_Operators_CahnHilliard
+import Mesh_QuadTree_Creation
+import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
+import Mesh_Tensor_Creation
+import MT_1D_ForwardAndInversion
+import MT_3D_Foward
+import Utils_surface2ind_topo
 
-__examples__ = ["EM_FDEM_1D_Inversion", "Mesh_QuadTree_Creation", "EM_TDEM_1D_Inversion", "Mesh_QuadTree_FaceDiv", "Mesh_Tensor_Creation", "FLOW_Richards_1D_Celia1990", "DC_Forward_PseudoSection", "Mesh_Operators_CahnHilliard", "Mesh_Basic_Types", "Inversion_IRLS", "Inversion_Linear", "EM_Schenkel_Morrison_Casing", "MT_3D_Foward", "Mesh_Basic_ForwardDC", "MT_1D_ForwardAndInversion", "Utils_surface2ind_topo", "MT_1D_analytic_nlayer_Earth", "EM_FDEM_Analytic_MagDipoleWholespace", "Mesh_Basic_PlotImage", "DC_Analytic_Dipole", "Mesh_QuadTree_HangingNodes"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_ForwardDC", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "Utils_surface2ind_topo"]
 
 ##### AUTOIMPORTS #####
 

SimPEG/Examples/simpegEMpaper_example_1.py

@@ -0,0 +1,157 @@
+from SimPEG import np, Mesh, Maps, Utils, DataMisfit, Regularization, Optimization, Inversion, InvProblem, Directives
+from SimPEG import SolverLU
+from SimPEG.EM import FDEM, TDEM, mu_0
+import matplotlib.pyplot as plt
+import matplotlib
+matplotlib.rcParams['font.size'] = 14
+
+def run(plotIt=True):
+    # Set up cylindrically symmeric mesh
+    cs, ncx, ncz, npad = 10., 15, 25, 13  # padded cyl mesh
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hz], '00C')
+
+    # Conductivity model
+    layerz = np.r_[-200., -100.]
+    layer = (mesh.vectorCCz>=layerz[0]) & (mesh.vectorCCz<=layerz[1])
+    active = mesh.vectorCCz<0.
+    sig_half  = 1e-2  # Half-space conductivity
+    sig_air   = 1e-8  # Air conductivity
+    sig_layer = 5e-2  # Layer conductivity
+    sigma = np.ones(mesh.nCz)*sig_air
+    sigma[active] = sig_half
+    sigma[layer] = sig_layer
+
+    # Mapping
+    actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
+    mtrue = np.log(sigma[active])
+
+    # FDEM problem & survey
+    rxlocs = Utils.ndgrid([np.r_[50.], np.r_[0], np.r_[0.]])
+    bzi = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'real')
+    bzr = FDEM.Rx.Point_bSecondary(rxlocs, 'z', 'imag')
+
+    freqs = np.logspace(2, 3, 5)
+    srcLoc = np.array([0., 0., 0.])
+
+    print 'min skin depth = ', 500./np.sqrt(freqs.max() * sig_half), 'max skin depth = ', 500./np.sqrt(freqs.min() * sig_half)
+    print 'max x ', mesh.vectorCCx.max(), 'min z ', mesh.vectorCCz.min(), 'max z ', mesh.vectorCCz.max()
+
+    srcList = []
+    [srcList.append(FDEM.Src.MagDipole([bzr, bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
+
+    surveyFD = FDEM.Survey(srcList)
+    prbFD = FDEM.Problem3D_b(mesh, mapping=mapping)
+    prbFD.pair(surveyFD)
+    std = 0.03
+    surveyFD.makeSyntheticData(mtrue, std)
+    surveyFD.eps = np.linalg.norm(surveyFD.dtrue)*1e-5
+
+    # FDEM inversion
+    np.random.seed(1)
+    dmisfit = DataMisfit.l2_DataMisfit(surveyFD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=10.)
+    invProb.beta = 1.
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta,target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbFD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.tolG = 1e-10
+    opt.eps = 1e-10
+    opt.remember('xc')
+    moptFD = inv.run(m0)
+
+    # TDEM problem
+    times = np.logspace(-4, np.log10(2e-3), 10)
+    print 'min diffusion distance ', 1.28*np.sqrt(times.min()/(sig_half*mu_0)), 'max diffusion distance ', 1.28*np.sqrt(times.max()/(sig_half*mu_0))
+    rx = TDEM.Rx(rxlocs, times, 'bz')
+    src = TDEM.Src.MagDipole([rx], waveform=TDEM.Src.StepOffWaveform(), loc=srcLoc) # same src location as FDEM problem
+
+    surveyTD = TDEM.Survey([src])
+    prbTD = TDEM.Problem_b(mesh, mapping=mapping)
+    prbTD.timeSteps = [(5e-5, 10),(1e-4, 10),(5e-4, 10)]
+    prbTD.pair(surveyTD)
+    prbTD.Solver = SolverLU
+
+    std = 0.03
+    surveyTD.makeSyntheticData(mtrue, std)
+    surveyTD.std = std
+    surveyTD.eps = np.linalg.norm(surveyTD.dtrue)*1e-5
+
+    # TDEM inversion
+    dmisfit = DataMisfit.l2_DataMisfit(surveyTD)
+    regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
+    reg = Regularization.Simple(regMesh)
+    opt = Optimization.InexactGaussNewton(maxIterCG=10, maxIter=4)
+    invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
+
+    # Inversion Directives
+    beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=3)
+    invProb.beta = 1.
+    # betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1.)
+    target = Directives.TargetMisfit()
+
+    inv = Inversion.BaseInversion(invProb, directiveList=[beta, target])
+    m0 = np.log(np.ones(mtrue.size)*sig_half)
+    reg.alpha_s = 5e-1
+    reg.alpha_x = 1.
+    prbTD.counter = opt.counter = Utils.Counter()
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+    moptTD = inv.run(m0)
+
+    if plotIt:
+        fig, ax = plt.subplots(1,1, figsize = (4, 6))
+        plt.semilogx(sigma[active], mesh.vectorCCz[active], 'k-', lw=2)
+        plt.semilogx(np.exp(moptFD), mesh.vectorCCz[active], 'ko', ms=3)
+        plt.semilogx(np.exp(moptTD), mesh.vectorCCz[active], 'k*')
+        ax.set_ylim(-1000, 0)
+        ax.set_xlim(5e-3, 1e-1)
+
+        ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
+        ax.set_ylabel('Depth (m)', fontsize = 14)
+        ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
+        plt.legend(['True', 'Pred (FD)', 'Pred (TD)'], fontsize=13, loc=4)
+        plt.show()
+
+        fig = plt.figure(figsize = (10*1.3, 5*1.3))
+        ax2 = plt.subplot(122)
+        ax2.plot(times, surveyTD.dobs, 'k-', lw=2)
+        ax2.plot(times, surveyTD.dpred(moptTD), 'ko', ms=4)
+        ax2.set_xscale('log')
+        ax2.set_yscale('log')
+        ax2.set_xlim(times.min(), times.max())
+        ax1 = plt.subplot(121)
+        ax1.plot(freqs, -surveyFD.dobs[::2], 'k-', lw=2)
+        ax1.plot(freqs, -surveyFD.dobs[1::2], 'k--', lw=2)
+        dpredFD = surveyFD.dpred(moptTD)
+        ax1.plot(freqs, -dpredFD[::2], 'ko', ms=4)
+        ax1.plot(freqs, -dpredFD[1::2], 'k+', markeredgewidth=2., ms=10)
+        ax1.set_xscale('log')
+        ax1.set_yscale('log')
+        ax2.set_xlabel('Time (s)', fontsize = 14)
+        ax1.set_xlabel('Frequency (Hz)', fontsize = 14)
+        ax1.set_ylabel('Vertical magnetic field (T)', fontsize = 14)
+        ax2.grid(True,which='minor')
+        ax1.grid(True,which='minor')
+        ax2.set_title("(b) TD observed vs. predicted", fontsize = 14)
+        ax1.set_title("(a) FD observed vs. predicted", fontsize = 14)
+        ax2.legend(("Obs", "Pred"), fontsize = 12)
+        ax1.legend(("Obs", "Pred (real)", "Pred (imag)"), fontsize = 12, loc=3)
+        ax1.set_xlim(freqs.max(), freqs.min())
+        plt.show()
+
+if __name__ == '__main__':
+    run()

SimPEG/MT/BaseMT.py

@@ -0,0 +1,132 @@
+from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
+from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
+from SurveyMT import Survey, Data
+from FieldsMT import BaseMTFields
+
+
+class BaseMTProblem(BaseFDEMProblem):
+    """
+        Base class for all Natural source problems.
+    """
+
+    def __init__(self, mesh, **kwargs):
+        BaseFDEMProblem.__init__(self, mesh, **kwargs)
+        Utils.setKwargs(self, **kwargs)
+    # Set the default pairs of the problem
+    surveyPair = Survey
+    dataPair = Data
+    fieldsPair = BaseMTFields
+
+    # Set the solver
+    Solver = SimpegSolver
+    solverOpts = {}
+
+    verbose = False
+    # Notes:
+    # Use the forward and devs from BaseFDEMProblem
+    # Might need to add more stuff here.
+
+    ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
+    def Jvec(self, m, v, f=None):
+        """
+        Function to calculate the data sensitivities dD/dm times a vector.
+
+            :param numpy.ndarray m (nC, 1) - conductive model
+            :param numpy.ndarray v (nC, 1) - random vector
+            :param MTfields object (optional) - MT fields object, if not given it is calculated
+            :rtype: MTdata object
+            :return: Data sensitivities wrt m
+        """
+
+        # Calculate the fields
+        if f is None:
+           f= self.fields(m)
+        # Set current model
+        self.curModel = m
+        # Initiate the Jv object
+        Jv = self.dataPair(self.survey)
+
+        # Loop all the frequenies
+        for freq in self.survey.freqs:
+            dA_du = self.getA(freq) #
+
+            dA_duI = self.Solver(dA_du, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                # 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,:]
+                # 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.
+                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 )
+                else:
+                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
+                # Calculate the projection derivatives
+                for rx in src.rxList:
+                    # Get the projection derivative
+                    # v should be of size 2*nE (for 2 polarizations)
+                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
+                    Jv[src, rx] = PDeriv_u(mkvc(du_dm))
+            dA_duI.clean()
+        # Return the vectorized sensitivities
+        return mkvc(Jv)
+
+    def Jtvec(self, m, v, f=None):
+        """
+        Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
+
+            :param numpy.ndarray m (nC, 1) - conductive model
+            :param numpy.ndarray v (nD, 1) - vector
+            :param MTfields object u (optional) - MT fields object, if not given it is calculated
+            :rtype: MTdata object
+            :return: Data sensitivities wrt m
+        """
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        # Ensure v is a data object.
+        if not isinstance(v, self.dataPair):
+            v = self.dataPair(self.survey, v)
+
+        Jtv = np.zeros(m.size)
+
+        for freq in self.survey.freqs:
+            AT = self.getA(freq).T
+
+            ATinv = self.Solver(AT, **self.solverOpts)
+
+            for src in self.survey.getSrcByFreq(freq):
+                ftype = self._fieldType + 'Solution'
+                f_src = f[src, :]
+
+                for rx in src.rxList:
+                    # Get the adjoint evalDeriv
+                    # PTv needs to be nE,
+                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
+                    # Get the
+                    dA_duIT = ATinv * PTv
+                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
+                    dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
+                    # Make du_dmT
+                    if dRHS_dmT is None:
+                        du_dmT = -dA_dmT
+                    else:
+                        du_dmT = -dA_dmT + dRHS_dmT
+                    # Select the correct component
+                    # du_dmT needs to be of size nC,
+                    real_or_imag = rx.projComp
+                    if real_or_imag == 'real':
+                        Jtv +=  du_dmT.real
+                    elif real_or_imag == 'imag':
+                        Jtv +=  -du_dmT.real
+                    else:
+                        raise Exception('Must be real or imag')
+            # Clean the factorization, clear memory.
+            ATinv.clean()
+        return Jtv

SimPEG/MT/FieldsMT.py

@@ -4,21 +4,18 @@ import sys
 from numpy.lib import recfunctions as recFunc
 from SimPEG.EM.Utils import omega
 
-
 ##############
 ### Fields ###
 ##############
-class BaseNSEMFields(Problem.Fields):
-    """Field Storage for a NSEM survey."""
+class BaseMTFields(Problem.Fields):
+    """Field Storage for a MT survey."""
     knownFields = {}
     dtype = complex
 
-###########
-# 1D Fields
-###########
-class Fields1D_ePrimSec(BaseNSEMFields):
+
+class Fields1D_e(BaseMTFields):
     """
-    Fields storage for the 1D NSEM solution.
+    Fields storage for the 1D MT solution.
     """
     knownFields = {'e_1dSolution':'F'}
     aliasFields = {
@@ -31,119 +28,7 @@ class Fields1D_ePrimSec(BaseNSEMFields):
                   }
 
     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, 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)
+        BaseMTFields.__init__(self,mesh,survey,**kwargs)
 
     def _ePrimary(self, eSolution, srcList):
         ePrimary = np.zeros_like(eSolution)
@@ -214,7 +99,7 @@ class Fields1D_eTotal(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -235,18 +120,9 @@ class Fields1D_eTotal(BaseNSEMFields):
         """
         return None
 
-
-###########
-# 2D Fields
-###########
-
-
-###########
-# 3D Fields
-###########
-class Fields3D_ePrimSec(BaseNSEMFields):
+class Fields3D_e(BaseMTFields):
     """
-    Fields storage for the 3D NSEM solution. Labels polarizations by px and py.
+    Fields storage for the 3D MT solution. Labels polarizations by px and py.
 
         :param SimPEG object mesh: The solution mesh
         :param SimPEG object survey: A survey object
@@ -271,7 +147,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
                   }
 
     def __init__(self,mesh,survey,**kwargs):
-        BaseNSEMFields.__init__(self,mesh,survey,**kwargs)
+        BaseMTFields.__init__(self,mesh,survey,**kwargs)
 
     def _e_pxPrimary(self, e_pxSolution, srcList):
         e_pxPrimary = np.zeros_like(e_pxSolution)
@@ -352,7 +228,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         b = (C * e_pxSolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the NSEM problem
+            # 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
@@ -363,7 +239,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         b = (C * e_pySolution)
         for i, src in enumerate(srcList):
             b[:,i] *= - 1./(1j*omega(src.freq))
-            # There is no magnetic source in the NSEM problem
+            # 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
@@ -426,7 +302,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 
@@ -443,7 +319,7 @@ class Fields3D_ePrimSec(BaseNSEMFields):
         """
         Derivative of the fields object wrt u.
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param numpy.ndarray v: random vector of f_sol.size
         This function stacks the fields derivatives appropriately
 

SimPEG/MT/Problem1D/Probs.py

@@ -0,0 +1,291 @@
+from SimPEG.EM.Utils import omega
+from SimPEG import mkvc
+from scipy.constants import mu_0
+from SimPEG.MT.BaseMT import BaseMTProblem
+from SimPEG.MT.SurveyMT import Survey, Data
+from SimPEG.MT.FieldsMT import Fields1D_e
+from SimPEG.MT.Utils.MT1Danalytic import getEHfields
+import numpy as np
+import multiprocessing, sys, time
+
+
+class eForm_psField(BaseMTProblem):
+    """
+    A MT 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 MTproblem.
+    _fieldType = 'e_1d'
+    _eqLocs    = 'EF'
+    _sigmaPrimary = None
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__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(BaseMTProblem):
+    """
+    A MT 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 MTproblem.
+    _fieldType = 'e'
+    _eqLocs    = 'EF'
+
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__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

SimPEG/MT/Problem1D/__init__.py

@@ -0,0 +1 @@
+from Probs import eForm_TotalField, eForm_psField

SimPEG/MT/Problem2D/__init__.py

@@ -0,0 +1 @@
+pass

SimPEG/MT/Problem3D/Probs.py

@@ -0,0 +1,138 @@
+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.MT.BaseMT import BaseMTProblem
+from SimPEG.MT.SurveyMT import Survey, Data
+from SimPEG.MT.FieldsMT import Fields3D_e
+import multiprocessing, sys, time
+
+
+
+class eForm_ps(BaseMTProblem):
+    """
+    A MT 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 MTproblem.
+    _fieldType = 'e'
+    _eqLocs    = 'FE'
+    fieldsPair = Fields3D_e
+    _sigmaPrimary = None
+
+    def __init__(self, mesh, **kwargs):
+        BaseMTProblem.__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
+

SimPEG/MT/Problem3D/__init__.py

@@ -0,0 +1 @@
+from Probs import eForm_ps

SimPEG/MT/SrcMT.py

@@ -11,9 +11,9 @@ import sys
 ###   Sources ###
 #################
 
-class BaseNSEMSrc(FDEMBaseSrc):
+class BaseMTSrc(FDEMBaseSrc):
     '''
-    Sources for the NSEM problem.
+    Sources for the MT problem.
     Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
 
     :param float freq: The frequency of the source
@@ -29,28 +29,28 @@ class BaseNSEMSrc(FDEMBaseSrc):
         FDEMBaseSrc.__init__(self, rxList)
 
 # 1D sources
-class polxy_1DhomotD(BaseNSEMSrc):
+class polxy_1DhomotD(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) for the total Domain.
+    MT source for both polarizations (x and y) for the total Domain.
 
     It calculates fields calculated based on conditions on the boundary of the domain.
     """
     def __init__(self, rxList, freq):
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
 
 
     # TODO: need to add the  primary fields calc and source terms into the problem.
 
 # Need to implement such that it works for all dims.
-class polxy_1Dprimary(BaseNSEMSrc):
+class polxy_1Dprimary(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) given a 1D primary models.
+    MT source for both polarizations (x and y) given a 1D primary models.
     It assigns fields calculated from the 1D model as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigma1d = None
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -86,7 +86,7 @@ class polxy_1Dprimary(BaseNSEMSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
+        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)
@@ -128,15 +128,15 @@ class polxy_1Dprimary(BaseNSEMSrc):
             # v should be nC size
             return MsigmaDeriv * v
 
-class polxy_3Dprimary(BaseNSEMSrc):
+class polxy_3Dprimary(BaseMTSrc):
     """
-    NSEM source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
+    MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
     as fields in the full space of the problem.
     """
     def __init__(self, rxList, freq):
         # assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
         self.sigmaPrimary = None
-        BaseNSEMSrc.__init__(self, rxList, freq)
+        BaseMTSrc.__init__(self, rxList, freq)
         # Hidden property of the ePrimary
         self._ePrimary = None
 
@@ -163,7 +163,7 @@ class polxy_3Dprimary(BaseNSEMSrc):
         Get the electrical field source
         """
         e_p = self.ePrimary(problem)
-        Map_sigma_p = Maps.SurjectVertical1D(problem.mesh)
+        Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
         sigma_p = Map_sigma_p._transform(self.sigma1d)
         # Make mass matrix
         # Note: M(sig) - M(sig_p) = M(sig - sig_p)

SimPEG/MT/SurveyMT.py

@@ -4,7 +4,7 @@ from SimPEG.EM.Utils import omega
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
 from Utils import rec2ndarr
-import SrcNSEM
+import SrcMT
 import sys
 
 #################
@@ -63,9 +63,9 @@ class Rx(SimPEGsurvey.BaseRx):
         '''
         Project the fields to natural source data.
 
-            :param SrcNSEM src: The source of the fields to project
+            :param SrcMT src: The source of the fields to project
             :param SimPEG.Mesh mesh:
-            :param FieldsNSEM f: Natural source fields object to project
+            :param FieldsMT f: Natural source fields object to project
         '''
 
         ## NOTE: Assumes that e is on t
@@ -143,9 +143,9 @@ class Rx(SimPEGsurvey.BaseRx):
         """
         The derivative of the projection wrt u
 
-        :param NSEMsrc src: NSEM source
+        :param MTsrc src: MT source
         :param TensorMesh mesh: Mesh defining the topology of the problem
-        :param NSEMfields f: NSEM fields object of the source
+        :param MTfields f: MT fields object of the source
         :param numpy.ndarray v: Random vector of size
         """
 
@@ -390,12 +390,12 @@ class Rx(SimPEGsurvey.BaseRx):
 #################
 class Survey(SimPEGsurvey.BaseSurvey):
     """
-        Survey class for NSEM. Contains all the sources associated with the survey.
+        Survey class for MT. Contains all the sources associated with the survey.
 
         :param list srcList: List of sources associated with the survey
 
     """
-    srcPair = SrcNSEM.BaseNSEMSrc
+    srcPair = SrcMT.BaseMTSrc
 
     def __init__(self, srcList, **kwargs):
         # Sort these by frequency
@@ -443,7 +443,7 @@ class Survey(SimPEGsurvey.BaseSurvey):
 #################
 class Data(SimPEGsurvey.Data):
     '''
-    Data class for NSEMdata. Stores the data vector indexed by the survey.
+    Data class for MTdata. Stores the data vector indexed by the survey.
 
     :param SimPEG survey object survey:
     :param v vector of the data in order matching of the survey
@@ -461,7 +461,7 @@ class Data(SimPEGsurvey.Data):
 
     def toRecArray(self,returnType='RealImag'):
         '''
-        Function that returns a numpy.recarray for a SimpegNSEM impedance data object.
+        Function that returns a numpy.recarray for a SimpegMT impedance data object.
 
         :param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
 
@@ -483,7 +483,7 @@ class Data(SimPEGsurvey.Data):
                 locs = np.hstack((np.array([[0.0]]),locs))
             tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
             # np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
-            # Get the type and the value for the DataNSEM object as a list
+            # Get the type and the value for the DataMT object as a list
             typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
             # Insert the values to the temp array
             for nr,(key,val) in enumerate(typeList):
@@ -517,17 +517,17 @@ class Data(SimPEGsurvey.Data):
     @classmethod
     def fromRecArray(cls, recArray, srcType='primary'):
         """
-        Class method that reads in a numpy record array to NSEMdata object.
+        Class method that reads in a numpy record array to MTdata object.
 
         Only imports the impedance data.
 
         """
         if srcType=='primary':
-            src = SrcNSEM.polxy_1Dprimary
+            src = SrcMT.polxy_1Dprimary
         elif srcType=='total':
-            src = SrcNSEM.polxy_1DhomotD
+            src = SrcMT.polxy_1DhomotD
         else:
-            raise NotImplementedError('{:s} is not a valid source type for NSEMdata')
+            raise NotImplementedError('{:s} is not a valid source type for MTdata')
 
         # Find all the frequencies in recArray
         uniFreq = np.unique(recArray['freq'])

SimPEG/MT/Utils/MT1Danalytic.py

@@ -3,7 +3,7 @@
 import numpy as np, SimPEG as simpeg
 from scipy.constants import mu_0, epsilon_0 as eps_0
 
-def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
+def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
     '''Analytic solution for MT 1D layered earth. Returns E and H fields.
 
     :param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
@@ -12,7 +12,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
     :param numpy array, vector zd: location to calculate EH fields at
     :param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
 
-    Assumes a halfspace with the same conductive as the deepest cell.
+    Assumes a halfspace with the same conductive as the last cell below.
 
     '''
     # Note add an error check for the mesh and sigma are the same size.
@@ -29,7 +29,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
 
     # Initiate the propagation matrix, in the order down up.
     UDp = np.zeros((2,m1d.nC+1),dtype=complex)
-    UDp[1,0] = scaleValue # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
+    UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
     # Loop over all the layers, starting at the bottom layer
     for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
         # Calculate
@@ -38,9 +38,9 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
         # Build the propagation matrix
 
         # Convert fields to down/up going components in layer below current layer
-        Pj1 = np.array([[1,1],[yp1,-yp1]],dtype=complex)
+        Pj1 = np.array([[1,1],[yp1,-yp1]])
         # Convert fields to down/up going components in current layer
-        Pjinv = 1./2*np.array([[1,zp],[1,-zp]],dtype=complex)
+        Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
         # Propagate down and up components through the current layer
         elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
 
@@ -48,14 +48,7 @@ def getEHfields(m1d,sigma,freq,zd,scaleUD=True,scaleValue=1):
         UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
 
         if scaleUD:
-            # Scale the values such that 1 at the top
-            scaleVal = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
-            if np.any(np.isnan(scaleVal)):
-                # If there is a nan (thickness very great), rebuild the move up cell
-                scaleVal = np.zeros_like(UDp[:,lnr+1::-1],dtype=complex)
-                scaleVal[1,0] = scaleValue
-
-            UDp[:,lnr+1::-1] = scaleVal
+             UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
 
     # Calculate the fields
     Ed = np.empty((zd.size,),dtype=complex)

SimPEG/MT/Utils/__init__.py

@@ -0,0 +1,4 @@
+from MT1Dsolutions import *  # Add the names of the functions
+from MT1Danalytic import *
+from dataUtils import *
+from ediFilesUtils import *

SimPEG/MT/Utils/dataUtils.py

@@ -5,25 +5,25 @@ import numpy.lib.recfunctions as recFunc
 from scipy.constants import mu_0
 from scipy import interpolate as sciint
 
-def getAppRes(NSEMdata):
+def getAppRes(MTdata):
     # Make impedance
     zList = []
-    for src in NSEMdata.survey.srcList:
+    for src in MTdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*NSEMdata[src,rx])
+            zc.append(m*MTdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
-def rotateData(NSEMdata,rotAngle):
+def rotateData(MTdata,rotAngle):
     '''
     Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
     '''
-    recData = NSEMdata.toRecArray('Complex')
+    recData = MTdata.toRecArray('Complex')
     impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
     # Make the rotation matrix
     # c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
@@ -40,8 +40,8 @@ def rotateData(NSEMdata,rotAngle):
     for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
         outRec[comp] = rotData[:,nr]
 
-    from SimPEG import NSEM
-    return NSEM.Data.fromRecArray(outRec)
+    from SimPEG import MT
+    return MT.Data.fromRecArray(outRec)
 
 
 def appResPhs(freq,z):
@@ -57,10 +57,10 @@ def rec2ndarr(x,dt=float):
     return x.view((dt, len(x.dtype.names)))
 
 def makeAnalyticSolution(mesh,model,elev,freqs):
-    from SimPEG import NSEM
+    from SimPEG import MT
     data1D = []
     for freq in freqs:
-        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
+        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
 
@@ -71,7 +71,7 @@ def makeAnalyticSolution(mesh,model,elev,freqs):
     return dataRec
 
 def plotMT1DModelData(problem,models,symList=None):
-    from SimPEG import NSEM
+    from SimPEG import MT
     # Setup the figure
     fontSize = 15
 
@@ -79,7 +79,7 @@ def plotMT1DModelData(problem,models,symList=None):
     axM = fig.add_axes([0.075,.1,.25,.875])
     axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
     axM.set_xlim(1e-1,1e5)
-    # axM.set_ylim(-10000,5000)
+    axM.set_ylim(-10000,5000)
     axM.set_ylabel('Depth [km]',fontsize=fontSize)
     axR = fig.add_axes([0.42,.575,.5,.4])
     axR.set_xscale('log')
@@ -132,94 +132,38 @@ def plotMT1DModelData(problem,models,symList=None):
     freq = simpeg.mkvc(data1D['freq'],2)
     res, phs = appResPhs(freq,allData)
 
-    if False:
-        stdCol = 'gray'
-        axRtw = axR.twinx()
-        axRtw.set_ylabel('Std of log10',color=stdCol)
-        [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-        axPtw = axP.twinx()
-        axPtw.set_ylabel('Std ',color=stdCol)
-        [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
-        axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
-        axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
+    stdCol = 'gray'
+    axRtw = axR.twinx()
+    axRtw.set_ylabel('Std of log10',color=stdCol)
+    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    axPtw = axP.twinx()
+    axPtw.set_ylabel('Std ',color=stdCol)
+    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
+    axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
 
     # Fix labels and ticks
 
-    # yMtick = [l/1000 for l in axM.get_yticks().tolist()]
-    # axM.set_yticklabels(yMtick)
+    yMtick = [l/1000 for l in axM.get_yticks().tolist()]
+    axM.set_yticklabels(yMtick)
     [ l.set_rotation(90) for l in axM.get_yticklabels()]
     [ l.set_rotation(90) for l in axR.get_yticklabels()]
-    # [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
-    # [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
+    [(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
+    [t.set_color(stdCol) for t in axPtw.get_yticklabels()]
     for ax in [axM,axR,axP]:
         ax.xaxis.set_tick_params(labelsize=fontSize)
         ax.yaxis.set_tick_params(labelsize=fontSize)
     return fig
 
-def plotImpAppRes(dataArrays,plotLoc,textStr=[]):
-    ''' Plots amplitude impedance and phase'''
-    # fig = plt.figure(1,(7, 7))
-    import plotDataTypes as pDt
-    # axes = ImageGrid(fig, (0.05,0.05,0.875,0.875),nrows_ncols = (2, 2),axes_pad = 0.25,add_all=True,share_all=True,label_mode = "L")
-    # Make the figure and axes
-    fig,axT=plt.subplots(2,2,sharex=True)
-    axes = axT.ravel()
-    fig.set_size_inches((13.5,7.0))
-    fig.suptitle('{:s}\nStation at: {:.1f}x ; {:.1f}y'.format(textStr,plotLoc[0],plotLoc[1]))
-    # Have to deal with axes
-    # Set log
-    for ax in axes.ravel():
-        ax.set_xscale('log')
-
-    axes[0].invert_xaxis()
-    axes[0].set_yscale('log')
-    axes[2].set_yscale('log')
-    # Set labels
-    axes[2].set_xlabel('Frequency [Hz]')
-    axes[3].set_xlabel('Frequency [Hz]')
-    axes[0].set_ylabel('Apperent resistivity [Ohm m]')
-    axes[1].set_ylabel('Apperent phase [degrees]')
-    axes[1].set_ylim(-180,180)
-    axes[2].set_ylabel('Impedance amplitude [V/A]')
-    axes[3].set_ylim(-180,180)
-    axes[3].set_ylabel('Impedance angle [degrees]')
-
-
-    # Plot the data
-    for nr,dataArray in enumerate(dataArrays):
-        if nr==1:
-            parSym = '*'
-        else:
-            parSym = 's'
-        # app res
-        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zxy',par='res',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[0],plotLoc,dataArray,'zyx',par='res',pSym=parSym)
-        # app phs
-        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[1],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
-        # imp abs
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxx',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zxy',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyx',par='abs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[2],plotLoc,dataArray,'zyy',par='abs',pSym=parSym)
-            # imp abs
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxx',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zxy',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyx',par='phs',pSym=parSym)
-        pDt.plotIsoStaImpedance(axes[3],plotLoc,dataArray,'zyy',par='phs',pSym=parSym)
-
-    return fig,axes
-
-
 def printTime():
     import time
     print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
 
-def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
-    from SimPEG import NSEM
+def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
+    from SimPEG import MT
     # Find the unique locations
     # Need to find the locations
-    recDataTemp = NSEMdata.toRecArray()
+    recDataTemp = MTdata.toRecArray()
     # Check if survey.std has been assigned.
     ## NEED TO: write this...
     # Calculte and add the DET of the tensor to the recArray
@@ -241,24 +185,24 @@ def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
         # Make the receiver list
         rx1DList = []
         for rxType in ['z1dr','z1di']:
-            rx1DList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
+            rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
         # Source list
         locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
         dat1DList = []
         src1DList = []
         for freq in locrecData['freq']:
-            src1DList.append(NSEM.SrcNSEM.src_polxy_1Dprimary(rx1DList,freq))
+            src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
             for comp  in ['r','i']:
                 dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
 
         # Make the survey
-        sur1D = NSEM.Survey(src1DList)
+        sur1D = MT.Survey(src1DList)
 
         # Make the data
         dataVec = np.hstack(dat1DList)
-        dat1D = NSEM.Data(sur1D,dataVec)
+        dat1D = MT.Data(sur1D,dataVec)
         sur1D.dobs = dataVec
-        # Need to take NSEMdata.survey.std and split it as well.
+        # Need to take MTdata.survey.std and split it as well.
         std=0.05
         sur1D.std =  np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
         mtData1DList.append(dat1D)
@@ -266,29 +210,29 @@ def convert3Dto1Dobject(NSEMdata,rxType3D='zyx'):
     # Return the the list of data.
     return mtData1DList
 
-def resampleNSEMdataAtFreq(NSEMdata,freqs):
+def resampleMTdataAtFreq(MTdata,freqs):
     """
-    Function to resample NSEMdata at set of frequencies
+    Function to resample MTdata at set of frequencies
 
     """
-    from SimPEG import NSEM
+    from SimPEG import MT
     # Make a rec array
-    NSEMrec = NSEMdata.toRecArray().data
+    MTrec = MTdata.toRecArray().data
 
     # Find unique locations
-    uniLoc = np.unique(NSEMrec[['x','y','z']])
-    uniFreq = NSEMdata.survey.freqs
+    uniLoc = np.unique(MTrec[['x','y','z']])
+    uniFreq = MTdata.survey.freqs
     # Get the comps
-    dNames = NSEMrec.dtype
+    dNames = MTrec.dtype
 
     # Loop over all the locations and interpolate
     for loc in uniLoc:
         # Find the index of the station
-        ind = np.sqrt(np.sum((rec2ndarr(NSEMrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
+        ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
         # Make a temporary recArray and interpolate all the components
         tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
         for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            int1d = sciint.interp1d(NSEMrec[ind]['freq'],NSEMrec[ind][comp],bounds_error=False)
+            int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
             tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
 
         # Join together
@@ -297,5 +241,5 @@ def resampleNSEMdataAtFreq(NSEMdata,freqs):
         except NameError as e:
             outRecArr = tArrRec
 
-    # Make the NSEMdata and return
-    return NSEM.Data.fromRecArray(outRecArr)
+    # Make the MTdata and return
+    return MT.Data.fromRecArray(outRecArr)

SimPEG/MT/Utils/ediFilesUtils.py

@@ -2,7 +2,7 @@
 from SimPEG import mkvc
 from scipy.constants import mu_0
 from numpy.lib import recfunctions as recFunc
-from SimPEG.NSEM.Utils.dataUtils import rec2ndarr
+from SimPEG.MT.Utils.dataUtils import rec2ndarr
 
 # Import modules
 import numpy as np

SimPEG/MT/Utils/sourceUtils.py

@@ -12,7 +12,7 @@ def homo1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.NSEM.Utils import get1DEfields
+    from SimPEG.MT.Utils import get1DEfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh
@@ -77,7 +77,7 @@ def analytic1DModelSource(mesh,freq,sigma_1d):
 
     '''
     # import
-    from SimPEG.NSEM.Utils import getEHfields
+    from SimPEG.MT.Utils import getEHfields
     # Get a 1d solution for a halfspace background
     if mesh.dim == 1:
         mesh1d = mesh

SimPEG/MT/Utils/srcUtils.py

@@ -0,0 +1,46 @@
+import SimPEG as simpeg, numpy as np
+
+def homo1DModelSource(mesh,freq,m_back):
+    '''
+        Function that calculates and return background fields for a 3D mesh and model.
+        The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
+        which is assigned at the fields everywhere for the respective polarizations.2
+
+        :param Simpeg mesh object mesh: Holds information on the discretization
+        :param float freq: The frequency to solve at
+        :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
+    from SimPEG.MT.Utils import get1DEfields
+    # Get a 1d solution for a halfspace background
+    mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
+    # Note: Everything is using e^iwt
+    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 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(ex_px,2),ey_px,ez_px))
+    # Setup y (north) polarization (_py)
+    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 np.arange(mesh.vnEy[0]):
+        for j in np.arange(mesh.vnEy[1]):
+            ey_py[i,j,:] = e0_1d
+    # ey_py[1:-1,1:-1,1:-1] = 0
+    eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
+
+    # Return the electric fields
+    eBG_bp = np.hstack((eBG_px,eBG_py))
+    return eBG_bp

SimPEG/MT/__init__.py

@@ -0,0 +1,5 @@
+import Utils
+from SurveyMT import Rx, Survey, Data
+from FieldsMT import Fields1D_e, Fields3D_e
+import Problem1D, Problem2D, Problem3D
+import SrcMT

SimPEG/Mesh/MeshIO.py

@@ -140,6 +140,7 @@ class TensorMeshIO(object):
                 vtkObj.GetCellData().AddArray(vtkDoubleArr)
             # Set the active scalar
             vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
+        # vtkObj.Update()
 
         # Check the extension of the fileName
         ext = os.path.splitext(fileName)[1]
@@ -156,50 +157,6 @@ class TensorMeshIO(object):
         vtrWriteFilter.SetFileName(fileName)
         vtrWriteFilter.Update()
 
-    def _toVTRObj(mesh,models=None):
-        """
-        Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model.
-
-        Input:
-        :param str, path to the output vtk file
-        :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK
-        :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells
-
-        """
-        # Import
-        from vtk import vtkRectilinearGrid as rectGrid, VTK_VERSION
-        from vtk.util.numpy_support import numpy_to_vtk
-
-        # Deal with dimensionalities
-        if mesh.dim >= 1:
-            vX = mesh.vectorNx
-            xD = mesh.nNx
-            yD,zD = 1,1
-            vY, vZ = np.array([0,0])
-        if mesh.dim >= 2:
-            vY = mesh.vectorNy
-            yD = mesh.nNy
-        if mesh.dim == 3:
-            vZ = mesh.vectorNz
-            zD = mesh.nNz
-        # Use rectilinear VTK grid.
-        # Assign the spatial information.
-        vtkObj = rectGrid()
-        vtkObj.SetDimensions(xD,yD,zD)
-        vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1))
-        vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1))
-        vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1))
-
-        # Assign the model('s) to the object
-        if models is not None:
-            for item in models.iteritems():
-                # Convert numpy array
-                vtkDoubleArr = numpy_to_vtk(item[1],deep=1)
-                vtkDoubleArr.SetName(item[0])
-                vtkObj.GetCellData().AddArray(vtkDoubleArr)
-            # Set the active scalar
-            vtkObj.GetCellData().SetActiveScalars(models.keys()[0])
-        return vtkObj
 
     def readModelUBC(mesh, fileName):
         """

SimPEG/NSEM/ProblemNSEM.py

@@ -1,560 +0,0 @@
-from SimPEG.EM.Utils.EMUtils import omega, mu_0
-from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
-from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem
-from SurveyNSEM import Survey, Data
-from FieldsNSEM import BaseNSEMFields, Fields1D_ePrimSec, Fields3D_ePrimSec
-from SimPEG.NSEM.Utils.MT1Danalytic import getEHfields
-import time, sys
-
-class BaseNSEMProblem(BaseFDEMProblem):
-    """
-        Base class for all Natural source problems.
-    """
-
-    def __init__(self, mesh, **kwargs):
-        BaseFDEMProblem.__init__(self, mesh, **kwargs)
-        Utils.setKwargs(self, **kwargs)
-    # Set the default pairs of the problem
-    surveyPair = Survey
-    dataPair = Data
-    fieldsPair = BaseNSEMFields
-
-    # Set the solver
-    Solver = SimpegSolver
-    solverOpts = {}
-
-    verbose = False
-    # Notes:
-    # Use the forward and devs from BaseFDEMProblem
-    # Might need to add more stuff here.
-
-    ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
-    def Jvec(self, m, v, f=None):
-        """
-        Function to calculate the data sensitivities dD/dm times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nC, 1) - random vector
-            :param NSEMfields object (optional) - NSEM fields object, if not given it is calculated
-            :rtype: NSEMdata object
-            :return: Data sensitivities wrt m
-        """
-
-        # Calculate the fields
-        if f is None:
-           f= self.fields(m)
-        # Set current model
-        self.curModel = m
-        # Initiate the Jv object
-        Jv = self.dataPair(self.survey)
-
-        # Loop all the frequenies
-        for freq in self.survey.freqs:
-            dA_du = self.getA(freq) #
-
-            dA_duI = self.Solver(dA_du, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                # 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.
-                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, 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 )
-                else:
-                    du_dm = dA_duI * ( -dA_dm + dRHS_dm )
-                # Calculate the projection derivatives
-                for rx in src.rxList:
-                    # Get the projection derivative
-                    # v should be of size 2*nE (for 2 polarizations)
-                    PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m
-                    Jv[src, rx] = PDeriv_u(mkvc(du_dm))
-            dA_duI.clean()
-        # Return the vectorized sensitivities
-        return mkvc(Jv)
-
-    def Jtvec(self, m, v, f=None):
-        """
-        Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
-
-            :param numpy.ndarray m (nC, 1) - conductive model
-            :param numpy.ndarray v (nD, 1) - vector
-            :param NSEMfields object f (optional) - NSEM fields object, if not given it is calculated
-            :rtype: NSEMdata object
-            :return: Data sensitivities wrt m
-        """
-
-        if f is None:
-            f = self.fields(m)
-
-        self.curModel = m
-
-        # Ensure v is a data object.
-        if not isinstance(v, self.dataPair):
-            v = self.dataPair(self.survey, v)
-
-        Jtv = np.zeros(m.size)
-
-        for freq in self.survey.freqs:
-            AT = self.getA(freq).T
-
-            ATinv = self.Solver(AT, **self.solverOpts)
-
-            for src in self.survey.getSrcByFreq(freq):
-                ftype = self._solutionType
-                f_src = f[src, :] # Need to fix this...
-
-                for rx in src.rxList:
-                    # Get the adjoint evalDeriv
-                    # PTv needs to be nE,
-                    PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
-                    # Get the
-                    dA_duIT = ATinv * PTv
-                    dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True)
-                    dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
-                    # Make du_dmT
-                    if dRHS_dmT is None:
-                        du_dmT = -dA_dmT
-                    else:
-                        du_dmT = -dA_dmT + dRHS_dmT
-                    # Select the correct component
-                    # du_dmT needs to be of size nC,
-                    real_or_imag = rx.projComp
-                    if real_or_imag == 'real':
-                        Jtv +=  du_dmT.real
-                    elif real_or_imag == 'imag':
-                        Jtv +=  -du_dmT.real
-                    else:
-                        raise Exception('Must be real or imag')
-            # 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._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

SimPEG/NSEM/Utils/__init__.py

@@ -1,5 +0,0 @@
-from MT1Dsolutions import get1DEfields  # Add the names of the functions
-from MT1Danalytic import getEHfields, getImpedance
-from dataUtils import *
-from ediFilesUtils import *
-from testUtils import *

SimPEG/NSEM/Utils/testUtils.py

@@ -1,198 +0,0 @@
-import unittest
-import sys
-from scipy.constants import mu_0
-import SimPEG as simpeg
-
-from SimPEG.Utils import meshTensor
-import numpy as np
-
-np.random.seed(1100)
-# Define the tolerances
-TOLr = 5e-2
-TOLp = 5e-2
-
-
-def getAppResPhs(NSEMdata):
-    # Make impedance
-    from SimPEG.NSEM.Utils import appResPhs
-    zList = []
-    for src in NSEMdata.survey.srcList:
-        zc = [src.freq]
-        for rx in src.rxList:
-            if 'i' in rx.rxType:
-                m=1j
-            else:
-                m = 1
-            zc.append(m*NSEMdata[src,rx])
-        zList.append(zc)
-    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
-
-
-def setup1DSurvey(sigmaHalf,tD=True,structure=False):
-    from SimPEG import NSEM
-    # Frequency
-    nFreq = 33
-    freqs = np.logspace(3,-3,nFreq)
-    # Make the mesh
-    ct = 5
-    air = meshTensor([(ct,25,1.3)])
-    # coreT0 = meshTensor([(ct,15,1.2)])
-    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
-    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
-    bot = meshTensor([(core[0],20,-1.3)])
-    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
-    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
-    # Make the model
-    sigma = np.zeros(m1d.nC) + sigmaHalf
-    sigma[m1d.gridCC > 0 ] = 1e-8
-    sigmaBack = sigma.copy()
-    # Add structure
-    if structure:
-        shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
-        deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
-        sigma[shallow] = 1
-        sigma[deep] = 0.1
-
-    rxList = []
-    for rxType in ['z1dr','z1di']:
-        rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
-    # Source list
-    srcList =[]
-    if tD:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1DhomotD(rxList,freq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-
-    survey = NSEM.Survey(srcList)
-    return survey, sigma, m1d
-
-
-def setupSimpegNSEM_ePrimSec(inputSetup,comp='Imp',singleFreq=False,expMap=True):
-    from SimPEG import NSEM
-
-    M,freqs,sig,sigBG,rx_loc = inputSetup
-    # Make a receiver list
-    rxList = []
-    if comp == 'All':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Imp':
-        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    elif comp == 'Tip':
-        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
-            rxList.append(NSEM.Rx(rx_loc,rxType))
-    else:
-        rxList.append(NSEM.Rx(rx_loc,comp))
-    # Source list
-    srcList =[]
-
-    if singleFreq:
-        srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,singleFreq))
-    else:
-        for freq in freqs:
-            srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
-    # Survey NSEM
-    survey = NSEM.Survey(srcList)
-
-    ## Setup the problem object
-    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
-    if expMap:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= np.log(sigma1d) )
-        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-        problem.curModel = np.log(sig)
-    else:
-        problem = NSEM.Problem3D_ePrimSec(M,sigmaPrimary= sigma1d)
-        problem.curModel = sig
-    problem.pair(survey)
-    problem.verbose = False
-    try:
-        from pymatsolver import MumpsSolver
-        problem.Solver = MumpsSolver
-    except:
-        pass
-
-    return (survey, problem)
-
-def getInputs():
-    """
-    Function that returns Mesh, freqs, rx_loc, elev.
-    """
-    # 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])
-    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
-    M = simpeg.Mesh.TensorMesh([[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,6,-1.5),(200.,4),(200,6,1.5)],[(200,8,-1.5),(200.,8),(200,8,1.5)]], x0=['C','C','C'])# Setup the model
-    # Set the frequencies
-    freqs = np.logspace(1,-3,5)
-    elev = 0
-
-    ## Setup the the survey object
-    # Receiver locations
-    rx_x, rx_y = np.meshgrid(np.arange(-350,350,200),np.arange(-350,350,200))
-    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
-
-    return M, freqs, rx_loc, elev
-
-def random(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Backround
-    sigBG = np.ones(M.nC)*conds
-    # Add randomness to the model (10% of the value).
-    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def halfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def blockInhalfSpace(conds):
-    ''' Returns a halfspace model based on the inputs'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    # conds = [1e-2]
-    groundInd = ccM[:,2] < elev
-    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
-    sig[~groundInd] = 1e-8
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-    return (M, freqs, sig, sigBG, rx_loc)
-
-def twoLayer(conds):
-    ''' Returns a 2 layer model based on the conductivity values given'''
-    M, freqs, rx_loc, elev = getInputs()
-
-    # Model
-    ccM = M.gridCC
-    groundInd = ccM[:,2] < elev
-    botInd = ccM[:,2] < -3000
-    sig = np.zeros(M.nC) + 1e-8
-    sig[groundInd] = conds[1]
-    sig[botInd] = conds[0]
-    # Set the background, not the same as the model
-    sigBG = np.zeros(M.nC) + 1e-8
-    sigBG[groundInd] = conds[1]
-
-
-    return (M, freqs, sig, sigBG, rx_loc)
-

SimPEG/NSEM/__init__.py

@@ -1,5 +0,0 @@
-import Utils
-from SurveyNSEM import Rx, Survey, Data
-from FieldsNSEM import Fields1D_ePrimSec, Fields3D_ePrimSec
-from ProblemNSEM import Problem1D_ePrimSec, Problem3D_ePrimSec
-import SrcNSEM

SimPEG/Survey.py

@@ -1,4 +1,5 @@
 import Utils, numpy as np, scipy.sparse as sp, uuid
+import gc
 
 class BaseRx(object):
     """SimPEG Receiver Object"""

docs/em/api_TDEM.rst

@@ -347,10 +347,10 @@ and
 
 
 
-TDEM - B formulation
-====================
+TDEM Problem
+============
 
-.. automodule:: SimPEG.EM.TDEM.TDEM_b
+.. automodule:: SimPEG.EM.TDEM.TDEM
     :show-inheritance:
     :members:
     :undoc-members:
@@ -359,7 +359,7 @@ TDEM - B formulation
 Field Storage
 =============
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.FieldsTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Fields
     :show-inheritance:
     :members:
     :undoc-members:
@@ -369,19 +369,19 @@ Field Storage
 TDEM Survey Classes
 ===================
 
-.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.SurveyTDEM
+.. autoclass:: SimPEG.EM.TDEM.SurveyTDEM.Survey
     :show-inheritance:
     :members:
     :undoc-members:
     :inherited-members:
 
 
-Base Classes
-============
+.. Base Classes
+.. ============
 
-.. automodule:: SimPEG.EM.TDEM.BaseTDEM
-    :show-inheritance:
-    :members:
-    :undoc-members:
-    :inherited-members:
+.. .. automodule:: SimPEG.EM.TDEM.BaseTDEM
+..     :show-inheritance:
+..     :members:
+..     :undoc-members:
+..     :inherited-members:
 

docs/examples/MT_1D_analytic_nlayer_Earth.rst

@@ -1,21 +0,0 @@
-.. _examples_MT_1D_analytic_nlayer_Earth:
-
-.. --------------------------------- ..
-..                                   ..
-..    THIS FILE IS AUTO GENEREATED   ..
-..                                   ..
-..    SimPEG/Examples/__init__.py    ..
-..                                   ..
-.. --------------------------------- ..
-
-MT 1D analytic nlayer Earth
-===========================
-
-.. plot::
-
-    from SimPEG import Examples
-    Examples.MT_1D_analytic_nlayer_Earth.run()
-
-.. literalinclude:: ../../SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
-    :language: python
-    :linenos:

tests/em/tdem/test_TDEM_b_DerivAdjoint.py

@@ -3,310 +3,218 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
-tol = 1e-6
 
-class TDEM_bDerivTests(unittest.TestCase):
+testDeriv   = True
+testAdjoint = True
 
-    def setUp(self):
+TOL = 1e-5
 
-        cs = 5.
-        ncx = 20
-        ncy = 6
-        npad = 20
-        hx = [(cs,ncx), (cs,npad,1.3)]
-        hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-        mesh = Mesh.CylMesh([hx,1,hy], '00C')
+def setUp(prbtype='b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
-        active = mesh.vectorCCz<0.
-        activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-        mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
 
-        rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    survey = EM.TDEM.Survey([src])
 
-        survey = EM.TDEM.SurveyTDEM([src])
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
 
-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
 
-        try:
-            from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
-        except ImportError, e:
-            self.prb.Solver = SolverLU
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
 
-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
 
-        self.prb.pair(survey)
-        self.mesh = mesh
+    prb.pair(survey)
+    mesh = mesh
 
-    def test_AhVec(self):
-        """
-            Test that fields and AhVec produce consistent results
-        """
+    return prb, m, mesh
 
-        prb = self.prb
-        sigma = self.sigma
 
-        u = prb.fields(sigma)
-        Ahu = prb._AhVec(sigma, u)
+class TDEM_DerivTests(unittest.TestCase):
 
-        V1 = Ahu[:,'b',1]
-        V2 = 1./prb.timeSteps[0]*prb.MfMui*u[:,'b',0]
-        self.assertLess(np.linalg.norm(V1-V2)/np.linalg.norm(V2), 1.e-6)
+# ====== TEST A ========== #
 
-        V1 = Ahu[:,'e',1]
-        return np.linalg.norm(V1) < 1.e-6
+    def AderivTest(self, prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2
+        if prbtype == 'b':
+            nu = mesh.nF
+        elif prbtype == 'e':
+            nu = mesh.nE
+        v = np.random.rand(nu)
 
-        for i in range(2,prb.nT):
+        def AderivFun(m):
+            prb.curModel = m
+            A = prb.getAdiag(tInd)
+            Av = A*v
+            prb.curModel = m0
+            ADeriv_dm = lambda dm: prb.getAdiagDeriv(tInd, v, dm)
 
-            dt = prb.timeSteps[i]
+            return Av, ADeriv_dm
 
-            V1 = Ahu[:,'b',i]
-            V2 = 1.0/dt*prb.MfMui*u[:,'b', i-1]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            self.assertLess(np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6)
+        print '\n Testing ADeriv %s'%(prbtype)
+        Tests.checkDerivative(AderivFun, m0, plotIt=False, num=4, eps=1e-20)
 
-            V1 = Ahu[:,'e',i]
-            V2 = prb.MeSigma*u[:,'e',i]
-            # print np.linalg.norm(V1), np.linalg.norm(V2)
-            return np.linalg.norm(V1)/np.linalg.norm(V2), 1.e-6
-
-    def test_AhVecVSMat_OneTS(self):
-
-        prb = self.prb
-        prb.timeSteps = [1e-05]
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        u1 = A*f.tovec()
-        u2 = prb._AhVec(sigma,f).tovec()
-
-        self.assertTrue(np.linalg.norm(u1-u2)/np.linalg.norm(u1)<1e-12)
-
-    def test_solveAhVSMat_OneTS(self):
-        prb = self.prb
-
-        prb.timeSteps = [1e-05]
-
-        sigma = self.sigma
-        prb.curModel = sigma
-
-        dt = prb.timeSteps[0]
-        a11 = 1.0/dt*prb.MfMui*sp.identity(prb.mesh.nF)
-        a12 = prb.MfMui*prb.mesh.edgeCurl
-        a21 = prb.mesh.edgeCurl.T*prb.MfMui
-        a22 = -prb.MeSigma
-        A = sp.bmat([[a11,a12],[a21,a22]])
-
-        f = prb.fields(sigma)
-        f[:,:,0] = {'b':0}
-        f[:,'b',1] = 0
-
-        self.assertTrue(np.all(np.r_[f[:,'b',1],f[:,'e',1]] == f.tovec()))
-
-        u1 = prb.solveAh(sigma,f).tovec().flatten()
-        u2 = sp.linalg.spsolve(A.tocsr(),f.tovec())
-
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_solveAhVsAhVec(self):
-
-        prb = self.prb
-        mesh = self.prb.mesh
-        sigma = self.sigma
-        self.prb.curModel = sigma
-
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        f[:,'b',:] = 0.0
-        for i in range(prb.nT):
-            f[:,'e', i] = np.random.rand(mesh.nE, 1)
-
-        Ahf = prb._AhVec(sigma, f)
-        f_test = prb.solveAh(sigma, Ahf)
-
-        u1 = f.tovec()
-        u2 = f_test.tovec()
-        self.assertTrue(np.linalg.norm(u1-u2)<1e-8)
-
-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
-
-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-        return passed
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
-
-
-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
-
-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        mesh = self.mesh
-
-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f[:,'e',i] = np.random.rand(mesh.nE, 1)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
-
-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
-
-        self.assertTrue((V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointAhVsAht(self):
-        prb = self.prb
-        mesh = self.mesh
-        sigma = self.sigma
-
-        f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-            f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = f2.tovec().dot(prb._AhVec(sigma, f1).tovec())
-        V2 = f1.tovec().dot(prb._AhtVec(sigma, f2).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    # def test_solveAhtVsAhtVec(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = np.random.rand(prb.mapping.nP)
-
-    #     f1 = EM.TDEM.FieldsTDEM(mesh,prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = prb.solveAht(sigma, f1)
-    #     f3 = prb._AhtVec(sigma, f2)
-
-    #     if True:
-    #         import matplotlib.pyplot as plt
-    #         plt.plot(f3.tovec(),'b')
-    #         plt.plot(f1.tovec(),'r')
-    #         plt.show()
-    #     V1 = np.linalg.norm(f3.tovec()-f1.tovec())
-    #     V2 = np.linalg.norm(f1.tovec())
-    #     print 'AhtVsAhtVec', V1, V2, f1.tovec()
-    #     print 'I am gunna fail this one: boo. :('
-    #     self.assertLess(V1/V2, 1e-6)
-
-    # def test_adjointsolveAhVssolveAht(self):
-    #     prb = self.prb
-    #     mesh = self.mesh
-    #     sigma = self.sigma
-
-    #     f1 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f1[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f1[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     f2 = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-    #     for i in range(1,prb.nT+1):
-    #         f2[:,'b',i] = np.random.rand(mesh.nF, 1)
-    #         f2[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-    #     V1 = f2.tovec().dot(prb.solveAh(sigma, f1).tovec())
-    #     V2 = f1.tovec().dot(prb.solveAht(sigma, f2).tovec())
-    #     print V1, V2
-    #     self.assertLess(np.abs(V1-V2)/np.abs(V1), 1e-6)
-
-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+    def A_adjointTest(self,prbtype):
+        prb, m0, mesh = setUp(prbtype)
+        tInd = 2
 
+        print '\n Testing A_adjoint'
         m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        if prbtype == 'b':
+            nu = prb.mesh.nF
+        elif prbtype == 'e':
+            nu = prb.mesh.nE
 
-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 1)
-            u[:,'e',i] = np.random.rand(mesh.nE, 1)
+        v = np.random.rand(nu)
+        u = np.random.rand(nu)
+        prb.curModel = m0
 
-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 1)
-            v[:,'e',i] = np.random.rand(mesh.nE, 1)
-
-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = v.tovec().dot(prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < tol)
-
-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
-
-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
-
-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        passed = np.abs(V1-V2)/np.abs(V1) < tol
-        print 'AdjointTest', V1, V2, passed
+        tInd = 2 # not actually used
+        V1 = v.dot(prb.getAdiagDeriv(tInd, u, m))
+        V2 = m.dot(prb.getAdiagDeriv(tInd, u, v, adjoint=True))
+        passed = np.abs(V1-V2) < TOL * (np.abs(V1) + np.abs(V2))/2.
+        print 'AdjointTest %s'%(prbtype), V1, V2, passed
         self.assertTrue(passed)
 
+    def test_Aderiv_b(self):
+        self.AderivTest('b')
+    def test_Aderiv_e(self):
+        self.AderivTest('e')
+
+    def test_Aadjoint_b(self):
+        self.A_adjointTest('b')
+    def test_Aadjoint_e(self):
+        self.A_adjointTest('e')
+
+# ====== TEST Fields Deriv Pieces ========== #
+
+    def test_eDeriv_m_adjoint(self):
+        prb, m0, mesh = setUp()
+        tInd = 0
+
+        v = np.random.rand(mesh.nF)
+
+        print '\n Testing eDeriv_m Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        m = np.random.rand(prb.mapping.nP)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_m(1, prb.survey.srcList[0], m))
+        V2 = m.dot(f._eDeriv_m(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+    def test_eDeriv_u_adjoint(self):
+        print '\n Testing eDeriv_u Adjoint'
+
+        prb, m0, mesh = setUp()
+        f = prb.fields(m0)
+
+        b = np.random.rand(prb.mesh.nF)
+        e = np.random.randn(prb.mesh.nE)
+        V1 = e.dot(f._eDeriv_u(1, prb.survey.srcList[0], b))
+        V2 = b.dot(f._eDeriv_u(1, prb.survey.srcList[0], e, adjoint=True))
+        tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+        passed = np.abs(V1-V2) < tol
+
+        print '    ', V1, V2, np.abs(V1-V2), tol, passed
+        self.assertTrue(passed)
+
+
+# ====== TEST Jvec ========== #
+
+    if testDeriv:
+
+        def JvecTest(self, prbtype, rxcomp):
+            prb, m, mesh = setUp(prbtype, rxcomp)
+
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m, mx)]
+            print '\n'
+            print 'test_Jvec_%s_%s' %(prbtype, rxcomp)
+            Tests.checkDerivative(derChk, m, plotIt=False, num=2, eps=1e-20)
+
+        def test_Jvec_b_bx(self):
+            self.JvecTest('b','bx')
+
+        def test_Jvec_b_bz(self):
+            self.JvecTest('b','bz')
+
+        def test_Jvec_b_dbxdt(self):
+            self.JvecTest('b','dbxdt')
+
+        def test_Jvec_b_dbzdt(self):
+            self.JvecTest('b','dbzdt')
+
+        def test_Jvec_b_ey(self):
+            self.JvecTest('b','ey')
+
+        def test_Jvec_e_ey(self):
+            self.JvecTest('e','ey')
+
+
+# ====== TEST Jtvec ========== #
+
+    if testAdjoint:
+
+        def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'):
+
+            print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
+
+            prb, m0, mesh = setUp(prbtype, rxcomp)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.randn(prb.survey.nD)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)
+
+        def test_Jvec_adjoint_b_bx(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_bz(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_dbxdt(self):
+            self.JvecVsJtvecTest('b', 'bx')
+
+        def test_Jvec_adjoint_b_dbzdt(self):
+            self.JvecVsJtvecTest('b', 'bz')
+
+        def test_Jvec_adjoint_b_ey(self):
+            self.JvecVsJtvecTest('b', 'ey')
+
+        # This is not working because Problem_e has not done
+        # def test_Jvec_adjoint_e_ey(self):
+        #     self.JvecVsJtvecTest('e', 'ey')
 
 
 

tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py

@@ -3,10 +3,12 @@ from SimPEG import *
 from SimPEG import EM
 
 plotIt = False
+testDeriv = True
+testAdjoint = True
 
-class TDEM_bDerivTests(unittest.TestCase):
+TOL = 1e-5
 
-    def setUp(self):
+def setUp(self, rxcomp='bz'):
 
         cs = 5.
         ncx = 20
@@ -21,131 +23,78 @@ class TDEM_bDerivTests(unittest.TestCase):
         mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
 
         rxOffset = 40.
-        rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
-        src = EM.TDEM.SrcTDEM_VMD_MVP( [rx], loc=np.array([0., 0., 0.]))
-        rx2 = EM.TDEM.RxTDEM(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), 'bz')
-        src2 = EM.TDEM.SrcTDEM_VMD_MVP( [rx2], loc=np.array([0., 0., 0.]))
+        rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), rxcomp)
+        src = EM.TDEM.Src.MagDipole( [rx], loc=np.array([0., 0., 0.]))
+        rx2 = EM.TDEM.Rx(np.array([[rxOffset-10, 0., 0.]]), np.logspace(-5,-4, 25), rxcomp)
+        src2 = EM.TDEM.Src.MagDipole( [rx2], loc=np.array([0., 0., 0.]))
 
-        survey = EM.TDEM.SurveyTDEM([src,src2])
+        survey = EM.TDEM.Survey([src,src2])
 
-        self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-        # self.prb.timeSteps = [1e-5]
-        self.prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-        # self.prb.timeSteps = [(1e-05, 100)]
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+        # prb.timeSteps = [1e-5]
+        prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+        # prb.timeSteps = [(1e-05, 100)]
 
         try:
             from pymatsolver import MumpsSolver
-            self.prb.Solver = MumpsSolver
+            prb.Solver = MumpsSolver
         except ImportError, e:
-            self.prb.Solver  = SolverLU
+            prb.Solver  = SolverLU
 
-        self.sigma = np.ones(mesh.nCz)*1e-8
-        self.sigma[mesh.vectorCCz<0] = 1e-1
-        self.sigma = np.log(self.sigma[active])
+        m = np.log(1e-1)*np.ones(prb.mapping.nP) + 1e-2*np.random.randn(prb.mapping.nP)
 
-        self.prb.pair(survey)
-        self.mesh = mesh
+        prb.pair(survey)
 
-    def test_DerivG(self):
-        """
-            Test the derivative of c with respect to sigma
-        """
+        return mesh, prb, m
 
-        # Random model and perturbation
-        sigma = np.random.rand(self.prb.mapping.nP)
-
-        f = self.prb.fields(sigma)
-        dm = 1000*np.random.rand(self.prb.mapping.nP)
-        h = 0.01
-
-        derChk = lambda m: [self.prb._AhVec(m, f).tovec(), lambda mx: self.prb.Gvec(sigma, mx, u=f).tovec()]
-        print '\ntest_DerivG'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_dUdM(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        dm = 10*np.random.rand(prb.mapping.nP)
-        f = prb.fields(sigma)
-
-        derChk = lambda m: [self.prb.fields(m).tovec(), lambda mx: -prb.solveAh(sigma, prb.Gvec(sigma, mx, u=f)).tovec()]
-        print '\n'
-        print 'test_Deriv_dUdM'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=dm, num=4, eps=1e-20)
-
-    def test_Deriv_J(self):
-
-        prb = self.prb
-        prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-        mesh = self.mesh
-        sigma = self.sigma
-
-        # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-        d_sig = 10*np.random.rand(prb.mapping.nP)
+class TDEM_bDerivTests(unittest.TestCase):
 
 
-        derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-        print '\n'
-        print 'test_Deriv_J'
-        Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=4, eps=1e-20)
+    if testDeriv:
+        def Deriv_J(self, rxcomp='bz'):
 
-    def test_projectAdjoint(self):
-        prb = self.prb
-        survey = prb.survey
-        nSrc = survey.nSrc
-        mesh = self.mesh
+            mesh, prb, m0 = setUp(rxcomp)
 
-        # Generate random fields and data
-        f = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(prb.nT):
-            f[:,'b',i] = np.random.rand(mesh.nF, nSrc)
-            f[:,'e',i] = np.random.rand(mesh.nE, nSrc)
-        d_vec = np.random.rand(survey.nD)
-        d = Survey.Data(survey,v=d_vec)
+            prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
 
-        # Check that d.T*Q*f = f.T*Q.T*d
-        V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
-        V2 = np.sum((f.tovec())*(survey.evalDeriv(None, v=d, adjoint=True).tovec()))
+            derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(m0, mx)]
+            print '\n'
+            print 'test_Deriv_J %s'%rxcomp
+            Tests.checkDerivative(derChk, m0, plotIt=False, num=3, eps=1e-20)
 
-        self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
+        def test_Jvec_bx(self):
+            self.Deriv_J('bx')
 
-    def test_adjointGvecVsGtvec(self):
-        mesh = self.mesh
-        prb = self.prb
+        def test_Jvec_bz(self):
+            self.Deriv_J('bz')
 
-        m = np.random.rand(prb.mapping.nP)
-        sigma = np.random.rand(prb.mapping.nP)
+        def test_Jvec_ey(self):
+            self.Deriv_J('ey')
 
-        u = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            u[:,'b',i] = np.random.rand(mesh.nF, 2)
-            u[:,'e',i] = np.random.rand(mesh.nE, 2)
+    if testAdjoint:
+        def adjointJvecVsJtvec(self, rxcomp='bz'):
+            print ' \n Testing Adjoint %s' %rxcomp
+            mesh, prb, m0 = setUp(rxcomp)
 
-        v = EM.TDEM.FieldsTDEM(prb.mesh, prb.survey)
-        for i in range(1,prb.nT+1):
-            v[:,'b',i] = np.random.rand(mesh.nF, 2)
-            v[:,'e',i] = np.random.rand(mesh.nE, 2)
+            m = np.random.rand(prb.mapping.nP)
+            d = np.random.rand(prb.survey.nD)
 
-        V1 = m.dot(prb.Gtvec(sigma, v, u))
-        V2 = np.sum(v.tovec()*prb.Gvec(sigma, m, u).tovec())
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) <1e-6)
+            V1 = d.dot(prb.Jvec(m0, m))
+            V2 = m.dot(prb.Jtvec(m0, d))
 
-    def test_adjointJvecVsJtvec(self):
-        mesh = self.mesh
-        prb = self.prb
-        sigma = self.sigma
+            tol = TOL * (np.abs(V1) + np.abs(V2)) / 2.
+            passed = np.abs(V1-V2) < tol
+            print '    ', V1, V2, np.abs(V1-V2), tol, passed
+            self.assertTrue(passed)
 
-        m = np.random.rand(prb.mapping.nP)
-        d = np.random.rand(prb.survey.nD)
+        def test_JvecVsJtvec_bx(self):
+            self.adjointJvecVsJtvec('bx')
 
-        V1 = d.dot(prb.Jvec(sigma, m))
-        V2 = m.dot(prb.Jtvec(sigma, d))
-        print 'AdjointTest', V1, V2
-        self.assertTrue(np.abs(V1-V2)/np.abs(V1) < 1e-6)
+        def test_JvecVsJtvec_bz(self):
+            self.adjointJvecVsJtvec('bz')
+
+        def test_JvecVsJtvec_ey(self):
+            self.adjointJvecVsJtvec('ey')
 
 
 

tests/em/tdem/test_TDEM_combos.py

@@ -1,94 +0,0 @@
-import unittest
-from SimPEG import *
-from SimPEG import EM
-
-plotIt = False
-
-def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
-    cs = 5.
-    ncx = 20
-    ncy = 6
-    npad = 20
-    hx = [(cs,ncx), (cs,npad,1.3)]
-    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
-    mesh = Mesh.CylMesh([hx,1,hy], '00C')
-
-    active = mesh.vectorCCz<0.
-    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
-    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
-
-    rxOffset = 40.
-
-    srcs = []
-    for ii in range(nSrc):
-        rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
-        srcs += [EM.TDEM.SrcTDEM_VMD_MVP(rxs,np.array([0., 0., 0.]))]
-
-    survey = EM.TDEM.SurveyTDEM(srcs)
-
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
-    # prb.timeSteps = [1e-5]
-    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
-    # prb.timeSteps = [(1e-05, 100)]
-
-    try:
-        from pymatsolver import MumpsSolver
-        prb.Solver = MumpsSolver
-    except ImportError, e:
-        prb.Solver  = SolverLU
-
-    sigma = np.ones(mesh.nCz)*1e-8
-    sigma[mesh.vectorCCz<0] = 1e-1
-    sigma = np.log(sigma[active])
-
-    prb.pair(survey)
-    return prb, mesh, sigma
-
-def dotestJvec(prb, mesh, sigma):
-    prb.timeSteps = [(1e-05, 10), (0.0001, 10), (0.001, 10)]
-    # d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
-    d_sig = 10*np.random.rand(prb.mapping.nP)
-    derChk = lambda m: [prb.survey.dpred(m), lambda mx: prb.Jvec(sigma, mx)]
-    return Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
-
-def dotestAdjoint(prb, mesh, sigma):
-    m = np.random.rand(prb.mapping.nP)
-    d = np.random.rand(prb.survey.nD)
-
-    V1 = d.dot(prb.Jvec(sigma, m))
-    V2 = m.dot(prb.Jtvec(sigma, d))
-    print 'AdjointTest', V1, V2
-    return np.abs(V1-V2)/np.abs(V1), 1e-6
-
-class TDEM_bDerivTests(unittest.TestCase):
-
-    def test_Jvec_bx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx')))
-    def test_Adjoint_bx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx')))
-
-    def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
-    def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
-
-    def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
-    def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
-
-    def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
-    def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
-
-    def test_Jvec_dbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt')))
-    def test_Adjoint_dbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt')))
-
-    def test_Jvec_dbzdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbzdt')))
-    def test_Adjoint_dbzdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbzdt')))
-
-    def test_Jvec_dbxdtbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='dbxdt,bz')))
-    def test_Adjoint_dbxdtbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='dbxdt,bz')))
-
-    def test_Jvec_ey(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey')))
-    def test_Adjoint_ey(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey')))
-
-    def test_Jvec_eybzdbxdt(self): self.assertTrue(dotestJvec(*getProb(rxTypes='ey,bz,dbxdt')))
-    def test_Adjoint_eybzdbxdt(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='ey,bz,dbxdt')))
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/em/tdem/test_TDEM_crosscheck.py

@@ -0,0 +1,76 @@
+import unittest
+from SimPEG import *
+from SimPEG import EM
+
+TOL = 1e-5
+FLR = 1e-20
+
+np.random.seed(seed=25) # set a seed so that the same conductivity model is used for all runs
+
+def setUp(prbtype = 'b', rxcomp='bz'):
+    cs = 5.
+    ncx = 20
+    ncy = 15
+    npad = 20
+    hx = [(cs,ncx), (cs,npad,1.3)]
+    hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
+    mesh = Mesh.CylMesh([hx,1,hy], '00C')
+#
+    active = mesh.vectorCCz<0.
+    activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
+    mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
+
+    rxOffset = 10.
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,]
+    src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.]))
+
+    survey = EM.TDEM.Survey([src])
+
+    if prbtype == 'b':
+        prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    elif prbtype == 'e':
+        prb = EM.TDEM.Problem_e(mesh, mapping=mapping)
+
+    prb.timeSteps = [(1e-05, 10), (5e-05, 10), (2.5e-4, 10)]
+    # prb.timeSteps = [(1e-05, 10), (1e-05, 50), (1e-05, 50) ] #, (2.5e-4, 10)]
+
+    try:
+        from pymatsolver import MumpsSolver
+        prb.Solver = MumpsSolver
+    except ImportError, e:
+        prb.Solver = SolverLU
+
+    m = np.log(1e-1)*np.ones(prb.mapping.nP) #+ 1e-2*np.random.randn(prb.mapping.nP)
+
+    prb.pair(survey)
+    mesh = mesh
+
+    return prb, m, mesh
+
+def CrossCheck(prbtype1='b', prbtype2='e', rxcomp='bz'):
+
+    prb1,m1,mesh1 = setUp(prbtype1, rxcomp)
+    prb2,m2,mesh2 = setUp(prbtype2, rxcomp)
+
+    assert (m1 == m2).all(), 'Models for two formulations are different'
+
+    d1 = prb1.survey.dpred(m1)
+    d2 = prb2.survey.dpred(m2)
+
+
+    check = np.linalg.norm(d1 - d2)
+    tol   = 0.5 * (np.linalg.norm(d1) + np.linalg.norm(d2)) * TOL
+    passed = check < tol
+
+    print 'Checking %s, %s for %s data'%(prbtype1, prbtype2, rxcomp)
+    print '  ', np.linalg.norm(d1), np.linalg.norm(d2), np.linalg.norm(check), tol, passed
+
+    assert passed
+
+class TDEM_cross_check_EB(unittest.TestCase):
+    def test_EB_ey(self):
+        CrossCheck('b','e','ey')
+
+if __name__ == '__main__':
+    unittest.main()
+

tests/em/tdem/test_TDEM_forward_Analytic.py

@@ -29,12 +29,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
     mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
 
-    rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
-    src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
+    rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
+    src = EM.TDEM.Src.MagDipole([rx], waveform= EM.TDEM.Src.StepOffWaveform(), loc=np.array([0., 0., 0.]))
     # src = EM.TDEM.SrcTDEM([rx], loc=np.array([0., 0., 0.]))
 
-    survey = EM.TDEM.SurveyTDEM([src])
-    prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
+    survey = EM.TDEM.Survey([src])
+    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
     prb.Solver = MumpsSolver
 
     prb.timeSteps = [(1e-06, 40), (5e-06, 40), (1e-05, 40), (5e-05, 40), (0.0001, 40), (0.0005, 40)]
@@ -50,6 +50,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
 
     ind = np.logical_and(rx.times > bounds[0],rx.times < bounds[1])
     log10diff = np.linalg.norm(np.log10(np.abs(bz_calc[ind])) - np.log10(np.abs(bz_ana[ind])))/np.linalg.norm(np.log10(np.abs(bz_ana[ind])))
+
+    print ' |bz_ana| = ',np.linalg.norm(bz_ana), ' |bz_num| = ', np.linalg.norm(bz_calc), ' |bz_ana - bz_num| =', np.linalg.norm(bz_ana-bz_calc)
     print 'Difference: ', log10diff
 
     if showIt == True:
@@ -61,6 +63,12 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=None,
     return log10diff
 
 
+class TDEM_SimpleSrcTests(unittest.TestCase):
+    def test_source(self):
+        waveform = EM.TDEM.Src.StepOffWaveform()
+        assert waveform.eval(0.) == 0.
+
+
 class TDEM_bTests(unittest.TestCase):
 
     def test_analytic_p2_CYL_50m(self):

tests/mt/forward/test_Problem1D_VsAnalyticHalfspace.py

@@ -1,102 +0,0 @@
-import unittest
-import SimPEG as simpeg
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-import numpy as np
-# Define the tolerances
-TOLr = 5e-1
-TOLp = 5e-1
-
-
-
-def appRes_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
-    problem = NSEM.Problem1D_eTotal(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app res and phs
-    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
-
-    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
-
-def appPhs_TotalFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
-    problem = NSEM.Problem1D_eTotal(mesh)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-def appRes_psFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app res and phs
-    app_r = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,0]
-
-    return np.linalg.norm(np.abs(np.log(app_r) - np.log(np.ones(survey.nFreq)/sigmaHalf))*np.log(sigmaHalf))
-
-def appPhs_psFieldNorm(sigmaHalf):
-
-    # Make the survey
-    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False)
-    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
-
-    # Get the fields
-    fields = problem.fields(sigma)
-
-    # Project the data
-    data = survey.eval(fields)
-
-    # Calculate the app  phs
-    app_p = np.array(NSEM.Utils.testUtils.getAppResPhs(data))[:,1]
-
-    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        pass
-    # Total Fields
-    # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
-    # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
-
-    # Primary/secondary
-    def test_appRes1en0_ps(self):self.assertLess(appRes_psFieldNorm(1e-0), TOLr)
-    def test_appPhs1en0_ps(self):self.assertLess(appPhs_psFieldNorm(1e-0), TOLp)
-    def test_appRes2en1_ps(self):self.assertLess(appRes_psFieldNorm(2e-1), TOLr)
-    def test_appPhs2en1_ps(self):self.assertLess(appPhs_psFieldNorm(2e-1), TOLp)
-    def test_appRes2en3_ps(self):self.assertLess(appRes_psFieldNorm(2e-3), TOLr)
-    def test_appPhs2en3_ps(self):self.assertLess(appPhs_psFieldNorm(2e-3), TOLp)
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/mt/forward/test_Problem3D_VsAnalyticSolution.py

@@ -1,54 +0,0 @@
-# Test functions
-from glob import glob
-import numpy as np, sys, os, time, scipy, subprocess
-import SimPEG as simpeg
-import unittest
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-from scipy.constants import mu_0
-
-np.random.seed(1100)
-
-TOLr = 1
-TOLp = 2
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = [1e-1, 2e-1]
-addrandoms = True
-
-def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
-    if appR:
-        label = 'resistivity'
-    else:
-        label = 'phase'
-    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
-
-    # Calculate the app  phs
-    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(NSEM.Utils.testUtils.halfSpace(sigmaHalf),expMap=expMap)
-    data = problem.dataPair(survey,survey.dpred(problem.curModel))
-    recData = data.toRecArray('Complex')
-    app_rpxy, app_rpyx = NSEM.Utils.appResPhs(recData['freq'],recData['zxy'])[0], NSEM.Utils.appResPhs(recData['freq'],recData['zyx'])[0]
-    if appR:
-        return np.linalg.norm( np.abs(np.log10(app_rpxy[0]) - np.log10(1./sigmaHalf)) * np.log10(sigmaHalf ))
-    else:
-        return np.linalg.norm( np.abs(app_rpxy[1] + 135) / 135 )
-
-
-class TestAnalytics(unittest.TestCase):
-
-    def setUp(self):
-        # Make the survey and the problem
-        pass
-
-    # # Test apparent resistivity and phase
-    def test_appRes1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2),TOLr)
-    def test_appPhs1en2(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False),TOLp)
-
-    def test_appRes1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1),TOLr)
-    def test_appPhs1en1(self):self.assertLess(appResPhsHalfspace_eFrom_ps_Norm(1e-1,False),TOLp)
-
-
-
-if __name__ == '__main__':
-    unittest.main()

tests/mt/inversion/__init__.py

@@ -1,12 +0,0 @@
-import os
-import glob
-import unittest
-
-if __name__ == '__main__':
-    test_file_strings = glob.glob('test_*.py')
-    module_strings = [str[0:len(str)-3] for str in test_file_strings]
-    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
-              in module_strings]
-    testSuite = unittest.TestSuite(suites)
-
-    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/mt/inversion/test_Problem3D_Adjoint.py

@@ -1,58 +0,0 @@
-# Test functions
-from glob import glob
-import numpy as np, sys, os, time, scipy, subprocess
-import SimPEG as simpeg
-import unittest
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-from scipy.constants import mu_0
-
-
-TOLr = 5e-2
-TOL = 1e-4
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = [1e-1, 2e-1]
-addrandoms = True
-
-
-
-def JvecAdjointTest(inputSetup,comp='All',freq=False):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp='All',singleFreq=freq)
-    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
-
-    m  = sig
-    u = problem.fields(m)
-
-    v = np.random.rand(survey.nD,)
-    # print problem.PropMap.PropModel.nP
-    w = np.random.rand(problem.mesh.nC,)
-
-    vJw = v.ravel().dot(problem.Jvec(m, w, u))
-    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
-    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
-    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
-    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
-    return np.abs(vJw - wJtv) < tol
-
-
-class NSEM_AdjointTests(unittest.TestCase):
-
-    def setUp(self):
-        pass
-
-    # Test the adjoint of Jvec and Jtvec
-    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
-    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
-    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
-    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
-    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
-    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
-    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
-    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
-    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
-
-if __name__ == '__main__':
-    unittest.main()

tests/mt/inversion/test_Problem3D_Derivs.py

@@ -1,83 +0,0 @@
-# Test functions
-from glob import glob
-import numpy as np, sys, os, time, scipy, subprocess
-import SimPEG as simpeg
-import unittest
-from SimPEG import NSEM
-from SimPEG.Utils import meshTensor
-from scipy.constants import mu_0
-
-np.random.seed(1100)
-
-TOLr = 5e-2
-TOL = 1e-4
-FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
-CONDUCTIVITY = 1e1
-MU = mu_0
-freq = [1e-1, 2e-1]
-addrandoms = True
-
-
-# Test the Jvec derivative
-def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
-    (M, freqs, sig, sigBG, rx_loc) = inputSetup
-    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
-    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
-    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-    # problem.sigmaPrimary = np.log(sigBG)
-    x0 = np.log(sigBG)
-    # cond = sig[0]
-    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
-    # problem.sigmaPrimary = x0
-    # if True:
-    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
-    survey = problem.survey
-    def fun(x):
-        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
-    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
-
-def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
-
-    survey, problem = NSEM.Utils.testUtils.setupSimpegNSEM_ePrimSec(inputSetup,comp,freq)
-    print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
-    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
-    # Initate things for the derivs Test
-    src = survey.srcList[0]
-    rx = src.rxList[0]
-
-    u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
-    u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
-    u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
-    f0 = problem.fieldsPair(survey.mesh,survey)
-    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
-    f0[src,'e_pxSolution'] =  u0[:len(u0)/2]#u0x
-    f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
-
-    def fun(u):
-        f = problem.fieldsPair(survey.mesh,survey)
-        f[src,'e_pxSolution'] = u[:len(u)/2]
-        f[src,'e_pySolution'] = u[len(u)/2::]
-        return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
-
-    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
-
-
-
-class NSEM_DerivTests(unittest.TestCase):
-
-    def setUp(self):
-        pass
-
-    # Do a derivative test of Jvec
-    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
-    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
-    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
-    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
-    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
-    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
-    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
-    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
-    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(NSEM.Utils.testUtils.random(1e-2),'All',.1))
-
-if __name__ == '__main__':
-    unittest.main()

tests/mt/test_ApparentResistivityAnalytic.py

@@ -1,10 +1,13 @@
 import unittest
 from SimPEG import *
-from SimPEG import NSEM
+from SimPEG import MT
 
 TOL = 1e-6
 
-
+def appResPhs(freq,z):
+    app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
+    app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
+    return app_res, app_phs
 
 def appResNorm(sigmaHalf):
     nFreq = 26
@@ -17,12 +20,12 @@ def appResNorm(sigmaHalf):
     freqs = np.logspace(4,-4,nFreq)
     Z = []
     for freq in freqs:
-        Ed, Eu, Hd, Hu = NSEM.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
+        Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
         Z.append((Ed + Eu)/(Hd + Hu))
 
     Zarr = np.concatenate(Z)
 
-    app_r, app_p = NSEM.Utils.appResPhs(freqs,Zarr)
+    app_r, app_p = appResPhs(freqs,Zarr)
 
     return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
 

tests/mt/test_Problem1D_againstAnalyticHalfspace.py

@@ -0,0 +1,162 @@
+import unittest
+import SimPEG as simpeg
+from SimPEG import MT
+from SimPEG.Utils import meshTensor
+import numpy as np
+# Define the tolerances
+TOLr = 5e-2
+TOLp = 5e-2
+
+
+def setupSurvey(sigmaHalf,tD=True):
+
+    # Frequency
+    nFreq = 33
+    freqs = np.logspace(3,-3,nFreq)
+    # Make the mesh
+    ct = 5
+    air = meshTensor([(ct,25,1.3)])
+    # coreT0 = meshTensor([(ct,15,1.2)])
+    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
+    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
+    bot = meshTensor([(core[0],10,-1.3)])
+    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
+    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
+    # Make the model
+    sigma = np.zeros(m1d.nC) + sigmaHalf
+    sigma[m1d.gridCC > 0 ] = 1e-8
+
+    rxList = []
+    for rxType in ['z1dr','z1di']:
+        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+    # Source list
+    srcList =[]
+    if tD:
+        for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
+    else:
+        for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+
+    survey = MT.Survey(srcList)
+    return survey, sigma, m1d
+
+def getAppResPhs(MTdata):
+    # Make impedance
+    def appResPhs(freq,z):
+        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
+        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
+        return app_res, app_phs
+    zList = []
+    for src in MTdata.survey.srcList:
+        zc = [src.freq]
+        for rx in src.rxList:
+            if 'i' in rx.rxType:
+                m=1j
+            else:
+                m = 1
+            zc.append(m*MTdata[src,rx])
+        zList.append(zc)
+    return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
+
+def appRes_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = setupSurvey(sigmaHalf)
+    problem = MT.Problem1D.eForm_TotalField(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
+
+def appPhs_TotalFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = setupSurvey(sigmaHalf)
+    problem = MT.Problem1D.eForm_TotalField(mesh)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+def appRes_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = setupSurvey(sigmaHalf,False)
+    problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app res and phs
+    app_r = np.array(getAppResPhs(data))[:,0]
+
+    return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
+
+def appPhs_psFieldNorm(sigmaHalf):
+
+    # Make the survey
+    survey, sigma, mesh = setupSurvey(sigmaHalf,False)
+    problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
+    problem.pair(survey)
+
+    # Get the fields
+    fields = problem.fields(sigma)
+
+    # Project the data
+    data = survey.eval(fields)
+
+    # Calculate the app  phs
+    app_p = np.array(getAppResPhs(data))[:,1]
+
+    return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        pass
+    # Total Fields
+    # def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
+    # def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
+
+    # def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr)
+    # def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp)
+
+    # def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
+    # def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
+
+    # def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
+    # def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
+
+    # def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
+    # def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
+
+    # def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
+    # def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
+
+    # Primary/secondary
+    def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
+    def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
+
+if __name__ == '__main__':
+    unittest.main()

tests/mt/test_Problem1D_totalDvsPSvsAnalytic.py

@@ -1,6 +1,6 @@
 import unittest
 import SimPEG as simpeg
-from SimPEG import NSEM
+from SimPEG import MT
 from SimPEG.Utils import meshTensor
 import numpy as np
 # Define the tolerances
@@ -8,30 +8,69 @@ TOLr = 5e-2
 TOLp = 5e-2
 
 
-def getAppResPhs(NSEMdata):
+def setupSurvey(sigmaHalf,tD=True):
+
+    # Frequency
+    nFreq = 33
+    freqs = np.logspace(3,-3,nFreq)
+    # Make the mesh
+    ct = 5
+    air = meshTensor([(ct,25,1.3)])
+    # coreT0 = meshTensor([(ct,15,1.2)])
+    # coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
+    core = np.concatenate( (  np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
+    bot = meshTensor([(core[0],15,-1.3)])
+    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
+    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
+    # Make the model
+    sigma = np.zeros(m1d.nC) + sigmaHalf
+    sigma[m1d.gridCC > 0 ] = 1e-8
+    sigmaBack = sigma.copy()
+    # Add structure
+    shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
+    deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
+    sigma[shallow] = 1
+    sigma[deep] = 0.1
+
+    rxList = []
+    for rxType in ['z1dr','z1di']:
+        rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
+    # Source list
+    srcList =[]
+    if tD:
+        for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
+    else:
+        for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+
+    survey = MT.Survey(srcList)
+    return survey, sigma, m1d
+
+def getAppResPhs(MTdata):
     # Make impedance
     def appResPhs(freq,z):
         app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
         app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
         return app_res, app_phs
     zList = []
-    for src in NSEMdata.survey.srcList:
+    for src in MTdata.survey.srcList:
         zc = [src.freq]
         for rx in src.rxList:
             if 'i' in rx.rxType:
                 m=1j
             else:
                 m = 1
-            zc.append(m*NSEMdata[src,rx])
+            zc.append(m*MTdata[src,rx])
         zList.append(zc)
     return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
 
 def calculateAnalyticSolution(srcList,mesh,model):
-    surveyAna = NSEM.Survey(srcList)
-    data1D = NSEM.Data(surveyAna)
+    surveyAna = MT.Survey(srcList)
+    data1D = MT.Data(surveyAna)
     for src in surveyAna.srcList:
         elev = src.rxList[0].locs[0]
-        anaEd, anaEu, anaHd, anaHu = NSEM.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
+        anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
         anaE = anaEd+anaEu
         anaH = anaHd+anaHu
         # Scale the solution
@@ -47,12 +86,12 @@ def dataMis_AnalyticTotalDomain(sigmaHalf):
     # Make the survey
 
     # Total domain solution
-    surveyTD, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf)
-    problemTD = NSEM.Problem1D_eTotal(mesh) # This not fully implemented
+    surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
+    problemTD = MT.Problem1D.eForm_TotalField(mesh)
     problemTD.pair(surveyTD)
     # Analytic data
     dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
-    # dataTDObj = NSEM.DataNSEM.DataNSEM(surveyTD, surveyTD.dpred(sigma))
+    # dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
     dataTD = surveyTD.dpred(sigma)
     dataAna = simpeg.mkvc(dataAnaObj)
     return np.all((dataTD - dataAna)/dataAna < 2.)
@@ -69,16 +108,16 @@ def dataMis_AnalyticPrimarySecondary(sigmaHalf):
 
     # Make the survey
     # Primary secondary
-    survey, sigma, mesh = NSEM.Utils.testUtils.setup1DSurvey(sigmaHalf,False,structure=True)
+    surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
+    problemPS = MT.Problem1D.eForm_psField(mesh)
+    problemPS.sigmaPrimary = sigmaPS
+    problemPS.pair(surveyPS)
     # Analytic data
-    problem = NSEM.Problem1D_ePrimSec(mesh, sigmaPrimary = sigma)
-    problem.pair(survey)
+    dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
 
-    dataAnaObj = calculateAnalyticSolution(survey.srcList,mesh,sigma)
-
-    data = survey.dpred(sigma)
+    dataPS = surveyPS.dpred(sigmaPS)
     dataAna = simpeg.mkvc(dataAnaObj)
-    return np.all((data - dataAna)/dataAna < 2.)
+    return np.all((dataPS - dataAna)/dataAna < 2.)
 
 
 

tests/mt/test_Problem3D_againstAnalytic.py

@@ -0,0 +1,268 @@
+# Test functions
+from glob import glob
+import numpy as np, sys, os, time, scipy, subprocess
+import SimPEG as simpeg
+import unittest
+from SimPEG import MT
+from SimPEG.Utils import meshTensor
+from scipy.constants import mu_0
+
+TOLr = 5e-2
+TOL = 1e-4
+FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
+CONDUCTIVITY = 1e1
+MU = mu_0
+freq = [1e-1, 2e-1]
+addrandoms = True
+
+
+def getInputs():
+    """
+    Function that returns Mesh, freqs, rx_loc, elev.
+    """
+    # 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])
+    # M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
+    M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model
+    # Set the frequencies
+    freqs = np.logspace(1,-3,5)
+    elev = 0
+
+    ## Setup the the survey object
+    # Receiver locations
+    rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500))
+    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
+
+    return M, freqs, rx_loc, elev
+
+def random(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Backround
+    sigBG = np.ones(M.nC)*conds
+    # Add randomness to the model (10% of the value).
+    sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def halfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def blockInhalfSpace(conds):
+    ''' Returns a halfspace model based on the inputs'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    # conds = [1e-2]
+    groundInd = ccM[:,2] < elev
+    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
+    sig[~groundInd] = 1e-8
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+def twoLayer(conds):
+    ''' Returns a 2 layer model based on the conductivity values given'''
+    M, freqs, rx_loc, elev = getInputs()
+
+    # Model
+    ccM = M.gridCC
+    groundInd = ccM[:,2] < elev
+    botInd = ccM[:,2] < -3000
+    sig = np.zeros(M.nC) + 1e-8
+    sig[groundInd] = conds[1]
+    sig[botInd] = conds[0]
+    # Set the background, not the same as the model
+    sigBG = np.zeros(M.nC) + 1e-8
+    sigBG[groundInd] = conds[1]
+
+
+    return (M, freqs, sig, sigBG, rx_loc)
+
+
+
+def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
+    M,freqs,sig,sigBG,rx_loc = inputSetup
+    # Make a receiver list
+    rxList = []
+    if comp == 'All':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(MT.Rx(rx_loc,rxType))
+    elif comp == 'Imp':
+        for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
+            rxList.append(MT.Rx(rx_loc,rxType))
+    elif comp == 'Tip':
+        for rxType in ['tzxr','tzxi','tzyr','tzyi']:
+            rxList.append(MT.Rx(rx_loc,rxType))
+    else:
+        rxList.append(MT.Rx(rx_loc,comp))
+    # Source list
+    srcList =[]
+
+    if singleFreq:
+        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
+    else:
+        for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+    # Survey MT
+    survey = MT.Survey(srcList)
+
+    ## Setup the problem object
+    sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
+    if expMap:
+        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
+        problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+        problem.curModel = np.log(sig)
+    else:
+        problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
+        problem.curModel = sig
+    problem.pair(survey)
+    problem.verbose = False
+    try:
+        from pymatsolver import MumpsSolver
+        problem.Solver = MumpsSolver
+    except:
+        pass
+
+    return (survey, problem)
+
+def getAppResPhs(MTdata):
+    # Make impedance
+    def appResPhs(freq,z):
+        app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
+        app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
+        return app_res, app_phs
+    recData = MTdata.toRecArray('Complex')
+    return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
+
+def JvecAdjointTest(inputSetup,comp='All',freq=False):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
+    print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
+
+    m  = sig
+    u = problem.fields(m)
+
+    v = np.random.rand(survey.nD,)
+    # print problem.PropMap.PropModel.nP
+    w = np.random.rand(problem.mesh.nC,)
+
+    vJw = v.ravel().dot(problem.Jvec(m, w, u))
+    wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
+    tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
+    print ' vJw   wJtv  vJw - wJtv     tol    abs(vJw - wJtv) < tol'
+    print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
+    return np.abs(vJw - wJtv) < tol
+
+# Test the Jvec derivative
+def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
+    (M, freqs, sig, sigBG, rx_loc) = inputSetup
+    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
+    print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
+    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+    # problem.sigmaPrimary = np.log(sigBG)
+    x0 = np.log(sigBG)
+    # cond = sig[0]
+    # x0 = np.log(np.ones(problem.mesh.nC)*cond)
+    # problem.sigmaPrimary = x0
+    # if True:
+    #     x0  = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
+    survey = problem.survey
+    def fun(x):
+        return survey.dpred(x), lambda x: problem.Jvec(x0, x)
+    return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
+
+def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
+
+    survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
+    print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
+    # problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
+    # Initate things for the derivs Test
+    src = survey.srcList[0]
+    rx = src.rxList[0]
+
+    u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
+    u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
+    u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
+    f0 = problem.fieldsPair(survey.mesh,survey)
+    # u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
+    f0[src,'e_pxSolution'] =  u0[:len(u0)/2]#u0x
+    f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
+
+    def fun(u):
+        f = problem.fieldsPair(survey.mesh,survey)
+        f[src,'e_pxSolution'] = u[:len(u)/2]
+        f[src,'e_pySolution'] = u[len(u)/2::]
+        return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
+
+    return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)
+
+def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
+    if appR:
+        label = 'resistivity'
+    else:
+        label = 'phase'
+    # Make the survey and the problem
+    survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
+    print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
+
+    data = problem.dataPair(survey,survey.dpred(problem.curModel))
+    # Calculate the app  phs
+    app_rpxy, app_rpyx = np.array(getAppResPhs(data))
+    if appR:
+        return np.all(np.abs(app_rpxy[0,:] - 1./sigmaHalf) * sigmaHalf < .4)
+    else:
+        return np.all(np.abs(app_rpxy[1,:] + 135) / 135 < .4)
+
+class TestAnalytics(unittest.TestCase):
+
+    def setUp(self):
+        pass
+    # # Test apparent resistivity and phase
+    def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
+    def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
+
+    def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
+    def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
+
+    # Do a derivative test of Jvec
+    # def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
+    # def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
+    # def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
+    # def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
+    # def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
+    # def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
+    # def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
+    # def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
+    def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1))
+
+    # Test the adjoint of Jvec and Jtvec
+    # def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
+    # def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
+    # def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
+    # def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
+    # def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
+    # def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
+    # def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
+    # def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
+    def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1))
+
+if __name__ == '__main__':
+    unittest.main()