start of Problem_e (e Jvec working)

SimPEG/EM/TDEM/FieldsTDEM.py

@@ -31,7 +31,15 @@ class Fields(SimPEG.Problem.TimeFields):
     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 = {
@@ -65,10 +73,10 @@ class Fields_b(Fields):
     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 _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 ) )
@@ -93,7 +101,42 @@ class Fields_b(Fields):
 
         return self.MeSigmaIDeriv(-S_e + self.edgeCurl.T * ( self.MfMui * bSolution)) * v - self.MeSigmaI * S_eDeriv(v)
 
-    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)
+
+
+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/SurveyTDEM.py

@@ -144,13 +144,18 @@ class BaseSrc(SimPEG.Survey.BaseSrc):
         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)
@@ -218,13 +223,41 @@ class MagDipole(BaseSrc):
 
 
     def bInitial(self, prob):
-        eqLocs = prob._eqLocs
 
         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)
+
+        if adjoint:
+            raise NotImplementedError
+
+        return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v
+
+
     def S_m(self, prob, time):
         if self.waveform.hasInitialFields is False:
             raise NotImplementedError

SimPEG/EM/TDEM/TDEM.py

@@ -116,9 +116,13 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
 
                 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 = 0
 
-                JRHS = dRHS_dm_v - dA_dm_v  # - dAsubdiag_dm_v (which is zero)
+                # if tInd >= 1:
+                dAsubdiag_dm_v = self.getAsubdiagDeriv(tInd, u[src,ftype,tInd], v)
+                # else:
+                #     dAsubdiag_dm_v = Zero()
+
+                JRHS = dRHS_dm_v - dAsubdiag_dm_v - dA_dm_v
 
                 # step in time and overwrite
                 if tInd != len(self.timeSteps+1):
@@ -180,8 +184,6 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
                     df_duT_v[src, '%sDeriv'%self._fieldType, tInd] = df_duT_v[src, '%sDeriv'%self._fieldType, tInd] + Utils.mkvc(cur[0],2)
                     JTv = JTv + cur[1]
 
-
-
         del PT_v # no longer need this
 
 
@@ -309,6 +311,8 @@ class BaseTDEMProblem(Problem.BaseTimeProblem, BaseEMProblem):
 ################################ E-B Formulation #########################################
 ##########################################################################################
 
+# ------------------------------- Problem_b -------------------------------------------- #
+
 class Problem_b(BaseTDEMProblem):
     """
     Starting from the quasi-static E-B formulation of Maxwell's equations (semi-discretized)
@@ -404,6 +408,9 @@ class Problem_b(BaseTDEMProblem):
 
         return Asubdiag
 
+    def getAsubdiagDeriv(self, tInd, u, v, adjoint=False):
+        return Zero() * v
+
 
 
     def getRHS(self, tInd):
@@ -450,6 +457,68 @@ class Problem_b(BaseTDEMProblem):
         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):
+
+        dt = self.timeSteps[tInd]
+        C = self.mesh.edgeCurl
+        MfMui = self.MfMui
+        MeSigmaDeriv = self.MeSigmaDeriv(u)
+
+        if adjoint:
+            raise 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 TDEM import BaseTDEMProblem, Problem_b
+from TDEM import BaseTDEMProblem, Problem_b, Problem_e
 from FieldsTDEM import Fields, Fields_b
-from SurveyTDEM import Survey, BaseSrc, Rx
+from SurveyTDEM import Survey, BaseSrc, Rx

tests/em/tdem/test_TDEM_b_DerivAdjoint.py

@@ -5,11 +5,11 @@ from SimPEG import EM
 plotIt = False
 
 testDeriv   = True
-testAdjoint = True
+testAdjoint = False
 
 TOL = 1e-5
 
-def setUp(rxcomp='bz'):
+def setUp(prbtype='b', rxcomp='bz'):
     cs = 5.
     ncx = 20
     ncy = 15
@@ -28,7 +28,11 @@ def setUp(rxcomp='bz'):
 
     survey = EM.TDEM.Survey([src])
 
-    prb = EM.TDEM.Problem_b(mesh, mapping=mapping)
+    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)]
 
@@ -129,40 +133,42 @@ class TDEM_DerivTests(unittest.TestCase):
 
     if testDeriv:
 
-        def JvecTest(self, rxcomp):
-            prb, m, mesh = setUp(rxcomp)
+        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' %(rxcomp)
+            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('bx')
+            self.JvecTest('b','bx')
 
         def test_Jvec_b_bz(self):
-            self.JvecTest('bz')
+            self.JvecTest('b','bz')
 
         def test_Jvec_b_dbxdt(self):
-            self.JvecTest('dbxdt')
+            self.JvecTest('b','dbxdt')
 
         def test_Jvec_b_dbzdt(self):
-            self.JvecTest('dbzdt')
+            self.JvecTest('b','dbzdt')
 
         def test_Jvec_b_ey(self):
-            self.JvecTest('ey')
+            self.JvecTest('b','ey')
 
+        def test_Jvec_e_ey(self):
+            self.JvecTest('e','ey')
 
 
 # ====== TEST Jtvec ========== #
 
     if testAdjoint:
 
-        def JvecVsJtvecTest(self, rxcomp='bz'):
+        def JvecVsJtvecTest(self, prbtype='b', rxcomp='bz'):
 
             print '\nAdjoint Testing Jvec, Jtvec %s' %(rxcomp)
 
-            prb, m0, mesh = setUp(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))
@@ -174,19 +180,19 @@ class TDEM_DerivTests(unittest.TestCase):
             self.assertTrue(passed)
 
         def test_Jvec_adjoint_b_bx(self):
-            self.JvecVsJtvecTest('bx')
+            self.JvecVsJtvecTest('b', 'bx')
 
         def test_Jvec_adjoint_b_bz(self):
-            self.JvecVsJtvecTest('bz')
+            self.JvecVsJtvecTest('b', 'bz')
 
         def test_Jvec_adjoint_b_dbxdt(self):
-            self.JvecVsJtvecTest('bx')
+            self.JvecVsJtvecTest('b', 'bx')
 
         def test_Jvec_adjoint_b_dbzdt(self):
-            self.JvecVsJtvecTest('bz')
+            self.JvecVsJtvecTest('b', 'bz')
 
         def test_Jvec_adjoint_b_ey(self):
-            self.JvecVsJtvecTest('ey')
+            self.JvecVsJtvecTest('b', 'ey')
 
 
 

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.SurveyTDEM.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()
+