Added AhVec method and test.

simpegEM/TDEM/BaseTDEM.py

@@ -132,12 +132,15 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
     @property
     def MfMui(self): return self._MfMui
 
+    @property
+    def MeSigma(self): return self._MeSigma
+
     @property
     def MeSigmaI(self): return self._MeSigmaI
 
     def makeMassMatrices(self, m):
-        MeSigma = self.mesh.getMass(m, loc='e')
-        self._MeSigmaI = sdiag(1/MeSigma.diagonal())
+        self._MeSigma = self.mesh.getMass(m, loc='e')
+        self._MeSigmaI = sdiag(1/self.MeSigma.diagonal())
         self._MfMui = self.mesh.getMass(1/mu_0, loc='f')
 
 
@@ -208,6 +211,12 @@ class FieldsTDEM(object):
         self.set_b(newFields['b'], tInd)
         self.set_e(newFields['e'], tInd)
 
+    def fieldVec(self):
+        u = np.ndarray((0,self.nTx))
+        for i in range(self.nTimes):
+            u = np.r_[u, self.get_b(i), self.get_e(i)]
+        return u
+
     ####################################################
     # Get Methods
     ####################################################

simpegEM/TDEM/TDEM_b.py

@@ -64,6 +64,26 @@ class ProblemTDEM_b(ProblemBaseTDEM):
         Y = self.fields(m, useThisRhs=AhRHS, useThisCalcFields=AhCalcFields)
         return Y
 
+    ####################################################
+    # Functions for tests
+    ####################################################
+
+    def AhVec(self, m, u):
+        self.makeMassMatrices(m)
+        dt = self.getDt(0)
+        b = 1/dt*u.get_b(0) + self.mesh.edgeCurl*u.get_e(0)
+        e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(0) - self.MeSigma*u.get_e(0)
+        f = FieldsTDEM(self.mesh, 1, self.times.size, 'b')
+        f.set_b(b, 0)
+        f.set_e(e, 0)
+        for i in range(1,self.times.size):
+            dt = self.getDt(i)
+            b = 1/dt*u.get_b(i) + self.mesh.edgeCurl*u.get_e(i) - 1/dt*u.get_b(i-1)
+            e = self.mesh.edgeCurl.T*self.MfMui*u.get_b(i) - self.MeSigma*u.get_e(i)
+            f.set_b(b, i)
+            f.set_e(e, i)
+        return f
+
 if __name__ == '__main__':
     from SimPEG import *
     import simpegEM as EM

simpegEM/Tests/test_forward_EMproblem.py

@@ -5,7 +5,7 @@ from scipy.constants import mu_0
 from simpegEM.Utils.Ana import hzAnalyticDipoleT
 
 
-class EMProblemTests(unittest.TestCase):
+class TDEM_bTests(unittest.TestCase):
 
     def setUp(self):
 
@@ -32,17 +32,27 @@ class EMProblemTests(unittest.TestCase):
        self.sigma[mesh.vectorCCz<0] = 0.1
        self.prb.pair(self.dat)
 
-
     def test_analitic_b(self):
        bz_calc = self.dat.dpred(self.sigma)
        bz_ana = mu_0*hzAnalyticDipoleT(self.dat.rxLoc[0], self.prb.times, self.sigma[0])
 
        diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten())
-       print diff
        self.assertTrue(diff<0.05)
 
-    def test_awesome(self):
-        self.assertTrue(True)
+    def test_AhVec(self):
+        """
+            Test that fields and AhVec produce consistent results
+        """
+
+        sigma = np.ones(self.prb.mesh.nCz)*1e-8
+        sigma[self.prb.mesh.vectorCCz<0] = 0.1
+        u = self.prb.fields(sigma)
+        Ahu = self.prb.AhVec(sigma, u)
+        self.assertTrue(np.linalg.norm(Ahu.get_b(0)-1/self.prb.getDt(0)*u.get_b(-1))/np.linalg.norm(u.get_b(0)) < 1.e-2)
+        self.assertTrue(np.linalg.norm(Ahu.get_e(0))/np.linalg.norm(u.get_e(0)) < 1.e-2)
+        for i in range(1,u.nTimes):
+            self.assertTrue(np.linalg.norm(Ahu.get_b(i))/np.linalg.norm(u.get_b(i)) < 1.e-2)
+            self.assertTrue(np.linalg.norm(Ahu.get_e(i))/np.linalg.norm(u.get_e(i)) < 1.e-2)