diff --git a/SimPEG/EM/Static/DC/FieldsDC_2D.py b/SimPEG/EM/Static/DC/FieldsDC_2D.py index 77e3199e..5b75031d 100644 --- a/SimPEG/EM/Static/DC/FieldsDC_2D.py +++ b/SimPEG/EM/Static/DC/FieldsDC_2D.py @@ -1,5 +1,6 @@ import SimPEG -import Utils, numpy as np, scipy.sparse as sp +from SimPEG.Utils import Identity, Zero +import numpy as np class Fields_ky(SimPEG.Problem.TimeFields): diff --git a/SimPEG/EM/Static/DC/ProblemDC_2D.py b/SimPEG/EM/Static/DC/ProblemDC_2D.py index ee04a560..8ebb9c67 100644 --- a/SimPEG/EM/Static/DC/ProblemDC_2D.py +++ b/SimPEG/EM/Static/DC/ProblemDC_2D.py @@ -12,7 +12,7 @@ class BaseDCProblem_2D(BaseEMProblem): surveyPair = Survey_ky fieldsPair = Fields_ky nky = 15 - ky = np.logspace(-4, 1, nky) + kys = np.logspace(-4, 1, nky) Ainv = [None for i in range(nky)] nT = nky # Only for using TimeFields @@ -26,7 +26,7 @@ class BaseDCProblem_2D(BaseEMProblem): f = self.fieldsPair(self.mesh, self.survey) Srcs = self.survey.srcList for iky in range(self.nky): - ky = self.ky[iky] + ky = self.kys[iky] A = self.getA(ky) self.Ainv[iky] = self.Solver(A, **self.solverOpts) RHS = self.getRHS(ky) @@ -34,28 +34,44 @@ class BaseDCProblem_2D(BaseEMProblem): f[Srcs, self._solutionType, iky] = u return f - # def Jvec(self, m, v, f=None): + def Jvec(self, m, v, f=None): - # if f is None: - # f = self.fields(m) + if f is None: + f = self.fields(m) - # self.curModel = m + self.curModel = m - # Jv = self.dataPair(self.survey) #same size as the data + Jv = self.dataPair(self.survey) #same size as the data + Jv0 = self.dataPair(self.survey) - # A = self.getA() + # Assume y=0. + # This needs some thoughts to implement in general when src is dipole + dky = np.diff(self.kys) + dky = np.r_[dky[0], dky] + y = 0. - # 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) - # du_dm_v = self.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) - # return Utils.mkvc(Jv) + for iky in range(self.nky): + ky = self.kys[iky] + A = self.getA(ky) + for src in self.survey.srcList: + u_src = f[src, self._solutionType, iky] # solution vector + dA_dm_v = self.getADeriv(ky, u_src, v) + dRHS_dm_v = self.getRHSDeriv(ky, src, v) + du_dm_v = self.Ainv[iky] * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: + df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dm_v = df_dmFun(iky, src, du_dm_v, v, adjoint=False) + # Trapezoidal intergration + Jv1_temp = 1./np.pi*rx.evalDeriv(ky, src, self.mesh, f, df_dm_v) + if iky==0: + #First assigment + Jv[src, rx] = Jv1_temp*dky[iky]*np.cos(ky*y) + else: + Jv[src, rx] += Jv1_temp*dky[iky] /2.*np.cos(ky*y) + Jv[src, rx] += Jv0[src, rx]*dky[iky]/2.*np.cos(ky*y) + Jv0[src, rx] = Jv1_temp.copy() + JV[iky,isrc,:] = Jv1_temp.copy() + return Utils.mkvc(Jv) # def Jtvec(self, m, v, f=None): # if f is None: @@ -146,11 +162,13 @@ class Problem2D_CC(BaseDCProblem_2D): D = self.Div G = self.Grad + vol = self.mesh.vol MfRhoIDeriv = self.MfRhoIDeriv - + rho = self.curModel.rho if adjoint: - return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + Utils.sdiag(ky**2*mesh.vol)*v - return D * ((MfRhoIDeriv( G * u )) * v) + Utils.sdiag(ky**2*mesh.vol)*v + return(MfRhoIDeriv( G * u ).T) * ( D.T * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v + + return D * ((MfRhoIDeriv( G * u )) * v) + ky**2*Utils.sdiag(u.flatten()*vol*(-1./rho**2))*v def getRHS(self, ky): """ diff --git a/SimPEG/EM/Static/DC/RxDC.py b/SimPEG/EM/Static/DC/RxDC.py index f7c7d352..d2ec6098 100644 --- a/SimPEG/EM/Static/DC/RxDC.py +++ b/SimPEG/EM/Static/DC/RxDC.py @@ -102,10 +102,10 @@ class Dipole_ky(BaseRx): self._Ps[mesh] = P return P - def eval(self, ky, src, mesh, f): + def eval(self, kys, src, mesh, f): P = self.getP(mesh, self.projGLoc(f)) Pf = P*f[src, self.projField,:] - return self.IntTrapezoidal(ky, Pf, y=0.) + return self.IntTrapezoidal(kys, Pf, y=0.) def evalDeriv(self, ky, src, mesh, f, v, adjoint=False): P = self.getP(mesh, self.projGLoc(f)) @@ -114,16 +114,16 @@ class Dipole_ky(BaseRx): elif adjoint: return P.T*v - def IntTrapezoidal(self, ky, Pf, y=0.): + def IntTrapezoidal(self, kys, Pf, y=0.): phi = np.zeros(Pf.shape[0]) - nky = ky.size - dky = np.diff(ky) + nky = kys.size + dky = np.diff(kys) dky = np.r_[dky[0], dky] - phi0 = Pf[:,0] + phi0 = 1./np.pi*Pf[:,0] for iky in range(nky): - phi1 = 2./np.pi*Pf[:,iky]/2. - phi += phi1*dky[iky]/2.*np.cos(ky[iky]*y) - phi += phi0*dky[iky]/2.*np.cos(ky[iky]*y) + phi1 = 1./np.pi*Pf[:,iky] + phi += phi1*dky[iky]/2.*np.cos(kys[iky]*y) + phi += phi0*dky[iky]/2.*np.cos(kys[iky]*y) phi0 = phi1.copy() return phi diff --git a/SimPEG/EM/Static/DC/SrcDC.py b/SimPEG/EM/Static/DC/SrcDC.py index 60a53dff..5f4ac0e2 100644 --- a/SimPEG/EM/Static/DC/SrcDC.py +++ b/SimPEG/EM/Static/DC/SrcDC.py @@ -36,11 +36,6 @@ class Dipole(BaseSrc): q = self.current * mkvc(qa+qb) return q - # def bc_contribution - - -# How to treat boundary conditions here - class Pole(BaseSrc): def __init__(self, rxList, loc, **kwargs): diff --git a/SimPEG/EM/Static/DC/SurveyDC.py b/SimPEG/EM/Static/DC/SurveyDC.py index 62e2922a..fb3d49a7 100644 --- a/SimPEG/EM/Static/DC/SurveyDC.py +++ b/SimPEG/EM/Static/DC/SurveyDC.py @@ -29,8 +29,8 @@ class Survey_ky(BaseEMSurvey): :return: data """ data = SimPEG.Survey.Data(self) - ky = self.prob.ky + kys = self.prob.kys for src in self.srcList: for rx in src.rxList: - data[src, rx] = rx.eval(ky, src, self.mesh, f) + data[src, rx] = rx.eval(kys, src, self.mesh, f) return data diff --git a/tests/em/static/test_DC_2D_jvecjtvecadj.py b/tests/em/static/test_DC_2D_jvecjtvecadj.py new file mode 100644 index 00000000..ad7198e9 --- /dev/null +++ b/tests/em/static/test_DC_2D_jvecjtvecadj.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +import SimPEG.EM.Static.DC as DC + + +class DCProblem_2DTestsCC(unittest.TestCase): + + def setUp(self): + + cs = 12.5 + hx = [(cs,7, -1.3),(cs,61),(cs,7, 1.3)] + hy = [(cs,7, -1.3),(cs,20)] + mesh = Mesh.TensorMesh([hx, hy],x0="CN") + x = np.linspace(-135, 250., 20) + M = Utils.ndgrid(x-12.5, np.r_[0.]) + N = Utils.ndgrid(x+12.5, np.r_[0.]) + A0loc = np.r_[-150, 0.] + A1loc = np.r_[-130, 0.] + rxloc = [np.c_[M, np.zeros(20)], np.c_[N, np.zeros(20)]] + rx = DC.Rx.Dipole_ky(M, N) + src0 = DC.Src.Pole([rx], A0loc) + src1 = DC.Src.Pole([rx], A1loc) + survey = DC.Survey_ky([src0, src1]) + problem = DC.Problem2D_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, 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-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, num=3) + # self.assertTrue(passed) + +# class DCProblemTestsN(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_N(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-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) +# self.assertTrue(passed) + +if __name__ == '__main__': + unittest.main()