From 3d11431f2fc88da5553179072f57a84f6e0cde69 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 11:47:57 -0700 Subject: [PATCH 1/7] use f where we are talking about fields --- SimPEG/EM/FDEM/FDEM.py | 112 ++++++++++++++++++++--------------------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index cf3dee7f..834c0187 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -18,9 +18,9 @@ class BaseFDEMProblem(BaseEMProblem): {\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`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + 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 + If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: @@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem): \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`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. + 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}\\\) """ @@ -39,73 +39,73 @@ class BaseFDEMProblem(BaseEMProblem): def fields(self, m=None): """ Solve the forward problem for the fields. - + :param numpy.array m: inversion model (nP,) :rtype numpy.array: - :return F: forward solution + :return f: forward solution """ self.curModel = m - F = self.fieldsPair(self.mesh, self.survey) + f = self.fieldsPair(self.mesh, self.survey) for freq in self.survey.freqs: A = self.getA(freq) rhs = self.getRHS(freq) Ainv = self.Solver(A, **self.solverOpts) - sol = Ainv * rhs + u = Ainv * rhs Srcs = self.survey.getSrcByFreq(freq) - F[Srcs, self._solutionType] = sol + f[Srcs, self._solutionType] = u Ainv.clean() - return F + return f - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ Sensitivity times a vector. :param numpy.array m: inversion model (nP,) :param numpy.array v: vector which we take sensitivity product with (nP,) - :param SimPEG.EM.FDEM.Fields u: fields object + :param SimPEG.EM.FDEM.Fields u: fields object :rtype numpy.array: - :return: Jv (ndata,) + :return: Jv (ndata,) """ - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) self.curModel = m Jv = self.dataPair(self.survey) for freq in self.survey.freqs: - A = self.getA(freq) - Ainv = self.Solver(A, **self.solverOpts) + A = self.getA(freq) + Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve) for src in self.survey.getSrcByFreq(freq): - u_src = u[src, self._solutionType] - dA_dm_v = self.getADeriv(freq, u_src, v) - dRHS_dm_v = self.getRHSDeriv(freq, src, v) + f_src = f[src, self._solutionType] + dA_dm_v = self.getADeriv(freq, f_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_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) - Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v) + Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) Ainv.clean() return Utils.mkvc(Jv) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): """ Sensitivity transpose times a vector :param numpy.array m: inversion model (nP,) :param numpy.array v: vector which we take adjoint product with (nP,) - :param SimPEG.EM.FDEM.Fields u: fields object + :param SimPEG.EM.FDEM.Fields u: fields object :rtype numpy.array: - :return: Jv (ndata,) + :return: Jv (ndata,) """ - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) self.curModel = m @@ -120,17 +120,17 @@ class BaseFDEMProblem(BaseEMProblem): ATinv = self.Solver(AT, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): - u_src = u[src, self._solutionType] + f_src = f[src, self._solutionType] for rx in src.rxList: - PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m + PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None) + df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) ATinvdf_duT = ATinv * df_duT - dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True) + dA_dmT = self.getADeriv(freq, f_src, ATinvdf_duT, adjoint=True) dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT @@ -144,7 +144,7 @@ class BaseFDEMProblem(BaseEMProblem): Jtv += - np.array(df_dmT, dtype=complex).real else: raise Exception('Must be real or imag') - + ATinv.clean() return Utils.mkvc(Jtv) @@ -154,7 +154,7 @@ class BaseFDEMProblem(BaseEMProblem): Evaluates the sources for a given frequency and puts them in matrix form :param float freq: Frequency - :rtype: (numpy.ndarray, numpy.ndarray) + :rtype: (numpy.ndarray, numpy.ndarray) :return: S_m, S_e (nE or nF, nSrc) """ Srcs = self.survey.getSrcByFreq(freq) @@ -207,7 +207,7 @@ class Problem_e(BaseFDEMProblem): def getA(self, freq): """ System matrix - + .. math :: \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} @@ -230,12 +230,12 @@ class Problem_e(BaseFDEMProblem): .. math :: \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = i \omega \\frac{d \mathbf{M^e_{\sigma}}\mathbf{v} }{d\mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nE,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nE,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ dsig_dm = self.curModel.sigmaDeriv @@ -248,13 +248,13 @@ class Problem_e(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \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 + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -266,7 +266,7 @@ class Problem_e(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -346,12 +346,12 @@ class Problem_b(BaseFDEMProblem): .. math :: \\frac{\mathbf{A}(\mathbf{m}) \mathbf{v}}{d \mathbf{m}} = \mathbf{C} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nF,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nF,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MfMui = self.MfMui @@ -373,13 +373,13 @@ class Problem_b(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e} :param float freq: Frequency - :rtype: numpy.ndarray + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -497,12 +497,12 @@ class Problem_j(BaseFDEMProblem): \\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,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nF,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MeMuI = self.MeMuI @@ -522,7 +522,7 @@ class Problem_j(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: @@ -546,7 +546,7 @@ class Problem_j(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -626,12 +626,12 @@ class Problem_h(BaseFDEMProblem): .. 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}} - :param float freq: frequency - :param numpy.ndarray u: solution vector (nE,) + :param float freq: frequency + :param numpy.ndarray u: solution vector (nE,) :param numpy.ndarray v: vector to take prodct with (nP,) or (nD,) for adjoint :param bool adjoint: adjoint? :rtype: numpy.ndarray - :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) + :return: derivative of the system matrix times a vector (nP,) or adjoint (nD,) """ MeMu = self.MeMu @@ -644,14 +644,14 @@ class Problem_h(BaseFDEMProblem): def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: \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 + :rtype: numpy.ndarray :return: RHS (nE, nSrc) """ @@ -663,7 +663,7 @@ class Problem_h(BaseFDEMProblem): def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source From d5f73d0fd35ac9b671cb70ffad4f4418e5955c65 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 11:56:58 -0700 Subject: [PATCH 2/7] f_src is actually u_src --- SimPEG/EM/FDEM/FDEM.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 834c0187..5c2683b2 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -81,8 +81,8 @@ class BaseFDEMProblem(BaseEMProblem): Ainv = self.Solver(A, **self.solverOpts) # create the concept of Ainv (actually a solve) for src in self.survey.getSrcByFreq(freq): - f_src = f[src, self._solutionType] - dA_dm_v = self.getADeriv(freq, f_src, v) + u_src = f[src, self._solutionType] + 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 ) @@ -120,7 +120,7 @@ class BaseFDEMProblem(BaseEMProblem): ATinv = self.Solver(AT, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): - f_src = f[src, self._solutionType] + u_src = f[src, self._solutionType] for rx in src.rxList: PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m @@ -130,7 +130,7 @@ class BaseFDEMProblem(BaseEMProblem): ATinvdf_duT = ATinv * df_duT - dA_dmT = self.getADeriv(freq, f_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 From 055061ac3b34c04e03d3037006729a22e0d5a295 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 12:15:52 -0700 Subject: [PATCH 3/7] abstracted FDEM survey to BaseEMSurvey (with methods eval and eval deriv) as this should be common to FDEM and TDEM problems (only implemented on FDEM problem, TDEM inheritance will be taken care of on the TDEM refactor branch) --- SimPEG/EM/Base.py | 51 ++++++++++++++++++++++++++---------- SimPEG/EM/FDEM/SurveyFDEM.py | 27 +++++-------------- 2 files changed, 43 insertions(+), 35 deletions(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 32018f7e..12d1ae43 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -2,14 +2,14 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve from scipy.constants import mu_0 class EMPropMap(Maps.PropMap): - """ + """ Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) """ sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) - rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) + rho = Maps.Property("Electrical Resistivity", propertyLink=('sigma', Maps.ReciprocalMap)) mui = Maps.Property("Inverse Magnetic Permeability", defaultVal = 1./mu_0, propertyLink=('mu', Maps.ReciprocalMap)) @@ -21,7 +21,7 @@ class BaseEMProblem(Problem.BaseProblem): surveyPair = Survey.BaseSurvey dataPair = Survey.Data - + PropMap = EMPropMap Solver = SimpegSolver @@ -51,7 +51,7 @@ class BaseEMProblem(Problem.BaseProblem): if self.mapping.muMap is not None or self.mapping.muiMap is not None: toDelete += ['_MeMu', '_MeMuI','_MfMui','_MfMuiI'] return toDelete - + @property def Me(self): """ @@ -71,7 +71,7 @@ class BaseEMProblem(Problem.BaseProblem): return self._Mf - # ----- Magnetic Permeability ----- # + # ----- Magnetic Permeability ----- # @property def MfMui(self): """ @@ -109,7 +109,7 @@ class BaseEMProblem(Problem.BaseProblem): return self._MeMuI - # ----- Electrical Conductivity ----- # + # ----- Electrical Conductivity ----- # #TODO: hardcoded to sigma as the model @property def MeSigma(self): @@ -120,18 +120,18 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma - # TODO: This should take a vector + # TODO: This should take a vector def MeSigmaDeriv(self, u): """ Derivative of MeSigma with respect to the model - """ + """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv - + @property def MeSigmaI(self): """ - Inverse of the edge inner product matrix for \\(\\sigma\\). + Inverse of the edge inner product matrix for \\(\\sigma\\). """ if getattr(self, '_MeSigmaI', None) is None: self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) @@ -140,8 +140,8 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MeSigmaIDeriv(self, u): """ - Derivative of :code:`MeSigma` with respect to the model - """ + Derivative of :code:`MeSigma` with respect to the model + """ # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG dMeSigmaI_dI = -self.MeSigmaI**2 @@ -163,7 +163,7 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MfRhoDeriv(self,u): """ - Derivative of :code:`MfRho` with respect to the model. + Derivative of :code:`MfRho` with respect to the model. """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @@ -181,6 +181,29 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MfRhoIDeriv(self,u): """ - Derivative of :code:`MfRhoI` with respect to the model. + Derivative of :code:`MfRhoI` with respect to the model. """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv + +class BaseEMSurvey(Survey.BaseSurvey): + + def __init__(self, srcList, **kwargs): + # Sort these by frequency + self.srcList = srcList + Survey.BaseSurvey.__init__(self, **kwargs) + + def eval(self, u): + """ + Project fields to receiver locations + :param Fields u: fields object + :rtype: numpy.ndarray + :return: data + """ + data = SimPEG.Survey.Data(self) + for src in self.srcList: + for rx in src.rxList: + data[src, rx] = rx.eval(src, self.mesh, u) + return data + + def evalDeriv(self, u): + raise Exception('Use Receivers to project fields deriv.') diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index ce803ed1..163f6914 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,5 +1,6 @@ import SimPEG from SimPEG.EM.Utils import * +from SimPEG.EM.Base import BaseEMSurvey from scipy.constants import mu_0 from SimPEG.Utils import Zero, Identity import SrcFDEM as Src @@ -84,7 +85,7 @@ class Rx(SimPEG.Survey.BaseRx): # get the real or imag component real_or_imag = self.projComp u_part = getattr(u_part_complex, real_or_imag) - + return P*u_part def evalDeriv(self, src, mesh, u, v, adjoint=False): @@ -123,7 +124,7 @@ class Rx(SimPEG.Survey.BaseRx): # Survey #################################################### -class Survey(SimPEG.Survey.BaseSurvey): +class Survey(BaseEMSurvey): """ Frequency domain electromagnetic survey @@ -131,12 +132,12 @@ class Survey(SimPEG.Survey.BaseSurvey): """ srcPair = Src.BaseSrc - rxPaair = Rx + rxPair = Rx def __init__(self, srcList, **kwargs): # Sort these by frequency self.srcList = srcList - SimPEG.Survey.BaseSurvey.__init__(self, **kwargs) + BaseEMSurvey.__init__(self, srcList, **kwargs) _freqDict = {} for src in srcList: @@ -171,24 +172,8 @@ class Survey(SimPEG.Survey.BaseSurvey): Returns the sources associated with a specific frequency. :param float freq: frequency for which we look up sources :rtype: dictionary - :return: sources at the sepcified frequency + :return: sources at the sepcified frequency """ assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] - def eval(self, u): - """ - Project fields to receiver locations - :param Fields u: fields object - :rtype: numpy.ndarray - :return: data - """ - data = SimPEG.Survey.Data(self) - for src in self.srcList: - for rx in src.rxList: - data[src, rx] = rx.eval(src, self.mesh, u) - return data - - def evalDeriv(self, u): - raise Exception('Use Receivers to project fields deriv.') - From cc9d2e5ac79b142ac3f412fa530193d77cfdf116 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 12:18:07 -0700 Subject: [PATCH 4/7] we don't support m is none --- SimPEG/EM/FDEM/FDEM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 5c2683b2..3c2183b1 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -36,7 +36,7 @@ class BaseFDEMProblem(BaseEMProblem): surveyPair = SurveyFDEM fieldsPair = Fields - def fields(self, m=None): + def fields(self, m): """ Solve the forward problem for the fields. From c51afa4aad66487fd2c9f104385a93f9fd8229b5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 13:04:02 -0700 Subject: [PATCH 5/7] s_m, s_e are vectors (so they should not be capitalized) --- SimPEG/EM/FDEM/FDEM.py | 62 ++++----- SimPEG/EM/FDEM/FieldsFDEM.py | 254 +++++++++++++++++------------------ SimPEG/EM/FDEM/SrcFDEM.py | 96 ++++++------- 3 files changed, 206 insertions(+), 206 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 3c2183b1..caca7602 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -155,22 +155,22 @@ class BaseFDEMProblem(BaseEMProblem): :param float freq: Frequency :rtype: (numpy.ndarray, numpy.ndarray) - :return: S_m, S_e (nE or nF, nSrc) + :return: s_m, s_e (nE or nF, nSrc) """ Srcs = self.survey.getSrcByFreq(freq) if self._formulation is 'EB': - S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) - S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) + s_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) + s_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) elif self._formulation is 'HJ': - S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) - S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) + s_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) + s_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) for i, src in enumerate(Srcs): smi, sei = src.eval(self) - S_m[:,i] = S_m[:,i] + smi - S_e[:,i] = S_e[:,i] + sei + s_m[:,i] = s_m[:,i] + smi + s_e[:,i] = s_e[:,i] + sei - return S_m, S_e + return s_m, s_e ########################################################################################## @@ -258,11 +258,11 @@ class Problem_e(BaseFDEMProblem): :return: RHS (nE, nSrc) """ - S_m, S_e = self.getSourceTerm(freq) + s_m, s_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MfMui = self.MfMui - return C.T * (MfMui * S_m) -1j * omega(freq) * S_e + return C.T * (MfMui * s_m) -1j * omega(freq) * s_e def getRHSDeriv(self, freq, src, v, adjoint=False): """ @@ -278,14 +278,14 @@ class Problem_e(BaseFDEMProblem): C = self.mesh.edgeCurl MfMui = self.MfMui - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) + s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: dRHS = MfMui * (C * v) - return S_mDeriv(dRHS) - 1j * omega(freq) * S_eDeriv(v) + return s_mDeriv(dRHS) - 1j * omega(freq) * s_eDeriv(v) else: - return C.T * (MfMui * S_mDeriv(v)) -1j * omega(freq) * S_eDeriv(v) + return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v) class Problem_b(BaseFDEMProblem): @@ -383,11 +383,11 @@ class Problem_b(BaseFDEMProblem): :return: RHS (nE, nSrc) """ - S_m, S_e = self.getSourceTerm(freq) + s_m, s_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MeSigmaI = self.MeSigmaI - RHS = S_m + C * ( MeSigmaI * S_e ) + RHS = s_m + C * ( MeSigmaI * s_e ) if self._makeASymmetric is True: MfMui = self.MfMui @@ -408,21 +408,21 @@ class Problem_b(BaseFDEMProblem): """ C = self.mesh.edgeCurl - S_m, S_e = src.eval(self) + s_m, s_e = src.eval(self) MfMui = self.MfMui if self._makeASymmetric and adjoint: v = self.MfMui * v - MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) + MeSigmaIDeriv = self.MeSigmaIDeriv(s_e) + s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) - SrcDeriv = S_mDeriv(v) + C * (self.MeSigmaI * S_eDeriv(v)) + SrcDeriv = s_mDeriv(v) + C * (self.MeSigmaI * s_eDeriv(v)) elif adjoint: RHSderiv = MeSigmaIDeriv.T * (C.T * v) - SrcDeriv = S_mDeriv(v) + self.MeSigmaI.T * (C.T * S_eDeriv(v)) + SrcDeriv = s_mDeriv(v) + self.MeSigmaI.T * (C.T * s_eDeriv(v)) if self._makeASymmetric is True and not adjoint: return MfMui.T * (SrcDeriv + RHSderiv) @@ -533,11 +533,11 @@ class Problem_j(BaseFDEMProblem): :return: RHS """ - S_m, S_e = self.getSourceTerm(freq) + s_m, s_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MeMuI = self.MeMuI - RHS = C * (MeMuI * S_m) - 1j * omega(freq) * S_e + RHS = C * (MeMuI * s_m) - 1j * omega(freq) * s_e if self._makeASymmetric is True: MfRho = self.MfRho return MfRho.T*RHS @@ -558,16 +558,16 @@ class Problem_j(BaseFDEMProblem): C = self.mesh.edgeCurl MeMuI = self.MeMuI - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) + s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: if self._makeASymmetric: MfRho = self.MfRho v = MfRho*v - return S_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * S_eDeriv(v) + return s_mDeriv(MeMuI.T * (C.T * v)) - 1j * omega(freq) * s_eDeriv(v) else: - RHSDeriv = C * (MeMuI * S_mDeriv(v)) - 1j * omega(freq) * S_eDeriv(v) + RHSDeriv = C * (MeMuI * s_mDeriv(v)) - 1j * omega(freq) * s_eDeriv(v) if self._makeASymmetric: MfRho = self.MfRho @@ -655,11 +655,11 @@ class Problem_h(BaseFDEMProblem): :return: RHS (nE, nSrc) """ - S_m, S_e = self.getSourceTerm(freq) + s_m, s_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MfRho = self.MfRho - return S_m + C.T * ( MfRho * S_e ) + return s_m + C.T * ( MfRho * s_e ) def getRHSDeriv(self, freq, src, v, adjoint=False): """ @@ -673,17 +673,17 @@ class Problem_h(BaseFDEMProblem): :return: product of rhs deriv with a vector """ - _, S_e = src.eval(self) + _, s_e = src.eval(self) C = self.mesh.edgeCurl MfRho = self.MfRho - MfRhoDeriv = self.MfRhoDeriv(S_e) + MfRhoDeriv = self.MfRhoDeriv(s_e) if not adjoint: RHSDeriv = C.T * (MfRhoDeriv * v) elif adjoint: RHSDeriv = MfRhoDeriv.T * (C * v) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) + s_mDeriv, s_eDeriv = src.evalDeriv(self, adjoint=adjoint) - return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v)) + 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 0f0a4963..e2193973 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -8,7 +8,7 @@ from SimPEG.Utils import Zero, Identity, sdiag class Fields(SimPEG.Problem.Fields): """ - + Fancy Field Storage for a FDEM survey. Only one field type is stored for each problem, the rest are computed. The fields obejct acts like an array and is indexed by @@ -34,56 +34,56 @@ class Fields(SimPEG.Problem.Fields): def _e(self, solution, srcList): """ - Total electric field is sum of primary and secondary - + 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: + 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 - + 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 + :return: total magnetic flux density """ - if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: + 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 - + 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: + 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 - + 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 + :return: total current density """ - if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: + 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) @@ -99,7 +99,7 @@ class Fields(SimPEG.Problem.Fields): :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: + 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: @@ -117,12 +117,12 @@ class Fields(SimPEG.Problem.Fields): :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: + 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 np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex) + return np.array(self._bDeriv_u(src, du_dm_v, adjoint) + self._bDeriv_m(src, v, adjoint), dtype = complex) def _hDeriv(self, src, du_dm_v, v, adjoint = False): """ @@ -135,10 +135,10 @@ class Fields(SimPEG.Problem.Fields): :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: + 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: + if adjoint: return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint) return np.array(self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint), dtype = complex) @@ -153,7 +153,7 @@ class Fields(SimPEG.Problem.Fields): :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: + 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: @@ -162,10 +162,10 @@ class Fields(SimPEG.Problem.Fields): class Fields_e(Fields): """ - Fields object for Problem_e. + Fields object for Problem_e. :param Mesh mesh: mesh - :param Survey survey: survey + :param Survey survey: survey """ knownFields = {'eSolution':'E'} @@ -233,9 +233,9 @@ class Fields_e(Fields): def _eDeriv_u(self, src, v, adjoint = False): """ - Partial derivative of the total electric field with respect to the thing we + 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 bool adjoint: adjoint? @@ -247,8 +247,8 @@ class Fields_e(Fields): def _eDeriv_m(self, src, v, adjoint = False): """ - 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. - + 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 :param bool adjoint: adjoint? @@ -289,14 +289,14 @@ class Fields_e(Fields): b = (C * eSolution) for i, src in enumerate(srcList): b[:,i] *= - 1./(1j*omega(src.freq)) - S_m, _ = src.eval(self.prob) - b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * S_m + s_m, _ = src.eval(self.prob) + b[:,i] = b[:,i]+ 1./(1j*omega(src.freq)) * s_m return b def _bDeriv_u(self, src, du_dm_v, adjoint = False): """ Derivative of the magnetic flux density with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -312,8 +312,8 @@ class Fields_e(Fields): def _bDeriv_m(self, src, v, adjoint = False): """ - Derivative of the magnetic flux density with respect to the inversion model. - + Derivative of the 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 :param bool adjoint: adjoint? @@ -321,8 +321,8 @@ class Fields_e(Fields): :return: product of the magnetic flux density derivative with respect to the inversion model with a vector """ - S_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) - return 1./(1j * omega(src.freq)) * S_mDeriv + s_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) + return 1./(1j * omega(src.freq)) * s_mDeriv def _j(self, eSolution, srcList): """ @@ -341,7 +341,7 @@ class Fields_e(Fields): def _jDeriv_u(self, src, du_dm_v, adjoint = False): """ Derivative of the current density with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -351,15 +351,15 @@ class Fields_e(Fields): n = int(self._aveE2CCV.shape[0] / self._nC) # number of components (instead of checking if cyl or not) VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - if adjoint: + if adjoint: return self._eDeriv_u(src, self._MeSigma.T * (self._aveE2CCV.T * (VI.T * du_dm_v) ), adjoint = adjoint) return VI * (self._aveE2CCV * (self._MeSigma * (self._eDeriv_u(src, du_dm_v, adjoint=adjoint) ) ) ) - + def _jDeriv_m(self, src, v, adjoint = False): """ - Derivative of the current density with respect to the inversion model. - + Derivative of the current density with respect to the inversion model. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -373,7 +373,7 @@ class Fields_e(Fields): if adjoint: return self._MeSigmaDeriv(e).T * (self._aveE2CCV.T * (VI.T * v)) + self._eDeriv_m(src, self._aveE2CCV.T * (VI.T * v), adjoint=adjoint) return VI * (self._aveE2CCV * ( self._eDeriv_m(src, v, adjoint=adjoint) + self._MeSigmaDeriv(e) * v)) - + def _h(self, eSolution, srcList): @@ -393,7 +393,7 @@ class Fields_e(Fields): def _hDeriv_u(self, src, du_dm_v, adjoint = False): """ Derivative of the magnetic field with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -409,8 +409,8 @@ class Fields_e(Fields): def _hDeriv_m(self, src, v, adjoint = False): """ - Derivative of the magnetic field with respect to the inversion model. - + Derivative of the magnetic field with respect to the inversion model. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -428,10 +428,10 @@ class Fields_e(Fields): class Fields_b(Fields): """ - Fields object for Problem_b. + Fields object for Problem_b. :param Mesh mesh: mesh - :param Survey survey: survey + :param Survey survey: survey """ knownFields = {'bSolution':'F'} @@ -506,9 +506,9 @@ class Fields_b(Fields): def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the total magnetic flux density with respect to the thing we + 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_v: vector to take product with :param bool adjoint: adjoint? @@ -520,8 +520,8 @@ class Fields_b(Fields): def _bDeriv_m(self, src, v, adjoint=False): """ - 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. - + 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 :param bool adjoint: adjoint? @@ -560,15 +560,15 @@ class Fields_b(Fields): e = ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): - _,S_e = src.eval(self.prob) - e[:,i] = e[:,i] + - S_e + _,s_e = src.eval(self.prob) + e[:,i] = e[:,i] + - s_e return self._MeSigmaI * e def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the 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 bool adjoint: adjoint? @@ -583,8 +583,8 @@ class Fields_b(Fields): def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the electric field with respect to the inversion model - + Derivative of the electric field with respect to the inversion model + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -593,15 +593,15 @@ class Fields_b(Fields): """ bSolution = Utils.mkvc(self[src, 'bSolution']) - _,S_e = src.eval(self.prob) + _,s_e = src.eval(self.prob) - w = -S_e + self._edgeCurl.T * (self._MfMui * bSolution) - _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint) + w = -s_e + self._edgeCurl.T * (self._MfMui * bSolution) + _, s_eDeriv = src.evalDeriv(self.prob, v, adjoint) if adjoint: - return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * S_eDeriv - return self._MeSigmaIDeriv(w) * v - self._MeSigmaI * S_eDeriv + return self._MeSigmaIDeriv(w).T * v - self._MeSigmaI.T * s_eDeriv + return self._MeSigmaIDeriv(w) * v - self._MeSigmaI * s_eDeriv def _j(self, bSolution, srcList): """ @@ -617,13 +617,13 @@ class Fields_b(Fields): VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (self._aveE2CCV * ( self._MeSigma * self._e(bSolution,srcList ) ) ) - + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the current density with respect to the thing we + Partial derivative of the current density with respect to the thing we solved for. - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -639,8 +639,8 @@ class Fields_b(Fields): def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the current density with respect to the inversion model - + Derivative of the current density with respect to the inversion model + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -664,9 +664,9 @@ class Fields_b(Fields): def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the magnetic field with respect to the thing we + Partial derivative of the magnetic field with respect to the thing we solved for. - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -682,8 +682,8 @@ class Fields_b(Fields): def _hDeriv_m(self, src, v, adjoint=False): """ - Derivative of the magnetic field with respect to the inversion model - + Derivative of the magnetic field with respect to the inversion model + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -695,10 +695,10 @@ class Fields_b(Fields): class Fields_j(Fields): """ - Fields object for Problem_j. + Fields object for Problem_j. :param Mesh mesh: mesh - :param Survey survey: survey + :param Survey survey: survey """ knownFields = {'jSolution':'F'} @@ -769,12 +769,12 @@ class Fields_j(Fields): def _j(self, jSolution, srcList): """ - Total current density is sum of primary and secondary - + Total current density is sum of primary and secondary + :param numpy.ndarray jSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: total current density + :return: total current density """ return self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, srcList) @@ -782,9 +782,9 @@ class Fields_j(Fields): def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the total current density with respect to the thing we + 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 bool adjoint: adjoint? @@ -797,8 +797,8 @@ class Fields_j(Fields): def _jDeriv_m(self, src, v, adjoint=False): """ - 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. - + 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 :param bool adjoint: adjoint? @@ -837,15 +837,15 @@ class Fields_j(Fields): h = (self._edgeCurl.T * (self._MfRho * jSolution) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) - S_m,_ = src.eval(self.prob) - h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (S_m) + s_m,_ = src.eval(self.prob) + h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m) return self._MeMuI * h def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the magnetic field with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -856,13 +856,13 @@ class Fields_j(Fields): if adjoint: return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * du_dm_v)) return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * du_dm_v) ) - + def _hDeriv_m(self, src, v, adjoint=False): """ - Derivative of the magnetic field with respect to the inversion model - + Derivative of the magnetic field with respect to the inversion model + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -875,19 +875,19 @@ class Fields_j(Fields): C = self._edgeCurl MfRho = self._MfRho MfRhoDeriv = self._MfRhoDeriv - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) + s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) if not adjoint: hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) - S_mDeriv = S_mDeriv(v) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( S_mDeriv) + s_mDeriv = s_mDeriv(v) + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * ( s_mDeriv) elif adjoint: hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - S_mDeriv = S_mDeriv(MeMuI.T * v) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv - + s_mDeriv = s_mDeriv(MeMuI.T * v) + hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * s_mDeriv + return hDeriv_m def _e(self, jSolution, srcList): @@ -901,12 +901,12 @@ class Fields_j(Fields): """ n = int(self._aveF2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList))) + return VI * (self._aveF2CCV * (self._MfRho * self._j(jSolution, srcList))) def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the electric field with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -921,8 +921,8 @@ class Fields_j(Fields): def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the electric field with respect to the inversion model - + Derivative of the electric field with respect to the inversion model + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -943,17 +943,17 @@ class Fields_j(Fields): :param numpy.ndarray hSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: secondary magnetic flux density + :return: secondary magnetic flux density """ n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) + return VI * (self._aveE2CCV * ( self._MeMu * self._h(jSolution,srcList)) ) def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the magnetic flux density with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -969,8 +969,8 @@ class Fields_j(Fields): def _bDeriv_m(self, src, v, adjoint=False): """ - Derivative of the magnetic flux density with respect to the inversion model - + Derivative of the 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 :param bool adjoint: adjoint? @@ -980,20 +980,20 @@ class Fields_j(Fields): jSolution = self[src,'jSolution'] n = int(self._aveE2CCV.shape[0] / self._nC) # number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) + s_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) if adjoint: v = self._aveE2CCV.T * ( VI.T * v) - return 1./(1j * omega(src.freq)) * ( S_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v )) - return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( S_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) ) + return 1./(1j * omega(src.freq)) * ( s_mDeriv(v) - self._MfRhoDeriv(jSolution).T * (self._edgeCurl * v )) + return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) ) class Fields_h(Fields): """ - Fields object for Problem_h. + Fields object for Problem_h. :param Mesh mesh: mesh - :param Survey survey: survey + :param Survey survey: survey """ knownFields = {'hSolution':'E'} @@ -1065,9 +1065,9 @@ class Fields_h(Fields): def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ - Partial derivative of the total magnetic field with respect to the thing we + 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_v: vector to take product with :param bool adjoint: adjoint? @@ -1079,8 +1079,8 @@ class Fields_h(Fields): def _hDeriv_m(self, src, v, adjoint=False): """ - 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. - + 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 :param bool adjoint: adjoint? @@ -1119,14 +1119,14 @@ class Fields_h(Fields): j = self._edgeCurl*hSolution for i, src in enumerate(srcList): - _,S_e = src.eval(self.prob) - j[:,i] = j[:,i]+ -S_e + _,s_e = src.eval(self.prob) + j[:,i] = j[:,i]+ -s_e return j def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the current density with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -1142,8 +1142,8 @@ class Fields_h(Fields): def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the current density with respect to the inversion model. - + Derivative of the current density with respect to the inversion model. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -1151,9 +1151,9 @@ class Fields_h(Fields): :return: product of the current density derivative with respect to the inversion model with a vector """ - _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - return -S_eDeriv - + _,s_eDeriv = src.evalDeriv(self.prob, v, adjoint) + return -s_eDeriv + def _e(self, hSolution, srcList): """ Electric field from hSolution @@ -1165,12 +1165,12 @@ class Fields_h(Fields): """ n = int(self._aveF2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList))) + return VI * (self._aveF2CCV * (self._MfRho * self._j(hSolution, srcList))) def _eDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the electric field with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? @@ -1181,12 +1181,12 @@ class Fields_h(Fields): VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) if adjoint: return self._edgeCurl.T * ( self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) ) - return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v )) + return VI * (self._aveF2CCV * (self._MfRho * self._edgeCurl * du_dm_v )) def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the electric field with respect to the inversion model. - + Derivative of the electric field with respect to the inversion model. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -1196,7 +1196,7 @@ class Fields_h(Fields): hSolution = Utils.mkvc(self[src,'hSolution']) n = int(self._aveF2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) - if adjoint: + if adjoint: return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) ) return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v )) @@ -1207,10 +1207,10 @@ class Fields_h(Fields): :param numpy.ndarray hSolution: field we solved for :param list srcList: list of sources :rtype: numpy.ndarray - :return: magnetic flux density + :return: magnetic flux density """ h = self._h(hSolution, srcList) - n = int(self._aveE2CCV.shape[0] / self._nC) #number of components + n = int(self._aveE2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) return VI * (self._aveE2CCV * (self._MeMu * h)) @@ -1218,14 +1218,14 @@ class Fields_h(Fields): def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ Derivative of the magnetic flux density with respect to the thing we solved for - + :param SimPEG.EM.FDEM.Src src: source :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 """ - n = int(self._aveE2CCV.shape[0] / self._nC) #number of components + n = int(self._aveE2CCV.shape[0] / self._nC) #number of components VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) if adjoint: return self._MeMu.T * (self._aveE2CCV.T * ( VI.T * du_dm_v )) @@ -1233,8 +1233,8 @@ class Fields_h(Fields): def _bDeriv_m(self, src, v, adjoint=False): """ - Derivative of the magnetic flux density with respect to the inversion model. - + Derivative of the 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 :param bool adjoint: adjoint? diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 31e4224f..87967dd5 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -15,22 +15,22 @@ class BaseSrc(Survey.BaseSrc): def eval(self, prob): """ Evaluate the source terms. - - :math:`S_m` : magnetic source term - - :math:`S_e` : electric source term + - :math:`s_m` : magnetic source term + - :math:`s_e` : electric source term :param Problem prob: FDEM Problem :rtype: (numpy.ndarray, numpy.ndarray) :return: tuple with magnetic source term and electric source term """ - S_m = self.S_m(prob) - S_e = self.S_e(prob) - return S_m, S_e + s_m = self.s_m(prob) + s_e = self.s_e(prob) + return s_m, s_e def evalDeriv(self, prob, v=None, adjoint=False): """ Derivatives of the source terms with respect to the inversion model - - :code:`S_mDeriv` : derivative of the magnetic source term - - :code:`S_eDeriv` : derivative of the electric source term + - :code:`s_mDeriv` : derivative of the magnetic source term + - :code:`s_eDeriv` : derivative of the electric source term :param Problem prob: FDEM Problem :param numpy.ndarray v: vector to take product with @@ -39,9 +39,9 @@ class BaseSrc(Survey.BaseSrc): :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) + 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): """ @@ -83,7 +83,7 @@ class BaseSrc(Survey.BaseSrc): """ return Zero() - def S_m(self, prob): + def s_m(self, prob): """ Magnetic source term @@ -93,7 +93,7 @@ class BaseSrc(Survey.BaseSrc): """ return Zero() - def S_e(self, prob): + def s_e(self, prob): """ Electric source term @@ -103,7 +103,7 @@ class BaseSrc(Survey.BaseSrc): """ return Zero() - def S_mDeriv(self, prob, v, adjoint = False): + def s_mDeriv(self, prob, v, adjoint = False): """ Derivative of magnetic source term with respect to the inversion model @@ -116,7 +116,7 @@ class BaseSrc(Survey.BaseSrc): return Zero() - def S_eDeriv(self, prob, v, adjoint = False): + def s_eDeriv(self, prob, v, adjoint = False): """ Derivative of electric source term with respect to the inversion model @@ -131,22 +131,22 @@ class BaseSrc(Survey.BaseSrc): class RawVec_e(BaseSrc): """ - RawVec electric source. It is defined by the user provided vector S_e + RawVec electric source. It is defined by the user provided vector s_e :param list rxList: receiver list :param float freq: frequency - :param numpy.array S_e: electric 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_e, integrate=True): #, 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): + def s_e(self, prob): """ Electric source term @@ -155,28 +155,28 @@ class RawVec_e(BaseSrc): :return: electric source term on mesh """ if prob._formulation is 'EB' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e + return prob.Me * self._s_e + return self._s_e class RawVec_m(BaseSrc): """ - RawVec magnetic source. It is defined by the user provided vector S_m + RawVec magnetic source. It is defined by the user provided vector s_m :param float freq: frequency :param rxList: receiver list - :param numpy.array S_m: magnetic source term + :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 BaseSrc.__init__(self, rxList) - def S_m(self, prob): + def s_m(self, prob): """ Magnetic source term @@ -185,28 +185,28 @@ class RawVec_m(BaseSrc): :return: magnetic source term on mesh """ if prob._formulation is 'HJ' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m + return prob.Me * self._s_m + return self._s_m class RawVec(BaseSrc): """ - RawVec source. It is defined by the user provided vectors S_m, S_e + RawVec source. It is defined by the user provided vectors s_m, s_e :param rxList: receiver list :param float freq: frequency - :param numpy.array S_m: magnetic source term - :param numpy.array S_e: electric source term + :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): + def s_m(self, prob): """ Magnetic source term @@ -215,10 +215,10 @@ class RawVec(BaseSrc): :return: magnetic source term on mesh """ if prob._formulation is 'HJ' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m + return prob.Me * self._s_m + return self._s_m - def S_e(self, prob): + def s_e(self, prob): """ Electric source term @@ -227,8 +227,8 @@ class RawVec(BaseSrc): :return: electric source term on mesh """ if prob._formulation is 'EB' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e + return prob.Me * self._s_e + return self._s_e class MagDipole(BaseSrc): @@ -335,9 +335,9 @@ class MagDipole(BaseSrc): :return: primary magnetic field """ b = self.bPrimary(prob) - return 1./self.mu * b + return 1./self.mu * b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term @@ -348,10 +348,10 @@ class MagDipole(BaseSrc): b_p = self.bPrimary(prob) if prob._formulation is 'HJ': - b_p = prob.Me * b_p + b_p = prob.Me * b_p return -1j*omega(self.freq)*b_p - def S_e(self, prob): + def s_e(self, prob): """ The electric source term @@ -453,7 +453,7 @@ class MagDipole_Bfield(BaseSrc): b = self.bPrimary(prob) return 1/self.mu * b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term @@ -466,7 +466,7 @@ class MagDipole_Bfield(BaseSrc): b = prob.Me * b return -1j*omega(self.freq)*b - def S_e(self, prob): + def s_e(self, prob): """ The electric source term @@ -565,7 +565,7 @@ class CircularLoop(BaseSrc): b = self.bPrimary(prob) return 1./self.mu*b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term @@ -578,7 +578,7 @@ class CircularLoop(BaseSrc): b = prob.Me * b return -1j*omega(self.freq)*b - def S_e(self, prob): + def s_e(self, prob): """ The electric source term @@ -604,6 +604,6 @@ class CircularLoop(BaseSrc): return -C.T * (MMui_s * self.bPrimary(prob)) - + From 579f1d7a65193bdf81cf05966363f7ef053641eb Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 13:36:19 -0700 Subject: [PATCH 6/7] FDEM uses f (so the u kwarg breaks). replace with f across the entire codebase?? --- SimPEG/DataMisfit.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index 425fe4ce..13574771 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -89,7 +89,7 @@ class l2_DataMisfit(BaseDataMisfit): """ if getattr(self, '_Wd', None) is None: - + survey = self.survey if getattr(survey,'std', None) is None: @@ -112,7 +112,7 @@ class l2_DataMisfit(BaseDataMisfit): "eval(m, u=None)" prob = self.prob survey = self.survey - R = self.Wd * survey.residual(m, u=u) + R = self.Wd * survey.residual(m, u) return 0.5*np.vdot(R, R) @Utils.timeIt @@ -121,11 +121,11 @@ class l2_DataMisfit(BaseDataMisfit): prob = self.prob survey = self.survey if u is None: u = prob.fields(m) - return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u=u)), u=u) + return prob.Jtvec(m, self.Wd * (self.Wd * survey.residual(m, u)), u) @Utils.timeIt def eval2Deriv(self, m, v, u=None): "eval2Deriv(m, v, u=None)" prob = self.prob if u is None: u = prob.fields(m) - return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u=u) + return prob.Jtvec_approx(m, self.Wd * (self.Wd * prob.Jvec_approx(m, v, u=u)), u) From c66db805af651598a8a087d4cd2e6fb0768008d6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 20 Mar 2016 14:44:34 -0700 Subject: [PATCH 7/7] typo fix --- SimPEG/EM/Base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 12d1ae43..a16cdb91 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -199,7 +199,7 @@ class BaseEMSurvey(Survey.BaseSurvey): :rtype: numpy.ndarray :return: data """ - data = SimPEG.Survey.Data(self) + data = Survey.Data(self) for src in self.srcList: for rx in src.rxList: data[src, rx] = rx.eval(src, self.mesh, u)