diff --git a/.bumpversion.cfg b/.bumpversion.cfg index ea37cced..84469e3f 100644 --- a/.bumpversion.cfg +++ b/.bumpversion.cfg @@ -1,4 +1,4 @@ [bumpversion] -current_version = 0.1.9 +current_version = 0.1.10 files = setup.py SimPEG/__init__.py docs/conf.py 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 e5a63547..542525d5 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 3e378a6a..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,88 +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) - ftype = self._fieldType + 'Solution' - F[Srcs, ftype] = 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): - ftype = self._fieldType + 'Solution' - u_src = u[src, ftype] - dA_dm = self.getADeriv_m(freq, u_src, v) - dRHS_dm = self.getRHSDeriv_m(freq, src, v) - du_dm = Ainv * ( - dA_dm + dRHS_dm ) - + u_src = f[src, self._solutionType] + dA_dm_v = self.getADeriv(freq, u_src, v) + dRHS_dm_v = self.getRHSDeriv(freq, src, v) + du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) + for rx in src.rxList: - df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) - df_dudu_dm = df_duFun(src, du_dm, adjoint=False) - - df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) - df_dm = df_dmFun(src, v, adjoint=False) - - - Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex) - - P = lambda v: rx.evalDeriv(src, self.mesh, u, v) # wrt u, also have wrt m - - Jv[src, rx] = P(Df_Dm) - + 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, 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 @@ -132,35 +120,31 @@ class BaseFDEMProblem(BaseEMProblem): ATinv = self.Solver(AT, **self.solverOpts) for src in self.survey.getSrcByFreq(freq): - ftype = self._fieldType + 'Solution' - u_src = u[src, ftype] + 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(f, '_%sDeriv'%rx.projField, None) + df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) - df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None) - df_duT = df_duTFun(src, PTv, adjoint=True) - ATinvdf_duT = ATinv * df_duT - dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True) - dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True) + dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True) + dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) - dfT_dm = df_dmFun(src, PTv, adjoint=True) + df_dmT = df_dmT + du_dmT - du_dmT += dfT_dm - - # TODO: this should be taken care of by the reciever + # TODO: this should be taken care of by the reciever? real_or_imag = rx.projComp if real_or_imag is 'real': - Jtv += np.array(du_dmT,dtype=complex).real + Jtv += np.array(df_dmT, dtype=complex).real elif real_or_imag is 'imag': - Jtv += - np.array(du_dmT,dtype=complex).real + Jtv += - np.array(df_dmT, dtype=complex).real else: raise Exception('Must be real or imag') - + ATinv.clean() return Utils.mkvc(Jtv) @@ -170,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._eqLocs is 'FE': - S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) - S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) - elif self._eqLocs is 'EF': - S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex) - S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex) + 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) + 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) 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 ########################################################################################## @@ -213,9 +197,9 @@ class Problem_e(BaseFDEMProblem): :param SimPEG.Mesh mesh: mesh """ - _fieldType = 'e' - _eqLocs = 'FE' - fieldsPair = Fields_e + _solutionType = 'eSolution' + _formulation = 'EB' + fieldsPair = Fields_e def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -223,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}} @@ -239,19 +223,19 @@ class Problem_e(BaseFDEMProblem): return C.T*MfMui*C + 1j*omega(freq)*MeSigma - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector .. 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 @@ -264,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_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -294,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): @@ -324,9 +308,9 @@ class Problem_b(BaseFDEMProblem): :param SimPEG.Mesh mesh: mesh """ - _fieldType = 'b' - _eqLocs = 'FE' - fieldsPair = Fields_b + _solutionType = 'bSolution' + _formulation = 'EB' + fieldsPair = Fields_b def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -354,7 +338,7 @@ class Problem_b(BaseFDEMProblem): return MfMui.T*A return A - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -362,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 @@ -389,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 @@ -411,7 +395,7 @@ class Problem_b(BaseFDEMProblem): return RHS - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ Derivative of the right hand side with respect to the model @@ -424,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) @@ -472,9 +456,9 @@ class Problem_j(BaseFDEMProblem): :param SimPEG.Mesh mesh: mesh """ - _fieldType = 'j' - _eqLocs = 'EF' - fieldsPair = Fields_j + _solutionType = 'jSolution' + _formulation = 'HJ' + fieldsPair = Fields_j def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -503,7 +487,7 @@ class Problem_j(BaseFDEMProblem): return A - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector @@ -513,32 +497,32 @@ 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 MfRho = self.MfRho C = self.mesh.edgeCurl - MfRhoDeriv_m = self.MfRhoDeriv(u) + MfRhoDeriv = self.MfRhoDeriv(u) if adjoint: if self._makeASymmetric is True: v = MfRho * v - return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v))) + return MfRhoDeriv.T * (C * (MeMuI.T * (C.T * v))) if self._makeASymmetric is True: - return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) ))) - return C * (MeMuI * (C.T * (MfRhoDeriv_m * v))) + return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv * v) ))) + return C * (MeMuI * (C.T * (MfRhoDeriv * v))) def getRHS(self, freq): """ - Right hand side for the system + Right hand side for the system .. math :: @@ -549,20 +533,20 @@ 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 return RHS - def getRHSDeriv_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -574,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 @@ -610,9 +594,9 @@ class Problem_h(BaseFDEMProblem): :param SimPEG.Mesh mesh: mesh """ - _fieldType = 'h' - _eqLocs = 'EF' - fieldsPair = Fields_h + _solutionType = 'hSolution' + _formulation = 'HJ' + fieldsPair = Fields_h def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -635,51 +619,51 @@ class Problem_h(BaseFDEMProblem): return C.T * (MfRho * C) + 1j*omega(freq)*MeMu - def getADeriv_m(self, freq, u, v, adjoint=False): + def getADeriv(self, freq, u, v, adjoint=False): """ Product of the derivative of our system matrix with respect to the model and a vector .. 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 C = self.mesh.edgeCurl - MfRhoDeriv_m = self.MfRhoDeriv(C*u) + MfRhoDeriv = self.MfRhoDeriv(C*u) if adjoint: - return MfRhoDeriv_m.T * (C * v) - return C.T * (MfRhoDeriv_m * v) + return MfRhoDeriv.T * (C * v) + return C.T * (MfRhoDeriv * v) def getRHS(self, freq): """ - 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_m(self, freq, src, v, adjoint=False): + def getRHSDeriv(self, freq, src, v, adjoint=False): """ - Derivative of the right hand side with respect to the model + Derivative of the right hand side with respect to the model :param float freq: frequency :param SimPEG.EM.FDEM.Src src: FDEM source @@ -689,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 e171a5c5..e2193973 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -3,12 +3,12 @@ import scipy.sparse as sp import SimPEG from SimPEG import Utils from SimPEG.EM.Utils import omega -from SimPEG.Utils import Zero, Identity +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 @@ -32,12 +32,140 @@ class Fields(SimPEG.Problem.Fields): knownFields = {} dtype = complex + def _e(self, solution, srcList): + """ + Total electric field is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total electric field + """ + if getattr(self, '_ePrimary', None) is None or getattr(self, '_eSecondary', None) is None: + raise NotImplementedError ('Getting e from %s is not implemented' %self.knownFields.keys()[0]) + + return self._ePrimary(solution,srcList) + self._eSecondary(solution,srcList) + + def _b(self, solution, srcList): + """ + Total magnetic flux density is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total magnetic flux density + """ + if getattr(self, '_bPrimary', None) is None or getattr(self, '_bSecondary', None) is None: + raise NotImplementedError ('Getting b from %s is not implemented' %self.knownFields.keys()[0]) + + return self._bPrimary(solution, srcList) + self._bSecondary(solution, srcList) + + def _h(self, solution, srcList): + """ + Total magnetic field is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total magnetic field + """ + if getattr(self, '_hPrimary', None) is None or getattr(self, '_hSecondary', None) is None: + raise NotImplementedError ('Getting h from %s is not implemented' %self.knownFields.keys()[0]) + + return self._hPrimary(solution, srcList) + self._hSecondary(solution, srcList) + + def _j(self, solution, srcList): + """ + Total current density is sum of primary and secondary + + :param numpy.ndarray solution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: total current density + """ + if getattr(self, '_jPrimary', None) is None or getattr(self, '_jSecondary', None) is None: + raise NotImplementedError ('Getting j from %s is not implemented' %self.knownFields.keys()[0]) + + return self._jPrimary(solution, srcList) + self._jSecondary(solution, srcList) + + def _eDeriv(self, src, du_dm_v, v, adjoint = False): + """ + Total derivative of e with respect to the inversion model. Returns :math:`d\mathbf{e}/d\mathbf{m}` for forward and (:math:`d\mathbf{e}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ + if getattr(self, '_eDeriv_u', None) is None or getattr(self, '_eDeriv_m', None) is None: + raise NotImplementedError ('Getting eDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._eDeriv_u(src, v, adjoint), self._eDeriv_m(src, v, adjoint) + return np.array(self._eDeriv_u(src, du_dm_v, adjoint) + self._eDeriv_m(src, v, adjoint), dtype = complex) + + def _bDeriv(self, src, du_dm_v, v, adjoint = False): + """ + Total derivative of b with respect to the inversion model. Returns :math:`d\mathbf{b}/d\mathbf{m}` for forward and (:math:`d\mathbf{b}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ + if getattr(self, '_bDeriv_u', None) is None or getattr(self, '_bDeriv_m', None) is None: + raise NotImplementedError ('Getting bDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._bDeriv_u(src, v, adjoint), self._bDeriv_m(src, v, adjoint) + return 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): + """ + Total derivative of h with respect to the inversion model. Returns :math:`d\mathbf{h}/d\mathbf{m}` for forward and (:math:`d\mathbf{h}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ + if getattr(self, '_hDeriv_u', None) is None or getattr(self, '_hDeriv_m', None) is None: + raise NotImplementedError ('Getting hDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._hDeriv_u(src, v, adjoint), self._hDeriv_m(src, v, adjoint) + return np.array(self._hDeriv_u(src, du_dm_v, adjoint) + self._hDeriv_m(src, v, adjoint), dtype = complex) + + def _jDeriv(self, src, du_dm_v, v, adjoint = False): + """ + Total derivative of j with respect to the inversion model. Returns :math:`d\mathbf{j}/d\mathbf{m}` for forward and (:math:`d\mathbf{j}/d\mathbf{u}`, :math:`d\mathb{u}/d\mathbf{m}`) for the adjoint + + :param Src src: sorce + :param numpy.ndarray du_dm_v: derivative of the solution vector with respect to the model times a vector (is None for adjoint) + :param numpy.ndarray v: vector to take sensitivity product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: derivative times a vector (or tuple for adjoint) + """ + if getattr(self, '_jDeriv_u', None) is None or getattr(self, '_jDeriv_m', None) is None: + raise NotImplementedError ('Getting jDerivs from %s is not implemented' %self.knownFields.keys()[0]) + + if adjoint: + return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) + return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex) + 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'} @@ -47,15 +175,34 @@ class Fields_e(Fields): 'eSecondary' : ['eSolution','E','_eSecondary'], 'b' : ['eSolution','F','_b'], 'bPrimary' : ['eSolution','F','_bPrimary'], - 'bSecondary' : ['eSolution','F','_bSecondary'] + 'bSecondary' : ['eSolution','F','_bSecondary'], + 'j' : ['eSolution','CCV','_j'], + 'h' : ['eSolution','CCV','_h'], } - def __init__(self,mesh,survey,**kwargs): + def __init__(self, mesh, survey, **kwargs): Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._nC = self.survey.prob.mesh.nC + self._MeSigma = self.survey.prob.MeSigma + self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv + self._MfMui = self.survey.prob.MfMui + + def _GLoc(self, fieldType): + if fieldType == 'e': + return 'E' + elif fieldType == 'b': + return 'F' + elif (fieldType == 'h') or (fieldType == 'j'): + return 'CCV' + else: + raise Exception('Field type must be e, b, h, j') + def _ePrimary(self, eSolution, srcList): """ @@ -67,7 +214,7 @@ class Fields_e(Fields): :return: primary electric field as defined by the sources """ - ePrimary = np.zeros_like(eSolution) + ePrimary = np.zeros([self.prob.mesh.nE,len(srcList)], dtype = complex) for i, src in enumerate(srcList): ep = src.ePrimary(self.prob) ePrimary[:,i] = ePrimary[:,i] + ep @@ -82,26 +229,13 @@ class Fields_e(Fields): :rtype: numpy.ndarray :return: secondary electric field """ - return eSolution - def _e(self, eSolution, srcList): - """ - Total electric field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total electric field - """ - - return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList) - def _eDeriv_u(self, src, v, adjoint = False): """ - Derivative of the total electric field with respect to the thing we - solved for - + Partial derivative of the total electric field with respect to the thing we + solved for. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -113,8 +247,8 @@ class Fields_e(Fields): def _eDeriv_m(self, src, v, adjoint = False): """ - Derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. - + Partial derivative of the total electric field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -135,7 +269,7 @@ class Fields_e(Fields): :return: primary magnetic flux density as defined by the sources """ - bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) + bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]], dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) bPrimary[:,i] = bPrimary[:,i] + bp @@ -155,87 +289,149 @@ 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 _bSecondaryDeriv_u(self, src, v, adjoint = False): + def _bDeriv_u(self, src, du_dm_v, adjoint = False): """ - Derivative of the secondary magnetic flux density with respect to the thing we solved for - + Derivative of the magnetic flux density with respect to the thing we solved for + :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic flux density with respect to the field we solved for with a vector - """ - - C = self._edgeCurl - if adjoint: - return - 1./(1j*omega(src.freq)) * (C.T * v) - return - 1./(1j*omega(src.freq)) * (C * v) - - def _bSecondaryDeriv_m(self, src, v, adjoint = False): - """ - Derivative of the secondary 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? - :rtype: numpy.ndarray - :return: product of the secondary 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 - - def _b(self, eSolution, srcList): - """ - Total magnetic flux density is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic flux density - """ - - return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList) - - def _bDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total magnetic flux density with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm_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 """ - # Primary does not depend on u - return self._bSecondaryDeriv_u(src, v, adjoint) + C = self._edgeCurl + if adjoint: + return - 1./(1j*omega(src.freq)) * (C.T * du_dm_v) + return - 1./(1j*omega(src.freq)) * (C * du_dm_v) - def _bDeriv_m(self, src, v, adjoint=False): + + def _bDeriv_m(self, src, v, adjoint = False): """ - Derivative of the total 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? - :rtype: SimPEG.Utils.Zero + :rtype: numpy.ndarray :return: product of the magnetic flux density derivative with respect to the inversion model with a vector """ - # Assuming the primary does not depend on the model - return self._bSecondaryDeriv_m(src, v, adjoint) + s_mDeriv, _ = src.evalDeriv(self.prob, v, adjoint) + return 1./(1j * omega(src.freq)) * s_mDeriv + + def _j(self, eSolution, srcList): + """ + Current density from eSolution + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: current density + """ + aveE2CCV = self._aveE2CCV + n = int(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)) + return VI * (aveE2CCV * (self._MeSigma * self._e(eSolution,srcList) ) ) + + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the current density with respect to the field we solved for with a vector + """ + 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: + 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. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the current density derivative with respect to the inversion model with a vector + """ + e = self[src, 'e'] + 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._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): + """ + Magnetic field from eSolution + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: magnetic field + """ + 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._MfMui * self._b(eSolution, srcList))) + + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector + """ + 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: + v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * du_dm_v)) + return self._bDeriv_u(src, v, adjoint=adjoint) + return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_u(src, du_dm_v, adjoint = adjoint))) + + def _hDeriv_m(self, src, v, adjoint = False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the magnetic field derivative with respect to the inversion model with a vector + """ + 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: + v = self._MfMui.T * (self._aveF2CCV.T * (VI.T * v)) + return self._bDeriv_m(src, v, adjoint=adjoint) + return VI * (self._aveF2CCV * (self._MfMui * self._bDeriv_m(src, v, adjoint = adjoint))) + 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'} @@ -246,6 +442,8 @@ class Fields_b(Fields): 'e' : ['bSolution','E','_e'], 'ePrimary' : ['bSolution','E','_ePrimary'], 'eSecondary' : ['bSolution','E','_eSecondary'], + 'j' : ['bSolution','CCV','_j'], + 'h' : ['bSolution','CCV','_h'], } def __init__(self,mesh,survey,**kwargs): @@ -254,10 +452,29 @@ class Fields_b(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeSigma = self.survey.prob.MeSigma self._MeSigmaI = self.survey.prob.MeSigmaI self._MfMui = self.survey.prob.MfMui + self._MeSigmaDeriv = self.survey.prob.MeSigmaDeriv self._MeSigmaIDeriv = self.survey.prob.MeSigmaIDeriv self._Me = self.survey.prob.Me + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._sigma = self.survey.prob.curModel.sigma + self._mui = self.survey.prob.curModel.mui + self._nC = self.survey.prob.mesh.nC + + + + def _GLoc(self,fieldType): + if fieldType == 'e': + return 'E' + elif fieldType == 'b': + return 'F' + elif (fieldType == 'h') or (fieldType == 'j'): + return'CCV' + else: + raise Exception('Field type must be e, b, h, j') def _bPrimary(self, bSolution, srcList): """ @@ -269,7 +486,7 @@ class Fields_b(Fields): :return: primary electric field as defined by the sources """ - bPrimary = np.zeros_like(bSolution) + bPrimary = np.zeros([self.prob.mesh.nF,len(srcList)], dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) bPrimary[:,i] = bPrimary[:,i] + bp @@ -287,35 +504,24 @@ class Fields_b(Fields): return bSolution - def _b(self, bSolution, srcList): - """ - Total magnetic flux density is sum of primary and secondary - - :param numpy.ndarray bSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic flux density + def _bDeriv_u(self, src, du_dm_v, adjoint=False): """ + Partial derivative of the total magnetic flux density with respect to the thing we + solved for. - return self._bPrimary(bSolution, srcList) + self._bSecondary(bSolution, srcList) - - def _bDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total magnetic flux density with respect to the thing we - solved for - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic flux density with respect to the field we solved for with a vector """ - return Identity()*v + + return Identity()*du_dm_v def _bDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. - + Partial derivative of the total magnetic flux density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -336,7 +542,7 @@ class Fields_b(Fields): :return: primary electric field as defined by the sources """ - ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) + ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]], dtype = complex) for i,src in enumerate(srcList): ep = src.ePrimary(self.prob) ePrimary[:,i] = ePrimary[:,i] + ep @@ -352,76 +558,17 @@ class Fields_b(Fields): :return: secondary electric field """ - e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) + e = ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): - _,S_e = src.eval(self.prob) - e[:,i] = e[:,i]+ -self._MeSigmaI * S_e - return e + _,s_e = src.eval(self.prob) + e[:,i] = e[:,i] + - s_e - def _eSecondaryDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the secondary 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? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary electric field with respect to the field we solved for with a vector + 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 - if not adjoint: - return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * v) ) - else: - return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * v)) - - def _eSecondaryDeriv_m(self, src, v, adjoint=False): - """ - Derivative of the secondary 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? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary electric field with respect to the model with a vector - """ - - bSolution = self[[src],'bSolution'] - _,S_e = src.eval(self.prob) - Me = self._Me - - if adjoint: - Me = Me.T - - w = self._edgeCurl.T * (self._MfMui * bSolution) - w = w - Utils.mkvc(Me * S_e,2) - - if not adjoint: - de_dm = self._MeSigmaIDeriv(w) * v - elif adjoint: - de_dm = self._MeSigmaIDeriv(w).T * v - - _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - - de_dm = de_dm - self._MeSigmaI * S_eDeriv - - return de_dm - - def _e(self, bSolution, srcList): - """ - Total electric field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total electric field - """ - - return self._ePrimary(bSolution, srcList) + self._eSecondary(bSolution, srcList) - - def _eDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total electric field with respect to the thing we solved for - :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -429,29 +576,129 @@ class Fields_b(Fields): :return: product of the derivative of the electric field with respect to the field we solved for with a vector """ - return self._eSecondaryDeriv_u(src, v, adjoint) + if not adjoint: + return self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * du_dm_v) ) + return self._MfMui.T * (self._edgeCurl * (self._MeSigmaI.T * du_dm_v)) + def _eDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total electric field density 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? :rtype: numpy.ndarray - :return: product of the electric field derivative with respect to the inversion model with a vector + :return: product of the derivative of the electric field with respect to the model with a vector """ - # assuming primary doesn't depend on model - return self._eSecondaryDeriv_m(src, v, adjoint) + bSolution = Utils.mkvc(self[src, 'bSolution']) + _,s_e = src.eval(self.prob) + + 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 + + def _j(self, bSolution, srcList): + """ + Secondary current density from bSolution + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current 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._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 + 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 current density with respect to the field we solved for with a vector + """ + 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._MfMui.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) + return VI * (self._aveE2CCV * (self._edgeCurl.T * ( self._MfMui * du_dm_v ) ) ) + + + def _jDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the current density with respect to the model with a vector + """ + return Zero() + + def _h(self, bSolution, srcList): + """ + Magnetic field from bSolution + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: magnetic field + """ + 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._MfMui * self._b(bSolution, srcList))) + + def _hDeriv_u(self, src, du_dm_v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector + """ + 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: + return self._MfMui.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v) ) + return VI * (self._aveF2CCV * (self._MfMui * du_dm_v)) + + def _hDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic field with respect to the model with a vector + """ + return Zero() 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'} @@ -462,6 +709,8 @@ class Fields_j(Fields): 'h' : ['jSolution','E','_h'], 'hPrimary' : ['jSolution','E','_hPrimary'], 'hSecondary' : ['jSolution','E','_hSecondary'], + 'e' : ['jSolution','CCV','_e'], + 'b' : ['jSolution','CCV','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -470,10 +719,25 @@ class Fields_j(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMu = self.survey.prob.MeMu self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho self._MfRhoDeriv = self.survey.prob.MfRhoDeriv - self._Me = self.survey.prob.Me + self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): + if fieldType == 'h': + return 'E' + elif fieldType == 'j': + return 'F' + elif (fieldType == 'e') or (fieldType == 'b'): + return 'CCV' + else: + raise Exception('Field type must be e, b, h, j') def _jPrimary(self, jSolution, srcList): """ @@ -485,7 +749,7 @@ class Fields_j(Fields): :return: primary current density as defined by the sources """ - jPrimary = np.zeros_like(jSolution,dtype = complex) + jPrimary = np.zeros_like(jSolution, dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) jPrimary[:,i] = jPrimary[:,i] + jp @@ -505,21 +769,22 @@ 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) - def _jDeriv_u(self, src, v, adjoint=False): + + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ - Derivative of the total current density with respect to the thing we - solved for - + Partial derivative of the total current density with respect to the thing we + solved for. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -527,12 +792,13 @@ class Fields_j(Fields): :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return Identity()*v + return Identity()*du_dm_v + def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. - + Partial derivative of the total current density with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -568,109 +834,166 @@ class Fields_j(Fields): :return: secondary magnetic field """ - h = self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) ) + 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)) * self._MeMuI * (S_m) - return h + s_m,_ = src.eval(self.prob) + h[:,i] = h[:,i] + 1./(1j*omega(src.freq)) * (s_m) + return self._MeMuI * h - def _hSecondaryDeriv_u(self, src, v, adjoint=False): + + def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ - Derivative of the secondary magnetic field with respect to the thing we solved for - + Derivative of the magnetic field with respect to the thing we solved for + :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with - :param bool adjoint: adjoint? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic field with respect to the field we solved for with a vector - """ - - if not adjoint: - return -1./(1j*omega(src.freq)) * self._MeMuI * (self._edgeCurl.T * (self._MfRho * v) ) - elif adjoint: - return -1./(1j*omega(src.freq)) * self._MfRho.T * (self._edgeCurl * ( self._MeMuI.T * v)) - - def _hSecondaryDeriv_m(self, src, v, adjoint=False): - """ - Derivative of the secondary 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? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary magnetic field with respect to the model with a vector - """ - - jSolution = self[[src],'jSolution'] - MeMuI = self._MeMuI - C = self._edgeCurl - MfRho = self._MfRho - MfRhoDeriv = self._MfRhoDeriv - Me = self._Me - - if not adjoint: - hDeriv_m = -1./(1j*omega(src.freq)) * MeMuI * (C.T * (MfRhoDeriv(jSolution)*v ) ) - elif adjoint: - hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) - - if not adjoint: - S_mDeriv = S_mDeriv(v) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * MeMuI * (Me * S_mDeriv) - elif adjoint: - S_mDeriv = S_mDeriv(Me.T * (MeMuI.T * v)) - hDeriv_m = hDeriv_m + 1./(1j*omega(src.freq)) * S_mDeriv - return hDeriv_m - - - def _h(self, jSolution, srcList): - """ - Total magnetic field is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic field - """ - - return self._hPrimary(jSolution, srcList) + self._hSecondary(jSolution, srcList) - - def _hDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total magnetic field with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return self._hSecondaryDeriv_u(src, v, adjoint) + 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 total magnetic field density 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? :rtype: numpy.ndarray - :return: product of the magnetic field derivative with respect to the inversion model with a vector + :return: product of the derivative of the magnetic field with respect to the model with a vector """ - # assuming the primary doesn't depend on the model - return self._hSecondaryDeriv_m(src, v, adjoint) + jSolution = Utils.mkvc(self[[src],'jSolution']) + MeMuI = self._MeMuI + C = self._edgeCurl + MfRho = self._MfRho + MfRhoDeriv = self._MfRhoDeriv + 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) + + 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 + + return hDeriv_m + + def _e(self, jSolution, srcList): + """ + Electric field from jSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: electric field + """ + 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))) + + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the field we solved for with a vector + """ + 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: + return self._MfRho.T * ( self._aveF2CCV.T * ( VI.T * du_dm_v ) ) + return VI * (self._aveF2CCV * (self._MfRho * du_dm_v)) + + def _eDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the model with a vector + """ + jSolution = Utils.mkvc(self[src,'jSolution']) + 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: + return self._MfRhoDeriv(jSolution).T * ( self._aveF2CCV.T * ( VI.T * v ) ) + return VI * (self._aveF2CCV * (self._MfRhoDeriv(jSolution) * v)) + + def _b(self, jSolution, srcList): + """ + Secondary magnetic flux density from jSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :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)) ) + + 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._aveF2CCV.shape[0] / self._nC) # number of components + VI = sdiag(np.kron(np.ones(n), 1./self.prob.mesh.vol)) + + if adjoint: + return -1./(1j*omega(src.freq)) * self._MfRho.T * ( self._edgeCurl * ( self._aveE2CCV.T * (VI.T * du_dm_v) ) ) + return -1./(1j*omega(src.freq)) * VI * (self._aveE2CCV * (self._edgeCurl.T * (self._MfRho * du_dm_v))) + + def _bDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the magnetic flux density with respect to the model with a vector + """ + 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) + + 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 ) ) ) 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'} @@ -680,7 +1003,9 @@ class Fields_h(Fields): 'hSecondary' : ['hSolution','E','_hSecondary'], 'j' : ['hSolution','F','_j'], 'jPrimary' : ['hSolution','F','_jPrimary'], - 'jSecondary' : ['hSolution','F','_jSecondary'] + 'jSecondary' : ['hSolution','F','_jSecondary'], + 'e' : ['hSolution','CCV','_e'], + 'b' : ['hSolution','CCV','_b'], } def __init__(self,mesh,survey,**kwargs): @@ -689,8 +1014,25 @@ class Fields_h(Fields): def startup(self): self.prob = self.survey.prob self._edgeCurl = self.survey.prob.mesh.edgeCurl + self._MeMu = self.survey.prob.MeMu self._MeMuI = self.survey.prob.MeMuI self._MfRho = self.survey.prob.MfRho + self._MfRhoDeriv = self.survey.prob.MfRhoDeriv + self._rho = self.survey.prob.curModel.rho + self._mu = self.survey.prob.curModel.mui + self._aveF2CCV = self.survey.prob.mesh.aveF2CCV + self._aveE2CCV = self.survey.prob.mesh.aveE2CCV + self._nC = self.survey.prob.mesh.nC + + def _GLoc(self,fieldType): + if fieldType == 'h': + return 'E' + elif fieldType == 'j': + return 'F' + elif (fieldType == 'e') or (fieldType == 'b'): + return 'CCV' + else: + raise Exception('Field type must be e, b, h, j') def _hPrimary(self, hSolution, srcList): """ @@ -720,36 +1062,25 @@ class Fields_h(Fields): return hSolution - def _h(self, hSolution, srcList): - """ - Total magnetic field is sum of primary and secondary - - :param numpy.ndarray hSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total magnetic field - """ - return self._hPrimary(hSolution, srcList) + self._hSecondary(hSolution, srcList) - - def _hDeriv_u(self, src, v, adjoint=False): + def _hDeriv_u(self, src, du_dm_v, adjoint=False): """ - Derivative of the total magnetic field with respect to the thing we - solved for - + Partial derivative of the total magnetic field with respect to the thing we + solved for. + :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the magnetic field with respect to the field we solved for with a vector """ - return Identity()*v + return Identity()*du_dm_v def _hDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. - + Partial derivative of the total magnetic field with respect to the inversion model. Here, we assume that the primary does not depend on the model. Note that this also includes derivative contributions from the sources. + :param SimPEG.EM.FDEM.Src src: source :param numpy.ndarray v: vector to take product with :param bool adjoint: adjoint? @@ -778,7 +1109,7 @@ class Fields_h(Fields): def _jSecondary(self, hSolution, srcList): """ - Secondary current density from eSolution + Secondary current density from hSolution :param numpy.ndarray hSolution: field we solved for :param list srcList: list of sources @@ -788,74 +1119,128 @@ 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 _jSecondaryDeriv_u(self, src, v, adjoint=False): + def _jDeriv_u(self, src, du_dm_v, adjoint=False): """ - Derivative of the secondary current density with respect to the thing we solved for - + Derivative of the 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? - :rtype: numpy.ndarray - :return: product of the derivative of the secondary current density with respect to the field we solved for with a vector - """ - - if not adjoint: - return self._edgeCurl*v - elif adjoint: - return self._edgeCurl.T*v - - def _jSecondaryDeriv_m(self, src, v, adjoint=False): - """ - Derivative of the secondary 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? - :rtype: numpy.ndarray - :return: product of the secondary current density derivative with respect to the inversion model with a vector - """ - - _,S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - return -S_eDeriv - - def _j(self, hSolution, srcList): - """ - Total current density is sum of primary and secondary - - :param numpy.ndarray eSolution: field we solved for - :param list srcList: list of sources - :rtype: numpy.ndarray - :return: total current density - """ - - return self._jPrimary(hSolution, srcList) + self._jSecondary(hSolution, srcList) - - def _jDeriv_u(self, src, v, adjoint=False): - """ - Derivative of the total current density with respect to the thing we solved for - - :param SimPEG.EM.FDEM.Src src: source - :param numpy.ndarray v: vector to take product with + :param numpy.ndarray du_dm_v: vector to take product with :param bool adjoint: adjoint? :rtype: numpy.ndarray :return: product of the derivative of the current density with respect to the field we solved for with a vector """ - return self._jSecondaryDeriv_u(src,v,adjoint) + + if not adjoint: + return self._edgeCurl*du_dm_v + elif adjoint: + return self._edgeCurl.T*du_dm_v + def _jDeriv_m(self, src, v, adjoint=False): """ - Derivative of the total 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? - :rtype: SimPEG.Utils.Zero - :return: product of the current density with respect to the inversion model with a vector + :rtype: numpy.ndarray + :return: product of the current density derivative with respect to the inversion model with a vector """ - # assuming the primary does not depend on the model - return self._jSecondaryDeriv_m(src,v,adjoint) + _,s_eDeriv = src.evalDeriv(self.prob, v, adjoint) + return -s_eDeriv + + def _e(self, hSolution, srcList): + """ + Electric field from hSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: electric field + """ + 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))) + + 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? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the field we solved for with a vector + """ + 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: + 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 )) + + def _eDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the electric field derivative with respect to the inversion model with a vector + """ + 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: + return ( self._MfRhoDeriv(self._edgeCurl * hSolution).T * ( self._aveF2CCV.T * (VI.T * v) ) ) + return VI * (self._aveF2CCV * (self._MfRhoDeriv(self._edgeCurl * hSolution) * v )) + + def _b(self, hSolution, srcList): + """ + Magnetic flux density from hSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: magnetic flux density + """ + h = self._h(hSolution, srcList) + 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)) + + 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 + 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 )) + return VI * (self._aveE2CCV * (self._MeMu * du_dm_v)) + + def _bDeriv_m(self, src, v, adjoint=False): + """ + 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? + :rtype: numpy.ndarray + :return: product of the magnetic flux density derivative with respect to the inversion model with a vector + """ + return Zero() + + diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 1213cef3..87967dd5 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -1,7 +1,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from scipy.constants import mu_0 from SimPEG.EM.Utils import * -from SimPEG.Utils import Zero +from SimPEG.Utils import Zero class BaseSrc(Survey.BaseSrc): """ @@ -14,34 +14,34 @@ class BaseSrc(Survey.BaseSrc): def eval(self, prob): """ - Evaluate the source terms. - - :math:`S_m` : magnetic source term - - :math:`S_e` : electric source term + Evaluate the source terms. + - :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 :param bool adjoint: adjoint? :rtype: (numpy.ndarray, numpy.ndarray) - :return: tuple with magnetic source term and electric source term derivatives times a vector + :return: tuple with magnetic source term and electric source term derivatives times a vector """ - if v is not None: - return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint) + if v is not None: + return self.s_mDeriv(prob, v, adjoint), self.s_eDeriv(prob, v, adjoint) else: - return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) + return lambda v: self.s_mDeriv(prob, v, adjoint), lambda v: self.s_eDeriv(prob, v, adjoint) def bPrimary(self, prob): """ @@ -49,7 +49,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary magnetic flux density + :return: primary magnetic flux density """ return Zero() @@ -59,7 +59,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ return Zero() @@ -69,7 +69,7 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary electric field + :return: primary electric field """ return Zero() @@ -79,13 +79,13 @@ class BaseSrc(Survey.BaseSrc): :param Problem prob: FDEM Problem :rtype: numpy.ndarray - :return: primary current density + :return: primary current density """ return Zero() - def S_m(self, prob): + def s_m(self, prob): """ - Magnetic source term + Magnetic source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray @@ -93,9 +93,9 @@ class BaseSrc(Survey.BaseSrc): """ return Zero() - def S_e(self, prob): + def s_e(self, prob): """ - Electric source term + Electric source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray @@ -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,88 +131,117 @@ 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): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_e = np.array(S_e,dtype=complex) - self.freq = float(freq) - BaseSrc.__init__(self, rxList) - - def S_e(self, prob): - - return self._S_e - - -class RawVec_m(BaseSrc): - """ - 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 - """ - - 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_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_m(self, prob): + def s_e(self, prob): """ - Magnetic source term + Electric source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :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 + + +class RawVec_m(BaseSrc): + """ + 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 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) + self.freq = float(freq) + self.integrate = integrate + + BaseSrc.__init__(self, rxList) + + def s_m(self, prob): + """ + Magnetic source term :param Problem prob: FDEM Problem :rtype: numpy.ndarray :return: magnetic source term on mesh """ - return self._S_m + if prob._formulation is 'HJ' and self.integrate is True: + 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): - if prob._eqLocs is 'EF' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m + def s_m(self, prob): + """ + Magnetic source term - def S_e(self, prob): - if prob._eqLocs is 'FE' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :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 + + def s_e(self, prob): + """ + Electric source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :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 class MagDipole(BaseSrc): - """ + """ Point magnetic dipole source calculated by taking the curl of a magnetic vector potential. By taking the discrete curl, we ensure that the magnetic - flux density is divergence free (no magnetic monopoles!). + flux density is divergence free (no magnetic monopoles!). This approach uses a primary-secondary in frequency. Here we show the derivation for E-B formulation noting that similar steps are followed for the H-J formulation. - .. math:: + .. math:: \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} @@ -225,17 +254,17 @@ class MagDipole(BaseSrc): and define a zero-frequency primary problem, noting that the source is generated by a divergence free electric current - .. math:: + .. math:: \mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\ {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}} - Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is + Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is - .. math:: + .. math:: \mathbf{e^P} = 0 \\\\ {\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}} - Our secondary problem is then + Our secondary problem is then .. math:: \mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\ @@ -245,15 +274,15 @@ class MagDipole(BaseSrc): :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0): self.freq = float(freq) self.loc = loc self.orientation = orientation + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.moment = moment self.mu = mu self.integrate = False @@ -265,17 +294,17 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': gridX = prob.mesh.gridEx gridY = prob.mesh.gridEy gridZ = prob.mesh.gridEz C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + elif formulation is 'HJ': gridX = prob.mesh.gridFx gridY = prob.mesh.gridFy gridZ = prob.mesh.gridFz @@ -303,44 +332,46 @@ class MagDipole(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) - return h_from_b(prob,b) + return 1./self.mu * b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b_p = self.bPrimary(prob) + if prob._formulation is 'HJ': + 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 :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': mui_s = prob.curModel.mui - 1./self.mu MMui_s = prob.mesh.getFaceInnerProduct(mui_s) C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + elif formulation is 'HJ': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) @@ -353,21 +384,20 @@ class MagDipole_Bfield(BaseSrc): fields from a magnetic dipole. No discrete curl is taken, so the magnetic flux density may not be strictly divergence free. - This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. + This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. :param list rxList: receiver list :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - #TODO: neither does moment def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): self.freq = float(freq) self.loc = loc + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.orientation = orientation self.moment = moment self.mu = mu @@ -379,18 +409,18 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': gridX = prob.mesh.gridFx gridY = prob.mesh.gridFy gridZ = prob.mesh.gridFz C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + elif formulation is 'HJ': gridX = prob.mesh.gridEx gridY = prob.mesh.gridEy gridZ = prob.mesh.gridEz @@ -418,42 +448,44 @@ class MagDipole_Bfield(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) - return h_from_b(prob, b) + return 1/self.mu * b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) + if prob._formulation is 'HJ': + b = prob.Me * b return -1j*omega(self.freq)*b - def S_e(self, prob): + def s_e(self, prob): """ The electric source term :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': mui_s = prob.curModel.mui - 1./self.mu MMui_s = prob.mesh.getFaceInnerProduct(mui_s) C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + elif formulation is 'HJ': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) @@ -463,22 +495,22 @@ class CircularLoop(BaseSrc): """ Circular loop magnetic source calculated by taking the curl of a magnetic vector potential. By taking the discrete curl, we ensure that the magnetic - flux density is divergence free (no magnetic monopoles!). + flux density is divergence free (no magnetic monopoles!). - This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. + This approach uses a primary-secondary in frequency in the same fashion as the MagDipole. :param list rxList: receiver list :param float freq: frequency :param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`) :param string orientation: 'X', 'Y', 'Z' - :param float moment: magnetic dipole moment - :param float mu: background magnetic permeability + :param float moment: magnetic dipole moment + :param float mu: background magnetic permeability """ - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', radius=1., mu=mu_0): self.freq = float(freq) self.orientation = orientation + assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.radius = radius self.mu = mu self.loc = loc @@ -491,17 +523,17 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': gridX = prob.mesh.gridEx gridY = prob.mesh.gridEy gridZ = prob.mesh.gridEz C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + elif formulation is 'HJ': gridX = prob.mesh.gridFx gridY = prob.mesh.gridFy gridZ = prob.mesh.gridFz @@ -528,44 +560,50 @@ class CircularLoop(BaseSrc): :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) return 1./self.mu*b - def S_m(self, prob): + def s_m(self, prob): """ The magnetic source term :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ b = self.bPrimary(prob) + if prob._formulation is 'HJ': + b = prob.Me * b return -1j*omega(self.freq)*b - def S_e(self, prob): + def s_e(self, prob): """ The electric source term :param Problem prob: FDEM problem :rtype: numpy.ndarray - :return: primary magnetic field + :return: primary magnetic field """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: - eqLocs = prob._eqLocs + formulation = prob._formulation - if eqLocs is 'FE': + if formulation is 'EB': mui_s = prob.curModel.mui - 1./self.mu MMui_s = prob.mesh.getFaceInnerProduct(mui_s) C = prob.mesh.edgeCurl - elif eqLocs is 'EF': + + + elif formulation is 'HJ': mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True) C = prob.mesh.edgeCurl.T return -C.T * (MMui_s * self.bPrimary(prob)) + + diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 4d220259..1552a12c 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,8 +1,10 @@ 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 +from SimPEG import sp #################################################### @@ -18,33 +20,33 @@ class Rx(SimPEG.Survey.BaseRx): """ knownRxTypes = { - 'exr':['e', 'Ex', 'real'], - 'eyr':['e', 'Ey', 'real'], - 'ezr':['e', 'Ez', 'real'], - 'exi':['e', 'Ex', 'imag'], - 'eyi':['e', 'Ey', 'imag'], - 'ezi':['e', 'Ez', 'imag'], + 'exr':['e', 'x', 'real'], + 'eyr':['e', 'y', 'real'], + 'ezr':['e', 'z', 'real'], + 'exi':['e', 'x', 'imag'], + 'eyi':['e', 'y', 'imag'], + 'ezi':['e', 'z', 'imag'], - 'bxr':['b', 'Fx', 'real'], - 'byr':['b', 'Fy', 'real'], - 'bzr':['b', 'Fz', 'real'], - 'bxi':['b', 'Fx', 'imag'], - 'byi':['b', 'Fy', 'imag'], - 'bzi':['b', 'Fz', 'imag'], + 'bxr':['b', 'x', 'real'], + 'byr':['b', 'y', 'real'], + 'bzr':['b', 'z', 'real'], + 'bxi':['b', 'x', 'imag'], + 'byi':['b', 'y', 'imag'], + 'bzi':['b', 'z', 'imag'], - 'jxr':['j', 'Fx', 'real'], - 'jyr':['j', 'Fy', 'real'], - 'jzr':['j', 'Fz', 'real'], - 'jxi':['j', 'Fx', 'imag'], - 'jyi':['j', 'Fy', 'imag'], - 'jzi':['j', 'Fz', 'imag'], + 'jxr':['j', 'x', 'real'], + 'jyr':['j', 'y', 'real'], + 'jzr':['j', 'z', 'real'], + 'jxi':['j', 'x', 'imag'], + 'jyi':['j', 'y', 'imag'], + 'jzi':['j', 'z', 'imag'], - 'hxr':['h', 'Ex', 'real'], - 'hyr':['h', 'Ey', 'real'], - 'hzr':['h', 'Ez', 'real'], - 'hxi':['h', 'Ex', 'imag'], - 'hyi':['h', 'Ey', 'imag'], - 'hzi':['h', 'Ez', 'imag'], + 'hxr':['h', 'x', 'real'], + 'hyr':['h', 'y', 'real'], + 'hzr':['h', 'z', 'real'], + 'hxi':['h', 'x', 'imag'], + 'hyi':['h', 'y', 'imag'], + 'hzi':['h', 'z', 'imag'], } radius = None @@ -56,16 +58,15 @@ class Rx(SimPEG.Survey.BaseRx): """Field Type projection (e.g. e b ...)""" return self.knownRxTypes[self.rxType][0] - @property - def projGLoc(self): - """Grid Location projection (e.g. Ex Fy ...)""" - return self.knownRxTypes[self.rxType][1] - @property def projComp(self): """Component projection (real/imag)""" return self.knownRxTypes[self.rxType][2] + def projGLoc(self, u): + """Grid Location projection (e.g. Ex Fy ...)""" + return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] + def eval(self, src, mesh, f): """ Project fields to recievers to get data. @@ -76,11 +77,16 @@ class Rx(SimPEG.Survey.BaseRx): :rtype: numpy.ndarray :return: fields projected to recievers """ - P = self.getP(mesh) # get interpolation to recievers - u_part_complex = f[src, self.projField] - real_or_imag = self.projComp # get the real or imag component - u_part = getattr(u_part_complex, real_or_imag) - return P*u_part + # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) + # projGLoc += self.knownRxTypes[self.rxType][1] + + 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 + f_part = getattr(f_part_complex, real_or_imag) + + return P*f_part def evalDeriv(self, src, mesh, f, v, adjoint=False): """ @@ -93,7 +99,8 @@ class Rx(SimPEG.Survey.BaseRx): :rtype: numpy.ndarray :return: fields projected to recievers """ - P = self.getP(mesh) + + P = self.getP(mesh, self.projGLoc(f)) if not adjoint: Pv_complex = P * v @@ -117,7 +124,7 @@ class Rx(SimPEG.Survey.BaseRx): # Survey #################################################### -class Survey(SimPEG.Survey.BaseSurvey): +class Survey(BaseEMSurvey): """ Frequency domain electromagnetic survey @@ -125,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: @@ -165,23 +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/EM/Utils/EMUtils.py b/SimPEG/EM/Utils/EMUtils.py index 4a342acb..e7dbf441 100644 --- a/SimPEG/EM/Utils/EMUtils.py +++ b/SimPEG/EM/Utils/EMUtils.py @@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0): beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) ) return alp - 1j*beta -# Constitutive relations -def e_from_j(prob,j): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MSigmaI = prob.MeSigmaI - elif eqLocs is 'EF': - MSigmaI = prob.MfRho - return MSigmaI*j - -def j_from_e(prob,e): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MSigma = prob.MeSigma - elif eqLocs is 'EF': - MSigma = prob.MfRhoI - return MSigma*e - -def b_from_h(prob,h): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MMu = prob.MfMuiI - elif eqLocs is 'EF': - MMu = prob.MeMu - return MMu*h - -def h_from_b(prob,b): - eqLocs = prob._eqLocs - if eqLocs is 'FE': - MMuI = prob.MfMui - elif eqLocs is 'EF': - MMuI = prob.MeMuI - return MMuI*b - diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index 18dddde9..ef779042 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,2 @@ -# import Sources -# import Ana -# import Solver -from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from EMUtils import omega, k from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 8c703083..569f8e6d 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -4,19 +4,28 @@ from SimPEG import EM import sys from scipy.constants import mu_0 -def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): - cs = 5. - ncx, ncy, ncz = 6, 6, 6 - npad = 3 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 5e-1 + + +def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): + cs = 10. + ncx, ncy, ncz = 0, 0, 0 + npad = 8 hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) - mapping = Maps.ExpMap(mesh) + if useMu is True: + mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] + else: + mapping = Maps.ExpMap(mesh) - x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source - XYZ = Utils.ndgrid(x,x,np.r_[0.]) + x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid + XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5)) Rx0 = EM.FDEM.Rx(XYZ, comp) Src = [] @@ -32,15 +41,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): if fdemType is 'e' or fdemType is 'b': S_m = np.zeros(mesh.nF) S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) elif fdemType is 'h' or fdemType is 'j': S_m = np.zeros(mesh.nE) S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) if verbose: @@ -70,6 +79,48 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): from pymatsolver import MumpsSolver prb.Solver = MumpsSolver except ImportError, e: - pass + prb.Solver = SolverLU - return prb \ No newline at end of file + return prb + +def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False): + + l2norm = lambda r: np.sqrt(r.dot(r)) + + prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) + mesh = prb1.mesh + print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) + + logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) + mu = np.ones(mesh.nC)*MU + + if addrandoms is True: + logsig += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 + mu += np.random.randn(mesh.nC)*MU*1e-1 + + if useMu is True: + m = np.r_[logsig, mu] + else: + m = logsig + + survey1 = prb1.survey + d1 = survey1.dpred(m) + + if verbose: + print ' Problem 1 solved' + + + prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose) + + survey2 = prb2.survey + d2 = survey2.dpred(m) + + if verbose: + print ' Problem 2 solved' + + r = d2-d1 + l2r = l2norm(r) + + tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR]) + print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol + return l2r < tol diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index e76b2439..29f51ed4 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -48,8 +48,7 @@ def run(plotIt=True): freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) - srcList = [] - [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs] + srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs] survey = EM.FDEM.Survey(srcList) prb = EM.FDEM.Problem_b(mesh, mapping=mapping) 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/__init__.py b/SimPEG/Examples/__init__.py index cce22296..4647ad90 100644 --- a/SimPEG/Examples/__init__.py +++ b/SimPEG/Examples/__init__.py @@ -5,6 +5,7 @@ import DC_Analytic_Dipole import DC_Forward_PseudoSection import EM_FDEM_1D_Inversion import EM_FDEM_Analytic_MagDipoleWholespace +import EM_Schenkel_Morrison_Casing import EM_TDEM_1D_Inversion import FLOW_Richards_1D_Celia1990 import Forward_BasicDirectCurrent @@ -19,7 +20,7 @@ import Mesh_Tensor_Creation import MT_1D_ForwardAndInversion import MT_3D_Foward -__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"] +__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "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_ForwardAndInversion", "MT_3D_Foward"] ##### AUTOIMPORTS ##### 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/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index 508f015c..1650b5cb 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -234,6 +234,9 @@ class BaseTensorMesh(BaseMesh): 'Fz' -> z-component of field defined on faces 'N' -> scalar field defined on nodes 'CC' -> scalar field defined on cell centers + 'CCVx' -> x-component of vector field defined on cell centers + 'CCVy' -> y-component of vector field defined on cell centers + 'CCVz' -> z-component of vector field defined on cell centers """ if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']: raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType) @@ -257,6 +260,16 @@ class BaseTensorMesh(BaseMesh): Q = sp.hstack(components) elif locType in ['CC', 'N']: Q = Utils.interpmat(loc, *self.getTensor(locType)) + elif locType in ['CCVx', 'CCVy', 'CCVz']: + Q = Utils.interpmat(loc, *self.getTensor('CC')) + Z = Utils.spzeros(loc.shape[0],self.nC) + if locType == 'CCVx': + Q = sp.hstack([Q,Z,Z]) + elif locType == 'CCVy': + Q = sp.hstack([Z,Q,Z]) + elif locType == 'CCVz': + Q = sp.hstack([Z,Z,Q]) + else: raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) 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 9f307c3f..fbc88276 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -34,7 +34,7 @@ class BaseRx(object): """Number of data in the receiver.""" return self.locs.shape[0] - def getP(self, mesh): + def getP(self, mesh, projGLoc=None): """ Returns the projection matrices as a list for all components collected by @@ -47,7 +47,10 @@ class BaseRx(object): if mesh in self._Ps: return self._Ps[mesh] - P = mesh.getInterpolationMat(self.locs, self.projGLoc) + if projGLoc is None: + projGLoc = self.projGLoc + + P = mesh.getInterpolationMat(self.locs, projGLoc) if self.storeProjections: self._Ps[mesh] = P return P @@ -292,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. @@ -334,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 @@ -349,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. @@ -369,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/__init__.py b/SimPEG/__init__.py index cc51fd1f..a1e989b6 100644 --- a/SimPEG/__init__.py +++ b/SimPEG/__init__.py @@ -15,7 +15,7 @@ import Directives import Inversion import Tests -__version__ = '0.1.9' +__version__ = '0.1.10' __author__ = 'Rowan Cockett' __license__ = 'MIT' __copyright__ = 'Copyright 2014 Rowan Cockett' diff --git a/docs/conf.py b/docs/conf.py index fee262de..45407435 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -51,9 +51,9 @@ copyright = u'2013, SimPEG Developers' # built documents. # # The short X.Y version. -version = '0.1.9' +version = '0.1.10' # The full version, including alpha/beta/rc tags. -release = '0.1.9' +release = '0.1.10' # The language for content autogenerated by Sphinx. Refer to documentation # for a list of supported languages. 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/setup.py b/setup.py index bcb5b8e3..3383c7f7 100644 --- a/setup.py +++ b/setup.py @@ -77,7 +77,7 @@ with open("README.rst") as f: setup( name = "SimPEG", - version = "0.1.9", + version = "0.1.10", packages = find_packages(), install_requires = ['numpy>=1.7', 'scipy>=0.13', diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py index 437f3708..da446675 100644 --- a/tests/em/fdem/forward/test_FDEM_forward.py +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -3,125 +3,75 @@ from SimPEG import * from SimPEG import EM import sys from scipy.constants import mu_0 -from SimPEG.EM.Utils.testingUtils import getFDEMProblem +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest testEB = True testHJ = True - +testEJ = True +testBH = True verbose = False -TOL = 1e-5 -FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order -CONDUCTIVITY = 1e1 -MU = mu_0 -freq = 1e-1 -addrandoms = True +TOLEBHJ = 1e-5 +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] -def crossCheckTest(fdemType, comp): - - l2norm = lambda r: np.sqrt(r.dot(r)) - - prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose) - mesh = prb1.mesh - print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp) - m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) - mu = np.log(np.ones(mesh.nC)*MU) - - if addrandoms is True: - m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(mesh.nC)*MU*1e-1 - - # prb1.PropMap.PropModel.mu = mu - # prb1.PropMap.PropModel.mui = 1./mu - survey1 = prb1.survey - d1 = survey1.dpred(m) - - if verbose: - print ' Problem 1 solved' - - if fdemType == 'e': - prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose) - elif fdemType == 'b': - prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose) - elif fdemType == 'j': - prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose) - elif fdemType == 'h': - prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose) - else: - raise NotImplementedError() - - # prb2.mu = mu - survey2 = prb2.survey - d2 = survey2.dpred(m) - - if verbose: - print ' Problem 2 solved' - - r = d2-d1 - l2r = l2norm(r) - - tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) - print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol - return l2r < tol - - class FDEM_CrossCheck(unittest.TestCase): if testEB: def test_EB_CrossCheck_exr_Eform(self): - self.assertTrue(crossCheckTest('e', 'exr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exr', verbose=verbose)) def test_EB_CrossCheck_eyr_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyr', verbose=verbose)) def test_EB_CrossCheck_ezr_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezr', verbose=verbose)) def test_EB_CrossCheck_exi_Eform(self): - self.assertTrue(crossCheckTest('e', 'exi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'exi', verbose=verbose)) def test_EB_CrossCheck_eyi_Eform(self): - self.assertTrue(crossCheckTest('e', 'eyi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'eyi', verbose=verbose)) def test_EB_CrossCheck_ezi_Eform(self): - self.assertTrue(crossCheckTest('e', 'ezi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'ezi', verbose=verbose)) def test_EB_CrossCheck_bxr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxr', verbose=verbose)) def test_EB_CrossCheck_byr_Eform(self): - self.assertTrue(crossCheckTest('e', 'byr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byr', verbose=verbose)) def test_EB_CrossCheck_bzr_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzr')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzr', verbose=verbose)) def test_EB_CrossCheck_bxi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bxi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bxi', verbose=verbose)) def test_EB_CrossCheck_byi_Eform(self): - self.assertTrue(crossCheckTest('e', 'byi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'byi', verbose=verbose)) def test_EB_CrossCheck_bzi_Eform(self): - self.assertTrue(crossCheckTest('e', 'bzi')) + self.assertTrue(crossCheckTest(SrcList, 'e', 'b', 'bzi', verbose=verbose)) if testHJ: def test_HJ_CrossCheck_jxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxr', verbose=verbose)) def test_HJ_CrossCheck_jyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyr', verbose=verbose)) def test_HJ_CrossCheck_jzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzr', verbose=verbose)) def test_HJ_CrossCheck_jxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jxi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jxi', verbose=verbose)) def test_HJ_CrossCheck_jyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jyi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jyi', verbose=verbose)) def test_HJ_CrossCheck_jzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'jzi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'jzi', verbose=verbose)) def test_HJ_CrossCheck_hxr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxr', verbose=verbose)) def test_HJ_CrossCheck_hyr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyr', verbose=verbose)) def test_HJ_CrossCheck_hzr_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzr')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzr', verbose=verbose)) def test_HJ_CrossCheck_hxi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hxi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hxi', verbose=verbose)) def test_HJ_CrossCheck_hyi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hyi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hyi', verbose=verbose)) def test_HJ_CrossCheck_hzi_Jform(self): - self.assertTrue(crossCheckTest('j', 'hzi')) + self.assertTrue(crossCheckTest(SrcList, 'j', 'h', 'hzi', verbose=verbose)) if __name__ == '__main__': unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/forward/test_FDEM_forwardEJHB.py b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py new file mode 100644 index 00000000..e6319dfc --- /dev/null +++ b/tests/em/fdem/forward/test_FDEM_forwardEJHB.py @@ -0,0 +1,125 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest + +testEJ = True +testBH = True + +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this + +SrcList = ['RawVec', 'MagDipole', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] + + +class FDEM_CrossCheck(unittest.TestCase): + if testEJ: + def test_EJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'jzi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'exi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'eyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'ezi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'byi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'bzi', TOL=TOLEJHB)) + + def test_EJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzr', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hxi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hyi', TOL=TOLEJHB)) + def test_EJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'e', 'j', 'hzi', TOL=TOLEJHB)) + + if testBH: + def test_HB_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'jzi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_exr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exr', TOL=TOLEJHB)) + def test_HB_CrossCheck_eyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_ezr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezr', TOL=TOLEJHB)) + def test_HB_CrossCheck_exi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'exi', TOL=TOLEJHB)) + def test_HB_CrossCheck_eyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'eyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_ezi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'ezi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_bxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_byr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byr', TOL=TOLEJHB)) + def test_HB_CrossCheck_bzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_bxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_byi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'byi', TOL=TOLEJHB)) + def test_HB_CrossCheck_bzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'bzi', TOL=TOLEJHB)) + + def test_HB_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hzr_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzr', TOL=TOLEJHB)) + def test_HB_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hxi', TOL=TOLEJHB)) + def test_HB_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hyi', TOL=TOLEJHB)) + def test_HB_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest(SrcList, 'h', 'b', 'hzi', TOL=TOLEJHB)) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py new file mode 100644 index 00000000..545a5014 --- /dev/null +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -0,0 +1,128 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem, crossCheckTest + +testEB = True +testHJ = True +testEJ = True +testBH = True +verbose = False + +TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) +#TODO: choose better testing parameters to lower this + +SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] + + +class FDEM_CrossCheck(unittest.TestCase): + if testBH: + def test_BH_CrossCheck_jxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_exi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) + + if testBH: + def test_BH_CrossCheck_jxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_jzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'jzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_exr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_exi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'exi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_eyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'eyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_ezi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'ezi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_bxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_byi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'byi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_bzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'bzi', verbose=verbose, TOL=TOLEJHB)) + + def test_BH_CrossCheck_hxr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzr(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzr', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hxi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hxi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hyi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hyi', verbose=verbose, TOL=TOLEJHB)) + def test_BH_CrossCheck_hzi(self): + self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) + +if __name__ == '__main__': + unittest.main() \ No newline at end of file diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py similarity index 57% rename from tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py rename to tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py index f77f2131..25762368 100644 --- a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointEB.py @@ -5,8 +5,8 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem -testEB = True -testHJ = True +testE = True +testB = True verbose = False @@ -17,10 +17,10 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' def adjointTest(fdemType, comp): - prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + prb = getFDEMProblem(fdemType, comp, SrcList, freq) print 'Adjoint %s formulation - %s' % (fdemType, comp) m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) @@ -45,7 +45,7 @@ def adjointTest(fdemType, comp): return np.abs(vJw - wJtv) < tol class FDEM_AdjointTests(unittest.TestCase): - if testEB: + if testE: def test_Jtvec_adjointTest_exr_Eform(self): self.assertTrue(adjointTest('e', 'exr')) def test_Jtvec_adjointTest_eyr_Eform(self): @@ -72,6 +72,33 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_bzi_Eform(self): self.assertTrue(adjointTest('e', 'bzi')) + def test_Jtvec_adjointTest_jxr_Eform(self): + self.assertTrue(adjointTest('e', 'jxr')) + def test_Jtvec_adjointTest_jyr_Eform(self): + self.assertTrue(adjointTest('e', 'jyr')) + def test_Jtvec_adjointTest_jzr_Eform(self): + self.assertTrue(adjointTest('e', 'jzr')) + def test_Jtvec_adjointTest_jxi_Eform(self): + self.assertTrue(adjointTest('e', 'jxi')) + def test_Jtvec_adjointTest_jyi_Eform(self): + self.assertTrue(adjointTest('e', 'jyi')) + def test_Jtvec_adjointTest_jzi_Eform(self): + self.assertTrue(adjointTest('e', 'jzi')) + + def test_Jtvec_adjointTest_hxr_Eform(self): + self.assertTrue(adjointTest('e', 'hxr')) + def test_Jtvec_adjointTest_hyr_Eform(self): + self.assertTrue(adjointTest('e', 'hyr')) + def test_Jtvec_adjointTest_hzr_Eform(self): + self.assertTrue(adjointTest('e', 'hzr')) + def test_Jtvec_adjointTest_hxi_Eform(self): + self.assertTrue(adjointTest('e', 'hxi')) + def test_Jtvec_adjointTest_hyi_Eform(self): + self.assertTrue(adjointTest('e', 'hyi')) + def test_Jtvec_adjointTest_hzi_Eform(self): + self.assertTrue(adjointTest('e', 'hzi')) + + if testB: def test_Jtvec_adjointTest_exr_Bform(self): self.assertTrue(adjointTest('b', 'exr')) def test_Jtvec_adjointTest_eyr_Bform(self): @@ -84,6 +111,7 @@ class FDEM_AdjointTests(unittest.TestCase): self.assertTrue(adjointTest('b', 'eyi')) def test_Jtvec_adjointTest_ezi_Bform(self): self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_bxr_Bform(self): self.assertTrue(adjointTest('b', 'bxr')) def test_Jtvec_adjointTest_byr_Bform(self): @@ -97,59 +125,31 @@ class FDEM_AdjointTests(unittest.TestCase): def test_Jtvec_adjointTest_bzi_Bform(self): self.assertTrue(adjointTest('b', 'bzi')) + def test_Jtvec_adjointTest_jxr_Bform(self): + self.assertTrue(adjointTest('b', 'jxr')) + def test_Jtvec_adjointTest_jyr_Bform(self): + self.assertTrue(adjointTest('b', 'jyr')) + def test_Jtvec_adjointTest_jzr_Bform(self): + self.assertTrue(adjointTest('b', 'jzr')) + def test_Jtvec_adjointTest_jxi_Bform(self): + self.assertTrue(adjointTest('b', 'jxi')) + def test_Jtvec_adjointTest_jyi_Bform(self): + self.assertTrue(adjointTest('b', 'jyi')) + def test_Jtvec_adjointTest_jzi_Bform(self): + self.assertTrue(adjointTest('b', 'jzi')) - if testHJ: - def test_Jtvec_adjointTest_jxr_Jform(self): - self.assertTrue(adjointTest('j', 'jxr')) - def test_Jtvec_adjointTest_jyr_Jform(self): - self.assertTrue(adjointTest('j', 'jyr')) - def test_Jtvec_adjointTest_jzr_Jform(self): - self.assertTrue(adjointTest('j', 'jzr')) - def test_Jtvec_adjointTest_jxi_Jform(self): - self.assertTrue(adjointTest('j', 'jxi')) - def test_Jtvec_adjointTest_jyi_Jform(self): - self.assertTrue(adjointTest('j', 'jyi')) - def test_Jtvec_adjointTest_jzi_Jform(self): - self.assertTrue(adjointTest('j', 'jzi')) - - def test_Jtvec_adjointTest_hxr_Jform(self): - self.assertTrue(adjointTest('j', 'hxr')) - def test_Jtvec_adjointTest_hyr_Jform(self): - self.assertTrue(adjointTest('j', 'hyr')) - def test_Jtvec_adjointTest_hzr_Jform(self): - self.assertTrue(adjointTest('j', 'hzr')) - def test_Jtvec_adjointTest_hxi_Jform(self): - self.assertTrue(adjointTest('j', 'hxi')) - def test_Jtvec_adjointTest_hyi_Jform(self): - self.assertTrue(adjointTest('j', 'hyi')) - def test_Jtvec_adjointTest_hzi_Jform(self): - self.assertTrue(adjointTest('j', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'hxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'hyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'hzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'hxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'hyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'jxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'jyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'jzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'jxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'jyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'jzi')) + def test_Jtvec_adjointTest_hxr_Bform(self): + self.assertTrue(adjointTest('b', 'hxr')) + def test_Jtvec_adjointTest_hyr_Bform(self): + self.assertTrue(adjointTest('b', 'hyr')) + def test_Jtvec_adjointTest_hzr_Bform(self): + self.assertTrue(adjointTest('b', 'hzr')) + def test_Jtvec_adjointTest_hxi_Bform(self): + self.assertTrue(adjointTest('b', 'hxi')) + def test_Jtvec_adjointTest_hyi_Bform(self): + self.assertTrue(adjointTest('b', 'hyi')) + def test_Jtvec_adjointTest_hzi_Bform(self): + self.assertTrue(adjointTest('b', 'hzi')) if __name__ == '__main__': diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py new file mode 100644 index 00000000..c3fb3d37 --- /dev/null +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjointHJ.py @@ -0,0 +1,155 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem + +testJ = True +testH = True + +verbose = False + +TOL = 1e-5 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 1e-1 +addrandoms = True + +SrcList = ['RawVec', 'MagDipole'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + +def adjointTest(fdemType, comp): + prb = getFDEMProblem(fdemType, comp, SrcList, freq) + print 'Adjoint %s formulation - %s' % (fdemType, comp) + + m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) + mu = np.ones(prb.mesh.nC)*MU + + if addrandoms is True: + m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 + + survey = prb.survey + u = prb.fields(m) + + v = np.random.rand(survey.nD) + w = np.random.rand(prb.mesh.nC) + + vJw = v.dot(prb.Jvec(m, w, u)) + wJtv = w.dot(prb.Jtvec(m, v, u)) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) + print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol + return np.abs(vJw - wJtv) < tol + +class FDEM_AdjointTests(unittest.TestCase): + + if testJ: + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) + + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + self.assertTrue(adjointTest('j', 'hzi')) + + def test_Jtvec_adjointTest_exr_Jform(self): + self.assertTrue(adjointTest('j', 'exr')) + def test_Jtvec_adjointTest_eyr_Jform(self): + self.assertTrue(adjointTest('j', 'eyr')) + def test_Jtvec_adjointTest_ezr_Jform(self): + self.assertTrue(adjointTest('j', 'ezr')) + def test_Jtvec_adjointTest_exi_Jform(self): + self.assertTrue(adjointTest('j', 'exi')) + def test_Jtvec_adjointTest_eyi_Jform(self): + self.assertTrue(adjointTest('j', 'eyi')) + def test_Jtvec_adjointTest_ezi_Jform(self): + self.assertTrue(adjointTest('j', 'ezi')) + + def test_Jtvec_adjointTest_bxr_Jform(self): + self.assertTrue(adjointTest('j', 'bxr')) + def test_Jtvec_adjointTest_byr_Jform(self): + self.assertTrue(adjointTest('j', 'byr')) + def test_Jtvec_adjointTest_bzr_Jform(self): + self.assertTrue(adjointTest('j', 'bzr')) + def test_Jtvec_adjointTest_bxi_Jform(self): + self.assertTrue(adjointTest('j', 'bxi')) + def test_Jtvec_adjointTest_byi_Jform(self): + self.assertTrue(adjointTest('j', 'byi')) + def test_Jtvec_adjointTest_bzi_Jform(self): + self.assertTrue(adjointTest('j', 'bzi')) + + if testH: + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) + + def test_Jtvec_adjointTest_jxr_Hform(self): + self.assertTrue(adjointTest('h', 'jxr')) + def test_Jtvec_adjointTest_jyr_Hform(self): + self.assertTrue(adjointTest('h', 'jyr')) + def test_Jtvec_adjointTest_jzr_Hform(self): + self.assertTrue(adjointTest('h', 'jzr')) + def test_Jtvec_adjointTest_jxi_Hform(self): + self.assertTrue(adjointTest('h', 'jxi')) + def test_Jtvec_adjointTest_jyi_Hform(self): + self.assertTrue(adjointTest('h', 'jyi')) + def test_Jtvec_adjointTest_jzi_Hform(self): + self.assertTrue(adjointTest('h', 'jzi')) + + def test_Jtvec_adjointTest_exr_Hform(self): + self.assertTrue(adjointTest('h', 'exr')) + def test_Jtvec_adjointTest_eyr_Hform(self): + self.assertTrue(adjointTest('h', 'eyr')) + def test_Jtvec_adjointTest_ezr_Hform(self): + self.assertTrue(adjointTest('h', 'ezr')) + def test_Jtvec_adjointTest_exi_Hform(self): + self.assertTrue(adjointTest('h', 'exi')) + def test_Jtvec_adjointTest_eyi_Hform(self): + self.assertTrue(adjointTest('h', 'eyi')) + def test_Jtvec_adjointTest_ezi_Hform(self): + self.assertTrue(adjointTest('h', 'ezi')) + + def test_Jtvec_adjointTest_bxr_Hform(self): + self.assertTrue(adjointTest('h', 'bxr')) + def test_Jtvec_adjointTest_byr_Hform(self): + self.assertTrue(adjointTest('h', 'byr')) + def test_Jtvec_adjointTest_bzr_Hform(self): + self.assertTrue(adjointTest('h', 'bzr')) + def test_Jtvec_adjointTest_bxi_Hform(self): + self.assertTrue(adjointTest('h', 'bxi')) + def test_Jtvec_adjointTest_byi_Hform(self): + self.assertTrue(adjointTest('h', 'byi')) + def test_Jtvec_adjointTest_bzi_Hform(self): + self.assertTrue(adjointTest('h', 'bzi')) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index d3bcb218..0a2e8b82 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -5,9 +5,11 @@ import sys from scipy.constants import mu_0 from SimPEG.EM.Utils.testingUtils import getFDEMProblem -testDerivs = True -testEB = True -testHJ = True + +testE = True +testB = True +testH = True +testJ = True verbose = False @@ -18,12 +20,12 @@ MU = mu_0 freq = 1e-1 addrandoms = True -SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' +SrcType = ['MagDipole', 'RawVec'] #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' def derivTest(fdemType, comp): - prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + prb = getFDEMProblem(fdemType, comp, SrcType, freq) print '%s formulation - %s' % (fdemType, comp) x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) mu = np.log(np.ones(prb.mesh.nC)*MU) @@ -32,9 +34,6 @@ def derivTest(fdemType, comp): x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1 - # prb.PropMap.PropModel.mu = mu - # prb.PropMap.PropModel.mui = 1./mu - survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) @@ -43,7 +42,7 @@ def derivTest(fdemType, comp): class FDEM_DerivTests(unittest.TestCase): - if testEB: + if testE: def test_Jvec_exr_Eform(self): self.assertTrue(derivTest('e', 'exr')) def test_Jvec_eyr_Eform(self): @@ -70,6 +69,33 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_bzi_Eform(self): self.assertTrue(derivTest('e', 'bzi')) + def test_Jvec_exr_Eform(self): + self.assertTrue(derivTest('e', 'jxr')) + def test_Jvec_eyr_Eform(self): + self.assertTrue(derivTest('e', 'jyr')) + def test_Jvec_ezr_Eform(self): + self.assertTrue(derivTest('e', 'jzr')) + def test_Jvec_exi_Eform(self): + self.assertTrue(derivTest('e', 'jxi')) + def test_Jvec_eyi_Eform(self): + self.assertTrue(derivTest('e', 'jyi')) + def test_Jvec_ezi_Eform(self): + self.assertTrue(derivTest('e', 'jzi')) + + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'hxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'hyr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'hzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'hxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'hyi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'hzi')) + + if testB: def test_Jvec_exr_Bform(self): self.assertTrue(derivTest('b', 'exr')) def test_Jvec_eyr_Bform(self): @@ -96,7 +122,33 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_bzi_Bform(self): self.assertTrue(derivTest('b', 'bzi')) - if testHJ: + def test_Jvec_jxr_Bform(self): + self.assertTrue(derivTest('b', 'jxr')) + def test_Jvec_jyr_Bform(self): + self.assertTrue(derivTest('b', 'jyr')) + def test_Jvec_jzr_Bform(self): + self.assertTrue(derivTest('b', 'jzr')) + def test_Jvec_jxi_Bform(self): + self.assertTrue(derivTest('b', 'jxi')) + def test_Jvec_jyi_Bform(self): + self.assertTrue(derivTest('b', 'jyi')) + def test_Jvec_jzi_Bform(self): + self.assertTrue(derivTest('b', 'jzi')) + + def test_Jvec_hxr_Bform(self): + self.assertTrue(derivTest('b', 'hxr')) + def test_Jvec_hyr_Bform(self): + self.assertTrue(derivTest('b', 'hyr')) + def test_Jvec_hzr_Bform(self): + self.assertTrue(derivTest('b', 'hzr')) + def test_Jvec_hxi_Bform(self): + self.assertTrue(derivTest('b', 'hxi')) + def test_Jvec_hyi_Bform(self): + self.assertTrue(derivTest('b', 'hyi')) + def test_Jvec_hzi_Bform(self): + self.assertTrue(derivTest('b', 'hzi')) + + if testJ: def test_Jvec_jxr_Jform(self): self.assertTrue(derivTest('j', 'jxr')) def test_Jvec_jyr_Jform(self): @@ -123,6 +175,34 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Jform(self): self.assertTrue(derivTest('j', 'hzi')) + def test_Jvec_exr_Jform(self): + self.assertTrue(derivTest('j', 'exr')) + def test_Jvec_eyr_Jform(self): + self.assertTrue(derivTest('j', 'eyr')) + def test_Jvec_ezr_Jform(self): + self.assertTrue(derivTest('j', 'ezr')) + def test_Jvec_exi_Jform(self): + self.assertTrue(derivTest('j', 'exi')) + def test_Jvec_eyi_Jform(self): + self.assertTrue(derivTest('j', 'eyi')) + def test_Jvec_ezi_Jform(self): + self.assertTrue(derivTest('j', 'ezi')) + + def test_Jvec_bxr_Jform(self): + self.assertTrue(derivTest('j', 'bxr')) + def test_Jvec_byr_Jform(self): + self.assertTrue(derivTest('j', 'byr')) + def test_Jvec_bzr_Jform(self): + self.assertTrue(derivTest('j', 'bzr')) + def test_Jvec_bxi_Jform(self): + self.assertTrue(derivTest('j', 'bxi')) + def test_Jvec_byi_Jform(self): + self.assertTrue(derivTest('j', 'byi')) + def test_Jvec_bzi_Jform(self): + self.assertTrue(derivTest('j', 'bzi')) + + + if testH: def test_Jvec_hxr_Hform(self): self.assertTrue(derivTest('h', 'hxr')) def test_Jvec_hyr_Hform(self): @@ -149,6 +229,32 @@ class FDEM_DerivTests(unittest.TestCase): def test_Jvec_hzi_Hform(self): self.assertTrue(derivTest('h', 'jzi')) + def test_Jvec_exr_Hform(self): + self.assertTrue(derivTest('h', 'exr')) + def test_Jvec_eyr_Hform(self): + self.assertTrue(derivTest('h', 'eyr')) + def test_Jvec_ezr_Hform(self): + self.assertTrue(derivTest('h', 'ezr')) + def test_Jvec_exi_Hform(self): + self.assertTrue(derivTest('h', 'exi')) + def test_Jvec_eyi_Hform(self): + self.assertTrue(derivTest('h', 'eyi')) + def test_Jvec_ezi_Hform(self): + self.assertTrue(derivTest('h', 'ezi')) + + def test_Jvec_bxr_Hform(self): + self.assertTrue(derivTest('h', 'bxr')) + def test_Jvec_byr_Hform(self): + self.assertTrue(derivTest('h', 'byr')) + def test_Jvec_bzr_Hform(self): + self.assertTrue(derivTest('h', 'bzr')) + def test_Jvec_bxi_Hform(self): + self.assertTrue(derivTest('h', 'bxi')) + def test_Jvec_byi_Hform(self): + self.assertTrue(derivTest('h', 'byi')) + def test_Jvec_bzi_Hform(self): + self.assertTrue(derivTest('h', 'bzi')) + if __name__ == '__main__': unittest.main() 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'