diff --git a/SimPEG/EM/Static/SIP/ProblemSIP.py b/SimPEG/EM/Static/SIP/ProblemSIP.py index e4ab33d8..88b85a7b 100644 --- a/SimPEG/EM/Static/SIP/ProblemSIP.py +++ b/SimPEG/EM/Static/SIP/ProblemSIP.py @@ -34,13 +34,20 @@ class BaseSIPProblem(BaseEMProblem): peta = self.curModel.eta*np.exp(-self.curModel.taui*t) return peta - def EtaDeriv(self, t, v): + def EtaDeriv(self, t, v, adjoint=False): v = np.array(v, dtype=float) - return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) + if adjoint: + return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v) + else: + return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v) - def TauiDeriv(self, t, v): + + def TauiDeriv(self, t, v, adjoint=False): v = np.array(v, dtype=float) - return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) + if adjoint: + return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v) + else: + return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v) def fields(self, m): self.curModel = m @@ -152,7 +159,7 @@ class BaseSIPProblem(BaseEMProblem): dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT), self.TauiDeriv(self.survey.times[tind], du_dmT)] + Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)] # Conductivity ((d u / d log sigma).T) if self._formulation is 'EB': diff --git a/tests/em/static/test_SIP_jvecjtvecadj.py b/tests/em/static/test_SIP_jvecjtvecadj.py index 6e6f071d..7cdb6def 100644 --- a/tests/em/static/test_SIP_jvecjtvecadj.py +++ b/tests/em/static/test_SIP_jvecjtvecadj.py @@ -150,5 +150,83 @@ class IPProblemTestsN(unittest.TestCase): passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3) self.assertTrue(passed) +class IPProblemTestsN_air(unittest.TestCase): + + def setUp(self): + + cs = 25. + hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)] + hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)] + mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC") + blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC) + blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC) + sigma = np.ones(mesh.nC)*1e-2 + airind = mesh.gridCC[:,2]>0. + sigma[airind] = 1e-8 + eta = np.zeros(mesh.nC) + tau = np.ones_like(sigma)*1. + eta[blkind0] = 0.1 + eta[blkind1] = 0.1 + tau[blkind0] = 0.1 + tau[blkind1] = 0.01 + + actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.) + actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.) + + x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)] + y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)] + Aloc = np.r_[-200., 0., 0.] + Bloc = np.r_[200., 0., 0.] + M = Utils.ndgrid(x-25.,y, np.r_[0.]) + N = Utils.ndgrid(x+25.,y, np.r_[0.]) + + times = np.arange(10)*1e-3 + 1e-3 + rx = SIP.Rx.Dipole(M, N, times) + src = SIP.Src.Dipole([rx], Aloc, Bloc) + survey = SIP.Survey([src]) + colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)] + problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap) + problem.Solver = MumpsSolver + problem.pair(survey) + mSynth = np.r_[eta[~airind], 1./tau[~airind]] + survey.makeSyntheticData(mSynth) + # Now set up the problem to do some minimization + dmis = DataMisfit.l2_DataMisfit(survey) + regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2)) + reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind) + 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, num=3) + 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-8 + 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, num=3) + self.assertTrue(passed) + if __name__ == '__main__': unittest.main()