From ffb87aafa7017ad17fd749a445e08f34fa437801 Mon Sep 17 00:00:00 2001
From: Dave Marchant <dwfmarchant@gmail.com>
Date: Wed, 12 Feb 2014 11:36:49 -0800
Subject: [PATCH] Added tests for G.

---
 simpegEM/TDEM/TDEM_b.py                  | 24 +++++++-------
 simpegEM/Tests/test_forward_EMproblem.py | 42 ++++++++++++++++++++++++
 2 files changed, 55 insertions(+), 11 deletions(-)

diff --git a/simpegEM/TDEM/TDEM_b.py b/simpegEM/TDEM/TDEM_b.py
index db15c5a5..3ddefb1d 100644
--- a/simpegEM/TDEM/TDEM_b.py
+++ b/simpegEM/TDEM/TDEM_b.py
@@ -70,7 +70,9 @@ class ProblemTDEM_b(ProblemBaseTDEM):
     # Functions for tests
     ####################################################
 
-    def AhVec(self, m, u):
+    def AhVec(self, m, u=None):
+        if u is None:
+            u = self.fields(m)
         self.makeMassMatrices(m)
         dt = self.getDt(0)
         b = 1/dt*u.get_b(0) + self.mesh.edgeCurl*u.get_e(0)
@@ -124,26 +126,26 @@ if __name__ == '__main__':
     # Ahu = prb.AhVec(sigma, u)
 
     # Random fields
-    f = FieldsTDEM(prb.mesh, 1, prb.times.size, 'b')
-    for i in range(f.nTimes):
-        f.set_b(np.random.rand(mesh.nF, 1), i)
-        f.set_e(np.random.rand(mesh.nE, 1), i)
-
-
     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.random.rand(mesh.nF, 1), i)
+        # f.set_e(np.random.rand(mesh.nE, 1), i)
+    f = prb.fields(sigma)
+
     dm = np.random.rand(mesh.nCz)
 
-    for h in np.logspace(30, 10, 20):
+    for h in np.logspace(0, -10, 10):
         # print h
         a = np.linalg.norm(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec())
         b = np.linalg.norm(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec() - h*prb.G(sigma, dm, u=f).fieldVec())
         print a, b, b/a
     # print 
     # h = 1.
-    # plt.semilogy(-prb.AhVec(sigma+h*dm, f).fieldVec() + prb.AhVec(sigma, f).fieldVec(), 'ko')
-    # plt.semilogy(-prb.G(sigma, dm, u=f).fieldVec(), 'rx')
+    plt.semilogy(np.abs(prb.AhVec(sigma+h*dm,f).fieldVec() - prb.AhVec(sigma, f).fieldVec()), 'ko')
+    plt.semilogy(np.abs(h*prb.G(sigma, dm, u=f).fieldVec()), 'rx')
     # plt.semilogy(prb.AhVec(sigma+h*dm, f).fieldVec() - prb.AhVec(sigma, f).fieldVec() - h*prb.G(sigma, dm, u=f).fieldVec(),'ko')
-    # plt.show()
+    plt.show()
 
     # plt.show()
 
diff --git a/simpegEM/Tests/test_forward_EMproblem.py b/simpegEM/Tests/test_forward_EMproblem.py
index f40cbb52..c1b47dcb 100644
--- a/simpegEM/Tests/test_forward_EMproblem.py
+++ b/simpegEM/Tests/test_forward_EMproblem.py
@@ -39,6 +39,34 @@ class TDEM_bTests(unittest.TestCase):
        diff = np.linalg.norm(bz_calc.flatten() - bz_ana.flatten())/np.linalg.norm(bz_ana.flatten())
        self.assertTrue(diff<0.05)
 
+
+class TDEM_bDerivTests(unittest.TestCase):
+
+    def setUp(self):
+
+       cs = 5.
+       ncx = 20
+       ncy = 6
+       npad = 20
+       hx = Utils.meshTensors(((0,cs), (ncx,cs), (npad,cs)))
+       hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs)))
+       mesh = Mesh.Cyl1DMesh([hx,hy], -hy.sum()/2)
+       model = Model.Vertical1DModel(mesh)
+
+       opts = {'txLoc':0.,
+               'txType':'VMD_MVP',
+               'rxLoc':np.r_[150., 0.],
+               'rxType':'bz',
+               'timeCh':np.logspace(-4,-2,20),
+               }
+       self.dat = EM.TDEM.DataTDEM1D(**opts)
+
+       self.prb = EM.TDEM.ProblemTDEM_b(mesh, model)
+       self.prb.setTimes([1e-5, 5e-5, 2.5e-4], [10, 10, 10])
+       self.sigma = np.ones(mesh.nCz)*1e-8
+       self.sigma[mesh.vectorCCz<0] = 0.1
+       self.prb.pair(self.dat)
+
     def test_AhVec(self):
         """
             Test that fields and AhVec produce consistent results
@@ -55,7 +83,21 @@ class TDEM_bTests(unittest.TestCase):
             self.assertTrue(np.linalg.norm(Ahu.get_e(i))/np.linalg.norm(u.get_e(i)) < 1.e-2)
 
 
+    def test_DerivG(self):
+        """
+            Test the derivative of c with respect to sigma
+        """
 
+        # Random model and perturbation
+        sigma = np.random.rand(self.prb.mesh.nCz)
+        f = self.prb.fields(sigma)
+        dm = np.random.rand(self.prb.mesh.nCz)
+        h = 1.
+
+        a = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec())
+        b = np.linalg.norm(self.prb.AhVec(sigma+h*dm, f).fieldVec() - self.prb.AhVec(sigma, f).fieldVec() - h*self.prb.G(sigma, dm, u=f).fieldVec())
+        # Assuming that the gradient is exact to machine precision
+        self.assertTrue(b<1e-16)
 
 
 if __name__ == '__main__':