start of the sketch of Jvec (not to be trusted yet!)

SimPEG/EM/Static/DC/ProblemDC.py

@@ -3,6 +3,7 @@ from SimPEG.EM.Base import BaseEMProblem
 from SurveyDC import Survey
 from FieldsDC import Fields, Fields_CC
 import numpy as np
+from SimPEG.Utils import Zero
 
 class BaseDCProblem(BaseEMProblem):
 
@@ -21,7 +22,30 @@ class BaseDCProblem(BaseEMProblem):
         return f
 
     def Jvec(self, m, v, f=None):
-        raise NotImplementedError
+
+        if f is None:
+            f = self.fields(m)
+
+        self.curModel = m
+
+        Jv = self.dataPair(self.survey) #same size as the data
+
+        A = self.getA()
+        Ainv = self.Solver(A, **self.solverOpts)
+
+        for src in self.survey.srcList:
+            u_src = f[src, self._solutionType] # solution vector
+            dA_dm_v = self.getADeriv(u_src, v)
+            dRHS_dm_v = self.getRHSDeriv(src, v)
+            print type(dA_dm_v + dRHS_dm_v), (dA_dm_v + dRHS_dm_v).shape
+            du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
+
+            for rx in src.rxList:
+                df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
+                df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
+                Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
+            Ainv.clean()
+        return Utils.mkvc(Jv)
 
     def Jtvec(self, m, v, f=None):
         raise NotImplementedError
@@ -85,6 +109,18 @@ class Problem3D_CC(BaseDCProblem):
 
     def getADeriv(self, u, v, adjoint= False):
 
+        D = self.mesh.faceDiv
+        MfRhoIDeriv = self.MfRhoIDeriv
+        V = self.Vol
+
+        if adjoint:
+            if self._makeASymmetric is True:
+                v = V * v
+            return V.T * ( D * ( MfRhoIDeriv(D * v) ) )
+
+        if self._makeASymmetric is True:
+            return V.T * ( D * ( MfRhoIDeriv( * D.T * ( V * u ) ) * v ) )
+        return D * ( MfRhoIDeriv( D.T * ( V * v ) ) )
 
     def getRHS(self):
         """
@@ -98,8 +134,8 @@ class Problem3D_CC(BaseDCProblem):
             return self.Vol.T * RHS
         return RHS
 
-    def getRHSDeriv():
-        raise NotImplementedError
+    def getRHSDeriv(self, src, v, adjoint=False):
+        return Zero()
 
 
 

SimPEG/EM/Static/DC/RxDC.py

@@ -48,7 +48,7 @@ class Dipole(BaseRx):
     @property
     def nD(self):
         """Number of data in the receiver."""
-        return self.locs[0].shape[0]
+        return int(self.locs[0].size / 2)
 
     def getP(self, mesh, Gloc):
         if mesh in self._Ps:

SimPEG/EM/Static/DC/Utils.py

@@ -0,0 +1,38 @@
+import numpy as np
+
+def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
+
+    import SimPEG.EM.Static.DC as DC
+
+    elocs = np.arange(0,aSpacing*nElecs,aSpacing)
+    elocs -= (nElecs*aSpacing - aSpacing)/2
+    space = 1
+    WENNER = np.zeros((0,),dtype=int)
+    for ii in range(nElecs):
+        for jj in range(nElecs):
+            test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
+            if np.any(test >= nElecs):
+                break
+            WENNER = np.r_[WENNER, test]
+        space += 1
+    WENNER = WENNER.reshape((-1,4))
+
+
+    if plotIt:
+        for i, s in enumerate('rbkg'):
+            plt.plot(elocs[WENNER[:,i]],s+'.')
+        plt.show()
+
+    # Create sources and receivers
+    i = 0
+    if in2D:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
+    else:
+        getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
+    srcList = []
+    for i in range(WENNER.shape[0]):
+        rx = DC.Rx.Dipole(getLoc(i,1),getLoc(i,2))
+        src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3))
+        srcList += [src]
+
+    return srcList

SimPEG/EM/Static/DC/__init__.py

@@ -3,3 +3,4 @@ from SurveyDC import Survey
 import SrcDC as Src #Pole
 import RxDC as Rx
 from FieldsDC import Fields_CC
+import Utils

tests/em/static/__init__.py

@@ -0,0 +1,12 @@
+import os
+import glob
+import unittest
+
+if __name__ == '__main__':
+    test_file_strings = glob.glob('test_*.py')
+    module_strings = [str[0:len(str)-3] for str in test_file_strings]
+    suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
+              in module_strings]
+    testSuite = unittest.TestSuite(suites)
+
+    unittest.TextTestRunner(verbosity=2).run(testSuite)

tests/em/static/test_DC.py

@@ -0,0 +1,77 @@
+import unittest
+from SimPEG import *
+import SimPEG.EM.Static.DC as DC
+
+
+class DCProblemTests(unittest.TestCase):
+
+    def setUp(self):
+
+        aSpacing=2.5
+        nElecs=10
+
+        surveySize = nElecs*aSpacing - aSpacing
+        cs = surveySize/nElecs/4
+
+        mesh = Mesh.TensorMesh([
+                [(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
+                [(cs,3, -1.3),(cs,3,1.3)],
+        #         [(cs,5, -1.3),(cs,10)]
+            ],'CN')
+
+        srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
+        survey = DC.Survey(srcList)
+        problem = DC.Problem3D_CC(mesh, mapping=[('rho', Maps.IdentityMap(mesh))])
+        problem.pair(survey)
+
+        mSynth = np.ones(mesh.nC)
+        survey.makeSyntheticData(mSynth)
+
+        # Now set up the problem to do some minimization
+        dmis = DataMisfit.l2_DataMisfit(survey)
+        reg = Regularization.Tikhonov(mesh)
+        opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
+        invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
+        inv = Inversion.BaseInversion(invProb)
+
+        self.inv = inv
+        self.reg = reg
+        self.p =     problem
+        self.mesh = mesh
+        self.m0 = mSynth
+        self.survey = survey
+        self.dmis = dmis
+
+    def test_misfit(self):
+        derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+    def test_adjoint(self):
+        # Adjoint Test
+        u = np.random.rand(self.mesh.nC*self.survey.nSrc)
+        v = np.random.rand(self.mesh.nC)
+        w = np.random.rand(self.survey.dobs.shape[0])
+        wtJv = w.dot(self.p.Jvec(self.m0, v))
+        vtJtw = v.dot(self.p.Jtvec(self.m0, w))
+        passed = np.abs(wtJv - vtJtw) < 1e-10
+        print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
+        self.assertTrue(passed)
+
+    def test_dataObj(self):
+        derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
+    def test_massMatrices(self):
+        Gu = np.random.rand(self.mesh.nF)
+        def derChk(m):
+            self.p.curModel = m
+            return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
+        passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
+        self.assertTrue(passed)
+
+
+if __name__ == '__main__':
+    unittest.main()