diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 97f08fc4..fd510abe 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -20,6 +20,7 @@ class BaseProblemFDEM(Problem.BaseProblem): def __init__(self, model, **kwargs): Problem.BaseProblem.__init__(self, model, **kwargs) + solType = None storeTheseFields = ['e', 'b'] surveyPair = SurveyFDEM @@ -49,7 +50,7 @@ class BaseProblemFDEM(Problem.BaseProblem): def MeSigma(self): #TODO: hardcoded to sigma as the model if getattr(self, '_MeSigma', None) is None: - sigma = self.currentTransformedModel + sigma = self.curTModel self._MeSigma = self.mesh.getEdgeInnerProduct(sigma) return self._MeSigma @@ -60,10 +61,22 @@ class BaseProblemFDEM(Problem.BaseProblem): self._MeSigmaI = Utils.sdiag(1/self.MeSigma.diagonal()) return self._MeSigmaI - currentTransformedModel = Utils.dependentProperty('_currentTransformedModel', None, ['_MeSigma', '_MeSigmaI'], 'Sets the current model, and removes dependent mass matrices.') + curModel = Utils.dependentProperty('_curModel', None, ['_MeSigma', '_MeSigmaI', '_curTModel', '_curTModelDeriv'], 'Sets the current model, and removes dependent mass matrices.') + + @property + def curTModel(self): + if getattr(self, '_curTModel', None) is None: + self._curTModel = self.model.transform(self.curModel) + return self._curTModel + + @property + def curTModelDeriv(self): + if getattr(self, '_curTModelDeriv', None) is None: + self._curTModelDeriv = self.model.transformDeriv(self.curModel) + return self._curTModelDeriv def fields(self, m): - self.currentTransformedModel = self.model.transform(m) + self.curModel = m F = self.forward(m, self.getRHS, self.calcFields) return F @@ -80,10 +93,54 @@ class BaseProblemFDEM(Problem.BaseProblem): return F + def Jvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + + self.curModel = m + + Jv = self.dataPair(self.survey) + + for freq in self.survey.freqs: + A = self.getA(freq) + solver = self.Solver(A, **self.solverOpts) + + for tx in self.survey.getTransmitters(freq): + w = self.getADeriv(freq, u[tx, self.solType], v) + Ainvw = solver.solve(w) + P = tx.projectFieldsDeriv(self.mesh, u) + Jv[tx] = -P*Ainvw + + return Utils.mkvc(Jv) + + def Jtvec(self, m, v, u=None): + if u is None: + u = self.fields(m) + + self.curModel = m + + # Ensure v is a data object. + if not isinstance(v, self.dataPair): + v = self.dataPair(self.survey, v) + + Jtv = np.zeros(self.model.nP, dtype=complex) + + for freq in self.survey.freqs: + AT = self.getA(freq).T + solver = self.Solver(AT, **self.solverOpts) + + for tx in self.survey.getTransmitters(freq): + P = tx.projectFieldsDeriv(self.mesh, u) + w = solver.solve( - P.T * v[tx]) + Jtv += self.getADeriv(freq, u[tx, self.solType], w, adjoint=True) + + return Jtv + class ProblemFDEM_e(BaseProblemFDEM): """ Solving for e! """ + solType = 'e' def __init__(self, model, **kwargs): BaseProblemFDEM.__init__(self, model, **kwargs) @@ -100,6 +157,16 @@ class ProblemFDEM_e(BaseProblemFDEM): return C.T*mui*C + 1j*omega(freq)*sig + def getADeriv(self, freq, u, v, adjoint=False): + sig = self.curTModel + dsig_dm = self.curTModelDeriv + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=u) + + if adjoint: + return 1j * omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * v ) ) + + return 1j * omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) + def getRHS(self, freq): """ :param float freq: Frequency @@ -139,65 +206,12 @@ class ProblemFDEM_e(BaseProblemFDEM): return fDict - def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) - - sig = self.model.transform(m) - self.currentTransformedModel = sig - - Jv = self.dataPair(self.survey) - dsig_dm = self.model.transformDeriv(m) - - for i, freq in enumerate(self.survey.freqs): - e = u[freq, 'e'] - A = self.getA(freq) - solver = self.Solver(A, **self.solverOpts) - - for txi, tx in enumerate(self.survey.getTransmitters(freq)): - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) - - P = tx.projectFieldsDeriv(self.mesh, u) - b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) - Ainvb = solver.solve(b) - Jv[tx] = -P*Ainvb - - return Utils.mkvc(Jv) - - - def Jtvec(self, m, v, u=None): - if u is None: - u = self.fields(m) - - sig = self.model.transform(m) - self.currentTransformedModel = sig - - # Ensure v is a data object. - if not isinstance(v, self.dataPair): - v = self.dataPair(self.survey, v) - - Jtv = np.zeros(self.model.nP, dtype=complex) - - dsig_dm = self.model.transformDeriv(m) - - for i, freq in enumerate(self.survey.freqs): - e = u[freq, 'e'] - AT = self.getA(freq).T - solver = self.Solver(AT, **self.solverOpts) - - for txi, tx in enumerate(self.survey.getTransmitters(freq)): - dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) - - P = tx.projectFieldsDeriv(self.mesh, u) - w = solver.solve(P.T * v[tx]) - Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) ) - - return Jtv class ProblemFDEM_b(BaseProblemFDEM): """ Solving for b! """ + solType = 'b' def __init__(self, model, **kwargs): BaseProblemFDEM.__init__(self, model, **kwargs) @@ -214,6 +228,23 @@ class ProblemFDEM_b(BaseProblemFDEM): return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui + def getADeriv(self, freq, u, v, adjoint=False): + + mui = self.MfMui + C = self.mesh.edgeCurl + sig = self.curTModel + dsig_dm = self.curTModelDeriv + #TODO: This only works if diagonal (no tensors)... + dMeSigmaI_dI = - self.MeSigmaI**2 + + vec = (C.T*(mui*u)) + dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=vec) + + if adjoint: + return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) ) + + return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) + def getRHS(self, freq): """ :param float freq: Frequency @@ -253,61 +284,3 @@ class ProblemFDEM_b(BaseProblemFDEM): return fDict - def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) - - raise NotImplemented('') - - # sig = self.model.transform(m) - # self.currentTransformedModel = sig - - # Jv = self.dataPair(self.survey) - # dsig_dm = self.model.transformDeriv(m) - - # for i, freq in enumerate(self.survey.freqs): - # e = u[freq, 'e'] - # A = self.getA(freq) - # solver = self.Solver(A, **self.solverOpts) - - # for txi, tx in enumerate(self.survey.getTransmitters(freq)): - # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) - - # P = tx.projectFieldsDeriv(self.mesh, u) - # b = 1j*omega(freq) * ( dMe_dsig * ( dsig_dm * v ) ) - # Ainvb = solver.solve(b) - # Jv[tx] = -P*Ainvb - - # return Utils.mkvc(Jv) - - - def Jtvec(self, m, v, u=None): - if u is None: - u = self.fields(m) - - # Ensure v is a data object. - if not isinstance(v, self.dataPair): - v = self.dataPair(self.survey, v) - - raise NotImplemented('') - - # sig = self.model.transform(m) - # self.currentTransformedModel = sig - - # Jtv = np.zeros(self.model.nP, dtype=complex) - - # dsig_dm = self.model.transformDeriv(m) - - # for i, freq in enumerate(self.survey.freqs): - # e = u[freq, 'e'] - # AT = self.getA(freq).T - # solver = self.Solver(AT, **self.solverOpts) - - # for txi, tx in enumerate(self.survey.getTransmitters(freq)): - # dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig, v=e[:,txi]) - - # P = tx.projectFieldsDeriv(self.mesh, u) - # w = solver.solve(P.T * v[tx]) - # Jtv += - 1j*omega(freq) * ( dsig_dm.T * ( dMe_dsig.T * w ) ) - - # return Jtv diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 82404962..4bb4d1c7 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -2,7 +2,15 @@ from SimPEG import Survey, Utils, np, sp class RxListFDEM(Survey.BaseRxList): - knownRxTypes = ['Ex', 'Ey', 'Ez'] + knownRxTypes = { + 'Ex':'Ex', + 'Ey':'Ey', + 'Ez':'Ez', + + 'Bx':'Fx', + 'By':'Fy', + 'Bz':'Fz', + } def __init__(self, locs, rxType): Survey.BaseRxList.__init__(self, locs, rxType) @@ -11,7 +19,7 @@ class RxListFDEM(Survey.BaseRxList): def getP(self, mesh): if mesh not in self._Ps: - self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.rxType) + self._Ps[mesh] = mesh.getInterpolationMat(self.locs, self.knownRxTypes[self.rxType]) return self._Ps[mesh] @@ -37,6 +45,8 @@ class TxFDEM(Survey.BaseTx): if self.rxList.rxType in ['Ex', 'Ey', 'Ez']: u_part = u[self, 'e'] + elif self.rxList.rxType in ['Bx', 'By', 'Bz']: + u_part = u[self, 'b'] else: raise NotImplemented('Unknown receiver type.') diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 3726427c..c515d061 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -4,44 +4,90 @@ import simpegEM as EM TOL = 1e-10 -class FDEM_bDerivTests(unittest.TestCase): +def getProblem(fdemType): + cs = 5. + ncx, ncy, ncz = 2, 2, 2 + npad = 3 + hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs))) + hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) + hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) + mesh = Mesh.TensorMesh([hx,hy,hz]) + + model = Model.LogModel(mesh) + + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0]) + if fdemType == 'e': + rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') + elif fdemType == 'b': + rxList = EM.FDEM.RxListFDEM(XYZ, 'Bx') + else: + raise NotImplementedError() + + Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList) + + x = np.linspace(5,10,3) + XYZ = Utils.ndgrid(x,x,np.r_[0]) + if fdemType == 'e': + rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey') + elif fdemType == 'b': + rxList = EM.FDEM.RxListFDEM(XYZ, 'By') + else: + raise NotImplementedError() + + Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList) + + survey = EM.FDEM.SurveyFDEM([Tx0, Tx1]) + + if fdemType == 'e': + prb = EM.FDEM.ProblemFDEM_e(model) + elif fdemType == 'b': + prb = EM.FDEM.ProblemFDEM_b(model) + else: + raise NotImplementedError() + prb.pair(survey) + + return prb + +class FDEM_DerivTests_e(unittest.TestCase): def setUp(self): - - cs = 5. - ncx, ncy, ncz = 2, 2, 2 - npad = 3 - hx = Utils.meshTensors(((npad,cs), (ncx,cs), (npad,cs))) - hy = Utils.meshTensors(((npad,cs), (ncy,cs), (npad,cs))) - hz = Utils.meshTensors(((npad,cs), (ncz,cs), (npad,cs))) - mesh = Mesh.TensorMesh([hx,hy,hz]) - - model = Model.LogModel(mesh) - - x = np.linspace(5,10,3) - XYZ = Utils.ndgrid(x,x,np.r_[0]) - rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex') - Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, rxList) - - x = np.linspace(5,10,3) - XYZ = Utils.ndgrid(x,x,np.r_[0]) - rxList = EM.FDEM.RxListFDEM(XYZ, 'Ey') - Tx1 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-4, rxList) - - survey = EM.FDEM.SurveyFDEM([Tx0, Tx1]) - - prb = EM.FDEM.ProblemFDEM_e(model) - prb.pair(survey) - - self.sigma = np.log(np.ones(mesh.nC)*1e-3) + prb = getProblem('e') self.prb = prb - self.survey = survey + self.sigma = np.log(np.ones(prb.mesh.nC)*1e-3) + self.survey = prb.survey def test_Jvec(self): x0 = self.sigma def fun(x): return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) - passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-18) + passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25) + self.assertTrue(passed) + + def test_Jtvec_adjointTest(self): + v = np.random.rand(self.survey.nD) + w = np.random.rand(self.prb.model.nP) + + m = self.sigma + u = self.prb.fields(m) + vJw = v.dot(self.prb.Jvec(m, w, u=u)) + wJtv = w.dot(self.prb.Jtvec(m, v, u=u)) + self.assertTrue(vJw - wJtv < TOL) + + +class FDEM_DerivTests_b(unittest.TestCase): + + def setUp(self): + prb = getProblem('e') + self.prb = prb + self.sigma = np.log(np.ones(prb.mesh.nC)*1e-3) + self.survey = prb.survey + + def test_Jvec(self): + x0 = self.sigma + def fun(x): + return self.survey.dpred(x), lambda x: self.prb.Jvec(x0, x) + passed = Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25) self.assertTrue(passed) def test_Jtvec_adjointTest(self):