minor fix

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

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/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/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:
 

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):