From 81c13b12e3ea69bbe6bf6ae2e9acf96c9bf6e8ad Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 7 Feb 2016 13:04:13 -0800 Subject: [PATCH] cleanup of docs, docs for survey --- SimPEG/EM/FDEM/FDEM.py | 50 +++++++++++++++++++----------------- SimPEG/EM/FDEM/FieldsFDEM.py | 13 ++++------ SimPEG/EM/FDEM/SurveyFDEM.py | 44 +++++++++++++++++++++++++++++-- docs/em/api_FDEM.rst | 18 +++++++++++++ 4 files changed, 91 insertions(+), 34 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 843ab245..09bbd143 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -15,18 +15,20 @@ class BaseFDEMProblem(BaseEMProblem): .. 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{M^e} \mathbf{s_e}} + {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} if using the E-B formulation (:code:`Problem_e` - or :code:`Problem_b`) or the magnetic field + or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + + If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: - \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\ + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). + if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ @@ -204,7 +206,7 @@ class Problem_e(BaseFDEMProblem): .. math :: - \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e} + \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e} which we solve for :math:`\mathbf{e}`. @@ -223,7 +225,7 @@ class Problem_e(BaseFDEMProblem): System matrix .. math :: - \mathbf{A} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} + \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} :param float freq: Frequency :rtype: scipy.sparse.csr_matrix @@ -265,7 +267,7 @@ class Problem_e(BaseFDEMProblem): Right hand side for the system .. math :: - \mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e} + \mathbf{RHS} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e} :param float freq: Frequency :rtype: numpy.ndarray @@ -294,7 +296,7 @@ class Problem_e(BaseFDEMProblem): C = self.mesh.edgeCurl MfMui = self.MfMui - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: dRHS = MfMui * (C * v) @@ -310,13 +312,13 @@ class Problem_b(BaseFDEMProblem): .. math :: - \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right) + \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right) and solve for :math:`\mathbf{b}` using: .. math :: - \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e} + \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e} .. note :: The inverse problem will not work with full anisotropy @@ -336,7 +338,7 @@ class Problem_b(BaseFDEMProblem): System matrix .. math :: - \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega + \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega :param float freq: Frequency :rtype: scipy.sparse.csr_matrix @@ -431,7 +433,7 @@ class Problem_b(BaseFDEMProblem): v = self.MfMui * v MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) @@ -458,13 +460,13 @@ class Problem_j(BaseFDEMProblem): .. math :: - \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right) + \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right) and solve for \\\(\\\mathbf{j}\\\) using .. math :: - \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e} + \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e} .. note:: This implementation does not yet work with full anisotropy!! @@ -484,7 +486,7 @@ class Problem_j(BaseFDEMProblem): System matrix .. math :: - \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{\\mu^{-1}}} \\mathbf{C}^T \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega + \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{\\mu^{-1}}} \\mathbf{C}^{\\top} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega :param float freq: Frequency :rtype: scipy.sparse.csr_matrix @@ -511,7 +513,7 @@ class Problem_j(BaseFDEMProblem): .. math :: - \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^T} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}} + \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \mathbf{M^e_{mu^{-1}}} \mathbf{C^{\\top}} \\frac{d \mathbf{M^f_{\sigma^{-1}}}\mathbf{v} }{d \mathbf{m}} :param float freq: frequency :param numpy.ndarray u: solution vector (nF,) @@ -543,6 +545,7 @@ class Problem_j(BaseFDEMProblem): .. math :: \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e} + :param float freq: Frequency :rtype: numpy.ndarray (nE, nSrc) :return: RHS @@ -573,7 +576,7 @@ class Problem_j(BaseFDEMProblem): C = self.mesh.edgeCurl MeMuI = self.MeMuI - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: if self._makeASymmetric: @@ -604,7 +607,7 @@ class Problem_h(BaseFDEMProblem): .. math :: - \\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} + \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e} :param SimPEG.Mesh mesh: mesh """ @@ -620,9 +623,8 @@ class Problem_h(BaseFDEMProblem): """ System matrix - .. math :: - - \mathbf{A} = \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e} + .. math:: + \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e} :param float freq: Frequency :rtype: scipy.sparse.csr_matrix @@ -640,7 +642,7 @@ class Problem_h(BaseFDEMProblem): Product of the derivative of our system matrix with respect to the model and a vector .. math:: - \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top} \\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}} + \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C}^{\\top}\\frac{d \mathbf{M^f_{\\rho}}\mathbf{v} }{d\mathbf{m}} :param float freq: frequency :param numpy.ndarray u: solution vector (nE,) @@ -664,7 +666,7 @@ class Problem_h(BaseFDEMProblem): .. math :: - \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} + \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e} :param float freq: Frequency :rtype: numpy.ndarray @@ -701,7 +703,7 @@ class Problem_h(BaseFDEMProblem): elif adjoint: RHSDeriv = MfRhoDeriv.T * (C * v) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v)) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 9f54855c..e171a5c5 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -186,8 +186,7 @@ class Fields_e(Fields): :return: product of the secondary magnetic flux density derivative with respect to the inversion model with a vector """ - S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) - S_mDeriv = S_mDeriv(v) + S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) return 1./(1j * omega(src.freq)) * S_mDeriv def _b(self, eSolution, srcList): @@ -401,10 +400,9 @@ class Fields_b(Fields): elif adjoint: de_dm = self._MeSigmaIDeriv(w).T * v - _, S_eDeriv = src.evalDeriv(self.prob, adjoint) - Se_Deriv = S_eDeriv(v) + _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - de_dm = de_dm - self._MeSigmaI * Se_Deriv + de_dm = de_dm - self._MeSigmaI * S_eDeriv return de_dm @@ -616,7 +614,7 @@ class Fields_j(Fields): elif adjoint: hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint) + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) if not adjoint: S_mDeriv = S_mDeriv(v) @@ -821,8 +819,7 @@ class Fields_h(Fields): :return: product of the secondary current density derivative with respect to the inversion model with a vector """ - _,S_eDeriv = src.evalDeriv(self.prob, adjoint) - S_eDeriv = S_eDeriv(v) + _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) return -S_eDeriv def _j(self, hSolution, srcList): diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index f60cbfdf..dd49b3ed 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -10,6 +10,12 @@ import SrcFDEM as Src #################################################### class Rx(SimPEG.Survey.BaseRx): + """ + Frequency domain receivers + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string rxType: reciever type from knownRxTypes + """ knownRxTypes = { 'exr':['e', 'Ex', 'real'], @@ -61,6 +67,15 @@ class Rx(SimPEG.Survey.BaseRx): return self.knownRxTypes[self.rxType][2] def projectFields(self, src, mesh, u): + """ + Project fields to recievers to get data. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields u: fields object + :rtype: numpy.ndarray + :return: fields projected to recievers + """ P = self.getP(mesh) u_part_complex = u[src, self.projField] # get the real or imag component @@ -69,6 +84,16 @@ class Rx(SimPEG.Survey.BaseRx): return P*u_part def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False): + """ + Derivative of projected fields with respect to the inversion model times a vector. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields u: fields object + :param numpy.ndarray v: vector to multiply + :rtype: numpy.ndarray + :return: fields projected to recievers + """ P = self.getP(mesh) if not adjoint: @@ -95,10 +120,13 @@ class Rx(SimPEG.Survey.BaseRx): class Survey(SimPEG.Survey.BaseSurvey): """ - docstring for SurveyFDEM + Frequency domain electromagnetic survey + + :param list srcList: list of FDEM sources used in the survey """ srcPair = Src.BaseSrc + rxPaair = Rx def __init__(self, srcList, **kwargs): # Sort these by frequency @@ -126,6 +154,7 @@ class Survey(SimPEG.Survey.BaseSurvey): @property def nSrcByFreq(self): + """Number of sources at each frequency""" if getattr(self, '_nSrcByFreq', None) is None: self._nSrcByFreq = {} for freq in self.freqs: @@ -133,11 +162,22 @@ class Survey(SimPEG.Survey.BaseSurvey): return self._nSrcByFreq def getSrcByFreq(self, freq): - """Returns the sources associated with a specific frequency.""" + """ + Returns the sources associated with a specific frequency. + :param float freq: frequency for which we look up sources + :rtype: dictionary + :returns: sources at the sepcified frequency + """ assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] def projectFields(self, u): + """ + Project fields to receiver locations + :param Fields u: fields object + :rtype: numpy.ndarray + :returns: data + """ data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: diff --git a/docs/em/api_FDEM.rst b/docs/em/api_FDEM.rst index e60c1bbf..f9de0979 100644 --- a/docs/em/api_FDEM.rst +++ b/docs/em/api_FDEM.rst @@ -144,6 +144,10 @@ H-J Formulation API === + +FDEM Problem +------------ + .. automodule:: SimPEG.EM.FDEM.FDEM :show-inheritance: :members: @@ -157,3 +161,17 @@ FDEM Survey :show-inheritance: :members: :undoc-members: + +.. automodule:: SimPEG.EM.FDEM.SrcFDEM + :show-inheritance: + :members: + :undoc-members: + +FDEM Fields +----------- + +.. automodule:: SimPEG.EM.FDEM.FieldsFDEM + :show-inheritance: + :members: + :undoc-members: +