diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index dd349b43..ac28b3a1 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -190,7 +190,13 @@ class BaseEMProblem(Problem.BaseProblem): """ Derivative of :code:`MfRhoI` with respect to the model. """ - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv + + dMfRhoI_dI = -self.MfRhoI**2 + dMf_drho = self.mesh.getFaceInnerProduct(self.curModel.rho)(u) + drho_dm = self.curModel.rhoDeriv + return dMfRhoI_dI * ( dMf_drho * ( drho_dm)) + + # return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv class BaseEMSurvey(Survey.BaseSurvey): diff --git a/SimPEG/EM/Static/DC/FieldsDC.py b/SimPEG/EM/Static/DC/FieldsDC.py index 91b243c1..46a779f5 100644 --- a/SimPEG/EM/Static/DC/FieldsDC.py +++ b/SimPEG/EM/Static/DC/FieldsDC.py @@ -1,5 +1,6 @@ import SimPEG from SimPEG.Utils import Identity, Zero +import numpy as np class Fields(SimPEG.Problem.Fields): knownFields = {} @@ -10,8 +11,8 @@ class Fields(SimPEG.Problem.Fields): raise NotImplementedError ('Getting phiDerivs from %s is not implemented' %self.knownFields.keys()[0]) if adjoint: - return self._phiDeriv_u(src, v, adjoint), self._phiDeriv_m(src, v, adjoint) - return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = complex) + return self._phiDeriv_u(src, v, adjoint=adjoint), self._phiDeriv_m(src, v, adjoint=adjoint) + return np.array(self._phiDeriv_u(src, du_dm_v, adjoint) + self._phiDeriv_m(src, v, adjoint), dtype = float) def _eDeriv(self, src, du_dm_v, v, adjoint=False): if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: @@ -19,7 +20,7 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) - return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = complex) + return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = float) def _jDeriv(self, src, du_dm_v, v, adjoint=False): if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: @@ -27,7 +28,7 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) - return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex) + return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = float) class Fields_CC(Fields): @@ -57,10 +58,10 @@ class Fields_CC(Fields): def _phi(self, phiSolution, srcList): return phiSolution - def _phiDeriv_u(): + def _phiDeriv_u(self, src, v, adjoint = False): return Identity() - def _phiDeriv_m(): + def _phiDeriv_m(self, src, v, adjoint = False): return Zero() def _j(self, phiSolution, srcList): @@ -96,10 +97,10 @@ class Fields_N(Fields): def _phi(self, phiSolution, srcList): return phiSolution - def _phiDeriv_u(): + def _phiDeriv_u(self, src, v, adjoint = False): return Identity() - def _phiDeriv_m(): + def _phiDeriv_m(self, src, v, adjoint = False): return Zero() def _j(self, phiSolution, srcList): diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index c41c8f62..25ac7fcf 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -4,6 +4,7 @@ from SurveyDC import Survey from FieldsDC import Fields, Fields_CC, Fields_N from SimPEG.Utils import sdiag import numpy as np +from SimPEG.Utils import Zero class BaseDCProblem(BaseEMProblem): @@ -22,7 +23,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 @@ -126,23 +150,32 @@ class Problem3D_CC(BaseDCProblem): """ + V = self.Vol + D = V * self.mesh.faceDiv # TODO: this won't work for full anisotropy - # V = self.Vol - # Div = V*self.mesh.faceDiv MfRhoI = self.MfRhoI - A = self.Div * MfRhoI * self.Div.T + A = D * MfRhoI * D.T + + # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: # return V.T * A return A - def getADeriv(self, u, v, adoint=False): - """ + def getADeriv(self, u, v, adjoint= False): - Product of the derivative of our system matrix with respect to the model and a vector + V = self.Vol + D = V * self.mesh.faceDiv + MfRhoIDeriv = self.MfRhoIDeriv - """ - return Div*self.MfRhoIDeriv(Div.T*u) + if adjoint: + # if self._makeASymmetric is True: + # v = V * v + return D * MfRhoIDeriv(D * v) + # I think we should deprecate this for DC problem. + # if self._makeASymmetric is True: + # return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) ) + return D * (MfRhoIDeriv( D.T * u ) * v) def getRHS(self): """ @@ -152,16 +185,21 @@ class Problem3D_CC(BaseDCProblem): """ RHS = self.getSourceTerm() + + # I think we should deprecate this for DC problem. # if self._makeASymmetric is True: # return self.Vol.T * RHS + return RHS def getRHSDeriv(self, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model """ - qDeriv = src.evalDeriv(self, adjoint=adjoint) - return qDeriv + # TODO: add qDeriv for RHS depending on m + # qDeriv = src.evalDeriv(self, adjoint=adjoint) + # return qDeriv + return Zero() diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 9fbd05cb..7692c880 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -36,6 +36,10 @@ class BaseRx(SimPEG.Survey.BaseRx): P = self.getP(mesh, self.projGLoc(f)) return P*f[src, self.projField] + def evalDeriv(self, src, mesh, f, v, adjoint=False): + P = self.getP(mesh, self.projGLoc(f)) + return P*v + # DC.Rx.Dipole(locs) class Dipole(BaseRx): @@ -50,6 +54,10 @@ class Dipole(BaseRx): """Number of data in the receiver.""" return self.locs[0].shape[0] + # Not sure why ... + # return int(self.locs[0].size / 2) + + def getP(self, mesh, Gloc): if mesh in self._Ps: return self._Ps[mesh] @@ -63,7 +71,6 @@ class Dipole(BaseRx): return P - # class Pole(BaseRx): diff --git a/SimPEG/EM/Static/DC/Utils.py b/SimPEG/EM/Static/DC/Utils.py new file mode 100644 index 00000000..626ae368 --- /dev/null +++ b/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).reshape([1,-1]),getLoc(i,2).reshape([1,-1])) + src = DC.Src.Dipole([rx], getLoc(i,0),getLoc(i,3)) + srcList += [src] + + return srcList diff --git a/SimPEG/EM/Static/DC/__init__.py b/SimPEG/EM/Static/DC/__init__.py index a3e1eba5..dcc40764 100644 --- a/SimPEG/EM/Static/DC/__init__.py +++ b/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 diff --git a/tests/em/static/__init__.py b/tests/em/static/__init__.py new file mode 100644 index 00000000..420388ef --- /dev/null +++ b/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) diff --git a/tests/em/static/test_DC.py b/tests/em/static/test_DC.py new file mode 100644 index 00000000..7de2ea45 --- /dev/null +++ b/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() diff --git a/tests/em/static/test_DC_deriv.py b/tests/em/static/test_DC_deriv.py new file mode 100644 index 00000000..e69de29b