diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 38172087..acd3af6c 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -94,6 +94,7 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma + # TODO: This should take a vector def MeSigmaDeriv(self, u): """ Deriv of MeSigma wrt sigma @@ -107,6 +108,7 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI + # TODO: This should take a vector def MeSigmaIDeriv(self, u): """ Deriv of MeSigma wrt sigma @@ -126,8 +128,9 @@ class BaseEMProblem(Problem.BaseProblem): self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) return self._MfRho + # TODO: This should take a vector def MfRhoDeriv(self,u): - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * self.curModel.rhoDeriv @property def MfRhoI(self): @@ -135,5 +138,7 @@ class BaseEMProblem(Problem.BaseProblem): self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) return self._MfRhoI + # TODO: This isn't going to work yet + # TODO: This should take a vector def dMfRhoIDeriv(self,u): - return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True) + return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index 58d7fcba..8db4474f 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -442,6 +442,7 @@ class ProblemFDEM_j(BaseFDEMProblem): MeMuI = self.MeMuI MfRho = self.MfRho + # print MfRho.max C = self.mesh.edgeCurl iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) @@ -463,19 +464,25 @@ class ProblemFDEM_j(BaseFDEMProblem): MeMuI = self.MeMuI MfRho = self.MfRho C = self.mesh.edgeCurl - sigi = self.sigmai - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) + MfRhoDeriv_m = self.MfRhoDeriv(u) + + # sigi = self.curModel.rho + # dsig_dm = self.curModel.rhoDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) + # MfRhoDeriv_m = (dMf_dsigi * ( dsigi_dsig * ( dsig_dm))) if adjoint: if self._makeASymmetric is True: v = MfRho * v - return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) + # return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) + return MfRhoDeriv_m.T * (C * (MeMuI.T * v)) if self._makeASymmetric is True: - return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) - return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) + # return MfRho.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) + return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) ))) + # return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) + return C * (MeMuI * (C.T * (MfRhoDeriv_m * v))) def getRHS(self, freq): @@ -497,9 +504,40 @@ class ProblemFDEM_j(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, v, adjoint=False): - raise NotImplementedError('getRHSDeriv not implemented yet') - return None + def getRHSDeriv(self, src, v, adjoint=False): + C = self.mesh.edgeCurl + MeMuI = self.MeMuI + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + + if adjoint: + if self._makeASymmetric: + v = MfRho*v + S_mDerivv = S_mDeriv(MeMuI.T * (C.T * v)) + S_eDerivv = S_eDeriv(v) + if S_mDerivv is not None and S_eDerivv is not None: + return S_mDerivv - 1j*omega(freq)*S_eDerivv + elif S_mDerivv is not None: + return S_mDerivv + elif S_eDerivv is not None: + return - 1j*omega(freq)*S_eDerivv + else: + return None + else: + S_mDerivv, S_eDerivv = S_mDeriv(v), S_eDeriv(v) + + if S_mDerivv is not None and S_eDerivv is not None: + RHSDeriv = C * (MeMuI * S_mDerivv) - 1j * omega(freq) * S_eDerivv + elif S_mDerivv is not None: + RHSDeriv = C * (MeMuI * S_mDerivv) + elif S_eDerivv is not None: + RHSDeriv = - 1j * omega(freq) * S_eDerivv + else: + return None + + if self._makeASymmetric: + return MfRho.T * RHSDeriv + return RHSDeriv + @@ -552,16 +590,18 @@ class ProblemFDEM_h(BaseFDEMProblem): MeMu = self.MeMu C = self.mesh.edgeCurl - sigi = self.sigmai - dsig_dm = self.curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 + # sigi = self.sigmai + # dsig_dm = self.curModel.transformDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) + MfRhoDeriv_m = self.MfRhoDeriv(C*u) + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) if adjoint: - return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) - return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) - + # return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) + return MfRhoDeriv_m.T * (C * v) + # return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) + return C.T * (MfRhoDeriv_m * v) def getRHS(self, freq): """ @@ -578,7 +618,30 @@ class ProblemFDEM_h(BaseFDEMProblem): return RHS - def getRHSDeriv(self, freq, v, adjoint=False): - raise NotImplementedError('getRHSDeriv not implemented yet') - return None + def getRHSDeriv(self, src, v, adjoint=False): + # raise NotImplementedError('getRHSDeriv not implemented yet') + # return None + _, S_e = self.getSourceTerm(src.freq) + C = self.mesh.edgeCurl + MfRho = self.MfRho + MfRhoDeriv = self.MfRhoDeriv(S_e) + + # if not adjoint: + MfRhoDeriv_m = self.MfRhoDeriv(S_e) + # elif adjoint: + # MfRhoDeriv_m = self.MfRhoDeriv + + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + + RHSDeriv = C.T * (MfRhoDeriv * v) + # = S_m + C.T * ( MfRho * S_e ) + + S_mDeriv = S_mDeriv(v) + S_eDeriv = S_eDeriv(v) + if S_mDeriv is not None: + RHSDeriv += S_mDeriv(v) + if S_eDeriv is not None: + RHSDeriv += C.T * (MfRho * S_e) + + return RHSDeriv diff --git a/simpegEM/FDEM/FieldsFDEM.py b/simpegEM/FDEM/FieldsFDEM.py index bf1175e8..2b599cdd 100644 --- a/simpegEM/FDEM/FieldsFDEM.py +++ b/simpegEM/FDEM/FieldsFDEM.py @@ -157,10 +157,10 @@ class FieldsFDEM_b(FieldsFDEM): return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v)) def _eSecondaryDeriv_m(self, src, v, adjoint=False): - bsol = self[[src],'bSolution'] + bSolution = self[[src],'bSolution'] _,S_e = src.eval(self.survey.prob) - w = self._edgeCurl.T * (self._MfMui * bsol) + w = self._edgeCurl.T * (self._MfMui * bSolution) if S_e is not None: w += -S_e @@ -207,6 +207,7 @@ class FieldsFDEM_j(FieldsFDEM): self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho self._curModel = self.survey.prob.curModel + self._MfRhoDeriv = self.survey.prob.MfRhoDeriv def _jPrimary(self, jSolution, srcList): jPrimary = np.zeros_like(jSolution) @@ -222,6 +223,13 @@ class FieldsFDEM_j(FieldsFDEM): def _j(self, jSolution, srcList): return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) + def _jDeriv_u(self, src, v, adjoint=False): + return None + + def _jDeriv_m(self, src, v, adjoint=False): + # assuming primary does not depend on the model + return None + def _hPrimary(self, jSolution, srcList): hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): @@ -242,27 +250,67 @@ class FieldsFDEM_j(FieldsFDEM): h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m return h + def _hSecondaryDeriv_u(self, src, v, adjoint=False): + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + if not adjoint: + return -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRho * v) ) + elif adjoint: + return -1./(1j*omega(src.freq)) * MfRho * (C * ( MeMuI .T * v)) + + def _hSecondaryDeriv_m(self, src, v, adjoint=False): + jSolution = self[[src],'jSolution'] + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + MfRhoDeriv = self._MfRhoDeriv + + if not adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) + elif adjoint: + hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) + + S_mDeriv,_ = src.evalDeriv(self.survey.prob, adjoint) + + if not adjoint: + S_mDeriv = S_mDeriv(v) + if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * MeMuI * S_mDeriv + elif adjoint: + S_mDeriv = S_mDeriv(MeMuI.T * v) + if S_mDeriv is not None: + hDeriv_m += 1./(1j*omega(src.freq)) * S_mDeriv + return h + + def _h(self, jSolution, srcList): return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) - def _hDeriv(self, jSolution, srcList, v, adjoint=False): - raise NotImplementedError('Fields Derivs Not Implemented Yet') - sig = self._curModel.transform - sigi = 1/sig - dsig_dm = self._curModel.transformDeriv - dsigi_dsig = -Utils.sdiag(sigi)**2 - dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) - sigi = self._MfRho + # raise NotImplementedError('Fields Derivs Not Implemented Yet') + # sig = self._curModel.transform + # sigi = 1/sig + # dsig_dm = self._curModel.transformDeriv + # dsigi_dsig = -Utils.sdiag(sigi)**2 + # dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) + # sigi = self._MfRho - S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) + # S_mDeriv,_ = src.getSourceDeriv(self.survey.prob, v, adjoint) - if not adjoint: - h_Deriv= -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) - else: - h_Deriv= -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) + # if not adjoint: + # h_Deriv= -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) + # else: + # h_Deriv= -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) - if S_mDeriv is not None: - return 1./(1j*omega(src.freq)) * S_mDeriv + h_Deriv + # if S_mDeriv is not None: + # return 1./(1j*omega(src.freq)) * S_mDeriv + h_Deriv + + def _hDeriv_u(self, src, v, adjoint=False): + return _hSecondaryDeriv_u(self, src, v, adjoint) + + def _hDeriv_m(self, src, v, adjoint=False): + # assuming the primary doesn't depend on the model + return _hSecondaryDeriv_u(self, src, v, adjoint) class FieldsFDEM_h(FieldsFDEM): @@ -298,6 +346,13 @@ class FieldsFDEM_h(FieldsFDEM): def _h(self, hSolution, srcList): return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) + def _hDeriv_u(self, src, v, adjoint=False): + return None + + def _hDeriv_m(self, src, v, adjoint=False): + # assuming primary does not depend on the model + return None + def _jPrimary(self, hSolution, srcList): jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]]) for i, src in enumerate(srcList): @@ -326,8 +381,8 @@ class FieldsFDEM_h(FieldsFDEM): return - S_eDeriv - # def calcFields(self, sol, freq, fieldType, adjoint=False): - # j = sol + # def calcFields(self, Solution, freq, fieldType, adjoint=False): + # j = Solution # if fieldType == 'j': # return j # elif fieldType == 'h': @@ -341,8 +396,8 @@ class FieldsFDEM_h(FieldsFDEM): # return h # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): - # j = sol + # def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False): + # j = Solution # if fieldType == 'j': # return None # elif fieldType == 'h': @@ -361,8 +416,8 @@ class FieldsFDEM_h(FieldsFDEM): # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # def calcFields(self, sol, freq, fieldType, adjoint=False): - # h = sol + # def calcFields(self, Solution, freq, fieldType, adjoint=False): + # h = Solution # if fieldType == 'j': # C = self.mesh.edgeCurl # if adjoint: @@ -372,5 +427,5 @@ class FieldsFDEM_h(FieldsFDEM): # return h # raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) - # def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + # def calcFieldsDeriv(self, Solution, freq, fieldType, v, adjoint=False): # return None \ No newline at end of file