diff --git a/SimPEG/DCIP/BaseDC.py b/SimPEG/DCIP/BaseDC.py index e3d353d7..475c621e 100644 --- a/SimPEG/DCIP/BaseDC.py +++ b/SimPEG/DCIP/BaseDC.py @@ -200,11 +200,11 @@ class ProblemDC_CC(Problem.BaseProblem): return F - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: Jv @@ -225,11 +225,10 @@ class ProblemDC_CC(Problem.BaseProblem): # Set current model; clear dependent property $\mathbf{A(m)}$ self.curModel = m sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if u is None: + if f is None: # Run forward simulation if $u$ not provided - u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol'] - else: - u = u[self.survey.srcList, 'phi_sol'] + f = self.fields(self.curModel) + u = f[self.survey.srcList, 'phi_sol'] D = self.mesh.faceDiv G = self.mesh.cellGrad @@ -251,19 +250,18 @@ class ProblemDC_CC(Problem.BaseProblem): if self.Ainv is None: self.Ainv = self.Solver(dA_du, **self.solverOpts) - P = self.survey.getP(self.mesh) + P = self.survey.getP(self.mesh) Jv = - P * mkvc( self.Ainv * dCdm_x_v ) return Jv - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): self.curModel = m sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ - if u is None: - # Run forward simulation if $u$ not provided - u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol'] - else: - u = u[self.survey.srcList, 'phi_sol'] + if f is None: + # Run forward simulation if $f$ not provided + f = self.fields(self.curModel) + u = f[self.survey.srcList, 'phi_sol'] shp = u.shape P = self.survey.getP(self.mesh) diff --git a/SimPEG/DCIP/BaseIP.py b/SimPEG/DCIP/BaseIP.py index cec0ea2e..a18b5a47 100644 --- a/SimPEG/DCIP/BaseIP.py +++ b/SimPEG/DCIP/BaseIP.py @@ -14,12 +14,12 @@ class SurveyIP(SurveyDC): Survey.BaseSurvey.__init__(self, **kwargs) self._Ps = {} - def dpred(self, m, u=None): + def dpred(self, m, f=None): """ Predicted data. .. math:: - d_\\text{pred} = Pu(m) + d_\\text{pred} = Pf(m) """ return self.prob.forward(m) @@ -143,10 +143,10 @@ class ProblemIP(Problem.BaseProblem): J_x_v = - P * mkvc( self.Ainv * dCdm_x_v ) return -J_x_v - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): return self.forward(v) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): self.curModel = m # sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$ diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index 425fe4ce..53728c4e 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -22,11 +22,11 @@ class BaseDataMisfit(object): Utils.setKwargs(self,**kwargs) @Utils.timeIt - def eval(self, m, u=None): - """eval(m, u=None) + def eval(self, m, f=None): + """eval(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param Fields f: fields :rtype: float :return: data misfit @@ -34,11 +34,11 @@ class BaseDataMisfit(object): raise NotImplementedError('This method should be overwritten.') @Utils.timeIt - def evalDeriv(self, m, u=None): - """evalDeriv(m, u=None) + def evalDeriv(self, m, f=None): + """evalDeriv(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: data misfit derivative @@ -47,12 +47,12 @@ class BaseDataMisfit(object): @Utils.timeIt - def eval2Deriv(self, m, v, u=None): - """eval2Deriv(m, v, u=None) + def eval2Deriv(self, m, v, f=None): + """eval2Deriv(m, v, f=None) :param numpy.array m: geophysical model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: data misfit derivative @@ -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: @@ -108,24 +108,20 @@ class l2_DataMisfit(BaseDataMisfit): self._Wd = value @Utils.timeIt - def eval(self, m, u=None): - "eval(m, u=None)" - prob = self.prob - survey = self.survey - R = self.Wd * survey.residual(m, u=u) + def eval(self, m, f=None): + "eval(m, f=None)" + if f is None: f = self.prob.fields(m) + R = self.Wd * self.survey.residual(m, f) return 0.5*np.vdot(R, R) @Utils.timeIt - def evalDeriv(self, m, u=None): - "evalDeriv(m, u=None)" - 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) + def evalDeriv(self, m, f=None): + "evalDeriv(m, f=None)" + if f is None: f = self.prob.fields(m) + return self.prob.Jtvec(m, self.Wd * (self.Wd * self.survey.residual(m, f=f)), f=f) @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) + def eval2Deriv(self, m, v, f=None): + "eval2Deriv(m, v, f=None)" + if f is None: f = self.prob.fields(m) + return self.prob.Jtvec_approx(m, self.Wd * (self.Wd * self.prob.Jvec_approx(m, v, f=f)), f=f) diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 2ed27c20..3d726712 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -123,10 +123,10 @@ class BetaEstimate_ByEig(InversionDirective): if self.debug: print 'Calculating the beta0 parameter.' m = self.invProb.curModel - u = self.invProb.getFields(m, store=True, deleteWarmstart=False) + f = self.invProb.getFields(m, store=True, deleteWarmstart=False) x0 = np.random.rand(*m.shape) - t = x0.dot(self.dmisfit.eval2Deriv(m,x0,u=u)) + t = x0.dot(self.dmisfit.eval2Deriv(m,x0,f=f)) b = x0.dot(self.reg.eval2Deriv(m, v=x0)) self.beta0 = self.beta0_ratio*(t/b) diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index 32018f7e..a16cdb91 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 = 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/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index cf3dee7f..caca7602 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}\\\) """ @@ -36,76 +36,76 @@ class BaseFDEMProblem(BaseEMProblem): surveyPair = SurveyFDEM fieldsPair = Fields - def fields(self, m=None): + def fields(self, m): """ 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] + u_src = f[src, self._solutionType] dA_dm_v = self.getADeriv(freq, u_src, v) - dRHS_dm_v = self.getRHSDeriv(freq, 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,12 +120,12 @@ class BaseFDEMProblem(BaseEMProblem): ATinv = self.Solver(AT, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): - u_src = u[src, self._solutionType] + u_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 @@ -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,23 +154,23 @@ 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) - :return: S_m, S_e (nE or nF, nSrc) + :rtype: (numpy.ndarray, numpy.ndarray) + :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 ########################################################################################## @@ -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,25 +248,25 @@ 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) """ - 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): """ - 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 @@ -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): @@ -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,21 +373,21 @@ 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) """ - 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) @@ -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 :: @@ -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 @@ -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 @@ -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 @@ -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,26 +644,26 @@ 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) """ - 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): """ - 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 @@ -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)) - + diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index ce803ed1..1552a12c 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 @@ -66,7 +67,7 @@ class Rx(SimPEG.Survey.BaseRx): """Grid Location projection (e.g. Ex Fy ...)""" return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] - def eval(self, src, mesh, u): + def eval(self, src, mesh, f): """ Project fields to recievers to get data. @@ -79,27 +80,27 @@ class Rx(SimPEG.Survey.BaseRx): # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) # projGLoc += self.knownRxTypes[self.rxType][1] - P = self.getP(mesh, self.projGLoc(u)) - u_part_complex = u[src, self.projField] + P = self.getP(mesh, self.projGLoc(f)) + f_part_complex = f[src, self.projField] # get the real or imag component real_or_imag = self.projComp - u_part = getattr(u_part_complex, real_or_imag) - - return P*u_part + f_part = getattr(f_part_complex, real_or_imag) - def evalDeriv(self, src, mesh, u, v, adjoint=False): + return P*f_part + + def evalDeriv(self, src, mesh, f, 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 Fields f: fields object :param numpy.ndarray v: vector to multiply :rtype: numpy.ndarray :return: fields projected to recievers """ - P = self.getP(mesh, self.projGLoc(u)) + P = self.getP(mesh, self.projGLoc(f)) if not adjoint: Pv_complex = P * v @@ -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.') - diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 0da22072..15fc19e3 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -108,11 +108,11 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): Ainv.clean() return F - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ :param numpy.array m: Conductivity model :param numpy.ndarray v: vector (model object) - :param simpegEM.TDEM.FieldsTDEM u: Fields resulting from m + :param simpegEM.TDEM.FieldsTDEM f: Fields resulting from m :rtype: numpy.ndarray :return: w (data object) @@ -125,15 +125,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): """ if self.verbose: print '%s\nCalculating J(v)\n%s'%('*'*50,'*'*50) self.curModel = m - if u is None: - u = self.fields(m) - p = self.Gvec(m, v, u) + if f is None: + f = self.fields(m) + p = self.Gvec(m, v, f) y = self.solveAh(m, p) - Jv = self.survey.evalDeriv(u, v=y) + Jv = self.survey.evalDeriv(f, v=y) if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50) return - mkvc(Jv) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): """ :param numpy.array m: Conductivity model :param numpy.ndarray,SimPEG.Survey.Data v: vector (data object) @@ -150,15 +150,15 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): """ if self.verbose: print '%s\nCalculating J^T(v)\n%s'%('*'*50,'*'*50) self.curModel = m - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) if not isinstance(v, self.dataPair): v = self.dataPair(self.survey, v) - p = self.survey.evalDeriv(u, v=v, adjoint=True) + p = self.survey.evalDeriv(f, v=v, adjoint=True) y = self.solveAht(m, p) - w = self.Gtvec(m, y, u) + w = self.Gtvec(m, y, f) if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50) return - mkvc(w) diff --git a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py new file mode 100644 index 00000000..76af4a3d --- /dev/null +++ b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py @@ -0,0 +1,275 @@ +from SimPEG import * +from SimPEG.EM import FDEM, Analytics, mu_0 +import time + +try: + from pymatsolver import MumpsSolver + solver = MumpsSolver +except Exception: + solver = SolverLU + pass + +def run(plotIt=True): + """ + EM: Schenkel and Morrison Casing Model + ====================================== + + Here we create and run a FDEM forward simulation to calculate the vertical + current inside a steel-cased. The model is based on the Schenkel and + Morrison Casing Model, and the results are used in a 2016 SEG abstract by + Yang et al. + + - Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686. + + + The model consists of: + - Air: Conductivity 1e-8 S/m, above z = 0 + - Background: conductivity 1e-2 S/m, below z = 0 + - Casing: conductivity 1e6 S/m + - 300m long + - radius of 0.1m + - thickness of 6e-3m + + Inside the casing, we take the same conductivity as the background. + + We are using an EM code to simulate DC, so we use frequency low enough + that the skin depth inside the casing is longer than the casing length (f + = 1e-6 Hz). The plot produced is of the current inside the casing. + + These results are shown in the SEG abstract by Yang et al., 2016: 3D DC + resistivity modeling of steel casing for reservoir monitoring using + equivalent resistor network. The solver used to produce these results and + achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_ + + .. _pymatsolver: https://github.com/rowanc1/pymatsolver + + This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1 + + If you would use this example for a code comparison, or build upon it, a + citation would be much appreciated! + + """ + + if plotIt: + import matplotlib.pylab as plt + + # ------------------ MODEL ------------------ + sigmaair = 1e-8 # air + sigmaback = 1e-2 # background + sigmacasing = 1e6 # casing + sigmainside = sigmaback # inside the casing + + + casing_t = 0.006 # 1cm thickness + casing_l = 300 # length of the casing + + casing_r = 0.1 + casing_a = casing_r - casing_t/2. # inner radius + casing_b = casing_r + casing_t/2. # outer radius + casing_z = np.r_[-casing_l,0.] + + + # ------------------ SURVEY PARAMETERS ------------------ + freqs = np.r_[1e-6] #[1e-1, 1, 5] # frequencies + dsz = -300 # down-hole z source location + src_loc = np.r_[0.,0.,dsz] + inf_loc = np.r_[0.,0.,1e4] + + print 'Skin Depth: ', [(500./np.sqrt(sigmaback*_)) for _ in freqs] + + + # ------------------ MESH ------------------ + # fine cells near well bore + csx1, csx2 = 2e-3, 60. + pfx1, pfx2 = 1.3, 1.3 + ncx1 = np.ceil(casing_b/csx1+2) + + # pad nicely to second cell size + npadx1 = np.floor(np.log(csx2/csx1) / np.log(pfx1)) + hx1a,hx1b = Utils.meshTensor([(csx1,ncx1)]),Utils.meshTensor([(csx1,npadx1,pfx1)]) + dx1 = sum(hx1a)+sum(hx1b) + dx1 = np.floor(dx1/csx2) + hx1b *= (dx1*csx2 - sum(hx1a))/sum(hx1b) + + # second chunk of mesh + dx2 = 300. # uniform mesh out to here + ncx2 = np.ceil((dx2 - dx1)/csx2) + npadx2 = 45 + hx2a, hx2b = Utils.meshTensor([(csx2,ncx2)]), Utils.meshTensor([(csx2,npadx2,pfx2)]) + hx = np.hstack([hx1a,hx1b,hx2a,hx2b]) + + # z-direction + csz = 0.05 + nza = 10 + ncz, npadzu, npadzd = np.int(np.ceil(np.diff(casing_z)[0]/csz))+10, 68, 68 # cell size, number of core cells, number of padding cells in the x- direction + hz = Utils.meshTensor([(csz,npadzd,-1.3), (csz,ncz), (csz,npadzu,1.3)]) # vector of cell widths in the z-direction + + # Mesh + mesh = Mesh.CylMesh([hx,1.,hz], [0.,0.,-np.sum(hz[:npadzu+ncz-nza])]) + + print 'Mesh Extent xmax: %f,: zmin: %f, zmax: %f'%(mesh.vectorCCx.max(), mesh.vectorCCz.min(), mesh.vectorCCz.max()) + print 'Number of cells', mesh.nC + + if plotIt is True: + fig, ax = plt.subplots(1, 1, figsize=(6, 4)) + ax.set_title('Simulation Mesh') + mesh.plotGrid(ax=ax) + plt.show() + + # Put the model on the mesh + sigWholespace = sigmaback*np.ones((mesh.nC)) + + sigBack = sigWholespace.copy() + sigBack[mesh.gridCC[:,2] > 0.] = sigmaair + + sigCasing = sigBack.copy() + iCasingZ = (mesh.gridCC[:,2] <= casing_z[1]) & (mesh.gridCC[:,2] >= casing_z[0]) + iCasingX = (mesh.gridCC[:,0] >= casing_a) & (mesh.gridCC[:,0] <= casing_b) + iCasing = iCasingX & iCasingZ + sigCasing[iCasing] = sigmacasing + + + if plotIt is True: + + # plotting parameters + xlim = np.r_[0., 0.2] + zlim = np.r_[-350., 10.] + clim_sig = np.r_[-8,6] + + # plot models + fig, ax = plt.subplots(1,1,figsize=(4,4)) + + f = plt.colorbar(mesh.plotImage(np.log10(sigCasing),ax=ax)[0], ax=ax) + ax.grid(which='both') + ax.set_title('Log_10 (Sigma)') + ax.set_xlim(xlim) + ax.set_ylim(zlim) + f.set_clim(clim_sig) + + plt.show() + + + # -------------- Sources -------------------- + # Define Custom Current Sources + + # surface source + sg_x = np.zeros(mesh.vnF[0],dtype=complex) + sg_y = np.zeros(mesh.vnF[1],dtype=complex) + sg_z = np.zeros(mesh.vnF[2],dtype=complex) + + nza = 2 # put the wire two cells above the surface + ncin = 2 + + # vertically directed wire + sgv_indx = (mesh.gridFz[:,0] > casing_a) & (mesh.gridFz[:,0] < casing_a + csx1) # hook it up to casing at the surface + sgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2) + sgv_ind = sgv_indx & sgv_indz + sg_z[sgv_ind] = -1. + + # horizontally directed wire + sgh_indx = (mesh.gridFx[:,0] > casing_a) & (mesh.gridFx[:,0] <= inf_loc[2]) + sgh_indz = (mesh.gridFx[:,2] > csz*(nza-0.5)) & (mesh.gridFx[:,2] < csz*(nza+0.5)) + sgh_ind = sgh_indx & sgh_indz + sg_x[sgh_ind] = -1. + + sgv2_indx = (mesh.gridFz[:,0] >= mesh.gridFx[sgh_ind,0].max()) & (mesh.gridFz[:,0] <= inf_loc[2]*1.2) # hook it up to casing at the surface + sgv2_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] >= -csz*2) + sgv2_ind = sgv2_indx & sgv2_indz + sg_z[sgv2_ind] = 1. + + # assemble the source + sg = np.hstack([sg_x,sg_y,sg_z]) + sg_p = [FDEM.Src.RawVec_e([],_,sg/mesh.area) for _ in freqs] + + # downhole source + dg_x = np.zeros(mesh.vnF[0],dtype=complex) + dg_y = np.zeros(mesh.vnF[1],dtype=complex) + dg_z = np.zeros(mesh.vnF[2],dtype=complex) + + # vertically directed wire + dgv_indx = (mesh.gridFz[:,0] < csx1) # go through the center of the well + dgv_indz = (mesh.gridFz[:,2] <= +csz*nza) & (mesh.gridFz[:,2] > dsz + csz/2.) + dgv_ind = dgv_indx & dgv_indz + dg_z[dgv_ind] = -1. + + # couple to the casing downhole + dgh_indx = mesh.gridFx[:,0] < casing_a + csx1 + dgh_indz = (mesh.gridFx[:,2] < dsz + csz) & (mesh.gridFx[:,2] >= dsz) + dgh_ind = dgh_indx & dgh_indz + dg_x[dgh_ind] = 1. + + # horizontal part at surface + dgh2_indx = mesh.gridFx[:,0] <= inf_loc[2]*1.2 + dgh2_indz = sgh_indz.copy() + dgh2_ind = dgh2_indx & dgh2_indz + dg_x[dgh2_ind] = -1. + + # vertical part at surface + dgv2_ind = sgv2_ind.copy() + dg_z[dgv2_ind] = 1. + + # assemble the source + dg = np.hstack([dg_x,dg_y,dg_z]) + dg_p = [FDEM.Src.RawVec_e([],_,dg/mesh.area) for _ in freqs] + + # ------------ Problem and Survey --------------- + survey = FDEM.Survey(sg_p + dg_p) + mapping = [('sigma', Maps.IdentityMap(mesh))] + problem = FDEM.Problem_h(mesh, mapping=mapping) + problem.pair(survey) + + # ------------- Solve --------------------------- + t0 = time.time() + fieldsCasing = problem.fields(sigCasing) + print 'Time to solve 2 sources', time.time() - t0 + + # Plot current + + # current density + jn0 = fieldsCasing[dg_p,'j'] + jn1 = fieldsCasing[sg_p,'j'] + + # current + in0 = [mesh.area*fieldsCasing[dg_p,'j'][:,i] for i in range(len(freqs))] + in1 = [mesh.area*fieldsCasing[sg_p,'j'][:,i] for i in range(len(freqs))] + + in0 = np.vstack(in0).T + in1 = np.vstack(in1).T + + # integrate to get z-current inside casing + inds_inx = (mesh.gridFz[:,0] >= casing_a) & (mesh.gridFz[:,0] <= casing_b) + inds_inz = (mesh.gridFz[:,2] >= dsz ) & (mesh.gridFz[:,2] <= 0) + inds_fz = inds_inx & inds_inz + + indsx = [False]*mesh.nFx + inds = list(indsx) + list(inds_fz) + + in0_in = in0[np.r_[inds]] + in1_in = in1[np.r_[inds]] + z_in = mesh.gridFz[inds_fz,2] + + in0_in = in0_in.reshape([in0_in.shape[0]/3,3]) + in1_in = in1_in.reshape([in1_in.shape[0]/3,3]) + z_in = z_in.reshape([z_in.shape[0]/3,3]) + + I0 = in0_in.sum(1).real + I1 = in1_in.sum(1).real + z_in = z_in[:,0] + + if plotIt is True: + fig, ax = plt.subplots(1,2,figsize=(12,4)) + + ax[0].plot(z_in,np.absolute(I0), z_in,np.absolute(I1)) + ax[0].legend(['top casing', 'bottom casing'],loc='best') + ax[0].set_title('Magnitude of Vertical Current in Casing') + + ax[1].semilogy(z_in,np.absolute(I0), z_in,np.absolute(I1)) + ax[1].legend(['top casing', 'bottom casing'],loc='best') + ax[1].set_title('Magnitude of Vertical Current in Casing') + ax[1].set_ylim([1e-2, 1.]) + + plt.show() + +if __name__ == '__main__': + run() + diff --git a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py index 6de7542e..62ab2a68 100644 --- a/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py +++ b/SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py @@ -1,7 +1,8 @@ from scipy.constants import epsilon_0, mu_0 import matplotlib.pyplot as plt import numpy as np -#from SimPEG.EM.Utils import k, omega +from ipywidgets import * +from SimPEG.EM.Utils import k, omega """ MT1D: n layered earth problem @@ -15,51 +16,45 @@ This code compute the analytic response of a n-layered Earth to a plane wave (Ma We start by looking at Maxwell's equations in the electric field \\\(\\\mathbf{E}\\) and the magnetic flux -\\\(\\\mathbf{H}\\) to write the wave equations +\\\(\\\mathbf{H}\\) to write the wave equations \\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) & \\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\) -Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of +Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of \\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\) \\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\) With U and D the Up and Down components of the E-field. -The iteration from one layer to another is ensure by: +The iteration from one layer to another is ensure by: -\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = +\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) = P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\) -And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0) +And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0) and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability) The layer 0 is assumed to be the air layer. """ -#Frequency conversion -omega = lambda f: 2.*np.pi*f - -#Evaluate k wavenumber -k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f)) - #Define a frquency range for a survey frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.) #Functions to create random physical Perties for a n-layered earth -thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]), +thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]), np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1) ,decimals =1)) -sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]), +sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]), np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1) ,decimals=3)) -mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]), +mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]), np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1) ,decimals=1)) -eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]), +eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]), np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1) ,decimals=1)) @@ -69,17 +64,8 @@ ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k #Complex Cole-Cole Conductivity - EM utils PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c))) - #Converted thickness array into top of layer array -def top(thick): - topv= np.zeros(len(thick)+1) - - topv[0]=-thick[0] - - for i in range(1,len(topv),1): - topv[i] = topv[i-1] + thick[i-1] - - return topv +top = lambda thick: np.cumsum(thick) #Propagation Matrix and theirs inverses @@ -104,36 +90,36 @@ H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U) #Plot the configuration of the problem def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): - + topn = top(thick) widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.) - + ax.set_ylim([z.min(),z.max()]) ax.set_xlim([-widthg,widthg]) - + ax.set_ylabel("Depth (m)", fontsize=16.) ax.yaxis.tick_right() ax.yaxis.set_label_position("right") - + #define filling for the different layers - hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ] - + hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ] + #Write the physical properties of air ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)), xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data', xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data', fontsize=14.) - + ax.annotate(("$\epsilon_r$= %1i")%(eps[0]), xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data', fontsize=14.) - + ax.annotate(("$\mu_r$= %1i")%(mu[0]), xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data', xytext=(0, -np.abs(z.max())/3.), textcoords='data', fontsize=14.) - + #Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern for i in range(1,len(topn)-1,1): if topn[i] == topn[i+1]: @@ -143,39 +129,39 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data', xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data', fontsize=14.) - + ax.annotate(("$\epsilon_r$= %1i")%(eps[i]), xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data', xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data', fontsize=14.) - + ax.annotate(("$\mu_r$= %1.2f")%(mu[i]), xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data', xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data', fontsize=14.) - - ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black') + + ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black') ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10]) - + #Write the physical properties of the n-th layer and fill it with pattern - ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black') + ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black') ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10]) - + ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)), xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data', xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data', fontsize=14.) - + ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]), xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data', xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data', fontsize=14.) - + ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]), xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data', xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data', fontsize=14.) - + #plot Trees! ax.annotate("", xy=(widthg/2., -1.*z.max()/5.), xycoords='data', @@ -194,7 +180,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): xytext=(widthg/2., 0.), textcoords='data', arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) ) - + ax.annotate("", xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data', xytext=(1.2*widthg/2., 0.), textcoords='data', @@ -231,7 +217,7 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.) ) - + ax.invert_yaxis() return ax @@ -239,83 +225,83 @@ def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z): #Propagate Up and Down component for a certain frequency & evaluate E and H field def Propagate(f,H,sig,chg,taux,c,mu,eps,n): - + sigcm = np.zeros_like(sig,dtype='complex_') - + for j in range(1,len(sig)): sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f) - - K = k(mu,sigcm,eps,f) + + K = k(f, sigcm, mu, eps) Z = ImpZ(f,mu,K) - + EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_') UD[1,-1] = 1. - + for i in range(-2,-(n+2),-1): - + UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1] UD = UD/((np.abs(UD[0,:]+UD[1,:])).max()) - - for j in range(0,n+1): + + for j in range(0,n+1): EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j] return UD, EH, Z ,K - + #Evaluate the apparent resistivity and phase for a frequency range def appres(F,H,sig,chg,taux,c,mu,eps,n): - + Res = np.zeros_like(F) Phase = np.zeros_like(F) App_ImpZ= np.zeros_like(F,dtype='complex_') - + for i in range(0,len(F)): - + UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n) - + App_ImpZ[i] = EH[0,1]/EH[1,1] - + Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i])) Phase[i] = np.angle(App_ImpZ[i], deg = True) - + return Res,Phase -#Evaluate Up, Down components, E and H field, for a frequency range, -#a discretized depth range and a time range (use to calculate envelope) +#Evaluate Up, Down components, E and H field, for a frequency range, +#a discretized depth range and a time range (use to calculate envelope) def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample): - + topc = top(H) - + layer = np.zeros(len(zsample),dtype=np.int)-1 - + Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_') UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_') - + for i in range(0,n+1,1): layer = layer+(zsample>=topc[i])*1 - + for j in range(0,len(F)): - + UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n) - + for p in range(0,len(zsample)): - + UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]]) - + for q in range(0,len(tsample)): - + Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F) Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F) Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F) - + return Exzt,Hyzt,Uz,Dz,UDaux,layer - + #Function to Plot Apparent Resistivity and Phase def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): @@ -334,44 +320,44 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): ax[0].grid(which='major') ax0 = ax[0].twiny() - + ax0.set_xlim([0.,90.]) ax0.set_ylim([F.min(),F.max()]) ax0.scatter(Phase,F,color='purple') ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple") - + zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.) - + ax[0].tick_params(labelsize=16) ax[1].tick_params(labelsize=16) ax0.tick_params(labelsize=16) - + if PlotEnvelope: - + widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0) fenvelope1n=np.ones(100)*fenvelope ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black') - + tc=np.arange(0.,1./fenvelope,0.01/(fenvelope)) Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc) - + ax1=ax[1].twiny() - + ax[1].tick_params(labelsize=16) ax1.tick_params(labelsize=16) ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16) ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16) - + ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))), - np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))), + np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))), color='blue', alpha=0.1) ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))), - np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))), + np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))), color='red', alpha=0.1) - + ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc) ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()]) @@ -379,12 +365,12 @@ def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope): print 'No envelop (if True, might be slow)' ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc) ax[1].get_xaxis().set_ticks([]) - + plt.show() #Interactive MT for Notebook def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope): - + frangn=frange(-5,5,100.) sig3= np.array([0.,0.001,0.1, 0.001]) thick3 = np.array([120000.,50.,50.]) @@ -394,7 +380,7 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2 chg3_0=np.array([0.,0.1,0.,0.]) taux3=np.array([0.,0.1,0.,0.1]) c3=np.array([1.,1.,1.,1.]) - + sig3[1]=sigl1 sig3[1]=10.**sig3[1] sig3[2]=sigl2 @@ -409,11 +395,11 @@ def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2 eps3[3]=epsl3 thick3[1]=h1 thick3[2]=h2 - - PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope) - -def run(n=3,plotIt=True): + PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope) + + +def run(n,plotIt=True): # something to make a plot F = frange(-5.,5.,20) @@ -429,14 +415,14 @@ def run(n=3,plotIt=True): if plotIt: - PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True) + PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True) return Res, Phase if __name__ == '__main__': - run() - - - + run(3) + + + diff --git a/SimPEG/Examples/__init__.py b/SimPEG/Examples/__init__.py index 8cf7dea0..f8e189c8 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -3,8 +3,11 @@ ##### AUTOIMPORTS ##### import DC_Analytic_Dipole import DC_Forward_PseudoSection +import DC_PseudoSection_Simulation import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace +import EM_FDEM_SusEffects +import EM_Schenkel_Morrison_Casing import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent @@ -16,17 +19,12 @@ import Mesh_QuadTree_Creation import Mesh_QuadTree_FaceDiv import Mesh_QuadTree_HangingNodes import Mesh_Tensor_Creation -<<<<<<< HEAD import MT_1D_analytic_nlayer_Earth -import sphereElectrostatic_example - -__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "sphereElectrostatic_example"] -======= import MT_1D_ForwardAndInversion import MT_3D_Foward +import sphereElectrostatic_example -__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"] ->>>>>>> master +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "sphereElectrostatic_example"] ##### AUTOIMPORTS ##### diff --git a/SimPEG/Examples/sphereElectrostatic_example.py b/SimPEG/Examples/sphereElectrostatic_example.py index aa5e3d8d..7bd633c5 100644 --- a/SimPEG/Examples/sphereElectrostatic_example.py +++ b/SimPEG/Examples/sphereElectrostatic_example.py @@ -16,8 +16,6 @@ finally the charges accumulation. Several plotting functions are defined for data visualisation. -Please visit http://em.geosci.xyz/en/latest/content/maxwell2_steady_state/electrostatic_sphere.html -for more examples using this code. ''' @@ -34,7 +32,7 @@ sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0) def conductivity_log_wrapper(log_sig0,log_sig1): sig0 = 10.**log_sig0 sig1 = 10.**log_sig1 - + return sig0,sig1 # Examples @@ -56,7 +54,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): dx = xr[1]-xr[0] top = np.sqrt(R**2-xplt**2) bot = -np.sqrt(R**2-xplt**2) - + if R != 0: ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5) ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 ) @@ -87,7 +85,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): ax.set_xticklabels([]) ax.set_yticklabels([]) ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14) - ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14) + ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14) ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'), xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data', xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data', @@ -98,7 +96,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): fontsize=14.) ax.set_xlabel('x',fontsize=12) ax.set_ylabel('y',fontsize=12) - + else: if label: ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)), @@ -116,7 +114,7 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): else: ax.set_xticklabels([]) ax.set_yticklabels([]) - ax.text(-0.05,-10,'$\sigma_0$',fontsize=14) + ax.text(-0.05,-10,'$\sigma_0$',fontsize=14) ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14) ax.set_xlabel('x',fontsize=12) ax.set_ylabel('y',fontsize=12) @@ -130,8 +128,8 @@ def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere): ax.set_aspect('equal') - - + + return ax def get_Conductivity(XYZ,sig0,sig1,R): @@ -140,61 +138,61 @@ def get_Conductivity(XYZ,sig0,sig1,R): ''' x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2] r_view=r(x,y,z) - + ind0= (r_view>R) ind1= (r_view<=R) - + assert (ind0 + ind1).all(), 'Some indicies not included' - + Sigma = np.zeros_like(x) - + Sigma[ind0] = sig0 Sigma[ind1] = sig1 - + return Sigma -def get_Potential(XYZ,sig0,sig1,R,E0): +def get_Potential(XYZ,sig0,sig1,R,E0): ''' Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin :input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field ''' - + x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2] - + sig_cur = sigf(sig0,sig1) - + r_cur = r(x,y,z) # current radius - + ind0 = (r_cur > R) ind1 = (r_cur <= R) - + assert (ind0 + ind1).all(), 'Some indicies not included' - + Vt = np.zeros_like(x) Vp = np.zeros_like(x) Vs = np.zeros_like(x) - + Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0) # inside the sphere - - + + Vp = - E0*x # primary potential - + Vs = Vt - Vp # secondary potential - + return Vt,Vp,Vs #plot the primary potential on ax def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax): - + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) - + xcirc = xr[np.abs(xr) <= R] - + Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size)) ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') ax.set_title('Primary Potential',fontsize=ftsize_title) @@ -207,19 +205,19 @@ def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) ax.set_aspect('equal') ax.tick_params(labelsize=ftsize_axis) - + return ax #plot the total potential on ax def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax): - + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) - + xcirc = xr[np.abs(xr) <= R] - + Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size)) ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') ax.set_title('Total Potential',fontsize=ftsize_title) @@ -232,16 +230,16 @@ def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) ax.set_aspect('equal') ax.tick_params(labelsize=ftsize_axis) - + return ax #plot the secondary potential on ax def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax): - + Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) - + xcirc = xr[np.abs(xr) <= R] Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size)) @@ -256,30 +254,30 @@ def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label) ax.set_aspect('equal') ax.tick_params(labelsize=ftsize_axis) - + return ax def get_ElectricField(XYZ,sig0,sig1,R,E0): ''' - Function that returns the total, the primary and the secondary electric fields, + Function that returns the total, the primary and the secondary electric fields, input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field ''' - + x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] - + r_cur=r(x,y,z) # current radius - + ind0= (r_cur>R) ind1= (r_cur<=R) - + assert (ind0 + ind1).all(), 'Some indicies not included' - + Ep = np.zeros(shape=(len(x),3)) Ep[:,0] = E0 - + Et = np.zeros(shape=(len(x),3)) - + Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.); Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1); Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1); @@ -287,16 +285,16 @@ def get_ElectricField(XYZ,sig0,sig1,R,E0): Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0; Et[ind1,1] = 0.; Et[ind1,2] = 0.; - + Es = Et - Ep - + return Et, Ep, Es #plot the total electric field on ax def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): - + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) xcirc = xr[np.abs(xr) <= R] @@ -304,7 +302,7 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): EtXr = Et[:,0].reshape(xr.size, yr.size) EtYr = Et[:,1].reshape(xr.size, yr.size) EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size) - + ax.set_xlim([xr.min(),xr.max()]) ax.set_ylim([yr.min(),yr.max()]) ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) @@ -312,22 +310,22 @@ def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax): ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') ax.tick_params(labelsize=ftsize_axis) ax.set_aspect('equal') - + Eplot = ax.pcolor(xr,yr,EtAmp) cb = plt.colorbar(Eplot,ax=ax) cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold') cb.ax.tick_params(labelsize=ftsize_axis) ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05) ax.set_title('Total Field',fontsize=ftsize_title) - - + + return ax - -#plot the secondary electric field on ax + +#plot the secondary electric field on ax def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): - + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) xcirc = xr[np.abs(xr) <= R] @@ -335,7 +333,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): EsXr = Es[:,0].reshape(xr.size, yr.size) EsYr = Es[:,1].reshape(xr.size, yr.size) EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size) - + ax.set_xlim([xr.min(),xr.max()]) ax.set_ylim([yr.min(),yr.max()]) ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label) @@ -343,7 +341,7 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') ax.tick_params(labelsize=ftsize_axis) ax.set_aspect('equal') - + Eplot = ax.pcolor(xr,yr,EsAmp) cb = plt.colorbar(Eplot,ax=ax) cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold') @@ -351,53 +349,53 @@ def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax): ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05) ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') ax.set_title('Secondary Field',fontsize=ftsize_title) - + return ax def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es): ''' - Function that returns the total, the primary and the secondary current densities, + Function that returns the total, the primary and the secondary current densities, :input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields, ''' - + x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] - + r_cur=r(x,y,z) - + ind0= (r_cur>R) ind1= (r_cur<=R) - + assert (ind0 + ind1).all(), 'Some indicies not included' - + Jt = np.zeros(shape=(len(x),3)) J0 = np.zeros(shape=(len(x),3)) Js = np.zeros(shape=(len(x),3)) - + Jp = sig0*Ep - - Jt[ind0,:] = sig0*Et[ind0,:] + + Jt[ind0,:] = sig0*Et[ind0,:] Jt[ind1,:] = sig1*Et[ind1,:] Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:]) Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:] - + return Jt,Jp,Js #plot the total currents density on ax def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax): - + Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) xcirc = xr[np.abs(xr) <= R] JtXr = Jt[:,0].reshape(xr.size, yr.size) JtYr = Jt[:,1].reshape(xr.size, yr.size) JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size) - + ax.set_xlim([xr.min(),xr.max()]) ax.set_ylim([yr.min(),yr.max()]) ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') @@ -405,30 +403,30 @@ def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) ax.tick_params(labelsize=ftsize_axis) ax.set_aspect('equal') - + Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size)) cb = plt.colorbar(Jplot,ax=ax) cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold') cb.ax.tick_params(labelsize=ftsize_axis) ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1) ax.set_title('Total Current Density',fontsize=ftsize_title) - + return ax #plot the secondary currents density on ax def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax): - + Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0) Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) xcirc = xr[np.abs(xr) <= R] - + JsXr = Js[:,0].reshape(xr.size, yr.size) JsYr = Js[:,1].reshape(xr.size, yr.size) JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Jt[:,2]**2).reshape(xr.size,yr.size) - + ax.set_xlim([xr.min(),xr.max()]) ax.set_ylim([yr.min(),yr.max()]) ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k') @@ -436,52 +434,52 @@ def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) ax.tick_params(labelsize=ftsize_axis) ax.set_aspect('equal') - + Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size)) cb = plt.colorbar(Jplot,ax=ax) cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold') cb.ax.tick_params(labelsize=ftsize_axis) ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1) ax.set_title('Secondary Current Density',fontsize=ftsize_title) - + return ax def get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep): ''' - Function that returns the charges accumulation at the background/sphere interface, + Function that returns the charges accumulation at the background/sphere interface, :input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields, ''' x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2] - + dx = x[1]-x[0] - + r_cur=r(x,y,z) - + ind0 = (r_cur > R) ind1 = (r_cur < R) ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) ) - + assert (ind0 + ind1 + ind2).all(), 'Some indicies not included' - + rho = np.zeros_like(x) - + rho[ind0] = 0 rho[ind1] = 0 rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.)) - + return rho #Plot charges density on ax def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax): - + xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) xcirc = xr[np.abs(xr) <= R] - + Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep) - + ax.set_xlim([xr.min(),xr.max()]) ax.set_ylim([yr.min(),yr.max()]) ax.set_aspect('equal') @@ -494,11 +492,11 @@ def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax): ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label) ax.tick_params(labelsize=ftsize_axis) ax.set_title('Charges Density', fontsize=ftsize_title) - + return ax def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): - + ''' Function that return array of midpoints electrodes, electrodes positions, potentials differences for total and secondary potentials fields, unormalized and @@ -515,20 +513,20 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): #D: total distance from start to end D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.) - + #MP: dipoles'midpoint positions (x,y) - MP = np.zeros(shape=(nbmp,2)) + MP = np.zeros(shape=(nbmp,2)) MP[:,0] = np.linspace(start[0],end[0],nbmp) MP[:,1] = np.linspace(start[1],end[1],nbmp) - + #Dipoles'Electrodes positions around each midpoints - EL = np.zeros(shape=(2*nbmp,2)) + EL = np.zeros(shape=(2*nbmp,2)) for n in range(0,len(EL),2): EL[n,0] = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2. EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2. EL[n,1] = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2. EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2. - + VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL) VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL) dVtMP = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) @@ -537,109 +535,109 @@ def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn): dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n) dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n) - - #Computing VtEL + + #Computing VtEL for m in range(0,2*nbmp): if (r(EL[m,0],EL[m,1],0) > R): VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.) else: VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0) - + #Computing VsEL VsEL = VtEL + E0*EL[:,0] - + #Computing dVtMP, dVsMP for p in range(0,nbmp): dVtMP[p] = VtEL[2*p]-VtEL[2*p+1] dVtMPn[p] = dVtMP[p]/mn dVsMP[p] = VsEL[2*p]-VsEL[2*p+1] dVsMPn[p] = dVsMP[p]/mn - + return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn #Compare the DC response of two configurations -def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt,ax): - +def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor): + #Define the mesh xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2]) - + #Defining the Profile start = np.array([xstart,ystart]) end = np.array([xend,yend]) - + #Calculating the data from the defined survey line for Configuration 0 and 1 MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing) MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing) # Initializing the figure - #fig = plt.figure(figsize=(20,20)) - #ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) - #ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) - #ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) - #ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) - #ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) + fig = plt.figure(figsize=(20,20)) + ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) + ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) + ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) + ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) + ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) #Plotting the Configuration 0 - ax[0] = get_Setup(XYZ,sig0,sig1,R0,E0,ax[0],True,[0.6,0.1,0.1]) - + ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1]) + #Plotting the Configuration 1 - ax[1] = get_Setup(XYZ,sig0,sig2,R1,E0,ax[1],True,[0.1,0.1,0.6]) - + ax1 = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6]) + #Plotting the Data (Legends) - ax[2].set_title('Potential Differences',fontsize=ftsize_title) - ax[2].set_ylabel('Potential difference ($V$)',fontsize=ftsize_label) - ax[2].set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label) - ax[2].tick_params(labelsize=ftsize_axis) - ax[2].grid() + ax2.set_title('Potential Differences',fontsize=ftsize_title) + ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label) + ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label) + ax2.tick_params(labelsize=ftsize_axis) + ax2.grid() if PlotOpt == 'Total': - ax[3]= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax[3]) - ax[4]= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax[4]) - - #Plot the Data (from Configuration 0) - gphy0 = ax[2].plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0 + ax3= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax3) + ax4= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax4) + + #Plot the Data (from Configuration 0) + gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0 ,marker='o',color='blue',linewidth=3.,label ='Left Model Response' ) #Plot the Data (from Configuration 1) - gphy1 = ax[2].plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1 + gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1 ,marker='o',color='red',linewidth=2.,label ='Right Model Response' ) - ax[2].legend(('Left Model Response','Right Model Response'),loc=4) + ax2.legend(('Left Model Response','Right Model Response'),loc=4) elif PlotOpt == 'Secondary': #plot the secondary potentials - ax[3]= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax[3]) - ax[4]= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax[4]) - + ax3= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax3) + ax4= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax4) + #Plot the data(from configuration 0) - gphy0 = ax[2].plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue' + gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue' ,marker='o',linewidth=3.,label ='Left Model Response' ) - + #Plot the Data (from Configuration 1) - gphy1 = ax[2].plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1 + gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1 ,marker='o',color='red',linewidth=2.,label ='Right Model Response' ) - ax[2].legend(('Left Model Response','Right Model Response'),loc=4 ) - + ax2.legend(('Left Model Response','Right Model Response'),loc=4 ) + else: print('What dont you get? Total or Secondary?') - + #Legends - ax[3].plot(MP0[:,0],MP0[:,1],color='gray') - Dip_Midpoint0 = ax[3].scatter(MP0[:,0],MP0[:,1],color='black') - Electrodes0 = ax[3].scatter(EL0[:,0],EL0[:,1],color='red') - ax[3].legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) - - ax[4].plot(MP1[:,0],MP1[:,1],color='gray') - Dip_Midpoint1 = ax[4].scatter(MP1[:,0],MP1[:,1],color='black') - Electrodes1 = ax[4].scatter(EL1[:,0],EL1[:,1],color='red') - ax[4].legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) - - return ax + ax3.plot(MP0[:,0],MP0[:,1],color='gray') + Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black') + Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red') + ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) + + ax4.plot(MP1[:,0],MP1[:,1],color='gray') + Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black') + Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red') + ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1) + + return fig #Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b): - + sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1) E0 = 1. # inducing field strength in V/m n = 100 #level of discretisation @@ -650,45 +648,45 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig fig, ax = plt.subplots(1,2,figsize=(18,6)) - #Setup figure 1 with options Configuration, Total or Secondary, + #Setup figure 1 with options Configuration, Total or Secondary, #then Potential, ElectricField, Current Density or Charges Density if Figure1a == 'Configuration': ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6]) - + elif Figure1a == 'Total': - + if Figure1b == 'Potential': ax[0] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0]) elif Figure1b == 'ElectricField': ax[0] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1b == 'CurrentDensity': ax[0] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1b == 'ChargesDensity': ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1a == 'Secondary': - + if Figure1b == 'Potential': ax[0] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1b == 'ElectricField': ax[0] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1b == 'CurrentDensity': ax[0] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[0]) - + elif Figure1b == 'ChargesDensity': ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0]) - - + + if Figure1a== 'Configuration': ax[1] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) print 'While figure1 is plotting Configuration, figure2 plots the primary field' - - elif Figure2a == 'Total': + + elif Figure2a == 'Total': if Figure2b == 'Potential': ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1]) @@ -701,8 +699,8 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig elif Figure2b == 'ChargesDensity': ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1]) - - elif Figure2a == 'Secondary': + + elif Figure2a == 'Secondary': if Figure2b == 'Potential': ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1]) @@ -717,10 +715,10 @@ def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Fig plt.tight_layout(True) plt.show() - + #Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example): - + sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1) sig2 = 10.**log_sig2 E0 = 1. # inducing field strength in V/m @@ -731,32 +729,21 @@ def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,x XYZ = ndgrid(xr,yr,zr) # Space Definition PlotOpt = 'Total' - #Initializing the figure - fig = plt.figure(figsize=(20,20)) - ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6) #Configuration Conductive Sphere - ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6) #Configuration Resistive Sphere - ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4) # Data - ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6) #Potential Conductive Sphere - ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6) #Potential Resistive Potential - ax = [ax0,ax1,ax2,ax3,ax4] - if matching_spheres_example: - sig0 = 10.**(-3) - sig1 = 10.**(-2) + sig0 = 10.**(-3) + sig1 = 10.**(-2) sig2 = 1.310344828 * 10**(-3) - R0 = 20. + R0 = 20. R1 = 40. - two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt,ax) + two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt) else: - two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt,ax) + two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt) plt.tight_layout(True) plt.show() - - def run(plotIt=True): sig0 = -3. # conductivity of the wholespace sig1 = -1. # conductivity of the sphere @@ -769,11 +756,6 @@ def run(plotIt=True): zr = np.r_[0] # identical to saying `zr = np.array([0])` XYZ = ndgrid(xr,yr,zr) # Space Definition - Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0) - Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0) - Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es) - rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep) - if plotIt: fig, ax = plt.subplots(2,5,figsize=(50,10)) ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1]) @@ -786,11 +768,13 @@ def run(plotIt=True): ax[1,3] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1,3]) ax[0,4] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[0,4]) ax[1,4] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1,4]) + else: + get_Potential(XYZ,sig0,sig1,R,E0) # This is so travis tests it + + plt.show() - return Vt,Vp,Vs,Et,Ep,Es,Jt,Jp,Js,rho - - if __name__ == '__main__': - run() + + diff --git a/SimPEG/FLOW/Richards/RichardsProblem.py b/SimPEG/FLOW/Richards/RichardsProblem.py index 4dcabe60..2346f4da 100644 --- a/SimPEG/FLOW/Richards/RichardsProblem.py +++ b/SimPEG/FLOW/Richards/RichardsProblem.py @@ -45,19 +45,19 @@ class RichardsSurvey(Survey.BaseSurvey): @Utils.count @Utils.requires('prob') - def dpred(self, m, u=None): + def dpred(self, m, f=None): """ Create the projected data from a model. - The field, u, (if provided) will be used for the predicted data + The field, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!). .. math:: - d_\\text{pred} = P(u(m), m) + d_\\text{pred} = P(f(m), m) Where P is a projection of the fields onto the data space. """ - if u is None: u = self.prob.fields(m) - return Utils.mkvc(self.eval(u, m)) + if f is None: f = self.prob.fields(m) + return Utils.mkvc(self.eval(f, m)) @Utils.requires('prob') def eval(self, U, m): @@ -233,16 +233,16 @@ class RichardsProblem(Problem.BaseTimeProblem): return r, J @Utils.timeIt - def Jfull(self, m, u=None): - if u is None: - u = self.fields(m) + def Jfull(self, m, f=None): + if f is None: + f = self.fields(m) - nn = len(u)-1 + nn = len(f)-1 Asubs, Adiags, Bs = range(nn), range(nn), range(nn) for ii in range(nn): dt = self.timeSteps[ii] - bc = self.getBoundaryConditions(ii, u[ii]) - Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, u[ii], u[ii+1], dt, bc) + bc = self.getBoundaryConditions(ii, f[ii]) + Asubs[ii], Adiags[ii], Bs[ii] = self.diagsJacobian(m, f[ii], f[ii+1], dt, bc) Ad = sp.block_diag(Adiags) zRight = Utils.spzeros((len(Asubs)-1)*Asubs[0].shape[0],Adiags[0].shape[1]) zTop = Utils.spzeros(Adiags[0].shape[0], len(Adiags)*Adiags[0].shape[1]) @@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem): B = np.array(sp.vstack(Bs).todense()) Ainv = self.Solver(A, **self.solverOpts) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) AinvB = Ainv * B z = np.zeros((self.mesh.nC, B.shape[1])) zAinvB = np.vstack((z, AinvB)) @@ -259,41 +259,41 @@ class RichardsProblem(Problem.BaseTimeProblem): return J @Utils.timeIt - def Jvec(self, m, v, u=None): - if u is None: - u = self.fields(m) + def Jvec(self, m, v, f=None): + if f is None: + f = self.fields(m) - JvC = range(len(u)-1) # Cell to hold each row of the long vector. + JvC = range(len(f)-1) # Cell to hold each row of the long vector. # This is done via forward substitution. - bc = self.getBoundaryConditions(0, u[0]) - temp, Adiag, B = self.diagsJacobian(m, u[0], u[1], self.timeSteps[0], bc) + bc = self.getBoundaryConditions(0, f[0]) + temp, Adiag, B = self.diagsJacobian(m, f[0], f[1], self.timeSteps[0], bc) Adiaginv = self.Solver(Adiag, **self.solverOpts) JvC[0] = Adiaginv * (B*v) - for ii in range(1,len(u)-1): - bc = self.getBoundaryConditions(ii, u[ii]) - Asub, Adiag, B = self.diagsJacobian(m, u[ii], u[ii+1], self.timeSteps[ii], bc) + for ii in range(1,len(f)-1): + bc = self.getBoundaryConditions(ii, f[ii]) + Asub, Adiag, B = self.diagsJacobian(m, f[ii], f[ii+1], self.timeSteps[ii], bc) Adiaginv = self.Solver(Adiag, **self.solverOpts) JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1]) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC) @Utils.timeIt - def Jtvec(self, m, v, u=None): - if u is None: - u = self.field(m) + def Jtvec(self, m, v, f=None): + if f is None: + f = self.field(m) - P = self.survey.evalDeriv(u, m) + P = self.survey.evalDeriv(f, m) PTv = P.T*v # This is done via backward substitution. minus = 0 BJtv = 0 - for ii in range(len(u)-1,0,-1): - bc = self.getBoundaryConditions(ii-1, u[ii-1]) - Asub, Adiag, B = self.diagsJacobian(m, u[ii-1], u[ii], self.timeSteps[ii-1], bc) + for ii in range(len(f)-1,0,-1): + bc = self.getBoundaryConditions(ii-1, f[ii-1]) + Asub, Adiag, B = self.diagsJacobian(m, f[ii-1], f[ii], self.timeSteps[ii-1], bc) #select the correct part of v vpart = range((ii)*Adiag.shape[0], (ii+1)*Adiag.shape[0]) AdiaginvT = self.Solver(Adiag.T, **self.solverOpts) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index 0296bf4b..fd6c48c3 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -82,23 +82,23 @@ class BaseInvProblem(object): self._warmstart = value def getFields(self, m, store=False, deleteWarmstart=True): - u = None + f = None for mtest, u_ofmtest in self.warmstart: if m is mtest: - u = u_ofmtest + f = u_ofmtest if self.debug: print 'InvProb is Warm Starting!' break - if u is None: - u = self.prob.fields(m) + if f is None: + f = self.prob.fields(m) if deleteWarmstart: self.warmstart = [] if store: - self.warmstart += [(m,u)] + self.warmstart += [(m,f)] - return u + return f @Utils.timeIt def evalFunction(self, m, return_g=True, return_H=True): @@ -109,21 +109,21 @@ class BaseInvProblem(object): gc.collect() # Store fields if doing a line-search - u = self.getFields(m, store=(return_g==False and return_H==False)) + f = self.getFields(m, store=(return_g==False and return_H==False)) - phi_d = self.dmisfit.eval(m, u=u) + phi_d = self.dmisfit.eval(m, f=f) phi_m = self.reg.eval(m) - self.dpred = self.survey.dpred(m, u=u) # This is a cheap matrix vector calculation. + self.dpred = self.survey.dpred(m, f=f) # This is a cheap matrix vector calculation. self.phi_d, self.phi_d_last = phi_d, self.phi_d self.phi_m, self.phi_m_last = phi_m, self.phi_m - f = phi_d + self.beta * phi_m + phi = phi_d + self.beta * phi_m - out = (f,) + out = (phi,) if return_g: - phi_dDeriv = self.dmisfit.evalDeriv(m, u=u) + phi_dDeriv = self.dmisfit.evalDeriv(m, f=f) phi_mDeriv = self.reg.evalDeriv(m) g = phi_dDeriv + self.beta * phi_mDeriv @@ -131,7 +131,7 @@ class BaseInvProblem(object): if return_H: def H_fun(v): - phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, u=u) + phi_d2Deriv = self.dmisfit.eval2Deriv(m, v, f=f) phi_m2Deriv = self.reg.eval2Deriv(m, v=v) return phi_d2Deriv + self.beta * phi_m2Deriv diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py index 36389430..c201dfb0 100644 --- a/SimPEG/MT/BaseMT.py +++ b/SimPEG/MT/BaseMT.py @@ -27,7 +27,7 @@ class BaseMTProblem(BaseFDEMProblem): # Might need to add more stuff here. ## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components. - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): """ Function to calculate the data sensitivities dD/dm times a vector. @@ -39,8 +39,8 @@ class BaseMTProblem(BaseFDEMProblem): """ # Calculate the fields - if u is None: - u = self.fields(m) + if f is None: + f= self.fields(m) # Set current model self.curModel = m # Initiate the Jv object @@ -56,9 +56,9 @@ class BaseMTProblem(BaseFDEMProblem): # We need fDeriv_m = df/du*du/dm + df/dm # Construct du/dm, it requires a solve # NOTE: need to account for the 2 polarizations in the derivatives. - u_src = u[src,:] + f_src = f[src,:] # dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations. - dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns. + dA_dm = self.getADeriv_m(freq, f_src, v) # Size: nE,2 (u_px,u_py) in the columns. dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns. if dRHS_dm is None: du_dm = dA_duI * ( -dA_dm ) @@ -68,13 +68,13 @@ class BaseMTProblem(BaseFDEMProblem): for rx in src.rxList: # Get the projection derivative # v should be of size 2*nE (for 2 polarizations) - PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m + PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, f, t) # wrt u, we don't have have PDeriv wrt m Jv[src, rx] = PDeriv_u(mkvc(du_dm)) dA_duI.clean() # Return the vectorized sensitivities return mkvc(Jv) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): """ Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector. @@ -85,8 +85,8 @@ class BaseMTProblem(BaseFDEMProblem): :return: Data sensitivities wrt m """ - if u is None: - u = self.fields(m) + if f is None: + f = self.fields(m) self.curModel = m @@ -103,15 +103,15 @@ class BaseMTProblem(BaseFDEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = u[src, :] + f_src = f[src, :] for rx in src.rxList: # Get the adjoint evalDeriv # PTv needs to be nE, - PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m + PTv = rx.evalDeriv(src, self.mesh, f, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m # Get the dA_duIT = ATinv * PTv - dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True) + dA_dmT = self.getADeriv_m(freq, f_src, mkvc(dA_duIT), adjoint=True) dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True) # Make du_dmT if dRHS_dmT is None: @@ -129,4 +129,4 @@ class BaseMTProblem(BaseFDEMProblem): raise Exception('Must be real or imag') # Clean the factorization, clear memory. ATinv.clean() - return Jtv \ No newline at end of file + return Jtv diff --git a/SimPEG/MT/SurveyMT.py b/SimPEG/MT/SurveyMT.py index 4e4a8688..0ec91a0e 100644 --- a/SimPEG/MT/SurveyMT.py +++ b/SimPEG/MT/SurveyMT.py @@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey): assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] - def eval(self, u): + def eval(self, f): data = Data(self) for src in self.srcList: sys.stdout.flush() for rx in src.rxList: - data[src, rx] = rx.eval(src, self.mesh, u) + data[src, rx] = rx.eval(src, self.mesh, f) return data - def evalDeriv(self, u): + def evalDeriv(self, f): raise Exception('Use Transmitters to project fields deriv.') ################# diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 4f2cb062..0a241710 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -888,6 +888,8 @@ class ProjectedGNCG(BFGS, Minimize, Remember): maxIterCG = 5 tolCG = 1e-1 + stepOffBoundsFact = 0.1 # perturbation of the inactive set off the bounds + lower = -np.inf upper = np.inf @@ -990,4 +992,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember): cgFlag = 1 # End CG Iterations + # Take a gradient step on the active cells if exist + if temp != self.xc.size: + + rhs_a = (Active) * -self.g + + dm_i = max( abs( delx ) ) + dm_a = max( abs(rhs_a) ) + + # perturb inactive set off of bounds so that they are included in the step + delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a) + + # Only keep gradients going in the right direction on the active set + indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0)) + delx[indx] = 0. + return delx diff --git a/SimPEG/Problem.py b/SimPEG/Problem.py index cd8a4aaa..f8520c5c 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -88,28 +88,28 @@ class BaseProblem(object): return self.survey is not None @Utils.timeIt - def Jvec(self, m, v, u=None): - """Jvec(m, v, u=None) + def Jvec(self, m, v, f=None): + """Jvec(m, v, f=None) Effect of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: Jv """ raise NotImplementedError('J is not yet implemented.') @Utils.timeIt - def Jtvec(self, m, v, u=None): - """Jtvec(m, v, u=None) + def Jtvec(self, m, v, f=None): + """Jtvec(m, v, f=None) Effect of transpose of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: JTv """ @@ -117,32 +117,32 @@ class BaseProblem(object): @Utils.timeIt - def Jvec_approx(self, m, v, u=None): - """Jvec_approx(m, v, u=None) + def Jvec_approx(self, m, v, f=None): + """Jvec_approx(m, v, f=None) Approximate effect of J(m) on a vector v :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: approxJv """ - return self.Jvec(m, v, u) + return self.Jvec(m, v, f) @Utils.timeIt - def Jtvec_approx(self, m, v, u=None): - """Jtvec_approx(m, v, u=None) + def Jtvec_approx(self, m, v, f=None): + """Jtvec_approx(m, v, f=None) Approximate effect of transpose of J(m) on a vector v. :param numpy.array m: model :param numpy.array v: vector to multiply - :param numpy.array u: fields + :param Fields f: fields :rtype: numpy.array :return: JTv """ - return self.Jtvec(m, v, u) + return self.Jtvec(m, v, f) def fields(self, m): """ @@ -224,9 +224,9 @@ class LinearProblem(BaseProblem): def fields(self, m): return self.G.dot(m) - def Jvec(self, m, v, u=None): + def Jvec(self, m, v, f=None): return self.G.dot(v) - def Jtvec(self, m, v, u=None): + def Jtvec(self, m, v, f=None): return self.G.T.dot(v) diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 37024028..fbc88276 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -295,38 +295,38 @@ class BaseSurvey(object): @Utils.count @Utils.requires('prob') - def dpred(self, m, u=None): - """dpred(m, u=None) + def dpred(self, m, f=None): + """dpred(m, f=None) Create the projected data from a model. - The field, u, (if provided) will be used for the predicted data + The fields, f, (if provided) will be used for the predicted data instead of recalculating the fields (which may be expensive!). .. math:: - d_\\text{pred} = P(u(m)) + d_\\text{pred} = P(f(m)) Where P is a projection of the fields onto the data space. """ - if u is None: u = self.prob.fields(m) - return Utils.mkvc(self.eval(u)) + if f is None: f = self.prob.fields(m) + return Utils.mkvc(self.eval(f)) @Utils.count - def eval(self, u): - """eval(u) + def eval(self, f): + """eval(f) This function projects the fields onto the data space. .. math:: - d_\\text{pred} = \mathbf{P} u(m) + d_\\text{pred} = \mathbf{P} f(m) """ raise NotImplemented('eval is not yet implemented.') @Utils.count - def evalDeriv(self, u): - """evalDeriv(u) + def evalDeriv(self, f): + """evalDeriv(f) This function s the derivative of projects the fields onto the data space. @@ -337,11 +337,11 @@ class BaseSurvey(object): raise NotImplemented('eval is not yet implemented.') @Utils.count - def residual(self, m, u=None): - """residual(m, u=None) + def residual(self, m, f=None): + """residual(m, f=None) :param numpy.array m: geophysical model - :param numpy.array u: fields + :param numpy.array f: fields :rtype: numpy.array :return: data residual @@ -352,14 +352,14 @@ class BaseSurvey(object): \mu_\\text{data} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs} """ - return Utils.mkvc(self.dpred(m, u=u) - self.dobs) + return Utils.mkvc(self.dpred(m, f=f) - self.dobs) @property def isSynthetic(self): "Check if the data is synthetic." return self.mtrue is not None - def makeSyntheticData(self, m, std=0.05, u=None, force=False): + def makeSyntheticData(self, m, std=0.05, f=None, force=False): """ Make synthetic data given a model, and a standard deviation. @@ -372,16 +372,16 @@ class BaseSurvey(object): if getattr(self, 'dobs', None) is not None and not force: raise Exception('Survey already has dobs. You can use force=True to override this exception.') self.mtrue = m - self.dtrue = self.dpred(m, u=u) + self.dtrue = self.dpred(m, f=f) noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape) self.dobs = self.dtrue+noise self.std = self.dobs*0 + std return self.dobs class LinearSurvey(BaseSurvey): - def eval(self, u): - return u - + def eval(self, f): + return f + @property def nD(self): return self.prob.G.shape[0] diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 18c1994f..3f5d62c8 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,3 +7,4 @@ from CounterUtils import * import ModelBuilder import SolverUtils from coordutils import * +from plottingUtils import * diff --git a/SimPEG/Utils/plottingUtils.py b/SimPEG/Utils/plottingUtils.py new file mode 100644 index 00000000..72ab7654 --- /dev/null +++ b/SimPEG/Utils/plottingUtils.py @@ -0,0 +1,3 @@ +# Plot Tree! +# Plot SphereSetup +# Plot LayerEarth diff --git a/docs/examples/EM_Schenkel_Morrison_Casing.rst b/docs/examples/EM_Schenkel_Morrison_Casing.rst new file mode 100644 index 00000000..55f00168 --- /dev/null +++ b/docs/examples/EM_Schenkel_Morrison_Casing.rst @@ -0,0 +1,58 @@ +.. _examples_EM_Schenkel_Morrison_Casing: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: Schenkel and Morrison Casing Model +====================================== + +Here we create and run a FDEM forward simulation to calculate the vertical +current inside a steel-cased. The model is based on the Schenkel and +Morrison Casing Model, and the results are used in a 2016 SEG abstract by +Yang et al. + +- Schenkel, C.J., and H.F. Morrison, 1990, Effects of well casing on potential field measurements using downhole current sources: Geophysical prospecting, 38, 663-686. + + +The model consists of: +- Air: Conductivity 1e-8 S/m, above z = 0 +- Background: conductivity 1e-2 S/m, below z = 0 +- Casing: conductivity 1e6 S/m + - 300m long + - radius of 0.1m + - thickness of 6e-3m + +Inside the casing, we take the same conductivity as the background. + +We are using an EM code to simulate DC, so we use frequency low enough +that the skin depth inside the casing is longer than the casing length (f += 1e-6 Hz). The plot produced is of the current inside the casing. + +These results are shown in the SEG abstract by Yang et al., 2016: 3D DC +resistivity modeling of steel casing for reservoir monitoring using +equivalent resistor network. The solver used to produce these results and +achieve the CPU time of ~30s is Mumps, which was installed using pymatsolver_ + +.. _pymatsolver: https://github.com/rowanc1/pymatsolver + +This example is on figshare: https://dx.doi.org/10.6084/m9.figshare.3126961.v1 + +If you would use this example for a code comparison, or build upon it, a +citation would be much appreciated! + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_Schenkel_Morrison_Casing.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_Schenkel_Morrison_Casing.py + :language: python + :linenos: diff --git a/tests/flow/test_Richards.py b/tests/flow/test_Richards.py index d63a6210..f67ec71d 100644 --- a/tests/flow/test_Richards.py +++ b/tests/flow/test_Richards.py @@ -116,8 +116,8 @@ class RichardsTests1D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print 'Richards Adjoint Test - PressureHead' @@ -188,8 +188,8 @@ class RichardsTests2D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print '2D: Richards Adjoint Test - PressureHead' @@ -260,8 +260,8 @@ class RichardsTests3D(unittest.TestCase): v = np.random.rand(self.survey.nD) z = np.random.rand(self.M.nC) Hs = self.prob.fields(self.Ks) - vJz = v.dot(self.prob.Jvec(self.Ks,z,u=Hs)) - zJv = z.dot(self.prob.Jtvec(self.Ks,v,u=Hs)) + vJz = v.dot(self.prob.Jvec(self.Ks,z,f=Hs)) + zJv = z.dot(self.prob.Jtvec(self.Ks,v,f=Hs)) tol = TOL*(10**int(np.log10(np.abs(zJv)))) passed = np.abs(vJz - zJv) < tol print '3D: Richards Adjoint Test - PressureHead'