testing the sensitivities

simpegEM/TDEM/BaseTDEM.py

@@ -135,6 +135,8 @@ class ProblemBaseTDEM(MixinTimeStuff, MixinInitialFieldCalc, BaseProblem):
         self.makeMassMatrices(m)
 
         F = self.getInitialFields()
+        #TODO: Split next code to forward and adjoint.
+        # fields would call forward
         dtFact = None
         for tInd, t in enumerate(self.times):
             dt = self.getDt(tInd)

simpegEM/TDEM/TDEM_b.py

@@ -131,23 +131,6 @@ if __name__ == '__main__':
     prb.pair(dat)
     sigma = np.random.rand(mesh.nCz)
 
-    f = FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
-    for i in range(f.nTimes):
-        f.set_b(np.zeros((mesh.nF, 1)), i)
-        f.set_e(np.random.rand(mesh.nE, 1), i)
-
-    Ahf = prb.AhVec(sigma, f)
-    f_test = prb.solveAh(sigma, Ahf)
-
-    print np.linalg.norm(f.fieldVec() - f_test.fieldVec())
-
-    e0 = f.get_e(0)
-    e1 = f_test.get_e(0)
-    b0 = f.get_b(0)
-    b1 = f_test.get_b(0)
-    plt.semilogy(np.abs(e0))
-    plt.semilogy(np.abs(e1),'r')
-    plt.show()
 
 
 

simpegEM/Tests/test_forward_EMproblem.py

@@ -3,6 +3,7 @@ from SimPEG import *
 import simpegEM as EM
 from scipy.constants import mu_0
 from simpegEM.Utils.Ana import hzAnalyticDipoleT
+import matplotlib.pyplot as plt
 
 
 class TDEM_bTests(unittest.TestCase):
@@ -66,6 +67,8 @@ class TDEM_bDerivTests(unittest.TestCase):
        self.sigma = np.ones(mesh.nCz)*1e-8
        self.sigma[mesh.vectorCCz<0] = 0.1
        self.prb.pair(self.dat)
+       self.mesh = mesh
+
 
     def test_AhVec(self):
         """
@@ -83,7 +86,7 @@ class TDEM_bDerivTests(unittest.TestCase):
             self.assertTrue(np.linalg.norm(Ahu.get_e(i))/np.linalg.norm(u.get_e(i)) < 1.e-2)
 
     def test_AhVecVSMat_OneTS(self):
-        
+
         self.prb.setTimes([1e-5], [1])
 
         sigma = np.ones(self.prb.mesh.nCz)*1e-8
@@ -165,6 +168,54 @@ class TDEM_bDerivTests(unittest.TestCase):
         # Assuming that the gradient is exact to machine precision
         self.assertTrue(b<1e-16)
 
+    def test_Deriv_dUdM(self):
+
+        prb = self.prb
+        prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
+        mesh = self.mesh
+        sigma = self.sigma
+
+        d_sig = sigma.copy() #np.random.rand(mesh.nCz)
+        d_sig[d_sig==1e-8] = 0
+
+        num = 10
+        error = np.zeros(num)
+        order = 0
+        hv = np.logspace(-1.2,-3, num)
+        for i, h in enumerate(hv):
+            f = prb.fields(sigma)
+            fstep = prb.fields(sigma + h*d_sig)
+            dcdm = prb.G(sigma, h*d_sig, u=f) # TODO: make negative!?!?
+            dudm = prb.solveAh(sigma, dcdm)
+
+            linear = np.linalg.norm(f.fieldVec() - fstep.fieldVec())
+            quad = np.linalg.norm(f.fieldVec() - fstep.fieldVec() - dudm.fieldVec())
+            error[i] = quad
+            if i > 0:
+                order = np.log(error[i]/error[i-1])/np.log(hv[i]/hv[i-1])
+
+            # print np.log(linearB/quadB)/np.log(h)
+            print h, linear, quad, order
+
+        self.assertTrue(order > 1.8)
+
+    def test_Deriv_J(self):
+
+        prb = self.prb
+        prb.setTimes([1e-5, 1e-4, 1e-3], [10, 10, 10])
+        mesh = self.mesh
+        sigma = self.sigma
+
+        d_sig = 0.8*sigma #np.random.rand(mesh.nCz)
+        d_sig[d_sig==1e-8] = 0
+
+
+        derChk = lambda m: [prb.data.dpred(m), lambda mx: -prb.J(sigma, mx)]
+        passed = Tests.checkDerivative(derChk, sigma, plotIt=False, dx=d_sig, num=2, eps=1e-20)
+        self.assertTrue(passed)
+
+
+
 
 if __name__ == '__main__':
     unittest.main()