From 64d89cfb6cef6fdac58e08ed209dd8240dae6976 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 14 Feb 2016 11:56:39 -0800 Subject: [PATCH 1/4] - simplified sensitivity calculation for FDEM : only look at fields once (pass it two vectors) - moved common elements from fields object into base fields object --- SimPEG/EM/FDEM/FDEM.py | 17 ++- SimPEG/EM/FDEM/FieldsFDEM.py | 276 ++++++++++++++++++----------------- 2 files changed, 150 insertions(+), 143 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 4b137b2c..39477138 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -89,18 +89,21 @@ class BaseFDEMProblem(BaseEMProblem): du_dm = Ainv * ( - dA_dm + dRHS_dm ) for rx in src.rxList: - df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) - df_dudu_dm = df_duFun(src, du_dm, adjoint=False) + df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) + df_dm = df_dmFun(src, du_dm, v, adjoint=False) - df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) - df_dm = df_dmFun(src, v, adjoint=False) + # df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) + # df_dudu_dm = df_duFun(src, du_dm, adjoint=False) + + # df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) + # df_dm = df_dmFun(src, v, adjoint=False) - Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) + Df_Dm = np.array(df_dm,dtype=complex) - P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m + # P = lambda v: - Jv[src, rx] = P(Df_Dm) + Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, Df_Dm) Ainv.clean() return Utils.mkvc(Jv) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index e171a5c5..6dc06cd7 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -32,6 +32,86 @@ class Fields(SimPEG.Problem.Fields): knownFields = {} dtype = complex + def _e(self, solution, srcList): + """ + Total electric field is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total electric field + """ + if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None: + raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0]) + + return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList) + + def _b(self, solution, srcList): + """ + Total magnetic flux density is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total magnetic flux density + """ + if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: + raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0]) + + return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList) + + def _h(self, solution, srcList): + """ + Total magnetic field is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total magnetic field + """ + if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None: + raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0]) + + return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList) + + def _j(self, solution, srcList): + """ + Total current density is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total current density + """ + if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: + raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0]) + + return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList) + + def _eDeriv(self, src, du_dm, v, adjoint = False): + if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: + raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + return self._eDeriv_u(src, du_dm, adjoint) + self._eDeriv_m(src, v, adjoint) + + def _bDeriv(self, src, du_dm, v, adjoint = False): + if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: + raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0]) + return self._bDeriv_u(src, du_dm, adjoint) + self._bDeriv_m(src, v, adjoint) + + def _hDeriv(self, src, du_dm, v, adjoint = False): + if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: + raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + return self._hDeriv_u(src, du_dm, adjoint) + self._hDeriv_m(src, v, adjoint) + + def _jDeriv(self, src, du_dm, v, adjoint = False): + if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: + raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + return self._jDeriv_u(src, du_dm, adjoint) + self._jDeriv_m(src, v, adjoint) + + class Fields_e(Fields): """ Fields object for Problem_e. @@ -85,22 +165,11 @@ class Fields_e(Fields): return eSolution - def _e(self, eSolution, srcList): - """ - Total electric field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total electric field - """ - - return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) def _eDeriv_u(self, src, v, adjoint = False): """ - Derivative of the total electric field with respect to the thing we - solved for + Partial derivative of the total electric field with respect to the thing we + solved for. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -113,7 +182,7 @@ class Fields_e(Fields): def _eDeriv_m(self, src, v, adjoint = False): """ - Derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. + Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -159,12 +228,12 @@ class Fields_e(Fields): b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b - def _bSecondaryDeriv_u(self, src, v, adjoint = False): + def _bSecondaryDeriv_u(self, src, du_dm, adjoint = False): """ Derivative of the secondary magnetic flux density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector @@ -172,8 +241,8 @@ class Fields_e(Fields): C = self._edgeCurl if adjoint: - return - 1./(1j*omega(src.freq)) * (C.T * v) - return - 1./(1j*omega(src.freq)) * (C * v) + return - 1./(1j*omega(src.freq)) * (C.T * du_dm) + return - 1./(1j*omega(src.freq)) * (C * du_dm) def _bSecondaryDeriv_m(self, src, v, adjoint = False): """ @@ -189,35 +258,23 @@ class Fields_e(Fields): S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) return 1./(1j * omega(src.freq)) * S_mDeriv - def _b(self, eSolution, srcList): - """ - Total magnetic flux density is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic flux density - """ - return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) - - def _bDeriv_u(self, src, v, adjoint=False): + def _bDeriv_u(self, src, du_dm, adjoint=False): """ - Derivative of the total magnetic flux density with respect to the thing we solved for + Partial derivative of the total magnetic flux density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ - # Primary does not depend on u - return self._bSecondaryDeriv_u(src, v, adjoint) + return self._bSecondaryDeriv_u(src, du_dm, adjoint) def _bDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic flux density with respect to the inversion model. + Partial derivative of the total magnetic flux density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -287,34 +344,23 @@ class Fields_b(Fields): return bSolution - def _b(self, bSolution, srcList): + def _bDeriv_u(self, src, du_dm, adjoint=False): """ - Total magnetic flux density is sum of primary and secondary - - :param numpy.ndarray bSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic flux density - """ - - return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList) - - def _bDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total magnetic flux density with respect to the thing we - solved for + Partial derivative of the total magnetic flux density with respect to the thing we + solved for. :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ - return Identity()*v + + return Identity()*du_dm def _bDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. + Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -358,7 +404,7 @@ class Fields_b(Fields): e[:,i] = e[:,i]+ -self._MeSigmaI * S_e return e - def _eSecondaryDeriv_u(self, src, v, adjoint=False): + def _eSecondaryDeriv_u(self, src, du_dm, adjoint=False): """ Derivative of the secondary electric field with respect to the thing we solved for @@ -370,9 +416,9 @@ class Fields_b(Fields): """ if not adjoint: - return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) ) + return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm) ) else: - return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v)) + return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm)) def _eSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -406,34 +452,22 @@ class Fields_b(Fields): return de_dm - def _e(self, bSolution, srcList): + def _eDeriv_u(self, src, du_dm, adjoint=False): """ - Total electric field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total electric field - """ - - return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList) - - def _eDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total electric field with respect to the thing we solved for + Partial derivative of the total electric field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the electric field with respect to the field we solved for with a vector """ - return self._eSecondaryDeriv_u(src, v, adjoint) + return self._eSecondaryDeriv_u(src, du_dm, adjoint) def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total electric field density with respect to the inversion model. + Partial derivative of the total electric field density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -515,10 +549,11 @@ class Fields_j(Fields): return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) - def _jDeriv_u(self, src, v, adjoint=False): + + def _jDeriv_u(self, src, du_dm, adjoint=False): """ - Derivative of the total current density with respect to the thing we - solved for + Partial derivative of the total current density with respect to the thing we + solved for. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -527,11 +562,12 @@ class Fields_j(Fields): :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return Identity()*v + return Identity()*du_dm + def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. + Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -575,21 +611,22 @@ class Fields_j(Fields): h[:,i] = h[:,i]+ 1./(1j*omega(src.freq)) * self._MeMuI * (S_m) return h - def _hSecondaryDeriv_u(self, src, v, adjoint=False): + + def _hSecondaryDeriv_u(self, src, du_dm, adjoint=False): """ Derivative of the secondary magnetic field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector """ if not adjoint: - return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) ) + return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm) ) elif adjoint: - return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v)) + return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm)) def _hSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -625,34 +662,22 @@ class Fields_j(Fields): return hDeriv_m - def _h(self, jSolution, srcList): + def _hDeriv_u(self, src, du_dm, adjoint=False): """ - Total magnetic field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic field - """ - - return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) - - def _hDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total magnetic field with respect to the thing we solved for + Partial derivative of the total magnetic field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return self._hSecondaryDeriv_u(src, v, adjoint) + return self._hSecondaryDeriv_u(src, du_dm, adjoint) def _hDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic field density with respect to the inversion model. + Partial derivative of the total magnetic field density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -720,35 +745,24 @@ class Fields_h(Fields): return hSolution - def _h(self, hSolution, srcList): - """ - Total magnetic field is sum of primary and secondary - - :param numpy.ndarray hSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic field - """ - return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) - - def _hDeriv_u(self, src, v, adjoint=False): + def _hDeriv_u(self, src, du_dm, adjoint=False): """ - Derivative of the total magnetic field with respect to the thing we - solved for + Partial derivative of the total magnetic field with respect to the thing we + solved for. :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return Identity()*v + return Identity()*du_dm def _hDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. + Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with @@ -792,21 +806,21 @@ class Fields_h(Fields): j[:,i] = j[:,i]+ -S_e return j - def _jSecondaryDeriv_u(self, src, v, adjoint=False): + def _jSecondaryDeriv_u(self, src, du_dm, adjoint=False): """ Derivative of the secondary current density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector """ if not adjoint: - return self._edgeCurl*v + return self._edgeCurl*du_dm elif adjoint: - return self._edgeCurl.T*v + return self._edgeCurl.T*du_dm def _jSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -822,33 +836,23 @@ class Fields_h(Fields): _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) return -S_eDeriv - def _j(self, hSolution, srcList): - """ - Total current density is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total current density - """ - return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) - - def _jDeriv_u(self, src, v, adjoint=False): + def _jDeriv_u(self, src, du_dm, adjoint=False): """ - Derivative of the total current density with respect to the thing we solved for + Partial derivative of the total current density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return self._jSecondaryDeriv_u(src,v,adjoint) + + return self._jSecondaryDeriv_u(src,du_dm,adjoint) def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total current density with respect to the inversion model. + Partial derivative of the total current density with respect to the inversion model. :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with From 8b152f38901f6dd4c2bc268b74cb1e25b05481cc Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 14 Feb 2016 12:28:25 -0800 Subject: [PATCH 2/4] adjoint of fields deriv now returns a tuple (so you don't call derivs wrt _u, _m independently --- SimPEG/EM/FDEM/FDEM.py | 30 ++++++---------------- SimPEG/EM/FDEM/FieldsFDEM.py | 49 ++++++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+), 22 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 39477138..d026e4b5 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -91,18 +91,7 @@ class BaseFDEMProblem(BaseEMProblem): for rx in src.rxList: df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) df_dm = df_dmFun(src, du_dm, v, adjoint=False) - - # df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) - # df_dudu_dm = df_duFun(src, du_dm, adjoint=False) - - # df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) - # df_dm = df_dmFun(src, v, adjoint=False) - - Df_Dm = np.array(df_dm,dtype=complex) - - # P = lambda v: - Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, Df_Dm) Ainv.clean() @@ -141,26 +130,23 @@ class BaseFDEMProblem(BaseEMProblem): for rx in src.rxList: PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m - df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None) - df_duT = df_duTFun(src, PTv, adjoint=True) - + df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) + ATinvdf_duT = ATinv * df_duT dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv_m(freq, src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) - dfT_dm = df_dmFun(src, PTv, adjoint=True) + df_dmT += du_dmT - du_dmT += dfT_dm - - # TODO: this should be taken care of by the reciever + # TODO: this should be taken care of by the reciever? real_or_imag = rx.projComp if real_or_imag is 'real': - Jtv += np.array(du_dmT,dtype=complex).real + Jtv += np.array(df_dmT,dtype=complex).real elif real_or_imag is 'imag': - Jtv += - np.array(du_dmT,dtype=complex).real + Jtv += - np.array(df_dmT,dtype=complex).real else: raise Exception('Must be real or imag') diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 6dc06cd7..da4b1749 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -89,26 +89,75 @@ class Fields(SimPEG.Problem.Fields): return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList) def _eDeriv(self, src, du_dm, v, adjoint = False): + """ + Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + if adjoint: + return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) return self._eDeriv_u(src, du_dm, adjoint) + self._eDeriv_m(src, v, adjoint) def _bDeriv(self, src, du_dm, v, adjoint = False): + """ + Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint) return self._bDeriv_u(src, du_dm, adjoint) + self._bDeriv_m(src, v, adjoint) def _hDeriv(self, src, du_dm, v, adjoint = False): + """ + Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0]) + if adjoint: + return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint) return self._hDeriv_u(src, du_dm, adjoint) + self._hDeriv_m(src, v, adjoint) def _jDeriv(self, src, du_dm, v, adjoint = False): + """ + Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + if adjoint: + return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) return self._jDeriv_u(src, du_dm, adjoint) + self._jDeriv_m(src, v, adjoint) From ea6ec4cb5556c845d90b4a745daa55a21fde743d Mon Sep 17 00:00:00 2001 From: Rowan Cockett Date: Mon, 15 Feb 2016 10:56:09 -0800 Subject: [PATCH 3/4] Integration in sources and documentation. --- SimPEG/EM/FDEM/SrcFDEM.py | 142 +++++++++++++++++++++++--------------- 1 file changed, 85 insertions(+), 57 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 1213cef3..6193b2f4 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -1,7 +1,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from scipy.constants import mu_0 from SimPEG.EM.Utils import * -from SimPEG.Utils import Zero +from SimPEG.Utils import Zero class BaseSrc(Survey.BaseSrc): """ @@ -14,7 +14,7 @@ class BaseSrc(Survey.BaseSrc): def eval(self, prob): """ - Evaluate the source terms. + Evaluate the source terms. - :math:`S_m` : magnetic source term - :math:`S_e` : electric source term @@ -36,12 +36,12 @@ class BaseSrc(Survey.BaseSrc): :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? :rtype: (numpy.ndarray, numpy.ndarray) - :return: tuple with magnetic source term and electric source term derivatives times a vector + :return: tuple with magnetic source term and electric source term derivatives times a vector """ - if v is not None: - return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint) + if v is not None: + return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint) else: - return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) + return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint) def bPrimary(self, prob): """ @@ -49,7 +49,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary magnetic flux density + :return: primary magnetic flux density """ return Zero() @@ -59,7 +59,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ return Zero() @@ -69,7 +69,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary electric field + :return: primary electric field """ return Zero() @@ -79,13 +79,13 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary current density + :return: primary current density """ return Zero() def S_m(self, prob): """ - Magnetic source term + Magnetic source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray @@ -95,7 +95,7 @@ class BaseSrc(Survey.BaseSrc): def S_e(self, prob): """ - Electric source term + Electric source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray @@ -136,15 +136,26 @@ class RawVec_e(BaseSrc): :param list rxList: receiver list :param float freq: frequency :param numpy.array S_e: electric source term + :param bool integrate: Integrate the source term (multiply by Me) [True] """ - def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_e = np.array(S_e,dtype=complex) + def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_e = np.array(S_e, dtype=complex) self.freq = float(freq) + self.integrate = integrate + BaseSrc.__init__(self, rxList) def S_e(self, prob): + """ + Electric source term + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: electric source term on mesh + """ + if prob._eqLocs is 'FE' and self.integrate is True: + return prob.Me * self._S_e return self._S_e @@ -155,10 +166,11 @@ class RawVec_m(BaseSrc): :param float freq: frequency :param rxList: receiver list :param numpy.array S_m: magnetic source term + :param bool integrate: Integrate the source term (multiply by Me) [True] """ - def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): - self._S_m = np.array(S_m,dtype=complex) + def __init__(self, rxList, freq, S_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): + self._S_m = np.array(S_m, dtype=complex) self.freq = float(freq) self.integrate = integrate @@ -166,12 +178,14 @@ class RawVec_m(BaseSrc): def S_m(self, prob): """ - Magnetic source term + Magnetic source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray :return: magnetic source term on mesh """ + if prob._eqLocs is 'EF' and self.integrate is True: + return prob.Me * self._S_m return self._S_m @@ -183,36 +197,51 @@ class RawVec(BaseSrc): :param float freq: frequency :param numpy.array S_m: magnetic source term :param numpy.array S_e: electric source term + :param bool integrate: Integrate the source term (multiply by Me) [True] """ - def __init__(self, rxList, freq, S_m, S_e, integrate = True): - self._S_m = np.array(S_m,dtype=complex) - self._S_e = np.array(S_e,dtype=complex) + def __init__(self, rxList, freq, S_m, S_e, integrate=True): + self._S_m = np.array(S_m, dtype=complex) + self._S_e = np.array(S_e, dtype=complex) self.freq = float(freq) self.integrate = integrate BaseSrc.__init__(self, rxList) def S_m(self, prob): + """ + Magnetic source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: magnetic source term on mesh + """ if prob._eqLocs is 'EF' and self.integrate is True: return prob.Me * self._S_m return self._S_m def S_e(self, prob): + """ + Electric source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: electric source term on mesh + """ if prob._eqLocs is 'FE' and self.integrate is True: return prob.Me * self._S_e return self._S_e class MagDipole(BaseSrc): - """ + """ Point magnetic dipole source calculated by taking the curl of a magnetic vector potential. By taking the discrete curl, we ensure that the magnetic - flux density is divergence free (no magnetic monopoles!). + flux density is divergence free (no magnetic monopoles!). This approach uses a primary-secondary in frequency. Here we show the derivation for E-B formulation noting that similar steps are followed for the H-J formulation. - .. math:: + .. math:: \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} @@ -225,17 +254,17 @@ class MagDipole(BaseSrc): and define a zero-frequency primary problem, noting that the source is generated by a divergence free electric current - .. math:: + .. math:: \mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\ {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}} - Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is + Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is - .. math:: + .. math:: \mathbf{e^P} = 0 \\\\ {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}} - Our secondary problem is then + Our secondary problem is then .. math:: \mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\ @@ -245,15 +274,15 @@ class MagDipole(BaseSrc): :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0): self.freq = float(freq) self.loc = loc self.orientation = orientation + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.moment = moment self.mu = mu self.integrate = False @@ -265,7 +294,7 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ eqLocs = prob._eqLocs @@ -303,7 +332,7 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return h_from_b(prob,b) @@ -314,7 +343,7 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b_p = self.bPrimary(prob) @@ -326,7 +355,7 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): @@ -340,7 +369,7 @@ class MagDipole(BaseSrc): C = prob.mesh.edgeCurl elif eqLocs is 'EF': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) @@ -353,21 +382,20 @@ class MagDipole_Bfield(BaseSrc): fields from a magnetic dipole. No discrete curl is taken, so the magnetic flux density may not be strictly divergence free. - This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. + This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. :param list rxList: receiver list :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - #TODO: neither does moment def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): self.freq = float(freq) self.loc = loc + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.orientation = orientation self.moment = moment self.mu = mu @@ -379,7 +407,7 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ eqLocs = prob._eqLocs @@ -418,7 +446,7 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return h_from_b(prob, b) @@ -429,7 +457,7 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return -1j*omega(self.freq)*b @@ -440,7 +468,7 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() @@ -453,7 +481,7 @@ class MagDipole_Bfield(BaseSrc): C = prob.mesh.edgeCurl elif eqLocs is 'EF': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) @@ -463,22 +491,22 @@ class CircularLoop(BaseSrc): """ Circular loop magnetic source calculated by taking the curl of a magnetic vector potential. By taking the discrete curl, we ensure that the magnetic - flux density is divergence free (no magnetic monopoles!). + flux density is divergence free (no magnetic monopoles!). - This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. + This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. :param list rxList: receiver list :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', radius=1., mu=mu_0): self.freq = float(freq) self.orientation = orientation + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.radius = radius self.mu = mu self.loc = loc @@ -491,7 +519,7 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ eqLocs = prob._eqLocs @@ -528,7 +556,7 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return 1./self.mu*b @@ -539,7 +567,7 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return -1j*omega(self.freq)*b @@ -550,7 +578,7 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() @@ -563,7 +591,7 @@ class CircularLoop(BaseSrc): C = prob.mesh.edgeCurl elif eqLocs is 'EF': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) From 649525fa884b5ade974d160e45c7d44b6e17794d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 19 Feb 2016 17:31:43 -0800 Subject: [PATCH 4/4] keep track of _v with notation --- SimPEG/EM/FDEM/FDEM.py | 46 +++++++++--------- SimPEG/EM/FDEM/FieldsFDEM.py | 94 ++++++++++++++++++------------------ 2 files changed, 70 insertions(+), 70 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index d026e4b5..d19f2ad1 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -84,15 +84,15 @@ class BaseFDEMProblem(BaseEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' u_src = u[src, ftype] - dA_dm = self.getADeriv_m(freq, u_src, v) - dRHS_dm = self.getRHSDeriv_m(freq, src, v) - du_dm = Ainv * ( - dA_dm + dRHS_dm ) + dA_dm_v = self.getADeriv(freq, u_src, v) + dRHS_dm_v = self.getRHSDeriv(freq, src, v) + du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) for rx in src.rxList: df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None) - df_dm = df_dmFun(src, du_dm, v, adjoint=False) - Df_Dm = np.array(df_dm,dtype=complex) - Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, Df_Dm) + df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) + df_dm_v = np.array(df_dm_v,dtype=complex) + Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, df_dm_v) Ainv.clean() return Utils.mkvc(Jv) @@ -135,8 +135,8 @@ class BaseFDEMProblem(BaseEMProblem): ATinvdf_duT = ATinv * df_duT - dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(freq, src, ATinvdf_duT, adjoint=True) + dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT df_dmT += du_dmT @@ -228,7 +228,7 @@ class Problem_e(BaseFDEMProblem): return C.T*MfMui*C + 1j*omega(freq)*MeSigma - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -269,7 +269,7 @@ class Problem_e(BaseFDEMProblem): return C.T * (MfMui * S_m) -1j * omega(freq) * S_e - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model @@ -343,7 +343,7 @@ class Problem_b(BaseFDEMProblem): return MfMui.T*A return A - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -400,7 +400,7 @@ class Problem_b(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model @@ -492,7 +492,7 @@ class Problem_j(BaseFDEMProblem): return A - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -513,16 +513,16 @@ class Problem_j(BaseFDEMProblem): MeMuI = self.MeMuI MfRho = self.MfRho C = self.mesh.edgeCurl - MfRhoDeriv_m = self.MfRhoDeriv(u) + MfRhoDeriv = self.MfRhoDeriv(u) if adjoint: if self._makeASymmetric is True: v = MfRho * v - return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v))) + return MfRhoDeriv.T * (C * (MeMuI.T * (C.T * v))) if self._makeASymmetric is True: - return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) ))) - return C * (MeMuI * (C.T * (MfRhoDeriv_m * v))) + return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv * v) ))) + return C * (MeMuI * (C.T * (MfRhoDeriv * v))) def getRHS(self, freq): @@ -549,7 +549,7 @@ class Problem_j(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model @@ -624,7 +624,7 @@ class Problem_h(BaseFDEMProblem): return C.T * (MfRho * C) + 1j*omega(freq)*MeMu - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -641,11 +641,11 @@ class Problem_h(BaseFDEMProblem): MeMu = self.MeMu C = self.mesh.edgeCurl - MfRhoDeriv_m = self.MfRhoDeriv(C*u) + MfRhoDeriv = self.MfRhoDeriv(C*u) if adjoint: - return MfRhoDeriv_m.T * (C * v) - return C.T * (MfRhoDeriv_m * v) + return MfRhoDeriv.T * (C * v) + return C.T * (MfRhoDeriv * v) def getRHS(self, freq): """ @@ -666,7 +666,7 @@ class Problem_h(BaseFDEMProblem): return S_m + C.T * ( MfRho * S_e ) - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index da4b1749..8d8c16f4 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -88,12 +88,12 @@ class Fields(SimPEG.Problem.Fields): return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList) - def _eDeriv(self, src, du_dm, v, adjoint = False): + def _eDeriv(self, src, du_dm_v, v, adjoint = False): """ Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint :param Src src: sorce - :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) :param numpy.ndarray v: vector to take sensitivity product with :param bool adjoint: adjoint? :rtype: numpy.ndarray @@ -104,14 +104,14 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) - return self._eDeriv_u(src, du_dm, adjoint) + self._eDeriv_m(src, v, adjoint) + return self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint) - def _bDeriv(self, src, du_dm, v, adjoint = False): + def _bDeriv(self, src, du_dm_v, v, adjoint = False): """ Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint :param Src src: sorce - :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) :param numpy.ndarray v: vector to take sensitivity product with :param bool adjoint: adjoint? :rtype: numpy.ndarray @@ -122,14 +122,14 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint) - return self._bDeriv_u(src, du_dm, adjoint) + self._bDeriv_m(src, v, adjoint) + return self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint) - def _hDeriv(self, src, du_dm, v, adjoint = False): + def _hDeriv(self, src, du_dm_v, v, adjoint = False): """ Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint :param Src src: sorce - :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) :param numpy.ndarray v: vector to take sensitivity product with :param bool adjoint: adjoint? :rtype: numpy.ndarray @@ -140,14 +140,14 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint) - return self._hDeriv_u(src, du_dm, adjoint) + self._hDeriv_m(src, v, adjoint) + return self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint) - def _jDeriv(self, src, du_dm, v, adjoint = False): + def _jDeriv(self, src, du_dm_v, v, adjoint = False): """ Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint :param Src src: sorce - :param numpy.ndarray du_dm: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) :param numpy.ndarray v: vector to take sensitivity product with :param bool adjoint: adjoint? :rtype: numpy.ndarray @@ -158,7 +158,7 @@ class Fields(SimPEG.Problem.Fields): if adjoint: return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) - return self._jDeriv_u(src, du_dm, adjoint) + self._jDeriv_m(src, v, adjoint) + return self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint) class Fields_e(Fields): @@ -277,12 +277,12 @@ class Fields_e(Fields): b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m return b - def _bSecondaryDeriv_u(self, src, du_dm, adjoint = False): + def _bSecondaryDeriv_u(self, src, du_dm_v, adjoint = False): """ Derivative of the secondary magnetic flux density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector @@ -290,8 +290,8 @@ class Fields_e(Fields): C = self._edgeCurl if adjoint: - return - 1./(1j*omega(src.freq)) * (C.T * du_dm) - return - 1./(1j*omega(src.freq)) * (C * du_dm) + return - 1./(1j*omega(src.freq)) * (C.T * du_dm_v) + return - 1./(1j*omega(src.freq)) * (C * du_dm_v) def _bSecondaryDeriv_m(self, src, v, adjoint = False): """ @@ -308,18 +308,18 @@ class Fields_e(Fields): return 1./(1j * omega(src.freq)) * S_mDeriv - def _bDeriv_u(self, src, du_dm, adjoint=False): + def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total magnetic flux density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ - return self._bSecondaryDeriv_u(src, du_dm, adjoint) + return self._bSecondaryDeriv_u(src, du_dm_v, adjoint) def _bDeriv_m(self, src, v, adjoint=False): """ @@ -393,19 +393,19 @@ class Fields_b(Fields): return bSolution - def _bDeriv_u(self, src, du_dm, adjoint=False): + def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total magnetic flux density with respect to the thing we solved for. :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ - return Identity()*du_dm + return Identity()*du_dm_v def _bDeriv_m(self, src, v, adjoint=False): """ @@ -453,7 +453,7 @@ class Fields_b(Fields): e[:,i] = e[:,i]+ -self._MeSigmaI * S_e return e - def _eSecondaryDeriv_u(self, src, du_dm, adjoint=False): + def _eSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the secondary electric field with respect to the thing we solved for @@ -465,9 +465,9 @@ class Fields_b(Fields): """ if not adjoint: - return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm) ) + return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) ) else: - return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm)) + return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v)) def _eSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -501,18 +501,18 @@ class Fields_b(Fields): return de_dm - def _eDeriv_u(self, src, du_dm, adjoint=False): + def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total electric field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the electric field with respect to the field we solved for with a vector """ - return self._eSecondaryDeriv_u(src, du_dm, adjoint) + return self._eSecondaryDeriv_u(src, du_dm_v, adjoint) def _eDeriv_m(self, src, v, adjoint=False): """ @@ -599,7 +599,7 @@ class Fields_j(Fields): return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) - def _jDeriv_u(self, src, du_dm, adjoint=False): + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total current density with respect to the thing we solved for. @@ -611,7 +611,7 @@ class Fields_j(Fields): :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return Identity()*du_dm + return Identity()*du_dm_v def _jDeriv_m(self, src, v, adjoint=False): @@ -661,21 +661,21 @@ class Fields_j(Fields): return h - def _hSecondaryDeriv_u(self, src, du_dm, adjoint=False): + def _hSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the secondary magnetic field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector """ if not adjoint: - return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm) ) + return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) elif adjoint: - return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm)) + return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) def _hSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -711,18 +711,18 @@ class Fields_j(Fields): return hDeriv_m - def _hDeriv_u(self, src, du_dm, adjoint=False): + def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total magnetic field with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return self._hSecondaryDeriv_u(src, du_dm, adjoint) + return self._hSecondaryDeriv_u(src, du_dm_v, adjoint) def _hDeriv_m(self, src, v, adjoint=False): """ @@ -795,19 +795,19 @@ class Fields_h(Fields): return hSolution - def _hDeriv_u(self, src, du_dm, adjoint=False): + def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total magnetic field with respect to the thing we solved for. :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return Identity()*du_dm + return Identity()*du_dm_v def _hDeriv_m(self, src, v, adjoint=False): """ @@ -855,21 +855,21 @@ class Fields_h(Fields): j[:,i] = j[:,i]+ -S_e return j - def _jSecondaryDeriv_u(self, src, du_dm, adjoint=False): + def _jSecondaryDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the secondary current density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector """ if not adjoint: - return self._edgeCurl*du_dm + return self._edgeCurl*du_dm_v elif adjoint: - return self._edgeCurl.T*du_dm + return self._edgeCurl.T*du_dm_v def _jSecondaryDeriv_m(self, src, v, adjoint=False): """ @@ -886,18 +886,18 @@ class Fields_h(Fields): return -S_eDeriv - def _jDeriv_u(self, src, du_dm, adjoint=False): + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ Partial derivative of the total current density with respect to the thing we solved for :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray du_dm: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return self._jSecondaryDeriv_u(src,du_dm,adjoint) + return self._jSecondaryDeriv_u(src,du_dm_v,adjoint) def _jDeriv_m(self, src, v, adjoint=False): """