From cd5339322eb545c0ef0759e56e3de756a86506ea Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Thu, 21 Apr 2016 11:00:20 -0700 Subject: [PATCH 1/3] sketch of DC --- SimPEG/EM/Static/DC/ProblemDC.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index f712e036..cc9c1ebb 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/SimPEG/EM/Static/DC/ProblemDC.py @@ -83,8 +83,8 @@ class Problem3D_CC(BaseDCProblem): return V.T * A return A - def getADeriv(): - raise NotImplementedError + def getADeriv(self, u, v, adjoint= False): + def getRHS(self): """ From 5ec6e79a39ab85549f71f03394411953e41df8a4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 22 Apr 2016 17:47:36 -0700 Subject: [PATCH 2/3] MfRhoIDeriv --- SimPEG/EM/Base.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index dd349b43..0df18cdb 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.getEdgeInnerProductDeriv(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): From 8775364d8f4b56d92e801c2df99cca4614a3a74a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 22 Apr 2016 17:48:26 -0700 Subject: [PATCH 3/3] start of the sketch of Jvec (not to be trusted yet!) --- SimPEG/EM/Static/DC/ProblemDC.py | 42 +++++++++++++++-- SimPEG/EM/Static/DC/RxDC.py | 2 +- SimPEG/EM/Static/DC/Utils.py | 38 ++++++++++++++++ SimPEG/EM/Static/DC/__init__.py | 1 + tests/em/static/__init__.py | 12 +++++ tests/em/static/test_DC.py | 77 ++++++++++++++++++++++++++++++++ tests/em/static/test_DC_deriv.py | 0 7 files changed, 168 insertions(+), 4 deletions(-) create mode 100644 SimPEG/EM/Static/DC/Utils.py create mode 100644 tests/em/static/__init__.py create mode 100644 tests/em/static/test_DC.py create mode 100644 tests/em/static/test_DC_deriv.py diff --git a/SimPEG/EM/Static/DC/ProblemDC.py b/SimPEG/EM/Static/DC/ProblemDC.py index cc9c1ebb..602ce3d8 100644 --- a/SimPEG/EM/Static/DC/ProblemDC.py +++ b/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() diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index 9fbd05cb..779ffde0 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/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: diff --git a/SimPEG/EM/Static/DC/Utils.py b/SimPEG/EM/Static/DC/Utils.py new file mode 100644 index 00000000..647590cc --- /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),getLoc(i,2)) + 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..a25dcf67 --- /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