diff --git a/.travis.yml b/.travis.yml index 258c4b7b..bec55e9c 100644 --- a/.travis.yml +++ b/.travis.yml @@ -33,7 +33,7 @@ before_install: # Install packages install: - - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose + - conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython ipython nose vtk - pip install nose-cov python-coveralls - git clone https://github.com/rowanc1/pymatsolver.git diff --git a/SimPEG/DataMisfit.py b/SimPEG/DataMisfit.py index c1203d47..425fe4ce 100644 --- a/SimPEG/DataMisfit.py +++ b/SimPEG/DataMisfit.py @@ -59,20 +59,6 @@ class BaseDataMisfit(object): """ raise NotImplementedError('This method should be overwritten.') - # TODO: implement target misfit as a property, or possibly as an inversion directive. - - # def target(self, forward): - # """target(forward) - - # Target for data misfit. By default this is the number of data, - # which satisfies the Discrepancy Principle. - - # :rtype: float - # :return: data misfit target - - # """ - # prob, survey = self.splitForward(forward) - # return survey.nD class l2_DataMisfit(BaseDataMisfit): @@ -103,10 +89,18 @@ class l2_DataMisfit(BaseDataMisfit): """ if getattr(self, '_Wd', None) is None: - print 'SimPEG.l2_DataMisfit is creating default weightings for Wd.' + survey = self.survey - eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5 - self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+eps)) + + if getattr(survey,'std', None) is None: + print 'SimPEG.DataMisfit.l2_DataMisfit assigning default std of 5%' + survey.std = 0.05 + + if getattr(survey, 'eps', None) is None: + print 'SimPEG.DataMisfit.l2_DataMisfit assigning default eps of 1e-5 * ||dobs||' + survey.eps = np.linalg.norm(Utils.mkvc(survey.dobs),2)*1e-5 + + self._Wd = Utils.sdiag(1/(abs(survey.dobs)*survey.std+survey.eps)) return self._Wd @Wd.setter diff --git a/SimPEG/Directives.py b/SimPEG/Directives.py index 48d7abcf..46576df5 100644 --- a/SimPEG/Directives.py +++ b/SimPEG/Directives.py @@ -206,6 +206,36 @@ class SaveOutputEveryIteration(_SaveEveryIteration): f.write(' %3d %1.4e %1.4e %1.4e %1.4e\n'%(self.opt.iter, self.invProb.beta, self.invProb.phi_d, self.invProb.phi_m, self.opt.f)) f.close() +class SaveOutputDictEveryIteration(_SaveEveryIteration): + """SaveOutputDictEveryIteration""" + + def initialize(self): + print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '###-%s.npz'"%self.fileName + + def endIter(self): + # Save the data. + ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) ) + phi_ms = 0.5*ms.dot(ms) + if self.reg.smoothModel == True: + mref = self.reg.mref + else: + mref = 0 + mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) ) + phi_mx = 0.5 * mx.dot(mx) + if self.prob.mesh.dim==2: + my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) ) + phi_my = 0.5 * my.dot(my) + else: + phi_my = 'NaN' + if self.prob.mesh.dim==3: + mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) ) + phi_mz = 0.5 * mz.dot(mz) + else: + phi_mz = 'NaN' + + + # Save the file as a npz + np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index f2167fd8..4b137b2c 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -15,18 +15,20 @@ class BaseFDEMProblem(BaseEMProblem): .. math :: \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ - {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}} + {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} if using the E-B formulation (:code:`Problem_e` - or :code:`Problem_b`) or the magnetic field + or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + + If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: - \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\ + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). + if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ @@ -36,7 +38,11 @@ class BaseFDEMProblem(BaseEMProblem): def fields(self, m=None): """ - Solve the forward problem for the fields. + Solve the forward problem for the fields. + + :param numpy.array m: inversion model (nP,) + :rtype numpy.array: + :return F: forward solution """ self.curModel = m @@ -50,16 +56,22 @@ class BaseFDEMProblem(BaseEMProblem): Srcs = self.survey.getSrcByFreq(freq) ftype = self._fieldType + 'Solution' F[Srcs, ftype] = sol - + Ainv.clean() return F - def Jvec(self, m, v, f=None): + def Jvec(self, m, v, u=None): """ - Sensitivity times a vector + 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 + :rtype numpy.array: + :return: Jv (ndata,) """ - if f is None: - f = self.fields(m) + if u is None: + u = self.fields(m) self.curModel = m @@ -71,33 +83,41 @@ class BaseFDEMProblem(BaseEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = f[src, ftype] + 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 ) for rx in src.rxList: - df_duFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None) df_dudu_dm = df_duFun(src, du_dm, adjoint=False) - df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + 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.projectFieldsDeriv(src, self.mesh, f, v) # wrt u, also have wrt m + P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m Jv[src, rx] = P(Df_Dm) + Ainv.clean() return Utils.mkvc(Jv) - def Jtvec(self, m, v, f=None): + def Jtvec(self, m, v, u=None): """ - Sensitivity transpose times a vector + 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 + :rtype numpy.array: + :return: Jv (ndata,) """ - if f is None: - f = self.fields(m) + if u is None: + u = self.fields(m) self.curModel = m @@ -113,12 +133,12 @@ class BaseFDEMProblem(BaseEMProblem): for src in self.survey.getSrcByFreq(freq): ftype = self._fieldType + 'Solution' - u_src = f[src, ftype] + u_src = u[src, ftype] for rx in src.rxList: - PTv = rx.projectFieldsDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt u, need possibility wrt m + PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m - df_duTFun = getattr(f, '_%sDeriv_u'%rx.projField, None) + df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None) df_duT = df_duTFun(src, PTv, adjoint=True) ATinvdf_duT = ATinv * df_duT @@ -127,11 +147,12 @@ class BaseFDEMProblem(BaseEMProblem): dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True) du_dmT = -dA_dmT + dRHS_dmT - df_dmFun = getattr(f, '_%sDeriv_m'%rx.projField, None) + df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None) dfT_dm = df_dmFun(src, PTv, adjoint=True) du_dmT += dfT_dm + # 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 @@ -139,16 +160,18 @@ class BaseFDEMProblem(BaseEMProblem): Jtv += - np.array(du_dmT,dtype=complex).real else: raise Exception('Must be real or imag') + + ATinv.clean() - return Jtv + return Utils.mkvc(Jtv) def getSourceTerm(self, freq): """ - Evaluates the sources for a given frequency and puts them in matrix form + Evaluates the sources for a given frequency and puts them in matrix form - :param float freq: Frequency - :rtype: numpy.ndarray (nE or nF, nSrc) - :return: S_m, S_e + :param float freq: Frequency + :rtype: (numpy.ndarray, numpy.ndarray) + :return: S_m, S_e (nE or nF, nSrc) """ Srcs = self.survey.getSrcByFreq(freq) if self._eqLocs is 'FE': @@ -172,20 +195,22 @@ class BaseFDEMProblem(BaseEMProblem): class Problem_e(BaseFDEMProblem): """ - By eliminating the magnetic flux density using - - .. math :: - - \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right) - - - we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only: + By eliminating the magnetic flux density using .. math :: - \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e} + \mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} + \mathbf{s_m}\\right) - which we solve for \\\(\\\mathbf{e}\\\). + + we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only: + + .. math :: + + \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C}+ i \omega \mathbf{M^e_{\sigma}} \\right)\mathbf{e} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M^e}\mathbf{s_e} + + which we solve for :math:`\mathbf{e}`. + + :param SimPEG.Mesh mesh: mesh """ _fieldType = 'e' @@ -197,13 +222,16 @@ class Problem_e(BaseFDEMProblem): def getA(self, freq): """ - .. math :: - \mathbf{A} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} + System matrix + + .. math :: + \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{C} + i \omega \mathbf{M^e_{\sigma}} - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A """ + MfMui = self.MfMui MeSigma = self.MeSigma C = self.mesh.edgeCurl @@ -212,6 +240,20 @@ class Problem_e(BaseFDEMProblem): def getADeriv_m(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 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,) + """ + dsig_dm = self.curModel.sigmaDeriv dMe_dsig = self.MeSigmaDeriv(u) @@ -222,26 +264,37 @@ class Problem_e(BaseFDEMProblem): def getRHS(self, freq): """ - .. math :: - \mathbf{RHS} = \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f}\mathbf{s_m} -i\omega\mathbf{M_e}\mathbf{s_e} + Right hand side for the system - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nSrc) - :return: RHS + .. 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 + :return: RHS (nE, nSrc) """ S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MfMui = self.MfMui - RHS = C.T * (MfMui * S_m) -1j * omega(freq) * S_e - - return RHS + return C.T * (MfMui * S_m) -1j * omega(freq) * S_e def getRHSDeriv_m(self, freq, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + + :param float freq: frequency + :param SimPEG.EM.FDEM.Src src: FDEM source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of rhs deriv with a vector + """ + C = self.mesh.edgeCurl MfMui = self.MfMui - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: dRHS = MfMui * (C * v) @@ -253,20 +306,22 @@ class Problem_e(BaseFDEMProblem): class Problem_b(BaseFDEMProblem): """ - We eliminate \\\(\\\mathbf{e}\\\) using + We eliminate :math:`\mathbf{e}` using - .. math :: + .. math :: - \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right) + \mathbf{e} = \mathbf{M^e_{\sigma}}^{-1} \\left(\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{s_e}\\right) - and solve for \\\(\\\mathbf{b}\\\) using: + and solve for :math:`\mathbf{b}` using: - .. math :: + .. math :: - \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e} + \\left(\mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega \\right)\mathbf{b} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{M^e}\mathbf{s_e} - .. note :: - The inverse problem will not work with full anisotropy + .. note :: + The inverse problem will not work with full anisotropy + + :param SimPEG.Mesh mesh: mesh """ _fieldType = 'b' @@ -278,12 +333,14 @@ class Problem_b(BaseFDEMProblem): def getA(self, freq): """ - .. math :: - \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} + i \omega + System matrix - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A + .. math :: + \mathbf{A} = \mathbf{C} \mathbf{M^e_{\sigma}}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} + i \omega + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A """ MfMui = self.MfMui @@ -299,6 +356,20 @@ class Problem_b(BaseFDEMProblem): def getADeriv_m(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} \\frac{\mathbf{M^e_{\sigma}} \mathbf{v}}{d\mathbf{m}} + + :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,) + """ + MfMui = self.MfMui C = self.mesh.edgeCurl MeSigmaIDeriv = self.MeSigmaIDeriv @@ -318,12 +389,14 @@ class Problem_b(BaseFDEMProblem): def getRHS(self, freq): """ - .. math :: - \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e} + Right hand side for the system - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nSrc) - :return: RHS + .. math :: + \mathbf{RHS} = \mathbf{s_m} + \mathbf{M^e_{\sigma}}^{-1}\mathbf{s_e} + + :param float freq: Frequency + :rtype: numpy.ndarray + :return: RHS (nE, nSrc) """ S_m, S_e = self.getSourceTerm(freq) @@ -339,6 +412,17 @@ class Problem_b(BaseFDEMProblem): return RHS def getRHSDeriv_m(self, freq, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + + :param float freq: frequency + :param SimPEG.EM.FDEM.Src src: FDEM source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of rhs deriv with a vector + """ + C = self.mesh.edgeCurl S_m, S_e = src.eval(self) MfMui = self.MfMui @@ -347,7 +431,7 @@ class Problem_b(BaseFDEMProblem): v = self.MfMui * v MeSigmaIDeriv = self.MeSigmaIDeriv(S_e) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if not adjoint: RHSderiv = C * (MeSigmaIDeriv * v) @@ -370,21 +454,22 @@ class Problem_b(BaseFDEMProblem): class Problem_j(BaseFDEMProblem): """ - We eliminate \\\(\\\mathbf{h}\\\) using + We eliminate \\\(\\\mathbf{h}\\\) using - .. math :: + .. math :: - \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right) + \mathbf{h} = \\frac{1}{i \omega} \mathbf{M_{\mu}^e}^{-1} \\left(-\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + \mathbf{M^e} \mathbf{s_m} \\right) - and solve for \\\(\\\mathbf{j}\\\) using + and solve for \\\(\\\mathbf{j}\\\) using - .. math :: + .. math :: - \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^T \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e} + \\left(\mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} + i \omega\\right)\mathbf{j} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1} \mathbf{M^e} \mathbf{s_m} -i\omega\mathbf{s_e} - .. note:: - This implementation does not yet work with full anisotropy!! + .. note:: + This implementation does not yet work with full anisotropy!! + :param SimPEG.Mesh mesh: mesh """ _fieldType = 'j' @@ -396,12 +481,14 @@ class Problem_j(BaseFDEMProblem): def getA(self, freq): """ - .. math :: - \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C}^T \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega + System matrix - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A + .. math :: + \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{\\mu^{-1}}} \\mathbf{C}^{\\top} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A """ MeMuI = self.MeMuI @@ -418,12 +505,20 @@ class Problem_j(BaseFDEMProblem): def getADeriv_m(self, freq, u, v, adjoint=False): """ - In this case, we assume that electrical conductivity, \\\(\\\sigma\\\) is the physical property of interest (i.e. \\\(\\\sigma\\\) = model.transform). Then we want + Product of the derivative of our system matrix with respect to the model and a vector - .. math :: + In this case, we assume that electrical conductivity, :math:`\sigma` is the physical property of interest (i.e. :math:`\sigma` = model.transform). Then we want - \\frac{\mathbf{A(\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{m}} - &= \\mathbf{C} \\mathbf{M^e_{mu}^{-1}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{d \\mathbf{\\sigma^{-1}}} \\frac{d \\mathbf{\\sigma^{-1}}}{d \\mathbf{\\sigma}} \\frac{d \\mathbf{\\sigma}}{d \\mathbf{m}} + .. math :: + + \\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 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,) """ MeMuI = self.MeMuI @@ -443,12 +538,15 @@ class Problem_j(BaseFDEMProblem): def getRHS(self, freq): """ - .. math :: + Right hand side for the system - \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e} - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nSrc) - :return: RHS + .. math :: + + \mathbf{RHS} = \mathbf{C} \mathbf{M_{\mu}^e}^{-1}\mathbf{s_m} -i\omega \mathbf{s_e} + + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nSrc) + :return: RHS """ S_m, S_e = self.getSourceTerm(freq) @@ -463,9 +561,20 @@ class Problem_j(BaseFDEMProblem): return RHS def getRHSDeriv_m(self, freq, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + + :param float freq: frequency + :param SimPEG.EM.FDEM.Src src: FDEM source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of rhs deriv with a vector + """ + C = self.mesh.edgeCurl MeMuI = self.MeMuI - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) if adjoint: if self._makeASymmetric: @@ -486,18 +595,19 @@ class Problem_j(BaseFDEMProblem): class Problem_h(BaseFDEMProblem): """ - We eliminate \\\(\\\mathbf{j}\\\) using + We eliminate \\\(\\\mathbf{j}\\\) using - .. math :: + .. math :: - \mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e} + \mathbf{j} = \mathbf{C} \mathbf{h} - \mathbf{s_e} - and solve for \\\(\\\mathbf{h}\\\) using + and solve for \\\(\\\mathbf{h}\\\) using - .. math :: + .. math :: - \\left(\mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} + \\left(\mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e}\\right) \mathbf{h} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{s_e} + :param SimPEG.Mesh mesh: mesh """ _fieldType = 'h' @@ -509,13 +619,14 @@ class Problem_h(BaseFDEMProblem): def getA(self, freq): """ - .. math :: + System matrix - \mathbf{A} = \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e} + .. math:: + \mathbf{A} = \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{C} + i \omega \mathbf{M_{\mu}^e} - :param float freq: Frequency - :rtype: scipy.sparse.csr_matrix - :return: A + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A """ MeMu = self.MeMu @@ -525,6 +636,19 @@ class Problem_h(BaseFDEMProblem): return C.T * (MfRho * C) + 1j*omega(freq)*MeMu def getADeriv_m(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 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,) + """ MeMu = self.MeMu C = self.mesh.edgeCurl @@ -536,24 +660,35 @@ class Problem_h(BaseFDEMProblem): def getRHS(self, freq): """ - .. math :: + Right hand side for the system - \mathbf{RHS} = \mathbf{M^e} \mathbf{s_m} + \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{s_e} + .. math :: - :param float freq: Frequency - :rtype: numpy.ndarray (nE, nSrc) - :return: RHS + \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 + :return: RHS (nE, nSrc) """ S_m, S_e = self.getSourceTerm(freq) C = self.mesh.edgeCurl MfRho = self.MfRho - RHS = S_m + C.T * ( MfRho * S_e ) - - return RHS + return S_m + C.T * ( MfRho * S_e ) def getRHSDeriv_m(self, freq, src, v, adjoint=False): + """ + Derivative of the right hand side with respect to the model + + :param float freq: frequency + :param SimPEG.EM.FDEM.Src src: FDEM source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of rhs deriv with a vector + """ + _, S_e = src.eval(self) C = self.mesh.edgeCurl MfRho = self.MfRho @@ -564,7 +699,7 @@ class Problem_h(BaseFDEMProblem): elif adjoint: RHSDeriv = MfRhoDeriv.T * (C * v) - S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint) + S_mDeriv, S_eDeriv = src.evalDeriv(self, adjoint=adjoint) return RHSDeriv + S_mDeriv(v) + C.T * (MfRho * S_eDeriv(v)) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 8f6fafe9..e171a5c5 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -7,11 +7,39 @@ from SimPEG.Utils import Zero, Identity class Fields(SimPEG.Problem.Fields): - """Fancy Field Storage for a FDEM survey.""" + """ + + 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 + + .. code-block:: python + + f = problem.fields(m) + e = f[srcList,'e'] + b = f[srcList,'b'] + + If accessing all sources for a given field, use the :code:`:` + + .. code-block:: python + + f = problem.fields(m) + e = f[:,'e'] + b = f[:,'b'] + + The array returned will be size (nE or nF, nSrcs :math:`\\times` nFrequencies) + """ + knownFields = {} dtype = complex class Fields_e(Fields): + """ + Fields object for Problem_e. + + :param Mesh mesh: mesh + :param Survey survey: survey + """ + knownFields = {'eSolution':'E'} aliasFields = { 'e' : ['eSolution','E','_e'], @@ -30,6 +58,15 @@ class Fields_e(Fields): self._edgeCurl = self.survey.prob.mesh.edgeCurl def _ePrimary(self, eSolution, srcList): + """ + Primary electric field from source + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary electric field as defined by the sources + """ + ePrimary = np.zeros_like(eSolution) for i, src in enumerate(srcList): ep = src.ePrimary(self.prob) @@ -37,19 +74,67 @@ class Fields_e(Fields): return ePrimary def _eSecondary(self, eSolution, srcList): + """ + Secondary electric field is the thing we solved for + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :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 + + :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 field we solved for with a vector + """ + return Identity()*v 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. + + :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 electric field derivative with respect to the inversion model with a vector + """ + # assuming primary does not depend on the model return Zero() def _bPrimary(self, eSolution, srcList): + """ + Primary magnetic flux density from source + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic flux density as defined by the sources + """ + bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) @@ -57,6 +142,15 @@ class Fields_e(Fields): return bPrimary def _bSecondary(self, eSolution, srcList): + """ + Secondary magnetic flux density from eSolution + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic flux density + """ + C = self._edgeCurl b = (C * eSolution) for i, src in enumerate(srcList): @@ -66,29 +160,84 @@ class Fields_e(Fields): return b def _bSecondaryDeriv_u(self, src, v, adjoint = False): + """ + Derivative of the secondary magnetic flux density with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray 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): - S_mDeriv, _ = src.evalDeriv(self.prob, adjoint) - S_mDeriv = S_mDeriv(v) + """ + 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 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) def _bDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the total magnetic flux density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: SimPEG.Utils.Zero + :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) class Fields_b(Fields): + """ + Fields object for Problem_b. + + :param Mesh mesh: mesh + :param Survey survey: survey + """ + knownFields = {'bSolution':'F'} aliasFields = { 'b' : ['bSolution','F','_b'], @@ -111,6 +260,15 @@ class Fields_b(Fields): self._Me = self.survey.prob.Me def _bPrimary(self, bSolution, srcList): + """ + Primary magnetic flux density from source + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary electric field as defined by the sources + """ + bPrimary = np.zeros_like(bSolution) for i, src in enumerate(srcList): bp = src.bPrimary(self.prob) @@ -118,19 +276,66 @@ class Fields_b(Fields): return bPrimary def _bSecondary(self, bSolution, srcList): + """ + Secondary magnetic flux density is the thing we solved for + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic flux density + """ + 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 + """ + 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 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 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. + + :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 magnetic flux density derivative with respect to the inversion model with a vector + """ + # assuming primary does not depend on the model return Zero() def _ePrimary(self, bSolution, srcList): + """ + Primary electric field from source + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary electric field as defined by the sources + """ + ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex) for i,src in enumerate(srcList): ep = src.ePrimary(self.prob) @@ -138,6 +343,15 @@ class Fields_b(Fields): return ePrimary def _eSecondary(self, bSolution, srcList): + """ + Secondary electric field from bSolution + + :param numpy.ndarray bSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary electric field + """ + e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution)) for i,src in enumerate(srcList): _,S_e = src.eval(self.prob) @@ -145,12 +359,32 @@ class Fields_b(Fields): return 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 + """ + 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 @@ -166,25 +400,60 @@ class Fields_b(Fields): elif adjoint: de_dm = self._MeSigmaIDeriv(w).T * v - _, S_eDeriv = src.evalDeriv(self.prob, adjoint) - Se_Deriv = S_eDeriv(v) + _, S_eDeriv = src.evalDeriv(self.prob, v, adjoint) - de_dm = de_dm - self._MeSigmaI * Se_Deriv + de_dm = de_dm - self._MeSigmaI * S_eDeriv return de_dm 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? + :rtype: numpy.ndarray + :return: product of the derivative of the electric field with respect to the field we solved for with a vector + """ + return self._eSecondaryDeriv_u(src, v, adjoint) def _eDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the total electric field density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the electric field derivative with respect to the inversion model with a vector + """ + # assuming primary doesn't depend on model return self._eSecondaryDeriv_m(src, v, adjoint) class Fields_j(Fields): + """ + Fields object for Problem_j. + + :param Mesh mesh: mesh + :param Survey survey: survey + """ + knownFields = {'jSolution':'F'} aliasFields = { 'j' : ['jSolution','F','_j'], @@ -207,6 +476,15 @@ class Fields_j(Fields): self._Me = self.survey.prob.Me def _jPrimary(self, jSolution, srcList): + """ + Primary current density from source + + :param numpy.ndarray jSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current density as defined by the sources + """ + jPrimary = np.zeros_like(jSolution,dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) @@ -214,19 +492,66 @@ class Fields_j(Fields): return jPrimary def _jSecondary(self, jSolution, srcList): + """ + Secondary current density is the thing we solved for + + :param numpy.ndarray jSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary current density + """ + return jSolution def _j(self, jSolution, srcList): + """ + 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 self._jPrimary(jSolution, srcList) + self._jSecondary(jSolution, 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 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 Identity()*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. + + :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 derivative with respect to the inversion model with a vector + """ # assuming primary does not depend on the model return Zero() def _hPrimary(self, jSolution, srcList): + """ + Primary magnetic field from source + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic field as defined by the sources + """ + hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex) for i, src in enumerate(srcList): hp = src.hPrimary(self.prob) @@ -234,6 +559,15 @@ class Fields_j(Fields): return hPrimary def _hSecondary(self, jSolution, srcList): + """ + Secondary magnetic field from bSolution + + :param numpy.ndarray jSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic field + """ + h = self._MeMuI * (self._edgeCurl.T * (self._MfRho * jSolution) ) for i, src in enumerate(srcList): h[:,i] *= -1./(1j*omega(src.freq)) @@ -242,12 +576,32 @@ class Fields_j(Fields): return h def _hSecondaryDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the secondary magnetic field with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray 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 @@ -260,7 +614,7 @@ class Fields_j(Fields): elif adjoint: hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) ) - S_mDeriv,_ = src.evalDeriv(self.prob, adjoint) + S_mDeriv,_ = src.evalDeriv(self.prob, adjoint = adjoint) if not adjoint: S_mDeriv = S_mDeriv(v) @@ -272,17 +626,53 @@ class Fields_j(Fields): 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 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) def _hDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the total magnetic field density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of the magnetic field derivative with respect to the inversion model with a vector + """ + # assuming the primary doesn't depend on the model return self._hSecondaryDeriv_m(src, v, adjoint) class Fields_h(Fields): + """ + Fields object for Problem_h. + + :param Mesh mesh: mesh + :param Survey survey: survey + """ + knownFields = {'hSolution':'E'} aliasFields = { 'h' : ['hSolution','E','_h'], @@ -303,6 +693,15 @@ class Fields_h(Fields): self._MfRho = self.survey.prob.MfRho def _hPrimary(self, hSolution, srcList): + """ + Primary magnetic field from source + + :param numpy.ndarray eSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary magnetic field as defined by the sources + """ + hPrimary = np.zeros_like(hSolution,dtype = complex) for i, src in enumerate(srcList): hp = src.hPrimary(self.prob) @@ -310,19 +709,67 @@ class Fields_h(Fields): return hPrimary def _hSecondary(self, hSolution, srcList): + """ + Secondary magnetic field is the thing we solved for + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary magnetic field + """ + 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): + """ + 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 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 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. + + :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 magnetic field derivative with respect to the inversion model with a vector + """ + # assuming primary does not depend on the model return Zero() def _jPrimary(self, hSolution, srcList): + """ + Primary current density from source + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: primary current density as defined by the sources + """ + jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]], dtype = complex) for i, src in enumerate(srcList): jp = src.jPrimary(self.prob) @@ -330,6 +777,15 @@ class Fields_h(Fields): return jPrimary def _jSecondary(self, hSolution, srcList): + """ + Secondary current density from eSolution + + :param numpy.ndarray hSolution: field we solved for + :param list srcList: list of sources + :rtype: numpy.ndarray + :return: secondary current density + """ + j = self._edgeCurl*hSolution for i, src in enumerate(srcList): _,S_e = src.eval(self.prob) @@ -337,22 +793,69 @@ class Fields_h(Fields): return j def _jSecondaryDeriv_u(self, src, v, adjoint=False): + """ + Derivative of the secondary current density with respect to the thing we solved for + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray 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): - _,S_eDeriv = src.evalDeriv(self.prob, adjoint) - S_eDeriv = S_eDeriv(v) + """ + 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 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) def _jDeriv_m(self, src, v, adjoint=False): + """ + Derivative of the total current density with respect to the inversion model. + + :param SimPEG.EM.FDEM.Src src: source + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: SimPEG.Utils.Zero + :return: product of the current density 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) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index b29768ac..1213cef3 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -1,55 +1,141 @@ 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 SurveyFDEM import Rx - +from SimPEG.Utils import Zero class BaseSrc(Survey.BaseSrc): + """ + Base source class for FDEM Survey + """ + freq = None - # rxPair = Rx + # rxPair = RxFDEM integrate = True def eval(self, prob): + """ + 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 - def evalDeriv(self, prob, v, adjoint=False): - return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) + 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 + + :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 + """ + 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) def bPrimary(self, prob): + """ + Primary magnetic flux density + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: primary magnetic flux density + """ return Zero() def hPrimary(self, prob): + """ + Primary magnetic field + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ return Zero() def ePrimary(self, prob): + """ + Primary electric field + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: primary electric field + """ return Zero() def jPrimary(self, prob): + """ + Primary current density + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: primary current density + """ return Zero() def S_m(self, prob): + """ + Magnetic source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: magnetic source term on mesh + """ return Zero() def S_e(self, prob): + """ + Electric source term + + :param Problem prob: FDEM Problem + :rtype: numpy.ndarray + :return: electric source term on mesh + """ return Zero() def S_mDeriv(self, prob, v, adjoint = False): + """ + Derivative of magnetic source term with respect to the inversion model + + :param Problem prob: FDEM Problem + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of magnetic source term derivative with a vector + """ + return Zero() def S_eDeriv(self, prob, v, adjoint = False): + """ + Derivative of electric source term with respect to the inversion model + + :param Problem prob: FDEM Problem + :param numpy.ndarray v: vector to take product with + :param bool adjoint: adjoint? + :rtype: numpy.ndarray + :return: product of electric source term derivative with a vector + """ return Zero() 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 numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list + :param list rxList: receiver list + :param float freq: frequency + :param numpy.array S_e: electric source term """ def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): @@ -58,16 +144,17 @@ class RawVec_e(BaseSrc): 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 + RawVec magnetic source. It is defined by the user provided vector S_m - :param numpy.array S_m: magnetic source term - :param float freq: frequency - :param rxList: receiver list + :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()): @@ -78,17 +165,24 @@ class RawVec_m(BaseSrc): 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 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 numpy.array S_m: magnetic source term - :param numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list + :param rxList: receiver list + :param float freq: frequency + :param numpy.array S_m: magnetic source term + :param numpy.array S_e: electric source term """ def __init__(self, rxList, freq, S_m, S_e, integrate = True): self._S_m = np.array(S_m,dtype=complex) @@ -109,6 +203,51 @@ class RawVec(BaseSrc): 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!). + + 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:: + \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}} + + We split up the fields and :math:`\mu^{-1}` into primary (:math:`\mathbf{P}`) and secondary (:math:`\mathbf{S}`) components + + - :math:`\mathbf{e} = \mathbf{e^P} + \mathbf{e^S}` + - :math:`\mathbf{b} = \mathbf{b^P} + \mathbf{b^S}` + - :math:`\\boldsymbol{\mu}^{\mathbf{-1}} = \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{P}} + \\boldsymbol{\mu}^{\mathbf{-1}^\mathbf{S}}` + + and define a zero-frequency primary problem, noting that the source is + generated by a divergence free electric current + + .. 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 + + .. 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 + + .. math:: + \mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\ + {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b^S} - \mathbf{M_{\sigma}^e} \mathbf{e^S} = -\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^S} \mathbf{b^P}} + + :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 + """ #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): @@ -121,6 +260,13 @@ class MagDipole(BaseSrc): BaseSrc.__init__(self, rxList) def bPrimary(self, prob): + """ + The primary magnetic flux density from a magnetic vector potential + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -152,14 +298,37 @@ class MagDipole(BaseSrc): return C*a def hPrimary(self, prob): + """ + The primary magnetic field from a magnetic vector potential + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ b = self.bPrimary(prob) return h_from_b(prob,b) def S_m(self, prob): + """ + The magnetic source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ + b_p = self.bPrimary(prob) return -1j*omega(self.freq)*b_p def S_e(self, prob): + """ + The electric source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: @@ -179,6 +348,21 @@ class MagDipole(BaseSrc): class MagDipole_Bfield(BaseSrc): + """ + Point magnetic dipole source calculated with the analytic solution for the + 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. + + :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 + """ + #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): @@ -190,6 +374,14 @@ class MagDipole_Bfield(BaseSrc): BaseSrc.__init__(self, rxList) def bPrimary(self, prob): + """ + The primary magnetic flux density from the analytic solution for magnetic fields from a dipole + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ + eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -221,14 +413,35 @@ class MagDipole_Bfield(BaseSrc): return b def hPrimary(self, prob): + """ + The primary magnetic field from a magnetic vector potential + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ b = self.bPrimary(prob) return h_from_b(prob, b) def S_m(self, prob): + """ + The magnetic source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ b = self.bPrimary(prob) return -1j*omega(self.freq)*b def S_e(self, prob): + """ + The electric source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: @@ -247,6 +460,20 @@ class MagDipole_Bfield(BaseSrc): 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!). + + 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 + """ #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): @@ -259,6 +486,13 @@ class CircularLoop(BaseSrc): BaseSrc.__init__(self, rxList) def bPrimary(self, prob): + """ + The primary magnetic flux density from a magnetic vector potential + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ eqLocs = prob._eqLocs if eqLocs is 'FE': @@ -289,14 +523,35 @@ class CircularLoop(BaseSrc): return C*a def hPrimary(self, prob): + """ + The primary magnetic field from a magnetic vector potential + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ b = self.bPrimary(prob) return 1./self.mu*b def S_m(self, prob): + """ + The magnetic source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ b = self.bPrimary(prob) return -1j*omega(self.freq)*b def S_e(self, prob): + """ + The electric source term + + :param Problem prob: FDEM problem + :rtype: numpy.ndarray + :return: primary magnetic field + """ if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): return Zero() else: diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index f60cbfdf..444df88d 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -10,6 +10,12 @@ import SrcFDEM as Src #################################################### class Rx(SimPEG.Survey.BaseRx): + """ + Frequency domain receivers + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string rxType: reciever type from knownRxTypes + """ knownRxTypes = { 'exr':['e', 'Ex', 'real'], @@ -61,6 +67,15 @@ class Rx(SimPEG.Survey.BaseRx): return self.knownRxTypes[self.rxType][2] def projectFields(self, src, mesh, u): + """ + Project fields to recievers to get data. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields u: fields object + :rtype: numpy.ndarray + :return: fields projected to recievers + """ P = self.getP(mesh) u_part_complex = u[src, self.projField] # get the real or imag component @@ -69,6 +84,16 @@ class Rx(SimPEG.Survey.BaseRx): return P*u_part def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False): + """ + Derivative of projected fields with respect to the inversion model times a vector. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields u: fields object + :param numpy.ndarray v: vector to multiply + :rtype: numpy.ndarray + :return: fields projected to recievers + """ P = self.getP(mesh) if not adjoint: @@ -95,10 +120,13 @@ class Rx(SimPEG.Survey.BaseRx): class Survey(SimPEG.Survey.BaseSurvey): """ - docstring for SurveyFDEM + Frequency domain electromagnetic survey + + :param list srcList: list of FDEM sources used in the survey """ srcPair = Src.BaseSrc + rxPaair = Rx def __init__(self, srcList, **kwargs): # Sort these by frequency @@ -126,6 +154,7 @@ class Survey(SimPEG.Survey.BaseSurvey): @property def nSrcByFreq(self): + """Number of sources at each frequency""" if getattr(self, '_nSrcByFreq', None) is None: self._nSrcByFreq = {} for freq in self.freqs: @@ -133,11 +162,22 @@ class Survey(SimPEG.Survey.BaseSurvey): return self._nSrcByFreq def getSrcByFreq(self, freq): - """Returns the sources associated with a specific frequency.""" + """ + Returns the sources associated with a specific frequency. + :param float freq: frequency for which we look up sources + :rtype: dictionary + :return: sources at the sepcified frequency + """ assert freq in self._freqDict, "The requested frequency is not in this survey." return self._freqDict[freq] def projectFields(self, u): + """ + Project fields to receiver locations + :param Fields u: fields object + :rtype: numpy.ndarray + :return: data + """ data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 4d4b39a7..2efb10ec 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -37,13 +37,21 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem): _FieldsForward_pair = FieldsTDEM #: used for the forward calculation only + waveformType = "STEPOFF" + current = None + + def currentwaveform(self, wave): + self._timeSteps = np.diff(wave[:,0]) + self.current = wave[:,1] + self.waveformType = "GENERAL" + def fields(self, m): if self.verbose: print '%s\nCalculating fields(m)\n%s'%('*'*50,'*'*50) self.curModel = m # Create a fields storage object F = self._FieldsForward_pair(self.mesh, self.survey) for src in self.survey.srcList: - # Set the initial conditions + # Set the initial conditions F[src,:,0] = src.getInitialFields(self.mesh) F = self.forward(m, self.getRHS, F=F) if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50) diff --git a/SimPEG/EM/__init__.py b/SimPEG/EM/__init__.py index 6a1ca774..565f63a8 100644 --- a/SimPEG/EM/__init__.py +++ b/SimPEG/EM/__init__.py @@ -1,6 +1,6 @@ -# from EM import * import TDEM import FDEM import Base import Analytics import Utils +from scipy.constants import mu_0, epsilon_0 diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py new file mode 100644 index 00000000..ff87b6a6 --- /dev/null +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -0,0 +1,116 @@ +from SimPEG import * +import SimPEG.EM as EM +from SimPEG.EM import mu_0 + + +def run(plotIt=True): + """ + EM: FDEM: 1D: Inversion + ======================= + + Here we will create and run a FDEM 1D inversion. + + """ + + cs, ncx, ncz, npad = 5., 25, 15, 15 + hx = [(cs,ncx), (cs,npad,1.3)] + hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hz], '00C') + + layerz = -100. + + active = mesh.vectorCCz<0. + layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz) + actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz) + mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap + sig_half = 2e-2 + sig_air = 1e-8 + sig_layer = 1e-2 + sigma = np.ones(mesh.nCz)*sig_air + sigma[active] = sig_half + sigma[layer] = sig_layer + mtrue = np.log(sigma[active]) + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (3, 6)) + plt.semilogx(sigma[active], mesh.vectorCCz[active]) + ax.set_ylim(-500, 0) + ax.set_xlim(1e-3, 1e-1) + ax.set_xlabel('Conductivity (S/m)', fontsize = 14) + ax.set_ylabel('Depth (m)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + + + rxOffset=10. + bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') + + 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] + + survey = EM.FDEM.Survey(srcList) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + prb.Solver = SolverLU + + prb.pair(survey) + + std = 0.05 + survey.makeSyntheticData(mtrue, std) + + survey.std = std + survey.eps = np.linalg.norm(survey.dtrue)*1e-5 + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (6, 6)) + ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-') + ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-') + ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16) + ax.set_xlabel('Time (s)', fontsize = 14) + ax.set_ylabel('$B_z$ (T)', fontsize = 16) + ax.set_xlabel('Time (s)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + + dmisfit = DataMisfit.l2_DataMisfit(survey) + regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]]) + reg = Regularization.Tikhonov(regMesh) + opt = Optimization.InexactGaussNewton(maxIter = 6) + invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt) + + # Create an inversion object + beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2) + betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0) + inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest]) + m0 = np.log(np.ones(mtrue.size)*sig_half) + reg.alpha_s = 1e-3 + reg.alpha_x = 1. + prb.counter = opt.counter = Utils.Counter() + opt.LSshorten = 0.5 + opt.remember('xc') + + mopt = inv.run(m0) + + if plotIt: + import matplotlib.pyplot as plt + fig, ax = plt.subplots(1,1, figsize = (3, 6)) + plt.semilogx(sigma[active], mesh.vectorCCz[active]) + plt.semilogx(np.exp(mopt), mesh.vectorCCz[active]) + ax.set_ylim(-500, 0) + ax.set_xlim(1e-3, 1e-1) + ax.set_xlabel('Conductivity (S/m)', fontsize = 14) + ax.set_ylabel('Depth (m)', fontsize = 14) + ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) + plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best') + plt.show() + + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/EM_FDEM_SusEffects.py b/SimPEG/Examples/EM_FDEM_SusEffects.py new file mode 100644 index 00000000..5021e87f --- /dev/null +++ b/SimPEG/Examples/EM_FDEM_SusEffects.py @@ -0,0 +1,145 @@ +from SimPEG import * +from SimPEG import EM +from pymatsolver import MumpsSolver +from scipy.constants import mu_0 + +def run(plotIt=True): + """ + FDEM: Effects of susceptibility + =============================== + + When airborne freqeuncy domain EM (AFEM) survey is flown over + the earth including significantly susceptible bodies (magnetite-rich rocks), + negative data is often observed in the real part of the lowest frequency + (e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization + occurs due to a susceptible body when the magnetic field applied. + + To clarify what is happening in the earth when we are exciting the earth with + a loop source in the frequency domain we run three forward modelling: + + - F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility + - F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity + - F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field + + We plot vector magnetic fields in the earth. For secondary fields we provide + F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following + figure show only real part, since that is our interest. + + """ + # Generate Cylindrical mesh + cs, ncx, ncz, npad = 5, 25, 24, 20. + hx = [(cs,ncx), (cs,npad,1.3)] + hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] + mesh = Mesh.CylMesh([hx,1,hz], '00C') + sighalf = 1e-3 + sigma = np.ones(mesh.nC)*1e-8 + sigmahomo = sigma.copy() + mu = np.ones(mesh.nC)*mu_0 + sigma[mesh.gridCC[:,-1]<0.] = sighalf + blkind = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50)) + sigma[blkind] = 1e-1 + mu[blkind] = mu_0*1.1 + offset = 0. + frequency = np.r_[10., 100., 1000.] + rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr') + rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi') + srcLists = [] + nfreq = frequency.size + for ifreq in range(nfreq): + src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.) + srcLists.append(src) + survey = EM.FDEM.Survey(srcLists) + iMap = Maps.IdentityMap(nP=int(mesh.nC)) + # Use PhysPropMap + maps = [('sigma', iMap), ('mu', iMap)] + prob = EM.FDEM.Problem_b(mesh, mapping=maps) + prob.Solver = MumpsSolver + survey.pair(prob) + m = np.r_[sigma, mu] + survey0 = EM.FDEM.Survey(srcLists) + prob0 = EM.FDEM.Problem_b(mesh, mapping=maps) + prob0.Solver = MumpsSolver + survey0.pair(prob0) + m = np.r_[sigma, mu] + m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0] + m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0] + # Anomalous conductivity and susceptibility + F = prob.fields(m) + # Only anomalous conductivity + F0 = prob.fields(m0) + # Primary field + F00 = prob.fields(m00) + + if plotIt: + import matplotlib.pyplot as plt + def vizfields(ifreq=0, primsec="secondary",realimag="real"): + + titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"] + actind = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200)) + + if primsec=="secondary": + bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') + bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim + bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim + elif primsec=="primary": + bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') + bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F') + + XYZ = mesh.gridCC[actind,:] + X = XYZ[:,0].reshape((31,43), order='F') + Z = XYZ[:,2].reshape((31,43), order='F') + bx = bCC[actind,0].reshape((31,43), order='F') + bz = bCC[actind,1].reshape((31,43), order='F') + bx0 = bCC0[actind,0].reshape((31,43), order='F') + bz0 = bCC0[actind,1].reshape((31,43), order='F') + + bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F') + bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F') + + absbreal = np.sqrt(bx.real**2+bz.real**2) + absbimag = np.sqrt(bx.imag**2+bz.imag**2) + absb0real = np.sqrt(bx0.real**2+bz0.real**2) + absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2) + + absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2) + absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2) + + fig = plt.figure(figsize=(15,5)) + ax1 = plt.subplot(131) + ax2 = plt.subplot(132) + ax3 = plt.subplot(133) + typefield="real" + scale=20 + if realimag=="real": + ax1.contourf(X, Z,np.log10(absbreal), 100) + ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5) + ax2.contourf(X, Z,np.log10(absb0real), 100) + ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5) + ax3.contourf(X, Z,np.log10(absbrealsec), 100) + ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5) + elif realimag=="imag": + ax1.contourf(X, Z,np.log10(absbimag), 100) + ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5) + ax2.contourf(X, Z,np.log10(absb0imag), 100) + ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5) + ax3.contourf(X, Z,np.log10(absbimagsec), 100) + ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5) + + ax = [ax1, ax2, ax3] + ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18) + ax2.text(30, 50, primsec, color="k", fontsize=18) + for i, axtemp in enumerate(ax): + axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3) + + axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3) + axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3) + axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3) + axtemp.set_ylim(-200, 100.) + axtemp.set_xlim(10, 100.) + axtemp.set_title(titles[i]) + plt.show() + vizfields(1, primsec="primary", realimag="real") + vizfields(1, primsec="secondary", realimag="real") + +if __name__ == '__main__': + run() diff --git a/SimPEG/Examples/EM_TDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py index d4d80e55..65ae6669 100644 --- a/SimPEG/Examples/EM_TDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -1,6 +1,6 @@ from SimPEG import * import SimPEG.EM as EM -from scipy.constants import mu_0 +from SimPEG.EM import mu_0 def run(plotIt=True): @@ -50,20 +50,18 @@ def run(plotIt=True): prb.Solver = SolverLU prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)] prb.pair(survey) - dtrue = survey.dpred(mtrue) - - survey.dtrue = dtrue + # create observed data std = 0.05 - noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape) - survey.dobs = survey.dtrue+noise - survey.std = survey.dobs*0 + std - survey.Wd = 1/(abs(survey.dobs)*std) + + survey.dobs = survey.makeSyntheticData(mtrue,std) + survey.std = std + survey.eps = 1e-5*np.linalg.norm(survey.dobs) if plotIt: import matplotlib.pyplot as plt fig, ax = plt.subplots(1,1, figsize = (10, 6)) - ax.loglog(rx.times, dtrue, 'b.-') + ax.loglog(rx.times, survey.dtrue, 'b.-') ax.loglog(rx.times, survey.dobs, 'r.-') ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16) ax.set_xlabel('Time (s)', fontsize = 14) @@ -76,6 +74,7 @@ def run(plotIt=True): reg = Regularization.Tikhonov(regMesh) opt = Optimization.InexactGaussNewton(maxIter = 5) invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt) + # Create an inversion object beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2) betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0) diff --git a/SimPEG/InvProblem.py b/SimPEG/InvProblem.py index 32a2195c..0296bf4b 100644 --- a/SimPEG/InvProblem.py +++ b/SimPEG/InvProblem.py @@ -66,8 +66,8 @@ class BaseInvProblem(object): self.curModel = m0 print """SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. - ***Done using same solver as the problem***""" - self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel)) + ***Done using same Solver and solverOpts as the problem***""" + self.opt.bfgsH0 = self.prob.Solver(self.reg.eval2Deriv(self.curModel), **self.prob.solverOpts) @property def warmstart(self): diff --git a/SimPEG/Maps.py b/SimPEG/Maps.py index 5b4782ac..a3c76e6a 100644 --- a/SimPEG/Maps.py +++ b/SimPEG/Maps.py @@ -10,21 +10,25 @@ class IdentityMap(object): SimPEG Map """ - __metaclass__ = Utils.SimPEGMetaClass - mesh = None #: A SimPEG Mesh - - def __init__(self, mesh, **kwargs): + def __init__(self, mesh=None, nP=None, **kwargs): Utils.setKwargs(self, **kwargs) + + if nP is not None: + assert type(nP) in [int, long], ' Number of parameters must be an integer.' + self.mesh = mesh + self._nP = nP @property def nP(self): """ :rtype: int - :return: number of parameters in the model + :return: number of parameters that the mapping accepts """ + if self._nP is not None: + return self._nP if self.mesh is None: return '*' return self.mesh.nC @@ -32,11 +36,15 @@ class IdentityMap(object): @property def shape(self): """ - The default shape is (mesh.nC, nP). + The default shape is (mesh.nC, nP) if the mesh is defined. + If this is a meshless mapping (i.e. nP is defined independently) + the shape will be the the shape (nP,nP). :rtype: (int,int) :return: shape of the operator as a tuple """ + if self._nP is not None: + return (self.nP, self.nP) if self.mesh is None: return ('*', self.nP) return (self.mesh.nC, self.nP) @@ -118,6 +126,7 @@ class IdentityMap(object): def __str__(self): return "%s(%s,%s)" % (self.__class__.__name__, self.shape[0], self.shape[1]) + class ComboMap(IdentityMap): """Combination of various maps.""" @@ -475,7 +484,7 @@ class ActiveCells(IdentityMap): else: self.valInactive = valInactive.copy() self.valInactive[self.indActive] = 0 - + inds = np.nonzero(self.indActive)[0] self.P = sp.csr_matrix((np.ones(inds.size),(inds, range(inds.size))), shape=(self.nC, self.nP)) @@ -708,7 +717,7 @@ class PolyMap(IdentityMap): Parameterize the model space using a polynomials in a wholespace. ..math:: - + y = \mathbf{V} c Define the model as: @@ -752,10 +761,10 @@ class PolyMap(IdentityMap): else: raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] - Z = self.mesh.gridCC[:,2] + Y = self.mesh.gridCC[:,1] + Z = self.mesh.gridCC[:,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X elif self.normal =='Y': @@ -766,43 +775,43 @@ class PolyMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) else: raise(Exception("Only supports 2D")) - + return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5) - + def deriv(self, m): alpha = self.slope sig1,sig2, c = m[0],m[1],m[2:] if self.logSigma: sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D - if self.mesh.dim == 2: + if self.mesh.dim == 2: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] if self.normal =='X': f = polynomial.polyval(Y, c) - X - V = polynomial.polyvander(Y, len(c)-1) + V = polynomial.polyvander(Y, len(c)-1) elif self.normal =='Y': f = polynomial.polyval(X, c) - Y - V = polynomial.polyvander(X, len(c)-1) + V = polynomial.polyvander(X, len(c)-1) else: - raise(Exception("Input for normal = X or Y or Z")) + raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] if self.normal =='X': f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X - V = polynomial.polyvander2d(Y, Z, self.order) + V = polynomial.polyvander2d(Y, Z, self.order) elif self.normal =='Y': f = polynomial.polyval2d(X, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - Y - V = polynomial.polyvander2d(X, Z, self.order) + V = polynomial.polyvander2d(X, Z, self.order) elif self.normal =='Z': f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z - V = polynomial.polyvander2d(X, Y, self.order) + V = polynomial.polyvander2d(X, Y, self.order) else: raise(Exception("Input for normal = X or Y or Z")) @@ -815,16 +824,16 @@ class PolyMap(IdentityMap): g3 = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*V - return sp.csr_matrix(np.c_[g1,g2,g3]) + return sp.csr_matrix(np.c_[g1,g2,g3]) class SplineMap(IdentityMap): """SplineMap - Parameterize the boundary of two geological units using a spline interpolation + Parameterize the boundary of two geological units using a spline interpolation ..math:: - + g = f(x)-y Define the model as: @@ -849,7 +858,7 @@ class SplineMap(IdentityMap): def nP(self): if self.mesh.dim == 2: return np.size(self.pts)+2 - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: return np.size(self.pts)*2+2 else: raise(Exception("Only supports 2D and 3D")) @@ -866,28 +875,28 @@ class SplineMap(IdentityMap): X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] self.spl = UnivariateSpline(self.pts, c, k=self.order, s=0) - if self.normal =='X': + if self.normal =='X': f = self.spl(Y) - X elif self.normal =='Y': f = self.spl(X) - Y else: raise(Exception("Input for normal = X or Y or Z")) - # 3D: - # Comments: + # 3D: + # Comments: # Make two spline functions and link them using linear interpolation. # This is not quite direct extension of 2D to 3D case # Using 2D interpolation is possible - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] - Y = self.mesh.gridCC[:,1] + Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] - npts = np.size(self.pts) + npts = np.size(self.pts) if np.mod(c.size, 2): raise(Exception("Put even points!")) - + self.spl = {"splb":UnivariateSpline(self.pts, c[:npts], k=self.order, s=0), "splt":UnivariateSpline(self.pts, c[npts:], k=self.order, s=0)} @@ -902,7 +911,7 @@ class SplineMap(IdentityMap): raise(Exception("Input for normal = X or Y or Z")) else: raise(Exception("Only supports 2D and 3D")) - + return sig1+(sig2-sig1)*(np.arctan(alpha*f)/np.pi+0.5) @@ -912,7 +921,7 @@ class SplineMap(IdentityMap): if self.logSigma: sig1, sig2 = np.exp(sig1), np.exp(sig2) #2D - if self.mesh.dim == 2: + if self.mesh.dim == 2: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] @@ -921,9 +930,9 @@ class SplineMap(IdentityMap): elif self.normal =='Y': f = self.spl(X) - Y else: - raise(Exception("Input for normal = X or Y or Z")) + raise(Exception("Input for normal = X or Y or Z")) #3D - elif self.mesh.dim == 3: + elif self.mesh.dim == 3: X = self.mesh.gridCC[:,0] Y = self.mesh.gridCC[:,1] Z = self.mesh.gridCC[:,2] @@ -931,7 +940,7 @@ class SplineMap(IdentityMap): zb = self.ptsv[0] zt = self.ptsv[1] flines = (self.spl["splt"](Y)-self.spl["splb"](Y))*(Z-zb)/(zt-zb) + self.spl["splb"](Y) - f = flines - X + f = flines - X # elif self.normal =='Y': # elif self.normal =='Z': else: @@ -944,7 +953,7 @@ class SplineMap(IdentityMap): g1 = -(np.arctan(alpha*f)/np.pi + 0.5) + 1.0 g2 = (np.arctan(alpha*f)/np.pi + 0.5) - + if self.mesh.dim ==2: g3 = np.zeros((self.mesh.nC, self.npts)) if self.normal =='Y': @@ -958,7 +967,7 @@ class SplineMap(IdentityMap): cb = c.copy() dy = self.mesh.hy[ind]*1.5 ca[i] = ctemp+dy - cb[i] = ctemp-dy + cb[i] = ctemp-dy spla = UnivariateSpline(self.pts, ca, k=self.order, s=0) splb = UnivariateSpline(self.pts, cb, k=self.order, s=0) fderiv = (spla(X)-splb(X))/(2*dy) @@ -968,7 +977,7 @@ class SplineMap(IdentityMap): g3 = np.zeros((self.mesh.nC, self.npts*2)) if self.normal =='X': # Here we use perturbation to compute sensitivity - for i in range(self.npts*2): + for i in range(self.npts*2): ctemp = c[i] ind = np.argmin(abs(self.mesh.vectorCCy-ctemp)) ca = c.copy() @@ -982,20 +991,20 @@ class SplineMap(IdentityMap): splbb = UnivariateSpline(self.pts, cb[:self.npts], k=self.order, s=0) flinesa = (self.spl["splt"](Y)-splba(Y))*(Z-zb)/(zt-zb) + splba(Y) - X flinesb = (self.spl["splt"](Y)-splbb(Y))*(Z-zb)/(zt-zb) + splbb(Y) - X - #treat top boundary + #treat top boundary else: splta = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0) spltb = UnivariateSpline(self.pts, ca[self.npts:], k=self.order, s=0) flinesa = (self.spl["splt"](Y)-splta(Y))*(Z-zb)/(zt-zb) + splta(Y) - X - flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X - fderiv = (flinesa-flinesb)/(2*dy) + flinesb = (self.spl["splt"](Y)-spltb(Y))*(Z-zb)/(zt-zb) + spltb(Y) - X + fderiv = (flinesa-flinesb)/(2*dy) g3[:,i] = Utils.sdiag(alpha*(sig2-sig1)/(1.+(alpha*f)**2)/np.pi)*fderiv else : raise(Exception("Not Implemented for Y and Z, your turn :)")) - return sp.csr_matrix(np.c_[g1,g2,g3]) + return sp.csr_matrix(np.c_[g1,g2,g3]) + - \ No newline at end of file diff --git a/SimPEG/Mesh/MeshIO.py b/SimPEG/Mesh/MeshIO.py new file mode 100644 index 00000000..7501a66f --- /dev/null +++ b/SimPEG/Mesh/MeshIO.py @@ -0,0 +1,416 @@ +import numpy as np, os +from SimPEG import Utils + +class TensorMeshIO(object): + + @classmethod + def readUBC(TensorMesh, fileName): + """ + Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD + + Input: + :param fileName, path to the UBC GIF mesh file + + Output: + :param SimPEG TensorMesh object + """ + + # Interal function to read cell size lines for the UBC mesh files. + def readCellLine(line): + for seg in line.split(): + if '*' in seg: + st = seg + sp = seg.split('*') + re = np.array(sp[0],dtype=int)*(' ' + sp[1]) + line = line.replace(st,re.strip()) + return np.array(line.split(),dtype=float) + + # Read the file as line strings, remove lines with comment = ! + msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') + + # Fist line is the size of the model + sizeM = np.array(msh[0].split(),dtype=float) + # Second line is the South-West-Top corner coordinates. + x0 = np.array(msh[1].split(),dtype=float) + # Read the cell sizes + h1 = readCellLine(msh[2]) + h2 = readCellLine(msh[3]) + h3temp = readCellLine(msh[4]) + h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom. + # Adjust the reference point to the bottom south west corner + x0[2] = x0[2] - np.sum(h3) + # Make the mesh + tensMsh = TensorMesh([h1,h2,h3],x0) + return tensMsh + + @classmethod + def readVTK(TensorMesh, fileName): + """ + Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model + + Input: + :param vtrFileName, path to the vtr model file to write to + + Output: + :return SimPEG TensorMesh object + :return SimPEG model dictionary + + """ + # Import + from vtk import vtkXMLRectilinearGridReader as vtrFileReader + from vtk.util.numpy_support import vtk_to_numpy + + # Read the file + vtrReader = vtrFileReader() + vtrReader.SetFileName(fileName) + vtrReader.Update() + vtrGrid = vtrReader.GetOutput() + # Sort information + hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates()))) + xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0] + hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates()))) + yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0] + zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates())) + # Check the direction of hz + if np.all(zD < 0): + hz = np.abs(zD[::-1]) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] + else: + hz = np.abs(zD) + zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] + x0 = np.array([xR,yR,zR]) + + # Make the SimPEG object + tensMsh = TensorMesh([hx,hy,hz],x0) + + # Grap the models + models = {} + for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()): + modelName = vtrGrid.GetCellData().GetArrayName(i) + if np.all(zD < 0): + modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + tM = tensMsh.r(modFlip,'CC','CC','M') + modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V') + else: + modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) + models[modelName] = modArr + + # Return the data + return tensMsh, models + + def writeVTK(mesh, fileName, models=None): + """ + Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. + + Input: + :param str, path to the output vtk file + :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK + :param models, dictionary of numpy.array - Name('s) and array('s). Match number of cells + + """ + # Import + from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter, VTK_VERSION + from vtk.util.numpy_support import numpy_to_vtk + + # Deal with dimensionalities + if mesh.dim >= 1: + vX = mesh.vectorNx + xD = mesh.nNx + yD,zD = 1,1 + vY, vZ = np.array([0,0]) + if mesh.dim >= 2: + vY = mesh.vectorNy + yD = mesh.nNy + if mesh.dim == 3: + vZ = mesh.vectorNz + zD = mesh.nNz + # Use rectilinear VTK grid. + # Assign the spatial information. + vtkObj = rectGrid() + vtkObj.SetDimensions(xD,yD,zD) + vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1)) + vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1)) + vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) + + # Assign the model('s) to the object + if models is not None: + for item in models.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtkObj.GetCellData().AddArray(vtkDoubleArr) + # Set the active scalar + vtkObj.GetCellData().SetActiveScalars(models.keys()[0]) + # vtkObj.Update() + + # Check the extension of the fileName + ext = os.path.splitext(fileName)[1] + if ext is '': + fileName = fileName + '.vtr' + elif ext not in '.vtr': + raise IOError('{:s} is an incorrect extension, has to be .vtr') + # Write the file. + vtrWriteFilter = rectWriter() + if float(VTK_VERSION.split('.')[0]) >=6: + vtrWriteFilter.SetInputData(vtkObj) + else: + vtuWriteFilter.SetInput(vtuObj) + vtrWriteFilter.SetFileName(fileName) + vtrWriteFilter.Update() + + + def readModelUBC(mesh, fileName): + """ + Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg + + Input: + :param fileName, path to the UBC GIF mesh file to read + :param mesh, TensorMesh object, mesh that coresponds to the model + + Output: + :return numpy array, model with TensorMesh ordered + """ + f = open(fileName, 'r') + model = np.array(map(float, f.readlines())) + f.close() + model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') + model = model[::-1,:,:] + model = np.transpose(model, (1, 2, 0)) + model = Utils.mkvc(model) + return model + + def writeModelUBC(mesh, fileName, model): + """ + Writes a model associated with a SimPEG TensorMesh + to a UBC-GIF format model file. + + :param str fileName: File to write to + :param simpeg.Mesh.TensorMesh mesh: The mesh + :param numpy.ndarray model: The model + """ + + # Reshape model to a matrix + modelMat = mesh.r(model,'CC','CC','M') + # Transpose the axes + modelMatT = modelMat.transpose((2,0,1)) + # Flip z to positive down + modelMatTR = Utils.mkvc(modelMatT[::-1,:,:]) + + np.savetxt(fileName, modelMatTR.ravel()) + + def writeUBC(mesh, fileName, models=None): + """ + Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. + + :param str fileName: File to write to + :param simpeg.Mesh.TensorMesh mesh: The mesh + + """ + assert mesh.dim == 3 + s = '' + s += '%i %i %i\n' %tuple(mesh.vnC) + origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. + origin.dtype = float + + s += '%.2f %.2f %.2f\n' %tuple(origin) + s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx) + s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy) + s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1]) + f = open(fileName, 'w') + f.write(s) + f.close() + + if models is None: return + assert type(models) is dict, 'models must be a dict' + for key in models: + assert type(key) is str, 'The dict key is a file name' + mesh.writeModelUBC(key, models[key]) + +class TreeMeshIO(object): + + def writeUBC(mesh, fileName, models=None): + """ + Write UBC ocTree mesh and model files from a simpeg ocTree mesh and model. + + :param str fileName: File to write to + :param simpeg.Mesh.TreeMesh mesh: The mesh + :param dictionary models: The models in a dictionary, where the keys is the name of the of the model file + """ + + # Calculate information to write in the file. + # Number of cells in the underlying mesh + nCunderMesh = np.array([h.size for h in mesh.h],dtype=np.int64) + # The top-south-west most corner of the mesh + tswCorn = mesh.x0 + np.array([0,0,np.sum(mesh.h[2])]) + # Smallest cell size + smallCell = np.array([h.min() for h in mesh.h]) + # Number of cells + nrCells = mesh.nC + + ## Extract iformation about the cells. + # cell pointers + cellPointers = np.array([c._pointer for c in mesh]) + # cell with + cellW = np.array([ mesh._levelWidth(i) for i in cellPointers[:,-1] ]) + # Need to shift the pointers to work with UBC indexing + # UBC Octree indexes always the top-left-close (top-south-west) corner first and orders the cells in z(top-down),x,y vs x,y,z(bottom-up). + # Shift index up by 1 + ubcCellPt = cellPointers[:,0:-1].copy() + np.array([1.,1.,1.]) + # Need reindex the z index to be from the top-left-close corner and to be from the global top. + ubcCellPt[:,2] = ( nCunderMesh[-1] + 2) - (ubcCellPt[:,2] + cellW) + + # Reorder the ubcCellPt + ubcReorder = np.argsort(ubcCellPt.view(','.join(3*['float'])),axis=0,order=['f2','f1','f0'])[:,0] + # Make a array with the pointers and the withs, that are order in the ubc ordering + indArr = np.concatenate((ubcCellPt[ubcReorder,:],cellW[ubcReorder].reshape((-1,1)) ),axis=1) + + ## Write the UBC octree mesh file + with open(fileName,'w') as mshOut: + mshOut.write('{:.0f} {:.0f} {:.0f}\n'.format(nCunderMesh[0],nCunderMesh[1],nCunderMesh[2])) + mshOut.write('{:.4f} {:.4f} {:.4f}\n'.format(tswCorn[0],tswCorn[1],tswCorn[2])) + mshOut.write('{:.3f} {:.3f} {:.3f}\n'.format(smallCell[0],smallCell[1],smallCell[2])) + mshOut.write('{:.0f} \n'.format(nrCells)) + np.savetxt(mshOut,indArr,fmt='%i') + + ## Print the models + # Assign the model('s) to the object + if models is not None: + # indUBCvector = np.argsort(cX0[np.argsort(np.concatenate((cX0[:,0:2],cX0[:,2:3].max() - cX0[:,2:3]),axis=1).view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0]].view(','.join(3*['float'])),axis=0,order=('f2','f1','f0'))[:,0] + for item in models.iteritems(): + # Save the data + np.savetxt(item[0],item[1][ubcReorder],fmt='%3.5e') + + @classmethod + def readUBC(TreeMesh, meshFile): + """ + Read UBC 3D OcTree mesh and/or modelFiles + + Input: + :param str meshFile: path to the UBC GIF OcTree mesh file to read + + Output: + :return SimPEG.Mesh.TreeMesh mesh: The octree mesh + :return list of ndarray's: models as a list of numpy array's + """ + + ## Read the file lines + fileLines = np.genfromtxt(meshFile,dtype=str,delimiter='\n') + # Extract the data + nCunderMesh = np.array(fileLines[0].split(),dtype=float) + # I think this is the case? + if np.unique(nCunderMesh).size >1: + raise Exception('SimPEG TreeMeshes have the same number of cell in all directions') + tswCorn = np.array(fileLines[1].split(),dtype=float) + smallCell = np.array(fileLines[2].split(),dtype=float) + nrCells = np.array(fileLines[3].split(),dtype=float) + # Read the index array + indArr = np.genfromtxt(fileLines[4::],dtype=np.int) + + ## Calculate simpeg parameters + h1,h2,h3 = [np.ones(nr)*sz for nr,sz in zip(nCunderMesh,smallCell)] + x0 = tswCorn - np.array([0,0,np.sum(h3)]) + # Need to convert the index array to a points list that complies with SimPEG TreeMesh. + # Shift to start at 0 + simpegCellPt = indArr[:,0:-1].copy() + simpegCellPt[:,2] = ( nCunderMesh[-1] + 2) - (simpegCellPt[:,2] + indArr[:,3]) + # Need reindex the z index to be from the bottom-left-close corner and to be from the global bottom. + simpegCellPt = simpegCellPt - np.array([1.,1.,1.]) + + # Calculate the cell level + simpegLevel = np.log2(np.min(nCunderMesh)) - np.log2(indArr[:,3]) + # Make a pointer matrix + simpegPointers = np.concatenate((simpegCellPt,simpegLevel.reshape((-1,1))),axis=1) + + ## Make the tree mesh + mesh = TreeMesh([h1,h2,h3],x0) + mesh._cells = set([mesh._index(p) for p in simpegPointers.tolist()]) + + # Figure out the reordering + mesh._simpegReorderUBC = np.argsort(np.array([mesh._index(i) for i in simpegPointers.tolist()])) + # mesh._simpegReorderUBC = np.argsort((np.array([[1,1,1,-1]])*simpegPointers).view(','.join(4*['float'])),axis=0,order=['f3','f2','f1','f0'])[:,0] + + return mesh + + + def readModelUBC(mesh, fileName): + """ + Read UBC OcTree model and get vector + + Input: + :param fileName, path to the UBC GIF model file to read + + Output: + :return numpy array, OcTree model + """ + + if type(fileName) is list: + out = {} + for f in fileName: + out[f] = mesh.readModelUBC(f) + return out + + assert hasattr(mesh, '_simpegReorderUBC'), 'The file must have been loaded from a UBC format.' + assert mesh.dim == 3 + + modList = [] + modArr = np.loadtxt(fileName) + if len(modArr.shape) == 1: + modList.append(modArr[mesh._simpegReorderUBC]) + else: + modList.append(modArr[mesh._simpegReorderUBC,:]) + return modList + + def writeVTK(mesh, fileName, models=None): + """ + Function to write a VTU file from a SimPEG TreeMesh and model. + """ + import vtk + from vtk import vtkXMLUnstructuredGridWriter as Writer, VTK_VERSION + from vtk.util.numpy_support import numpy_to_vtk, numpy_to_vtkIdTypeArray + + if str(type(mesh)).split()[-1][1:-2] not in 'SimPEG.Mesh.TreeMesh.TreeMesh': + raise IOError('mesh is not a SimPEG TreeMesh.') + + # Make the data parts for the vtu object + # Points + mesh.number() + ptsMat = mesh._gridN + mesh.x0 + + vtkPts = vtk.vtkPoints() + vtkPts.SetData(numpy_to_vtk(ptsMat,deep=True)) + # Cells + cellConn = np.array([c.nodes for c in mesh],dtype=np.int64) + + cellsMat = np.concatenate((np.ones((cellConn.shape[0],1),dtype=np.int64)*cellConn.shape[1],cellConn),axis=1).ravel() + cellsArr = vtk.vtkCellArray() + cellsArr.SetNumberOfCells(cellConn.shape[0]) + cellsArr.SetCells(cellConn.shape[0],numpy_to_vtkIdTypeArray(cellsMat,deep=True)) + + # Make the object + vtuObj = vtk.vtkUnstructuredGrid() + vtuObj.SetPoints(vtkPts) + vtuObj.SetCells(vtk.VTK_VOXEL,cellsArr) + # Add the level of refinement as a cell array + cellSides = np.array([np.array(vtuObj.GetCell(i).GetBounds()).reshape((3,2)).dot(np.array([-1, 1])) for i in np.arange(vtuObj.GetNumberOfCells())]) + uniqueLevel, indLevel = np.unique(np.prod(cellSides,axis=1),return_inverse=True) + refineLevelArr = numpy_to_vtk(indLevel.max() - indLevel,deep=1) + refineLevelArr.SetName('octreeLevel') + vtuObj.GetCellData().AddArray(refineLevelArr) + # Assign the model('s) to the object + if models is not None: + for item in models.iteritems(): + # Convert numpy array + vtkDoubleArr = numpy_to_vtk(item[1],deep=1) + vtkDoubleArr.SetName(item[0]) + vtuObj.GetCellData().AddArray(vtkDoubleArr) + + # Make the writer + vtuWriteFilter = Writer() + if float(VTK_VERSION.split('.')[0]) >=6: + vtuWriteFilter.SetInputData(vtuObj) + else: + vtuWriteFilter.SetInput(vtuObj) + vtuWriteFilter.SetFileName(fileName) + # Write the file + vtuWriteFilter.Update() + diff --git a/SimPEG/Mesh/TensorMesh.py b/SimPEG/Mesh/TensorMesh.py index c76306fe..508f015c 100644 --- a/SimPEG/Mesh/TensorMesh.py +++ b/SimPEG/Mesh/TensorMesh.py @@ -1,558 +1,559 @@ -from SimPEG import Utils, np, sp -from BaseMesh import BaseMesh, BaseRectangularMesh -from View import TensorView -from DiffOperators import DiffOperators -from InnerProducts import InnerProducts - -class BaseTensorMesh(BaseMesh): - - __metaclass__ = Utils.SimPEGMetaClass - - _meshType = 'BASETENSOR' - - _unitDimensions = [1, 1, 1] - - def __init__(self, h_in, x0_in=None): - assert type(h_in) in [list, tuple], 'h_in must be a list' - assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3' - h = range(len(h_in)) - for i, h_i in enumerate(h_in): - if Utils.isScalar(h_i) and type(h_i) is not np.ndarray: - # This gives you something over the unit cube. - h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) - elif type(h_i) is list: - h_i = Utils.meshTensor(h_i) - assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i) - assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) - h[i] = h_i[:] # make a copy. - - x0 = np.zeros(len(h)) - if x0_in is not None: - assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" - for i in range(len(h)): - x_i, h_i = x0_in[i], h[i] - if Utils.isScalar(x_i): - x0[i] = x_i - elif x_i == '0': - x0[i] = 0.0 - elif x_i == 'C': - x0[i] = -h_i.sum()*0.5 - elif x_i == 'N': - x0[i] = -h_i.sum() - else: - raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i) - - if isinstance(self, BaseRectangularMesh): - BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) - else: - BaseMesh.__init__(self, np.array([x.size for x in h]), x0) - - # Ensure h contains 1D vectors - self._h = [Utils.mkvc(x.astype(float)) for x in h] - - @property - def h(self): - """h is a list containing the cell widths of the tensor mesh in each dimension.""" - return self._h - - @property - def hx(self): - "Width of cells in the x direction" - return self._h[0] - - @property - def hy(self): - "Width of cells in the y direction" - return None if self.dim < 2 else self._h[1] - - @property - def hz(self): - "Width of cells in the z direction" - return None if self.dim < 3 else self._h[2] - - @property - def vectorNx(self): - """Nodal grid vector (1D) in the x direction.""" - return np.r_[0., self.hx.cumsum()] + self.x0[0] - - @property - def vectorNy(self): - """Nodal grid vector (1D) in the y direction.""" - return None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] - - @property - def vectorNz(self): - """Nodal grid vector (1D) in the z direction.""" - return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] - - @property - def vectorCCx(self): - """Cell-centered grid vector (1D) in the x direction.""" - return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0] - - @property - def vectorCCy(self): - """Cell-centered grid vector (1D) in the y direction.""" - return None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] - - @property - def vectorCCz(self): - """Cell-centered grid vector (1D) in the z direction.""" - return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] - - @property - def gridCC(self): - """Cell-centered grid.""" - return self._getTensorGrid('CC') - - @property - def gridN(self): - """Nodal grid.""" - return self._getTensorGrid('N') - - @property - def gridFx(self): - """Face staggered grid in the x direction.""" - if self.nFx == 0: return - return self._getTensorGrid('Fx') - - @property - def gridFy(self): - """Face staggered grid in the y direction.""" - if self.nFy == 0 or self.dim < 2: return - return self._getTensorGrid('Fy') - - @property - def gridFz(self): - """Face staggered grid in the z direction.""" - if self.nFz == 0 or self.dim < 3: return - return self._getTensorGrid('Fz') - - @property - def gridEx(self): - """Edge staggered grid in the x direction.""" - if self.nEx == 0: return - return self._getTensorGrid('Ex') - - @property - def gridEy(self): - """Edge staggered grid in the y direction.""" - if self.nEy == 0 or self.dim < 2: return - return self._getTensorGrid('Ey') - - @property - def gridEz(self): - """Edge staggered grid in the z direction.""" - if self.nEz == 0 or self.dim < 3: return - return self._getTensorGrid('Ez') - - def _getTensorGrid(self, key): - if getattr(self, '_grid' + key, None) is None: - setattr(self, '_grid' + key, Utils.ndgrid(self.getTensor(key))) - return getattr(self, '_grid' + key) - - def getTensor(self, key): - """ Returns a tensor list. - - :param str key: What tensor (see below) - :rtype: list - :return: list of the tensors that make up the mesh. - - key can be:: - - 'CC' -> scalar field defined on cell centers - 'N' -> scalar field defined on nodes - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - - """ - - if key == 'Fx': - ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] - elif key == 'Fy': - ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] - elif key == 'Fz': - ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] - elif key == 'Ex': - ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] - elif key == 'Ey': - ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] - elif key == 'Ez': - ten = [self.vectorNx , self.vectorNy , self.vectorCCz] - elif key == 'CC': - ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] - elif key == 'N': - ten = [self.vectorNx , self.vectorNy , self.vectorNz ] - - return [t for t in ten if t is not None] - - # --------------- Methods --------------------- - - def isInside(self, pts, locType='N'): - """ - Determines if a set of points are inside a mesh. - - :param numpy.ndarray pts: Location of points to test - :rtype numpy.ndarray - :return inside, numpy array of booleans - """ - pts = Utils.asArray_N_x_Dim(pts, self.dim) - - tensors = self.getTensor(locType) - - if locType == 'N' and self._meshType == 'CYL': - #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! - tensors[0] = np.r_[0.,tensors[0]] - tensors[1] = np.r_[tensors[1], 2.0*np.pi] - - inside = np.ones(pts.shape[0],dtype=bool) - for i, tensor in enumerate(tensors): - TOL = np.diff(tensor).min() * 1.0e-10 - inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) - return inside - - def getInterpolationMat(self, loc, locType, zerosOutside=False): - """ Produces interpolation matrix - - :param numpy.ndarray loc: Location of points to interpolate to - :param str locType: What to interpolate (see below) - :rtype: scipy.sparse.csr.csr_matrix - :return: M, the interpolation matrix - - locType can be:: - - 'Ex' -> x-component of field defined on edges - 'Ey' -> y-component of field defined on edges - 'Ez' -> z-component of field defined on edges - 'Fx' -> x-component of field defined on faces - 'Fy' -> y-component of field defined on faces - 'Fz' -> z-component of field defined on faces - 'N' -> scalar field defined on nodes - 'CC' -> scalar 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) - - loc = Utils.asArray_N_x_Dim(loc, self.dim) - - if zerosOutside is False: - assert np.all(self.isInside(loc)), "Points outside of mesh" - else: - indZeros = np.logical_not(self.isInside(loc)) - loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) - - if locType in ['Fx','Fy','Fz','Ex','Ey','Ez']: - ind = {'x':0, 'y':1, 'z':2}[locType[1]] - assert self.dim >= ind, 'mesh is not high enough dimension.' - nF_nE = self.vnF if 'F' in locType else self.vnE - components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] - components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) - # remove any zero blocks (hstack complains) - components = [comp for comp in components if comp.shape[1] > 0] - Q = sp.hstack(components) - elif locType in ['CC', 'N']: - Q = Utils.interpmat(loc, *self.getTensor(locType)) - else: - raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) - - if zerosOutside: - Q[indZeros, :] = 0 - - return Q.tocsr() - - - def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): - """ - Fast version of getFaceInnerProduct. - This does not handle the case of a full tensor prop. - - :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) - :param str projType: 'E' or 'F' - :param bool returnP: returns the projection matrices - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: scipy.csr_matrix - :return: M, the inner product matrix (nF, nF) - """ - assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - - if prop is None: - prop = np.ones(self.nC) - - if invProp: - prop = 1./prop - - if Utils.isScalar(prop): - prop = prop*np.ones(self.nC) - - if prop.size == self.nC: - Av = getattr(self, 'ave'+projType+'2CC') - Vprop = self.vol * Utils.mkvc(prop) - M = self.dim * Utils.sdiag(Av.T * Vprop) - elif prop.size == self.nC*self.dim: - Av = getattr(self, 'ave'+projType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) - else: - return None - - if invMat: - return Utils.sdInv(M) - else: - return M - - def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False): - """ - :param str projType: 'E' or 'F' - :param TensorType tensorType: type of the tensor - :param bool invProp: inverts the material property - :param bool invMat: inverts the matrix - :rtype: function - :return: dMdmu, the derivative of the inner product matrix - """ - assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" - tensorType = Utils.TensorType(self, prop) - - dMdprop = None - - if invMat: - MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat) - - if tensorType == 0: - Av = getattr(self, 'ave'+projType+'2CC') - V = Utils.sdiag(self.vol) - ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) - if not invMat and not invProp: - dMdprop = self.dim * Av.T * V * ones - elif invMat and invProp: - dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2) - - if tensorType == 1: - Av = getattr(self, 'ave'+projType+'2CC') - V = Utils.sdiag(self.vol) - if not invMat and not invProp: - dMdprop = self.dim * Av.T * V - elif invMat and invProp: - dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) - - if tensorType == 2: # anisotropic - Av = getattr(self, 'ave'+projType+'2CCV') - V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) - if not invMat and not invProp: - dMdprop = Av.T * V - elif invMat and invProp: - dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) - - if dMdprop is not None: - def innerProductDeriv(v=None): - if v is None: - print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop' - return dMdprop - return Utils.sdiag(v) * dMdprop - return innerProductDeriv - else: - return None - - - -class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, InnerProducts): - """ - TensorMesh is a mesh class that deals with tensor product meshes. - - Any Mesh that has a constant width along the entire axis - such that it can defined by a single width vector, called 'h'. - - :: - - hx = np.array([1,1,1]) - hy = np.array([1,2]) - hz = np.array([1,1,1,1]) - - mesh = Mesh.TensorMesh([hx, hy, hz]) - - Example of a padded tensor mesh using :func:`SimPEG.Utils.meshutils.meshTensor`: - - .. plot:: - :include-source: - - from SimPEG import Mesh, Utils - M = Mesh.TensorMesh([[(10,10,-1.3),(10,40),(10,10,1.3)], [(10,10,-1.3),(10,20)]]) - M.plotGrid() - - For a quick tensor mesh on a (10x12x15) unit cube:: - - mesh = Mesh.TensorMesh([10, 12, 15]) - - """ - - __metaclass__ = Utils.SimPEGMetaClass - - _meshType = 'TENSOR' - - def __init__(self, h_in, x0=None): - BaseTensorMesh.__init__(self, h_in, x0) - - def __str__(self): - outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) - def printH(hx, outStr=''): - i = -1 - while True: - i = i + 1 - if i > hx.size: - break - elif i == hx.size: - break - h = hx[i] - n = 1 - for j in range(i+1, hx.size): - if hx[j] == h: - n = n + 1 - i = i + 1 - else: - break - - if n == 1: - outStr += ' {0:.2f},'.format(h) - else: - outStr += ' {0:d}*{1:.2f},'.format(n,h) - - return outStr[:-1] - - if self.dim == 1: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += printH(self.hx, outStr='\n hx:') - pass - elif self.dim == 2: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n y0: {0:.2f}'.format(self.x0[1]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += '\n nCy: {0:d}'.format(self.nCy) - outStr += printH(self.hx, outStr='\n hx:') - outStr += printH(self.hy, outStr='\n hy:') - elif self.dim == 3: - outStr += '\n x0: {0:.2f}'.format(self.x0[0]) - outStr += '\n y0: {0:.2f}'.format(self.x0[1]) - outStr += '\n z0: {0:.2f}'.format(self.x0[2]) - outStr += '\n nCx: {0:d}'.format(self.nCx) - outStr += '\n nCy: {0:d}'.format(self.nCy) - outStr += '\n nCz: {0:d}'.format(self.nCz) - outStr += printH(self.hx, outStr='\n hx:') - outStr += printH(self.hy, outStr='\n hy:') - outStr += printH(self.hz, outStr='\n hz:') - - return outStr - - - # --------------- Geometries --------------------- - @property - def vol(self): - """Construct cell volumes of the 3D model as 1d array.""" - if getattr(self, '_vol', None) is None: - vh = self.h - # Compute cell volumes - if self.dim == 1: - self._vol = Utils.mkvc(vh[0]) - elif self.dim == 2: - # Cell sizes in each direction - self._vol = Utils.mkvc(np.outer(vh[0], vh[1])) - elif self.dim == 3: - # Cell sizes in each direction - self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) - return self._vol - - @property - def area(self): - """Construct face areas of the 3D model as 1d array.""" - if getattr(self, '_area', None) is None: - # Ensure that we are working with column vectors - vh = self.h - # The number of cell centers in each direction - n = self.vnC - # Compute areas of cell faces - if(self.dim == 1): - self._area = np.ones(n[0]+1) - elif(self.dim == 2): - area1 = np.outer(np.ones(n[0]+1), vh[1]) - area2 = np.outer(vh[0], np.ones(n[1]+1)) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] - elif(self.dim == 3): - area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2]))) - area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] - return self._area - - @property - def edge(self): - """Construct edge legnths of the 3D model as 1d array.""" - if getattr(self, '_edge', None) is None: - # Ensure that we are working with column vectors - vh = self.h - # The number of cell centers in each direction - n = self.vnC - # Compute edge lengths - if(self.dim == 1): - self._edge = Utils.mkvc(vh[0]) - elif(self.dim == 2): - l1 = np.outer(vh[0], np.ones(n[1]+1)) - l2 = np.outer(np.ones(n[0]+1), vh[1]) - self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)] - elif(self.dim == 3): - l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) - l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) - l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) - self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] - return self._edge - - @property - def faceBoundaryInd(self): - """ - Find indices of boundary faces in each direction - """ - if self.dim==1: - indxd = (self.gridFx==min(self.gridFx)) - indxu = (self.gridFx==max(self.gridFx)) - return indxd, indxu - elif self.dim==2: - indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) - indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) - indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) - indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) - return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) - indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) - indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) - indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) - indzd = (self.gridFz[:,2]==min(self.gridFz[:,2])) - indzu = (self.gridFz[:,2]==max(self.gridFz[:,2])) - return indxd, indxu, indyd, indyu, indzd, indzu - - @property - def cellBoundaryInd(self): - """ - Find indices of boundary faces in each direction - """ - if self.dim==1: - indxd = (self.gridCC==min(self.gridCC)) - indxu = (self.gridCC==max(self.gridCC)) - return indxd, indxu - elif self.dim==2: - indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) - indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) - indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) - indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) - return indxd, indxu, indyd, indyu - elif self.dim==3: - indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) - indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) - indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) - indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) - indzd = (self.gridCC[:,2]==min(self.gridCC[:,2])) - indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) - return indxd, indxu, indyd, indyu, indzd, indzu +from SimPEG import Utils, np, sp +from BaseMesh import BaseMesh, BaseRectangularMesh +from View import TensorView +from DiffOperators import DiffOperators +from InnerProducts import InnerProducts +from MeshIO import TensorMeshIO + +class BaseTensorMesh(BaseMesh): + + __metaclass__ = Utils.SimPEGMetaClass + + _meshType = 'BASETENSOR' + + _unitDimensions = [1, 1, 1] + + def __init__(self, h_in, x0_in=None): + assert type(h_in) in [list, tuple], 'h_in must be a list' + assert len(h_in) in [1,2,3], 'h_in must be of dimension 1, 2, or 3' + h = range(len(h_in)) + for i, h_i in enumerate(h_in): + if Utils.isScalar(h_i) and type(h_i) is not np.ndarray: + # This gives you something over the unit cube. + h_i = self._unitDimensions[i] * np.ones(int(h_i))/int(h_i) + elif type(h_i) is list: + h_i = Utils.meshTensor(h_i) + assert isinstance(h_i, np.ndarray), ("h[%i] is not a numpy array." % i) + assert len(h_i.shape) == 1, ("h[%i] must be a 1D numpy array." % i) + h[i] = h_i[:] # make a copy. + + x0 = np.zeros(len(h)) + if x0_in is not None: + assert len(h) == len(x0_in), "Dimension mismatch. x0 != len(h)" + for i in range(len(h)): + x_i, h_i = x0_in[i], h[i] + if Utils.isScalar(x_i): + x0[i] = x_i + elif x_i == '0': + x0[i] = 0.0 + elif x_i == 'C': + x0[i] = -h_i.sum()*0.5 + elif x_i == 'N': + x0[i] = -h_i.sum() + else: + raise Exception("x0[%i] must be a scalar or '0' to be zero, 'C' to center, or 'N' to be negative." % i) + + if isinstance(self, BaseRectangularMesh): + BaseRectangularMesh.__init__(self, np.array([x.size for x in h]), x0) + else: + BaseMesh.__init__(self, np.array([x.size for x in h]), x0) + + # Ensure h contains 1D vectors + self._h = [Utils.mkvc(x.astype(float)) for x in h] + + @property + def h(self): + """h is a list containing the cell widths of the tensor mesh in each dimension.""" + return self._h + + @property + def hx(self): + "Width of cells in the x direction" + return self._h[0] + + @property + def hy(self): + "Width of cells in the y direction" + return None if self.dim < 2 else self._h[1] + + @property + def hz(self): + "Width of cells in the z direction" + return None if self.dim < 3 else self._h[2] + + @property + def vectorNx(self): + """Nodal grid vector (1D) in the x direction.""" + return np.r_[0., self.hx.cumsum()] + self.x0[0] + + @property + def vectorNy(self): + """Nodal grid vector (1D) in the y direction.""" + return None if self.dim < 2 else np.r_[0., self.hy.cumsum()] + self.x0[1] + + @property + def vectorNz(self): + """Nodal grid vector (1D) in the z direction.""" + return None if self.dim < 3 else np.r_[0., self.hz.cumsum()] + self.x0[2] + + @property + def vectorCCx(self): + """Cell-centered grid vector (1D) in the x direction.""" + return np.r_[0, self.hx[:-1].cumsum()] + self.hx*0.5 + self.x0[0] + + @property + def vectorCCy(self): + """Cell-centered grid vector (1D) in the y direction.""" + return None if self.dim < 2 else np.r_[0, self.hy[:-1].cumsum()] + self.hy*0.5 + self.x0[1] + + @property + def vectorCCz(self): + """Cell-centered grid vector (1D) in the z direction.""" + return None if self.dim < 3 else np.r_[0, self.hz[:-1].cumsum()] + self.hz*0.5 + self.x0[2] + + @property + def gridCC(self): + """Cell-centered grid.""" + return self._getTensorGrid('CC') + + @property + def gridN(self): + """Nodal grid.""" + return self._getTensorGrid('N') + + @property + def gridFx(self): + """Face staggered grid in the x direction.""" + if self.nFx == 0: return + return self._getTensorGrid('Fx') + + @property + def gridFy(self): + """Face staggered grid in the y direction.""" + if self.nFy == 0 or self.dim < 2: return + return self._getTensorGrid('Fy') + + @property + def gridFz(self): + """Face staggered grid in the z direction.""" + if self.nFz == 0 or self.dim < 3: return + return self._getTensorGrid('Fz') + + @property + def gridEx(self): + """Edge staggered grid in the x direction.""" + if self.nEx == 0: return + return self._getTensorGrid('Ex') + + @property + def gridEy(self): + """Edge staggered grid in the y direction.""" + if self.nEy == 0 or self.dim < 2: return + return self._getTensorGrid('Ey') + + @property + def gridEz(self): + """Edge staggered grid in the z direction.""" + if self.nEz == 0 or self.dim < 3: return + return self._getTensorGrid('Ez') + + def _getTensorGrid(self, key): + if getattr(self, '_grid' + key, None) is None: + setattr(self, '_grid' + key, Utils.ndgrid(self.getTensor(key))) + return getattr(self, '_grid' + key) + + def getTensor(self, key): + """ Returns a tensor list. + + :param str key: What tensor (see below) + :rtype: list + :return: list of the tensors that make up the mesh. + + key can be:: + + 'CC' -> scalar field defined on cell centers + 'N' -> scalar field defined on nodes + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + + """ + + if key == 'Fx': + ten = [self.vectorNx , self.vectorCCy, self.vectorCCz] + elif key == 'Fy': + ten = [self.vectorCCx, self.vectorNy , self.vectorCCz] + elif key == 'Fz': + ten = [self.vectorCCx, self.vectorCCy, self.vectorNz ] + elif key == 'Ex': + ten = [self.vectorCCx, self.vectorNy , self.vectorNz ] + elif key == 'Ey': + ten = [self.vectorNx , self.vectorCCy, self.vectorNz ] + elif key == 'Ez': + ten = [self.vectorNx , self.vectorNy , self.vectorCCz] + elif key == 'CC': + ten = [self.vectorCCx, self.vectorCCy, self.vectorCCz] + elif key == 'N': + ten = [self.vectorNx , self.vectorNy , self.vectorNz ] + + return [t for t in ten if t is not None] + + # --------------- Methods --------------------- + + def isInside(self, pts, locType='N'): + """ + Determines if a set of points are inside a mesh. + + :param numpy.ndarray pts: Location of points to test + :rtype numpy.ndarray + :return inside, numpy array of booleans + """ + pts = Utils.asArray_N_x_Dim(pts, self.dim) + + tensors = self.getTensor(locType) + + if locType == 'N' and self._meshType == 'CYL': + #NOTE: for a CYL mesh we add a node to check if we are inside in the radial direction! + tensors[0] = np.r_[0.,tensors[0]] + tensors[1] = np.r_[tensors[1], 2.0*np.pi] + + inside = np.ones(pts.shape[0],dtype=bool) + for i, tensor in enumerate(tensors): + TOL = np.diff(tensor).min() * 1.0e-10 + inside = inside & (pts[:,i] >= tensor.min()-TOL) & (pts[:,i] <= tensor.max()+TOL) + return inside + + def getInterpolationMat(self, loc, locType='CC', zerosOutside=False): + """ Produces interpolation matrix + + :param numpy.ndarray loc: Location of points to interpolate to + :param str locType: What to interpolate (see below) + :rtype: scipy.sparse.csr.csr_matrix + :return: M, the interpolation matrix + + locType can be:: + + 'Ex' -> x-component of field defined on edges + 'Ey' -> y-component of field defined on edges + 'Ez' -> z-component of field defined on edges + 'Fx' -> x-component of field defined on faces + 'Fy' -> y-component of field defined on faces + 'Fz' -> z-component of field defined on faces + 'N' -> scalar field defined on nodes + 'CC' -> scalar 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) + + loc = Utils.asArray_N_x_Dim(loc, self.dim) + + if zerosOutside is False: + assert np.all(self.isInside(loc)), "Points outside of mesh" + else: + indZeros = np.logical_not(self.isInside(loc)) + loc[indZeros, :] = np.array([v.mean() for v in self.getTensor('CC')]) + + if locType in ['Fx','Fy','Fz','Ex','Ey','Ez']: + ind = {'x':0, 'y':1, 'z':2}[locType[1]] + assert self.dim >= ind, 'mesh is not high enough dimension.' + nF_nE = self.vnF if 'F' in locType else self.vnE + components = [Utils.spzeros(loc.shape[0], n) for n in nF_nE] + components[ind] = Utils.interpmat(loc, *self.getTensor(locType)) + # remove any zero blocks (hstack complains) + components = [comp for comp in components if comp.shape[1] > 0] + Q = sp.hstack(components) + elif locType in ['CC', 'N']: + Q = Utils.interpmat(loc, *self.getTensor(locType)) + else: + raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim)) + + if zerosOutside: + Q[indZeros, :] = 0 + + return Q.tocsr() + + + def _fastInnerProduct(self, projType, prop=None, invProp=False, invMat=False): + """ + Fast version of getFaceInnerProduct. + This does not handle the case of a full tensor prop. + + :param numpy.array prop: material property (tensor properties are possible) at each cell center (nC, (1, 3, or 6)) + :param str projType: 'E' or 'F' + :param bool returnP: returns the projection matrices + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix + :rtype: scipy.csr_matrix + :return: M, the inner product matrix (nF, nF) + """ + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + + if prop is None: + prop = np.ones(self.nC) + + if invProp: + prop = 1./prop + + if Utils.isScalar(prop): + prop = prop*np.ones(self.nC) + + if prop.size == self.nC: + Av = getattr(self, 'ave'+projType+'2CC') + Vprop = self.vol * Utils.mkvc(prop) + M = self.dim * Utils.sdiag(Av.T * Vprop) + elif prop.size == self.nC*self.dim: + Av = getattr(self, 'ave'+projType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + M = Utils.sdiag(Av.T * V * Utils.mkvc(prop)) + else: + return None + + if invMat: + return Utils.sdInv(M) + else: + return M + + def _fastInnerProductDeriv(self, projType, prop, invProp=False, invMat=False): + """ + :param str projType: 'E' or 'F' + :param TensorType tensorType: type of the tensor + :param bool invProp: inverts the material property + :param bool invMat: inverts the matrix + :rtype: function + :return: dMdmu, the derivative of the inner product matrix + """ + assert projType in ['F', 'E'], "projType must be 'F' for faces or 'E' for edges" + tensorType = Utils.TensorType(self, prop) + + dMdprop = None + + if invMat: + MI = self._fastInnerProduct(projType, prop, invProp=invProp, invMat=invMat) + + if tensorType == 0: + Av = getattr(self, 'ave'+projType+'2CC') + V = Utils.sdiag(self.vol) + ones = sp.csr_matrix((np.ones(self.nC), (range(self.nC), np.zeros(self.nC))), shape=(self.nC,1)) + if not invMat and not invProp: + dMdprop = self.dim * Av.T * V * ones + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * ones * Utils.sdiag(1./prop**2) + + if tensorType == 1: + Av = getattr(self, 'ave'+projType+'2CC') + V = Utils.sdiag(self.vol) + if not invMat and not invProp: + dMdprop = self.dim * Av.T * V + elif invMat and invProp: + dMdprop = self.dim * Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) + + if tensorType == 2: # anisotropic + Av = getattr(self, 'ave'+projType+'2CCV') + V = sp.kron(sp.identity(self.dim), Utils.sdiag(self.vol)) + if not invMat and not invProp: + dMdprop = Av.T * V + elif invMat and invProp: + dMdprop = Utils.sdiag(MI.diagonal()**2) * Av.T * V * Utils.sdiag(1./prop**2) + + if dMdprop is not None: + def innerProductDeriv(v=None): + if v is None: + print 'Depreciation Warning: TensorMesh.innerProductDeriv. You should be supplying a vector. Use: sdiag(u)*dMdprop' + return dMdprop + return Utils.sdiag(v) * dMdprop + return innerProductDeriv + else: + return None + + + +class TensorMesh(BaseTensorMesh, BaseRectangularMesh, TensorView, DiffOperators, InnerProducts, TensorMeshIO): + """ + TensorMesh is a mesh class that deals with tensor product meshes. + + Any Mesh that has a constant width along the entire axis + such that it can defined by a single width vector, called 'h'. + + :: + + hx = np.array([1,1,1]) + hy = np.array([1,2]) + hz = np.array([1,1,1,1]) + + mesh = Mesh.TensorMesh([hx, hy, hz]) + + Example of a padded tensor mesh using :func:`SimPEG.Utils.meshutils.meshTensor`: + + .. plot:: + :include-source: + + from SimPEG import Mesh, Utils + M = Mesh.TensorMesh([[(10,10,-1.3),(10,40),(10,10,1.3)], [(10,10,-1.3),(10,20)]]) + M.plotGrid() + + For a quick tensor mesh on a (10x12x15) unit cube:: + + mesh = Mesh.TensorMesh([10, 12, 15]) + + """ + + __metaclass__ = Utils.SimPEGMetaClass + + _meshType = 'TENSOR' + + def __init__(self, h_in, x0=None): + BaseTensorMesh.__init__(self, h_in, x0) + + def __str__(self): + outStr = ' ---- {0:d}-D TensorMesh ---- '.format(self.dim) + def printH(hx, outStr=''): + i = -1 + while True: + i = i + 1 + if i > hx.size: + break + elif i == hx.size: + break + h = hx[i] + n = 1 + for j in range(i+1, hx.size): + if hx[j] == h: + n = n + 1 + i = i + 1 + else: + break + + if n == 1: + outStr += ' {0:.2f},'.format(h) + else: + outStr += ' {0:d}*{1:.2f},'.format(n,h) + + return outStr[:-1] + + if self.dim == 1: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += printH(self.hx, outStr='\n hx:') + pass + elif self.dim == 2: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n y0: {0:.2f}'.format(self.x0[1]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += '\n nCy: {0:d}'.format(self.nCy) + outStr += printH(self.hx, outStr='\n hx:') + outStr += printH(self.hy, outStr='\n hy:') + elif self.dim == 3: + outStr += '\n x0: {0:.2f}'.format(self.x0[0]) + outStr += '\n y0: {0:.2f}'.format(self.x0[1]) + outStr += '\n z0: {0:.2f}'.format(self.x0[2]) + outStr += '\n nCx: {0:d}'.format(self.nCx) + outStr += '\n nCy: {0:d}'.format(self.nCy) + outStr += '\n nCz: {0:d}'.format(self.nCz) + outStr += printH(self.hx, outStr='\n hx:') + outStr += printH(self.hy, outStr='\n hy:') + outStr += printH(self.hz, outStr='\n hz:') + + return outStr + + + # --------------- Geometries --------------------- + @property + def vol(self): + """Construct cell volumes of the 3D model as 1d array.""" + if getattr(self, '_vol', None) is None: + vh = self.h + # Compute cell volumes + if self.dim == 1: + self._vol = Utils.mkvc(vh[0]) + elif self.dim == 2: + # Cell sizes in each direction + self._vol = Utils.mkvc(np.outer(vh[0], vh[1])) + elif self.dim == 3: + # Cell sizes in each direction + self._vol = Utils.mkvc(np.outer(Utils.mkvc(np.outer(vh[0], vh[1])), vh[2])) + return self._vol + + @property + def area(self): + """Construct face areas of the 3D model as 1d array.""" + if getattr(self, '_area', None) is None: + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.vnC + # Compute areas of cell faces + if(self.dim == 1): + self._area = np.ones(n[0]+1) + elif(self.dim == 2): + area1 = np.outer(np.ones(n[0]+1), vh[1]) + area2 = np.outer(vh[0], np.ones(n[1]+1)) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2)] + elif(self.dim == 3): + area1 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], vh[2]))) + area2 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + area3 = np.outer(vh[0], Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + self._area = np.r_[Utils.mkvc(area1), Utils.mkvc(area2), Utils.mkvc(area3)] + return self._area + + @property + def edge(self): + """Construct edge legnths of the 3D model as 1d array.""" + if getattr(self, '_edge', None) is None: + # Ensure that we are working with column vectors + vh = self.h + # The number of cell centers in each direction + n = self.vnC + # Compute edge lengths + if(self.dim == 1): + self._edge = Utils.mkvc(vh[0]) + elif(self.dim == 2): + l1 = np.outer(vh[0], np.ones(n[1]+1)) + l2 = np.outer(np.ones(n[0]+1), vh[1]) + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2)] + elif(self.dim == 3): + l1 = np.outer(vh[0], Utils.mkvc(np.outer(np.ones(n[1]+1), np.ones(n[2]+1)))) + l2 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(vh[1], np.ones(n[2]+1)))) + l3 = np.outer(np.ones(n[0]+1), Utils.mkvc(np.outer(np.ones(n[1]+1), vh[2]))) + self._edge = np.r_[Utils.mkvc(l1), Utils.mkvc(l2), Utils.mkvc(l3)] + return self._edge + + @property + def faceBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridFx==min(self.gridFx)) + indxu = (self.gridFx==max(self.gridFx)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridFx[:,0]==min(self.gridFx[:,0])) + indxu = (self.gridFx[:,0]==max(self.gridFx[:,0])) + indyd = (self.gridFy[:,1]==min(self.gridFy[:,1])) + indyu = (self.gridFy[:,1]==max(self.gridFy[:,1])) + indzd = (self.gridFz[:,2]==min(self.gridFz[:,2])) + indzu = (self.gridFz[:,2]==max(self.gridFz[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu + + @property + def cellBoundaryInd(self): + """ + Find indices of boundary faces in each direction + """ + if self.dim==1: + indxd = (self.gridCC==min(self.gridCC)) + indxu = (self.gridCC==max(self.gridCC)) + return indxd, indxu + elif self.dim==2: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + return indxd, indxu, indyd, indyu + elif self.dim==3: + indxd = (self.gridCC[:,0]==min(self.gridCC[:,0])) + indxu = (self.gridCC[:,0]==max(self.gridCC[:,0])) + indyd = (self.gridCC[:,1]==min(self.gridCC[:,1])) + indyu = (self.gridCC[:,1]==max(self.gridCC[:,1])) + indzd = (self.gridCC[:,2]==min(self.gridCC[:,2])) + indzu = (self.gridCC[:,2]==max(self.gridCC[:,2])) + return indxd, indxu, indyd, indyu, indzd, indzu diff --git a/SimPEG/Mesh/TreeMesh.py b/SimPEG/Mesh/TreeMesh.py index 1d5c8c8c..02c23dee 100644 --- a/SimPEG/Mesh/TreeMesh.py +++ b/SimPEG/Mesh/TreeMesh.py @@ -100,11 +100,12 @@ except Exception, e: from InnerProducts import InnerProducts from TensorMesh import TensorMesh, BaseTensorMesh +from MeshIO import TreeMeshIO import time MAX_BITS = 20 -class TreeMesh(BaseTensorMesh, InnerProducts): +class TreeMesh(BaseTensorMesh, InnerProducts, TreeMeshIO): _meshType = 'TREE' @@ -564,15 +565,18 @@ class TreeMesh(BaseTensorMesh, InnerProducts): return [p - (p % mod) for p in pointer[:-1]] + [pointer[-1]-1] def _cellN(self, p): + """Node location [x,y(,z)] of a single cell, closest to origin, given a pointer.""" p = self._asPointer(p) return [hi[:p[ii]].sum() for ii, hi in enumerate(self.h)] def _cellH(self, p): + """Widths of a single cell given a pointer.""" p = self._asPointer(p) w = self._levelWidth(p[-1]) return [hi[p[ii]:p[ii]+w].sum() for ii, hi in enumerate(self.h)] def _cellC(self, p): + """Cell center of a single cell (without origin correction), given a pointer.""" return (np.array(self._cellH(p))/2.0 + self._cellN(p)).tolist() def _levelWidth(self, level): @@ -827,8 +831,10 @@ class TreeMesh(BaseTensorMesh, InnerProducts): def _numberCells(self, force=False): if not self.__dirtyCells__ and not force: return self._cc2i = dict() + self._i2cc = dict() for ii, c in enumerate(sorted(self._cells)): self._cc2i[c] = ii + self._i2cc[ii] = c self.__dirtyCells__ = False def _numberNodes(self, force=False): @@ -1704,9 +1710,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts): "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CC', None) is None: if self.dim == 2: - self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]) + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC]).tocsr() elif self.dim == 3: - self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) + self._aveF2CC = 1./self.dim*sp.hstack([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr() return self._aveF2CC @property @@ -1714,9 +1720,9 @@ class TreeMesh(BaseTensorMesh, InnerProducts): "Construct the averaging operator on cell faces to cell centers." if getattr(self, '_aveF2CCV', None) is None: if self.dim == 2: - self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]) + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC]).tocsr() elif self.dim == 3: - self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]) + self._aveF2CCV = sp.block_diag([self.aveFx2CC, self.aveFy2CC, self.aveFz2CC]).tocsr() return self._aveF2CCV @property @@ -2218,6 +2224,25 @@ class TreeMesh(BaseTensorMesh, InnerProducts): if showIt: plt.show() return tuple(out) + def __len__(self): return self.nC + + def __getitem__(self, key): + if isinstance( key, slice ) : + #Get the start, stop, and step from the slice + return [self[ii] for ii in xrange(*key.indices(len(self)))] + elif isinstance( key, int ) : + if key < 0 : #Handle negative indices + key += len( self ) + if key >= len( self ) : + raise IndexError, "The index (%d) is out of range."%key + + self._numberCells() # no-op if numbered + index = self._i2cc[key] + pointer = self._asPointer(index) + return Cell(self, index, pointer) + else: + raise TypeError, "Invalid argument type." + class Cell(object): def __init__(self, mesh, index, pointer): @@ -2225,6 +2250,35 @@ class Cell(object): self._index = index self._pointer = pointer + @property + def nodes(self): + """The node index in _gridN (this may include hanging nodes).""" + M = self.mesh + M._numberNodes() + p = self._pointer + i = self._index + w = M._levelWidth(p[-1]) + + if M.dim == 2: + n = [ + i, + M._index([ p[0] + w, p[1] , p[2]]), + M._index([ p[0] , p[1]+ w, p[2]]), + M._index([ p[0] + w, p[1]+ w, p[2]]), + ] + elif self.dim == 3: + n = [ + i, + M._index([ p[0] + w, p[1] , p[2] ,p[3]]), + M._index([ p[0] , p[1] + w, p[2] ,p[3]]), + M._index([ p[0] + w, p[1] + w, p[2] ,p[3]]), + M._index([ p[0] , p[1] , p[2] + w,p[3]]), + M._index([ p[0] + w, p[1] , p[2] + w,p[3]]), + M._index([ p[0] , p[1] + w, p[2] + w,p[3]]), + M._index([ p[0] + w, p[1] + w, p[2] + w,p[3]]), + ] + return [M._n2i[_] for _ in n] + @property def center(self): if getattr(self, '_center', None) is None: @@ -2282,121 +2336,3 @@ class NotBalancedException(TreeException): pass class CellLookUpException(TreeException): pass - -if __name__ == '__main__': - - - import matplotlib.pyplot as plt - import matplotlib - from mpl_toolkits.mplot3d import Axes3D - import matplotlib.colors as colors - import matplotlib.cm as cmx - - def topo(x): - return np.sin(x*(2.*np.pi))*0.3 + 0.5 - - def function(cell): - r = cell.center - np.array([0.5]*len(cell.center)) - dist = np.sqrt(r.dot(r)) - # dist2 = np.abs(cell.center[-1] - topo(cell.center[0])) - - # dist = min([dist1,dist2]) - # if dist < 0.05: - # return 5 - if dist < 0.1: - return 5 - if dist < 0.2: - return 4 - if dist < 0.4: - return 3 - return 2 - - # T = TreeMesh([[(1,128)],[(1,128)],[(1,128)]],levels=7) - # T = TreeMesh([128,128,128]) - # T = TreeMesh([64,64],levels=6) - T = TreeMesh([4,4,4]) - # T = TreeMesh([[(1,128)],[(1,128)]],levels=7) - # T.refine(lambda xc:2, balance=False) - # T._index([0,0,0]) - # T._pointer(0) - - - # tic = time.time() - T.refine(function)#, balance=False) - # print time.time() - tic - # print T.nC - T.plotSlice(np.log(T.vol))#np.random.rand(T.nC)) - - plt.show() - blah - - # T.plotImage(np.arange(len(T.vol)),showIt=True) - - # print T.getFaceInnerProduct() - # print T.gridFz - - - # T._refineCell([8,0,1]) - # T._refineCell([8,0,2]) - # T._refineCell([12,0,2]) - # T._refineCell([8,4,2]) - # T._refineCell([6,0,3]) - # T._refineCell([8,8,1]) - # T._refineCell([0,0,0,1]) - # T.__dirty__ = True - - - # print T.gridFx.shape[0], T.nFx - - - - ax = plt.subplot(211) - ax.spy(T.edgeCurl) - - # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - # print T.edgeCurl.todense() - # print Mesh.TensorMesh([2,2,2]).edgeCurl.todense() - T.edgeCurl.todense() - # print T.gridEy - Mesh.TensorMesh([2,2,2]).gridEy - - # print T.edge - # T.plotGrid(ax=ax) - - # R = deflationMatrix(T._facesX, T._hangingFx, T._fx2i) - # print R - - ax = plt.subplot(212)#, projection='3d') - ax.spy(Mesh.TensorMesh([2,2,2]).edgeCurl) - - # ax = plt.subplot(313) - # ax.spy(T.faceDiv[:,:T.nFx] * R) - - - # T.balance() - # T.plotGrid(ax=ax) - - # cx = T._getNextCell([0,0,1],direction=0,positive=True) - # print cx - # # print [T._asPointer(_) for _ in cx] - # cx = T._getNextCell([8,0,3],direction=0,positive=False) - # print T._asPointer(cx) - # cx = T._getNextCell([8,8,1],direction=1,positive=False) - # print cx, #[T._asPointer(_) for _ in cx] - # cm = T._getNextCell([64,80,4],direction=0,positive=False) - # cy = T._getNextCell([64,80,4],direction=1,positive=True) - # cp = T._getNextCell([64,80,4],direction=1,positive=False) - - # ax.plot( T._cellN([4,0,1])[0],T._cellN([4,0,1])[1], 'yd') - # ax.plot( T._cellN(cx)[0],T._cellN(cx)[1], 'ys') - # ax.plot( T._cellN(cm)[0],T._cellN(cm)[1], 'ys') - # ax.plot( T._cellN(cy)[0],T._cellN(cy)[1], 'ys') - # ax.plot( T._cellN(cp[0])[0],T._cellN(cp[0])[1], 'ys') - # ax.plot( T._cellN(cp[1])[0],T._cellN(cp[1])[1], 'ys') - - - - - - # print T.nN - - plt.show() - diff --git a/SimPEG/Optimization.py b/SimPEG/Optimization.py index 4f2cb062..54f6c4ff 100644 --- a/SimPEG/Optimization.py +++ b/SimPEG/Optimization.py @@ -990,4 +990,18 @@ 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) ) + + delx = delx + rhs_a * dm_i / dm_a /10. + + # 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 be29ff01..f24b57de 100644 --- a/SimPEG/Problem.py +++ b/SimPEG/Problem.py @@ -32,8 +32,8 @@ class BaseProblem(object): val._assertMatchesPair(self.mapPair) self._mapping = val else: - self._mapping = self.PropMap(val) - + self._mapping = self.PropMap(val) + def __init__(self, mesh, mapping=None, **kwargs): Utils.setKwargs(self, **kwargs) assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." @@ -158,9 +158,6 @@ class BaseProblem(object): class BaseTimeProblem(BaseProblem): """Sets up that basic needs of a time domain problem.""" - - waveformType = "STEPOFF" - current = None @property def timeSteps(self): @@ -187,11 +184,6 @@ class BaseTimeProblem(BaseProblem): self._timeSteps = Utils.meshTensor(value) del self.timeMesh - def currentwaveform(self, wave): - self._timeSteps = np.diff(wave[:,0]) - self.current = wave[:,1] - self.waveformType = "GENERAL" - @property def nT(self): "Number of time steps." diff --git a/SimPEG/Regularization.py b/SimPEG/Regularization.py index 6c5b5f97..677fba5c 100644 --- a/SimPEG/Regularization.py +++ b/SimPEG/Regularization.py @@ -20,12 +20,13 @@ class BaseRegularization(object): mesh = None #: A SimPEG.Mesh instance. mref = None #: Reference model. - def __init__(self, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive=None, **kwargs): Utils.setKwargs(self, **kwargs) self.mesh = mesh assert isinstance(mesh, Mesh.BaseMesh), "mesh must be a SimPEG.Mesh object." self.mapping = mapping or self.mapPair(mesh) self.mapping._assertMatchesPair(self.mapPair) + self.indActive = indActive @property def parent(self): @@ -112,8 +113,6 @@ class BaseRegularization(object): return mD.T * ( self.W.T * ( self.W * ( mD * v) ) ) - - class Tikhonov(BaseRegularization): """ """ @@ -126,14 +125,18 @@ class Tikhonov(BaseRegularization): alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction") alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction") - def __init__(self, mesh, mapping=None, **kwargs): + def __init__(self, mesh, mapping=None, indActive = None, **kwargs): BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs) + self.indActive = indActive @property def Ws(self): """Regularization matrix Ws""" if getattr(self,'_Ws', None) is None: - self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5) + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Ws = Pac.T * self._Ws * Pac return self._Ws @property @@ -142,6 +145,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wx', None) is None: Ave_x_vol = self.mesh.aveF2CC[:,:self.mesh.nFx].T*self.mesh.vol self._Wx = Utils.sdiag((Ave_x_vol*self.alpha_x)**0.5)*self.mesh.cellGradx + + if self.indActive is not None: + indActive_Fx = (self.mesh.aveFx2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafx = Utils.speye(self.mesh.nFx)[:,indActive_Fx] + self._Wx = Pafx.T*self._Wx*Pac + return self._Wx @property @@ -150,6 +160,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wy', None) is None: Ave_y_vol = self.mesh.aveF2CC[:,self.mesh.nFx:np.sum(self.mesh.vnF[:2])].T*self.mesh.vol self._Wy = Utils.sdiag((Ave_y_vol*self.alpha_y)**0.5)*self.mesh.cellGrady + + if self.indActive is not None: + indActive_Fy = (self.mesh.aveFy2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafy = Utils.speye(self.mesh.nFy)[:,indActive_Fy] + self._Wy = Pafy.T*self._Wy*Pac + return self._Wy @property @@ -158,6 +175,13 @@ class Tikhonov(BaseRegularization): if getattr(self, '_Wz', None) is None: Ave_z_vol = self.mesh.aveF2CC[:,np.sum(self.mesh.vnF[:2]):].T*self.mesh.vol self._Wz = Utils.sdiag((Ave_z_vol*self.alpha_z)**0.5)*self.mesh.cellGradz + + if self.indActive is not None: + indActive_Fz = (self.mesh.aveFz2CC.T * self.indActive) == 1 + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + Pafz = Utils.speye(self.mesh.nFz)[:,indActive_Fz] + self._Wz = Pafz.T*self._Wz*Pac + return self._Wz @property @@ -165,6 +189,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wxx""" if getattr(self, '_Wxx', None) is None: self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wxx = Pac.T*self._Wxx*Pac + return self._Wxx @property @@ -172,6 +201,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wyy""" if getattr(self, '_Wyy', None) is None: self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wyy = Pac.T*self._Wyy*Pac + return self._Wyy @property @@ -179,6 +213,11 @@ class Tikhonov(BaseRegularization): """Regularization matrix Wzz""" if getattr(self, '_Wzz', None) is None: self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz + + if self.indActive is not None: + Pac = Utils.speye(self.mesh.nC)[:,self.indActive] + self._Wzz = Pac.T*self._Wzz*Pac + return self._Wzz @property diff --git a/SimPEG/Survey.py b/SimPEG/Survey.py index 88355df1..47a88ae2 100644 --- a/SimPEG/Survey.py +++ b/SimPEG/Survey.py @@ -205,6 +205,7 @@ class BaseSurvey(object): __metaclass__ = Utils.SimPEGMetaClass std = None #: Estimated Standard Deviations + eps = None #: Estimated Noise Floor dobs = None #: Observed data dtrue = None #: True data, if data is synthetic mtrue = None #: True model, if data is synthetic diff --git a/SimPEG/Utils/SolverUtils.py b/SimPEG/Utils/SolverUtils.py index 26ff3e2a..47c899d3 100644 --- a/SimPEG/Utils/SolverUtils.py +++ b/SimPEG/Utils/SolverUtils.py @@ -26,7 +26,14 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): def __init__(self, A, **kwargs): self.A = A.tocsc() + + self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy) + if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"] + self.accuracyTol = kwargs.get("accuracyTol", accuracyTol) + if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"] + self.kwargs = kwargs + if factorize: self.solver = fun(self.A, **kwargs) @@ -57,8 +64,8 @@ def SolverWrapD(fun, factorize=True, checkAccuracy=True, accuracyTol=1e-6): else: X[:,i] = fun(self.A, b[:,i], **self.kwargs) - if checkAccuracy: - _checkAccuracy(self.A, b, X, accuracyTol) + if self.checkAccuracy: + _checkAccuracy(self.A, b, X, self.accuracyTol) return X def clean(self): @@ -81,6 +88,12 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5): def __init__(self, A, **kwargs): self.A = A + + self.checkAccuracy = kwargs.get("checkAccuracy", checkAccuracy) + if kwargs.has_key("checkAccuracy"): del kwargs["checkAccuracy"] + self.accuracyTol = kwargs.get("accuracyTol", accuracyTol) + if kwargs.has_key("accuracyTol"): del kwargs["accuracyTol"] + self.kwargs = kwargs def __mul__(self, b): @@ -108,8 +121,8 @@ def SolverWrapI(fun, checkAccuracy=True, accuracyTol=1e-5): else: X[:,i] = out - if checkAccuracy: - _checkAccuracy(self.A, b, X, accuracyTol) + if self.checkAccuracy: + _checkAccuracy(self.A, b, X, self.accuracyTol) return X def clean(self): diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 250146c1..18c1994f 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -1,6 +1,6 @@ from matutils import * from codeutils import * -from meshutils import exampleLrmGrid, meshTensor, closestPoints, readUBCTensorMesh, writeUBCTensorMesh, writeUBCTensorModel, readVTRFile, writeVTRFile +from meshutils import * from curvutils import volTetra, faceInfo, indexCube from interputils import interpmat from CounterUtils import * diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 4a9a28a7..bfd00889 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -17,7 +17,7 @@ def memProfileWrapper(towrap, *funNames): For example:: - foo_mem = memProfile(foo,'my_func') + foo_mem = memProfileWrapper(foo,['my_func']) fooi = foo_mem() for i in range(5): fooi.my_func() diff --git a/SimPEG/Utils/matutils.py b/SimPEG/Utils/matutils.py index b38bb4a1..3a6030fa 100644 --- a/SimPEG/Utils/matutils.py +++ b/SimPEG/Utils/matutils.py @@ -2,7 +2,6 @@ import numpy as np import scipy.sparse as sp from codeutils import isScalar - def mkvc(x, numDims=1): """Creates a vector with the number of dimension specified @@ -26,6 +25,9 @@ def mkvc(x, numDims=1): if hasattr(x, 'tovec'): x = x.tovec() + if isinstance(x, Zero): + return x + assert isinstance(x, np.ndarray), "Vector must be a numpy array" if numDims == 1: @@ -37,6 +39,9 @@ def mkvc(x, numDims=1): def sdiag(h): """Sparse diagonal matrix""" + if isinstance(h, Zero): + return Zero() + return sp.spdiags(mkvc(h), 0, h.size, h.size, format="csr") def sdInv(M): @@ -417,6 +422,12 @@ class Zero(object): def __ge__(self, v):return 0 >= v def __gt__(self, v):return 0 > v + @property + def transpose(self): return Zero() + + @property + def T(self): return Zero() + class Identity(object): _positive = True def __init__(self, positive=True): diff --git a/SimPEG/Utils/meshutils.py b/SimPEG/Utils/meshutils.py index 585fcc9a..eb5d13a1 100644 --- a/SimPEG/Utils/meshutils.py +++ b/SimPEG/Utils/meshutils.py @@ -102,223 +102,6 @@ def closestPoints(mesh, pts, gridLoc='CC'): return nodeInds -def readUBCTensorMesh(fileName): - """ - Read UBC GIF 3DTensor mesh and generate 3D Tensor mesh in simpegTD - - Input: - :param fileName, path to the UBC GIF mesh file - - Output: - :param SimPEG TensorMesh object - :return - """ - - # Interal function to read cell size lines for the UBC mesh files. - def readCellLine(line): - for seg in line.split(): - if '*' in seg: - st = seg - sp = seg.split('*') - re = np.array(sp[0],dtype=int)*(' ' + sp[1]) - line = line.replace(st,re.strip()) - return np.array(line.split(),dtype=float) - - # Read the file as line strings, remove lines with comment = ! - msh = np.genfromtxt(fileName,delimiter='\n',dtype=np.str,comments='!') - - # Fist line is the size of the model - sizeM = np.array(msh[0].split(),dtype=float) - # Second line is the South-West-Top corner coordinates. - x0 = np.array(msh[1].split(),dtype=float) - # Read the cell sizes - h1 = readCellLine(msh[2]) - h2 = readCellLine(msh[3]) - h3temp = readCellLine(msh[4]) - h3 = h3temp[::-1] # Invert the indexing of the vector to start from the bottom. - # Adjust the reference point to the bottom south west corner - x0[2] = x0[2] - np.sum(h3) - # Make the mesh - from SimPEG import Mesh - tensMsh = Mesh.TensorMesh([h1,h2,h3],x0) - return tensMsh - -def readUBCTensorModel(fileName, mesh): - """ - Read UBC 3DTensor mesh model and generate 3D Tensor mesh model in simpeg - - Input: - :param fileName, path to the UBC GIF mesh file to read - :param mesh, TensorMesh object, mesh that coresponds to the model - - Output: - :return numpy array, model with TensorMesh ordered - """ - f = open(fileName, 'r') - model = np.array(map(float, f.readlines())) - f.close() - model = np.reshape(model, (mesh.nCz, mesh.nCx, mesh.nCy), order = 'F') - model = model[::-1,:,:] - model = np.transpose(model, (1, 2, 0)) - model = mkvc(model) - - return model - -def writeUBCTensorMesh(fileName, mesh): - """ - Writes a SimPEG TensorMesh to a UBC-GIF format mesh file. - - :param str fileName: File to write to - :param simpeg.Mesh.TensorMesh mesh: The mesh - - """ - assert mesh.dim == 3 - s = '' - s += '%i %i %i\n' %tuple(mesh.vnC) - origin = mesh.x0 + np.array([0,0,mesh.hz.sum()]) # Have to it in the same operation or use mesh.x0.copy(), otherwise the mesh.x0 is updated. - origin.dtype = float - - s += '%.2f %.2f %.2f\n' %tuple(origin) - s += ('%.2f '*mesh.nCx+'\n')%tuple(mesh.hx) - s += ('%.2f '*mesh.nCy+'\n')%tuple(mesh.hy) - s += ('%.2f '*mesh.nCz+'\n')%tuple(mesh.hz[::-1]) - f = open(fileName, 'w') - f.write(s) - f.close() - -def writeUBCTensorModel(fileName, mesh, model): - """ - Writes a model associated with a SimPEG TensorMesh - to a UBC-GIF format model file. - - :param str fileName: File to write to - :param simpeg.Mesh.TensorMesh mesh: The mesh - :param numpy.ndarray model: The model - """ - - # Reshape model to a matrix - modelMat = mesh.r(model,'CC','CC','M') - # Transpose the axes - modelMatT = modelMat.transpose((2,0,1)) - # Flip z to positive down - modelMatTR = mkvc(modelMatT[::-1,:,:]) - - np.savetxt(fileName, modelMatTR.ravel()) - - -def readVTRFile(fileName): - """ - Read VTK Rectilinear (vtr xml file) and return SimPEG Tensor mesh and model - - Input: - :param vtrFileName, path to the vtr model file to write to - - Output: - :return SimPEG TensorMesh object - :return SimPEG model dictionary - - """ - # Import - from vtk import vtkXMLRectilinearGridReader as vtrFileReader - from vtk.util.numpy_support import vtk_to_numpy - - # Read the file - vtrReader = vtrFileReader() - vtrReader.SetFileName(fileName) - vtrReader.Update() - vtrGrid = vtrReader.GetOutput() - # Sort information - hx = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetXCoordinates()))) - xR = vtk_to_numpy(vtrGrid.GetXCoordinates())[0] - hy = np.abs(np.diff(vtk_to_numpy(vtrGrid.GetYCoordinates()))) - yR = vtk_to_numpy(vtrGrid.GetYCoordinates())[0] - zD = np.diff(vtk_to_numpy(vtrGrid.GetZCoordinates())) - # Check the direction of hz - if np.all(zD < 0): - hz = np.abs(zD[::-1]) - zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[-1] - else: - hz = np.abs(zD) - zR = vtk_to_numpy(vtrGrid.GetZCoordinates())[0] - x0 = np.array([xR,yR,zR]) - - # Make the SimPEG object - from SimPEG import Mesh - tensMsh = Mesh.TensorMesh([hx,hy,hz],x0) - - # Grap the models - modelDict = {} - for i in np.arange(vtrGrid.GetCellData().GetNumberOfArrays()): - modelName = vtrGrid.GetCellData().GetArrayName(i) - if np.all(zD < 0): - modFlip = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) - tM = tensMsh.r(modFlip,'CC','CC','M') - modArr = tensMsh.r(tM[:,:,::-1],'CC','CC','V') - else: - modArr = vtk_to_numpy(vtrGrid.GetCellData().GetArray(i)) - modelDict[modelName] = modArr - - # Return the data - return tensMsh, modelDict - -def writeVTRFile(fileName,mesh,model=None): - """ - Makes and saves a VTK rectilinear file (vtr) for a simpeg Tensor mesh and model. - - Input: - :param str, path to the output vtk file - :param mesh, SimPEG TensorMesh object - mesh to be transfer to VTK - :param model, dictionary of numpy.array - Name('s) and array('s). Match number of cells - - """ - # Import - from vtk import vtkRectilinearGrid as rectGrid, vtkXMLRectilinearGridWriter as rectWriter - from vtk.util.numpy_support import numpy_to_vtk - - # Deal with dimensionalities - if mesh.dim >= 1: - vX = mesh.vectorNx - xD = mesh.nNx - yD,zD = 1,1 - vY, vZ = np.array([0,0]) - if mesh.dim >= 2: - vY = mesh.vectorNy - yD = mesh.nNy - if mesh.dim == 3: - vZ = mesh.vectorNz - zD = mesh.nNz - # Use rectilinear VTK grid. - # Assign the spatial information. - vtkObj = rectGrid() - vtkObj.SetDimensions(xD,yD,zD) - vtkObj.SetXCoordinates(numpy_to_vtk(vX,deep=1)) - vtkObj.SetYCoordinates(numpy_to_vtk(vY,deep=1)) - vtkObj.SetZCoordinates(numpy_to_vtk(vZ,deep=1)) - - # Assign the model('s) to the object - for item in model.iteritems(): - # Convert numpy array - vtkDoubleArr = numpy_to_vtk(item[1],deep=1) - vtkDoubleArr.SetName(item[0]) - vtkObj.GetCellData().AddArray(vtkDoubleArr) - # Set the active scalar - vtkObj.GetCellData().SetActiveScalars(model.keys()[0]) - vtkObj.Update() - - - # Check the extension of the fileName - ext = os.path.splitext(fileName)[1] - if ext is '': - fileName = fileName + '.vtr' - elif ext not in '.vtr': - raise IOError('{:s} is an incorrect extension, has to be .vtr') - # Write the file. - vtrWriteFilter = rectWriter() - vtrWriteFilter.SetInput(vtkObj) - vtrWriteFilter.SetFileName(fileName) - vtrWriteFilter.Update() - - def ExtractCoreMesh(xyzlim, mesh, meshType='tensor'): """ Extracts Core Mesh from Global mesh diff --git a/docs/em/api_FDEM.rst b/docs/em/api_FDEM.rst index bf5bdcb4..778e0b3a 100644 --- a/docs/em/api_FDEM.rst +++ b/docs/em/api_FDEM.rst @@ -19,14 +19,14 @@ Electromagnetic phenomena are governed by Maxwell's equations. They describe the Fourier Transform Convention ---------------------------- -In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the \\(e^{i \\omega t} \\) convention, so we define our Fourier Transform pair as +In order to examine Maxwell's equations in the frequency domain, we must first define our choice of harmonic time-dependence by choosing a Fourier transform convention. We use the :math:`e^{i \omega t}` convention, so we define our Fourier Transform pair as .. math :: - F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\ + F(\omega) = \int_{-\infty}^{\infty} f(t) e^{- i \omega t} dt \\ - f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega + f(t) = \frac{1}{2\pi}\int_{-\infty}^{\infty} F(\omega) e^{i \omega t} d \omega -where \\(\\omega\\) is angular frequency, \\(t\\) is time, \\(F(\\omega)\\) is the function defined in the frequency domain and \\(f(t)\\) is the function defined in the time domain. +where :math:`\omega` is angular frequency, :math:`t` is time, :math:`F(\omega)` is the function defined in the frequency domain and :math:`f(t)` is the function defined in the time domain. Maxwell's Equations @@ -34,44 +34,46 @@ Maxwell's Equations In the frequency domain, Maxwell's equations are given by .. math :: - \curl \vec{E} = - i \omega \vec{B} \\ + \curl \vec{E} + i \omega \vec{B} = \vec{S_m}\\ - \curl \vec{H} = \vec{J} + i \omega \vec{D} + \vec{S} \\ + \curl \vec{H} - \vec{J} - i \omega \vec{D} = \vec{S_e} \\ - \div \vec{B} = 0 \\ + \div \vec{B} = 0 \\ - \div \vec{D} = \rho_f + \div \vec{D} = \rho_f where: - - \\(\\vec{E}\\) : electric field (\\(V/m\\)) - - \\(\\vec{H}\\) : magnetic field (\\(A/m\\)) - - \\(\\vec{B}\\) : magnetic flux density (\\(Wb/m^2\\)) - - \\(\\vec{D}\\) : electric displacement / electric flux density (\\(C/m^2\\)) - - \\(\\vec{J}\\) : electric current density (\\(A/m^2\\)) - - \\(\\rho_f\\) : free charge density + - :math:`\vec{E}` : electric field (:math:`V/m` ) + - :math:`\vec{H}` : magnetic field (:math:`A/m` ) + - :math:`\vec{B}` : magnetic flux density (:math:`Wb/m^2` ) + - :math:`\vec{D}` : electric displacement / electric flux density (:math:`C/m^2` ) + - :math:`\vec{J}` : electric current density (:math:`A/m^2` ) + - :math:`\vec{S_m}` : magnetic source term (:math:`V/m^2` ) + - :math:`\vec{S_e}` : electric source term (:math:`A/m^2` ) + - :math:`\rho_f` : free charge density (:math:`\Omega m` ) -The source term is \\(\\vec{S}\\) Constitutive Relations ---------------------- + The fields and fluxes are related through the constitutive relations. At each frequency, they are given by .. math :: - \vec{J} = \sigma \vec{E} \\ + \vec{J} = \sigma \vec{E} \\ - \vec{B} = \mu \vec{H} \\ + \vec{B} = \mu \vec{H} \\ - \vec{D} = \varepsilon \vec{E} + \vec{D} = \varepsilon \vec{E} where: -- \\(\\sigma\\) : electrical conductivity \\(S/m\\) -- \\(\\mu\\) : magnetic permeability \\(H/m\\) -- \\(\\varepsilon\\) : dielectric permittivity \\(F/m\\) +- :math:`\sigma` : electrical conductivity (:math:`S/m`) +- :math:`\mu` : magnetic permeability (:math:`H/m`) +- :math:`\varepsilon` : dielectric permittivity (:math:`F/m`) -\\(\\sigma\\), \\(\\mu\\), \\(\\varepsilon\\) are physical properties which depend on the material. \\(\\sigma\\) describes how easily electric current passes through a material, \\(\\mu\\) describes how easily a material is magnetized, and \\(\\varepsilon\\) describes how easily a material is electrically polarized. In most geophysical applications of EM, \\(\\sigma\\) is the the primary physical property of interest, and \\(\\mu\\), \\(\\varepsilon\\) are assumed to have their free-space values \\(\\mu_0 = 4\\pi \\times 10^{-7} H/m \\), \\(\\varepsilon_0 = 8.85 \\times 10^{-12} F/m\\) +:math:`\sigma`, :math:`\mu`, :math:`\varepsilon` are physical properties which depend on the material. :math:`\sigma` describes how easily electric current passes through a material, :math:`\mu` describes how easily a material is magnetized, and :math:`\varepsilon` describes how easily a material is electrically polarized. In most geophysical applications of EM, :math:`\sigma` is the the primary physical property of interest, and :math:`\mu`, :math:`\varepsilon` are assumed to have their free-space values :math:`\mu_0 = 4\pi \times 10^{-7} H/m` , :math:`\varepsilon_0 = 8.85 \times 10^{-12} F/m` Quasi-static Approximation @@ -80,8 +82,8 @@ Quasi-static Approximation For the frequency range typical of most geophysical surveys, the contribution of the electric displacement is negligible compared to the electric current density. In this case, we use the Quasi-static approximation and assume that this term can be neglected, giving .. math :: - \nabla \times \vec{E} = -i \omega \vec{B} \\ - \nabla \times \vec{H} = \vec{J} + \vec{S} + \nabla \times \vec{E} + i \omega \vec{B} = \vec{S_m} \\ + \nabla \times \vec{H} - \vec{J} = \vec{S_e} Implementation in SimPEG.EM @@ -90,14 +92,14 @@ Implementation in SimPEG.EM We consider two formulations in SimPEG.EM, both first-order and both in terms of one field and one flux. We allow for the definition of magnetic and electric sources (see for example: Ward and Hohmann, starting on page 144). The E-B formulation is in terms of the electric field and the magnetic flux: .. math :: - \nabla \times \vec{E} + i \omega \vec{B} = \vec{S}_m \\ - \nabla \times \mu^{-1} \vec{B} - \sigma \vec{E} = \vec{S}_e + \nabla \times \vec{E} + i \omega \vec{B} = \vec{S}_m \\ + \nabla \times \mu^{-1} \vec{B} - \sigma \vec{E} = \vec{S}_e The H-J formulation is in terms of the current density and the magnetic field: .. math :: - \nabla \times \sigma^{-1} \vec{J} + i \omega \mu \vec{H} = \vec{S}_m \\ - \nabla \times \vec{H} - \vec{J} = \vec{S}_e + \nabla \times \sigma^{-1} \vec{J} + i \omega \mu \vec{H} = \vec{S}_m \\ + \nabla \times \vec{H} - \vec{J} = \vec{S}_e Discretizing @@ -106,34 +108,34 @@ For both formulations, we use a finite volume discretization and discretize fields on cell edges, fluxes on cell faces and physical properties in cell centers. This is particularly important when using symmetry to reduce the dimensionality of a problem -(for instance on a 2D CylMesh, there are \\(r\\), \\(z\\) faces and \\(\\theta\\) edges) +(for instance on a 2D CylMesh, there are :math:`r`, :math:`z` faces and :math:`\theta` edges) .. figure:: ../images/finitevolrealestate.png - :align: center - :scale: 60 % + :align: center + :scale: 60 % For the two formulations, the discretization of the physical properties, fields and fluxes are summarized below. .. figure:: ../images/ebjhdiscretizations.png - :align: center - :scale: 60 % + :align: center + :scale: 60 % -Note that resistivity is the inverse of conductivity, \\(\\rho = \\sigma^{-1}\\). +Note that resistivity is the inverse of conductivity, :math:`\rho = \sigma^{-1}`. -E-B Formulation: -**************** +E-B Formulation +--------------- .. math :: - \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\ - \mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{M^e} \mathbf{s_e} + \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\ + \mathbf{C^T} \mathbf{M^f_{\mu^{-1}}} \mathbf{b} - \mathbf{M^e_\sigma} \mathbf{e} = \mathbf{M^e} \mathbf{s_e} -H-J Formulation: -**************** +H-J Formulation +--------------- .. math :: - \mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\ - \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} + \mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\ + \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} .. Forward Problem @@ -144,6 +146,10 @@ H-J Formulation: API === + +FDEM Problem +------------ + .. automodule:: SimPEG.EM.FDEM.FDEM :show-inheritance: :members: @@ -157,3 +163,17 @@ FDEM Survey :show-inheritance: :members: :undoc-members: + +.. automodule:: SimPEG.EM.FDEM.SrcFDEM + :show-inheritance: + :members: + :undoc-members: + +FDEM Fields +----------- + +.. automodule:: SimPEG.EM.FDEM.FieldsFDEM + :show-inheritance: + :members: + :undoc-members: + diff --git a/docs/em/api_TDEM.rst b/docs/em/api_TDEM.rst index cbbc48b8..fe3dc613 100644 --- a/docs/em/api_TDEM.rst +++ b/docs/em/api_TDEM.rst @@ -48,6 +48,305 @@ \newcommand{\I}{\vec{I}} +Time Domain Electromagnetics +**************************** + +.. _api_TDEM_derivation: + +Time-Domain EM Derivation +========================= + +The following shows the derivation for the TDEM problem. We use the b-formulation below. +(More to come soon..!) + + +Sensitivity Calculation +----------------------- + +.. math:: + + \begin{align} + \dcurl \e^{(t+1)} + \frac{\b^{(t+1)} - \b^{(t)}}{\delta t} = 0 \\ + \dcurl^\top \MfMui \b^{(t+1)} - \MeSig \e^{(t+1)} = \Me \j_s^{(t+1)} + \end{align} + +Using Gauss-Newton to solve the inverse problem requires the ability to calculate the product of the +Jacobian and a vector, as well as the transpose of the Jacobian times a vector. +The above system can be rewritten as: + +.. math:: + + \begin{align} + \mathbf{A} \u^{(t+1)} + \mathbf{B} \u^{(t)}= \s^{(t+1)} + \end{align} + +where + +.. math:: + + \begin{align} + \mathbf{A} = + \left[ + \begin{array}{cc} + \frac{1}{\delta t} \MfMui & \MfMui\dcurl \\ + \dcurl^\top \MfMui & -\MeSig + \end{array} + \right] \\ + \mathbf{B} = + \left[ + \begin{array}{cc} + -\frac{1}{\delta t} \MfMui & 0 \\ + 0 & 0 + \end{array} + \right] \\ + \u^{(k)} = \left[ + \begin{array}{c} + \b^{(k)}\\ + \e^{(k)} + \end{array} + \right] \\ + \s^{(k)} = \left[ + \begin{array}{c} + 0\\ + \Me \j^{(k)}_s + \end{array} + \right] + \end{align} + +.. note:: + + Here we have multiplied through by \\(\\MfMui\\) to make A and B symmetric! + +The entire time dependent system can be written in a single matrix expression + +.. math:: + + \begin{align} + \hat{\mathbf{A}} \hat{u} = \hat{s} + \end{align} + +where + +.. math:: + + \begin{align} + \mathbf{\hat{A}} = \left[ + \begin{array}{cccc} + A & 0 & & \\ + B & A & & \\ + & \ddots & \ddots & \\ + & & B & A + \end{array} + \right] \\ + \hat{u} = \left[ + \begin{array}{c} + \u^{(1)} \\ + \u^{(2)} \\ + \vdots \\ + \u^{(N)} + \end{array} \right]\\ + \hat{s} = \left[ + \begin{array}{c} + \s^{(1)} - \mathbf{B} \u^{(0)} \\ + \s^{(2)} \\ + \vdots \\ + \s^{(N)} + \end{array} + \right] + \end{align} + +For the fields \\(\\u\\), the measured data is given by + +.. math:: + + \begin{align} + \vec{d} = \mathbf{Q} \u + \end{align} + +The sensitivity matrix **J** is then defined as + +.. math:: + + \begin{align} + \mathbf{J} = \mathbf{Q} \frac{\partial \u}{\partial \sigma} + \end{align} + + +Defining the function \\(\\c(m,\\u)\\) to be + +.. math:: + + \begin{align} + \vec{c}(m,\u) = \hat{\mathbf{A}} \vec{u} - \vec{q} = \vec{0} + \end{align} + +then + +.. math:: + + \begin{align} + \frac{\partial \vec{c}}{\partial m} \partial m + + \frac{\partial \vec{c}}{\partial \u} \partial \vec{u} = 0 + \end{align} + +or + +.. math:: + + \begin{align} + \frac{\partial \vec{u}}{\partial m} = -\left(\frac{\partial \vec{c}}{\partial \u} \right)^{-1} \frac{\partial \vec{c}}{\partial m} + \end{align} + + +Differentiating, we find that + +.. math:: + + \begin{align} + \frac{\partial \vec{c}}{\partial \hat{u}} = \hat{\mathbf{A}} + \end{align} + +and + +.. math:: + + \begin{align} + \frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma = + \left[ + \begin{array}{c} + g_\sigma^{(1)}\\ + g_\sigma^{(2)}\\ + \vdots \\ + g_\sigma^{(N)} + \end{array} + \right] + \end{align} + +with + +.. math:: + + \begin{align} + g_\sigma^{(n)} = + \left[ + \begin{array}{c} + \mathbf{0} \\ + - \diag{\e^{(n)}} \Ace \diag{\vec{V}} + \end{array} + \right] + \end{align} + + +Implementing **J** times a vector +--------------------------------- + +Multiplying **J** onto a vector can be broken into three steps + + +* Compute \\(\\vec{p} = \\mathbf{G}m\\) +* Solve \\(\\hat{\\mathbf{A}} \\vec{y} = \\vec{p}\\) +* Compute \\(\\vec{w} = -\\mathbf{Q} \\vec{y}\\) + +.. math:: + + \begin{align} + \vec{p}^{(n)} = \left[ + \begin{array}{c} + \vec{p}_b^{(n)} \\ + \vec{p}_e^{(n)} + \end{array} + \right] \\ + \vec{p}_b^{(n)} = 0 \\ + \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m + \end{align} + + +For all time steps: + +.. math:: + + \begin{align} + \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t+1)} + \MfMui\dcurl \vec{y}_{e}^{(t+1)} + - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)} + = \vec{p}_b^{(t+1)} \\ + \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig \vec{y}_e^{(t+1)} = \vec{p}_e^{(t+1)} + \end{align} + +and + +.. math:: + + \begin{align} + \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t+1)} = + \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} + + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} + \vec{p}_b^{(t+1)} \\ + \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)} + \end{align} + +.. note:: + + For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero. + + + + +Implementing **J** transpose times a vector +------------------------------------------- + +Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps + + +* Compute \\(\\vec{p} = \\mathbf{Q}^\\top \\vec{v}\\) +* Solve \\(\\hat{\\mathbf{A}}^\\top \\vec{y} = \\vec{p}\\) +* Compute \\(\\vec{w} = -\\mathbf{G}^\\top y\\) + + +.. math:: + + \mathbf{\hat{A}}^\top = \left[ + \begin{array}{cccc} + A & B & & \\ + & \ddots & \ddots & \\ + & & A & B \\ + & & 0 & A + \end{array} + \right] + +For the all time-steps (going backwards in time): + + +.. math:: + + A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)} + + +.. math:: + + \begin{align} + \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)} + - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)} + = \vec{p}_b^{(t)} \\ + \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)} + \end{align} + +and + +.. math:: + + \begin{align} + \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} = + \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)} + + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\ + \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)} + \end{align} + + +.. note:: + + For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero. + + + TDEM - B formulation ==================== diff --git a/docs/em/api_TDEM_derivation.rst b/docs/em/api_TDEM_derivation.rst deleted file mode 100644 index af3fc2fc..00000000 --- a/docs/em/api_TDEM_derivation.rst +++ /dev/null @@ -1,341 +0,0 @@ -.. _api_TDEM_derivation: - - -.. math:: - - \renewcommand{\div}{\nabla\cdot\,} - \newcommand{\grad}{\vec \nabla} - \newcommand{\curl}{{\vec \nabla}\times\,} - \newcommand {\J}{{\vec J}} - \renewcommand{\H}{{\vec H}} - \newcommand {\E}{{\vec E}} - \newcommand{\dcurl}{{\mathbf C}} - \newcommand{\dgrad}{{\mathbf G}} - \newcommand{\Acf}{{\mathbf A_c^f}} - \newcommand{\Ace}{{\mathbf A_c^e}} - \renewcommand{\S}{{\mathbf \Sigma}} - \newcommand{\St}{{\mathbf \Sigma_\tau}} - \newcommand{\T}{{\mathbf T}} - \newcommand{\Tt}{{\mathbf T_\tau}} - \newcommand{\diag}[1]{\,{\sf diag}\left( #1 \right)} - \newcommand{\M}{{\mathbf M}} - \newcommand{\MfMui}{{\M^f_{\mu^{-1}}}} - \newcommand{\MeSig}{{\M^e_\sigma}} - \newcommand{\MeSigInf}{{\M^e_{\sigma_\infty}}} - \newcommand{\MeSigO}{{\M^e_{\sigma_0}}} - \newcommand{\Me}{{\M^e}} - \newcommand{\Mes}[1]{{\M^e_{#1}}} - \newcommand{\Mee}{{\M^e_e}} - \newcommand{\Mej}{{\M^e_j}} - \newcommand{\BigO}[1]{\mathcal{O}\bigl(#1\bigr)} - \newcommand{\bE}{\mathbf{E}} - \newcommand{\bH}{\mathbf{H}} - \newcommand{\B}{\vec{B}} - \newcommand{\D}{\vec{D}} - \renewcommand{\H}{\vec{H}} - \newcommand{\s}{\vec{s}} - \newcommand{\bfJ}{\bf{J}} - \newcommand{\vecm}{\vec m} - \renewcommand{\Re}{\mathsf{Re}} - \renewcommand{\Im}{\mathsf{Im}} - \renewcommand {\j} { {\vec j} } - \newcommand {\h} { {\vec h} } - \renewcommand {\b} { {\vec b} } - \newcommand {\e} { {\vec e} } - \newcommand {\c} { {\vec c} } - \renewcommand {\d} { {\vec d} } - \renewcommand {\u} { {\vec u} } - \newcommand{\I}{\vec{I}} - - -Time-Domain EM Derivation -************************* - -The following shows the derivation for the TDEM problem. We use the b-formulation below. -(More to come soon..!) - - -Sensitivity Calculation -======================= - -.. math:: - - \begin{align} - \dcurl \e^{(t+1)} + \frac{\b^{(t+1)} - \b^{(t)}}{\delta t} = 0 \\ - \dcurl^\top \MfMui \b^{(t+1)} - \MeSig \e^{(t+1)} = \Me \j_s^{(t+1)} - \end{align} - -Using Gauss-Newton to solve the inverse problem requires the ability to calculate the product of the -Jacobian and a vector, as well as the transpose of the Jacobian times a vector. -The above system can be rewritten as: - -.. math:: - - \begin{align} - \mathbf{A} \u^{(t+1)} + \mathbf{B} \u^{(t)}= \s^{(t+1)} - \end{align} - -where - -.. math:: - - \begin{align} - \mathbf{A} = - \left[ - \begin{array}{cc} - \frac{1}{\delta t} \MfMui & \MfMui\dcurl \\ - \dcurl^\top \MfMui & -\MeSig - \end{array} - \right] \\ - \mathbf{B} = - \left[ - \begin{array}{cc} - -\frac{1}{\delta t} \MfMui & 0 \\ - 0 & 0 - \end{array} - \right] \\ - \u^{(k)} = \left[ - \begin{array}{c} - \b^{(k)}\\ - \e^{(k)} - \end{array} - \right] \\ - \s^{(k)} = \left[ - \begin{array}{c} - 0\\ - \Me \j^{(k)}_s - \end{array} - \right] - \end{align} - -.. note:: - - Here we have multiplied through by \\(\\MfMui\\) to make A and B symmetric! - -The entire time dependent system can be written in a single matrix expression - -.. math:: - - \begin{align} - \hat{\mathbf{A}} \hat{u} = \hat{s} - \end{align} - -where - -.. math:: - - \begin{align} - \mathbf{\hat{A}} = \left[ - \begin{array}{cccc} - A & 0 & & \\ - B & A & & \\ - & \ddots & \ddots & \\ - & & B & A - \end{array} - \right] \\ - \hat{u} = \left[ - \begin{array}{c} - \u^{(1)} \\ - \u^{(2)} \\ - \vdots \\ - \u^{(N)} - \end{array} \right]\\ - \hat{s} = \left[ - \begin{array}{c} - \s^{(1)} - \mathbf{B} \u^{(0)} \\ - \s^{(2)} \\ - \vdots \\ - \s^{(N)} - \end{array} - \right] - \end{align} - -For the fields \\(\\u\\), the measured data is given by - -.. math:: - - \begin{align} - \vec{d} = \mathbf{Q} \u - \end{align} - -The sensitivity matrix **J** is then defined as - -.. math:: - - \begin{align} - \mathbf{J} = \mathbf{Q} \frac{\partial \u}{\partial \sigma} - \end{align} - - -Defining the function \\(\\c(m,\\u)\\) to be - -.. math:: - - \begin{align} - \vec{c}(m,\u) = \hat{\mathbf{A}} \vec{u} - \vec{q} = \vec{0} - \end{align} - -then - -.. math:: - - \begin{align} - \frac{\partial \vec{c}}{\partial m} \partial m - + \frac{\partial \vec{c}}{\partial \u} \partial \vec{u} = 0 - \end{align} - -or - -.. math:: - - \begin{align} - \frac{\partial \vec{u}}{\partial m} = -\left(\frac{\partial \vec{c}}{\partial \u} \right)^{-1} \frac{\partial \vec{c}}{\partial m} - \end{align} - - -Differentiating, we find that - -.. math:: - - \begin{align} - \frac{\partial \vec{c}}{\partial \hat{u}} = \hat{\mathbf{A}} - \end{align} - -and - -.. math:: - - \begin{align} - \frac{\partial \vec{c}}{\partial \sigma} = \mathbf{G}_\sigma = - \left[ - \begin{array}{c} - g_\sigma^{(1)}\\ - g_\sigma^{(2)}\\ - \vdots \\ - g_\sigma^{(N)} - \end{array} - \right] - \end{align} - -with - -.. math:: - - \begin{align} - g_\sigma^{(n)} = - \left[ - \begin{array}{c} - \mathbf{0} \\ - - \diag{\e^{(n)}} \Ace \diag{\vec{V}} - \end{array} - \right] - \end{align} - - -Implementing **J** times a vector -================================= - -Multiplying **J** onto a vector can be broken into three steps - - -* Compute \\(\\vec{p} = \\mathbf{G}m\\) -* Solve \\(\\hat{\\mathbf{A}} \\vec{y} = \\vec{p}\\) -* Compute \\(\\vec{w} = -\\mathbf{Q} \\vec{y}\\) - -.. math:: - - \begin{align} - \vec{p}^{(n)} = \left[ - \begin{array}{c} - \vec{p}_b^{(n)} \\ - \vec{p}_e^{(n)} - \end{array} - \right] \\ - \vec{p}_b^{(n)} = 0 \\ - \vec{p}_e^{(n)} = - \diag{\e^{(n)}} \Ace \diag{V} m - \end{align} - - -For all time steps: - -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t+1)} + \MfMui\dcurl \vec{y}_{e}^{(t+1)} - - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t)} - = \vec{p}_b^{(t+1)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig \vec{y}_e^{(t+1)} = \vec{p}_e^{(t+1)} - \end{align} - -and - -.. math:: - - \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t+1)} = - \frac{1}{\delta t} \MfMui \vec{y}_b^{(t)} - + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t+1)} + \vec{p}_b^{(t+1)} \\ - \vec{y}_e^{(t+1)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t+1)} - \MeSig^{-1} \vec{p}_e^{(t+1)} - \end{align} - -.. note:: - - For the first time step, \\\(t=0\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(0)}\\\) is zero. - - - - -Implementing **J** transpose times a vector -=========================================== - -Multiplying \\(\\mathbf{J}^\\top\\) onto a vector can be broken into three steps - - -* Compute \\(\\vec{p} = \\mathbf{Q}^\\top \\vec{v}\\) -* Solve \\(\\hat{\\mathbf{A}}^\\top \\vec{y} = \\vec{p}\\) -* Compute \\(\\vec{w} = -\\mathbf{G}^\\top y\\) - - -.. math:: - - \mathbf{\hat{A}}^\top = \left[ - \begin{array}{cccc} - A & B & & \\ - & \ddots & \ddots & \\ - & & A & B \\ - & & 0 & A - \end{array} - \right] - -For the all time-steps (going backwards in time): - - -.. math:: - - A \vec{y}^{(t)} + B \vec{y}^{(t+1)} = \vec{p}^{(t)} - - -.. math:: - - \begin{align} - \frac{1}{\delta t} \MfMui\vec{y}_{b}^{(t)} + \MfMui\dcurl \vec{y}_{e}^{(t)} - - \frac{1}{\delta t} \MfMui \vec{y}_{b}^{(t+1)} - = \vec{p}_b^{(t)} \\ - \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig \vec{y}_e^{(t)} = \vec{p}_e^{(t)} - \end{align} - -and - -.. math:: - - \begin{align} - \left( \MfMui \dcurl \MeSig^{-1} \dcurl^\top \MfMui + \frac{1}{\delta t} \MfMui \right) \vec{y}_{b}^{(t)} = - \frac{1}{\delta t} \MfMui \vec{y}_b^{(t+1)} - + \MfMui \dcurl \MeSig^{-1} \vec{p}_e^{(t)} + \vec{p}_b^{(t)} \\ - \vec{y}_e^{(t)} = \MeSig^{-1} \dcurl^\top \MfMui \vec{y}_b^{(t)} - \MeSig^{-1} \vec{p}_e^{(t)} - \end{align} - - -.. note:: - - For the last time step, \\\(t=N\\\), the term: \\\(\\frac{1}{\\delta t} \\MfMui \\vec{y}_b^{(N+1)}\\\) is zero. diff --git a/docs/em/api_Utils.rst b/docs/em/api_Utils.rst index 8ae98855..ac8f9d34 100644 --- a/docs/em/api_Utils.rst +++ b/docs/em/api_Utils.rst @@ -4,6 +4,16 @@ simpegEM Utilities SimPEG for EM provides a few EM specific utility codes, sources, and analytic functions. +Utilities for Electromagnetics +============================== + +.. automodule:: SimPEG.EM.Utils + :show-inheritance: + :members: + :undoc-members: + :inherited-members: + + Analytic Functions - Time ========================= @@ -22,12 +32,3 @@ Analytic Functions - Frequency :members: :undoc-members: :inherited-members: - - -Sources -======= - -.. autoclass:: SimPEG.EM.FDEM.SrcFDEM.MagDipole - :show-inheritance: - :members: - :undoc-members: diff --git a/docs/em/index.rst b/docs/em/index.rst index fdf4dc19..a86ebb69 100644 --- a/docs/em/index.rst +++ b/docs/em/index.rst @@ -3,42 +3,24 @@ Electromagnetics ================ `SimPEG.EM` uses SimPEG as the framework for the forward and inverse -electromagnetics geophysical problems. +electromagnetics geophysical problems. -Time Domian Electromagnetics ----------------------------- - -.. toctree:: - :maxdepth: 2 - - api_TDEM_derivation +To solve for predicted data, we follow the framework shown below. The model is +what we invert for. This is mapped to a physical property on the simulation +mesh. A source which is used to excite the system is specified. Having a model +and a source, we can solve Maxwell's equations for fields. We sample these +fields with recievers to give us predicted data. -Code for Time Domian Electromagnetics -------------------------------------- +.. image:: ../images/simpegEM_noMath.png + :scale: 50% -.. toctree:: - :maxdepth: 2 - - api_TDEM - -Frequency Domian Electromagnetics ---------------------------------- .. toctree:: :maxdepth: 2 api_FDEM - - -Utility Codes -------------- - -.. toctree:: - :maxdepth: 2 - + api_TDEM api_Utils - - diff --git a/docs/examples/EM_FDEM_1D_Inversion.rst b/docs/examples/EM_FDEM_1D_Inversion.rst new file mode 100644 index 00000000..acbc8cdc --- /dev/null +++ b/docs/examples/EM_FDEM_1D_Inversion.rst @@ -0,0 +1,26 @@ +.. _examples_EM_FDEM_1D_Inversion: + +.. --------------------------------- .. +.. .. +.. THIS FILE IS AUTO GENEREATED .. +.. .. +.. SimPEG/Examples/__init__.py .. +.. .. +.. --------------------------------- .. + + +EM: FDEM: 1D: Inversion +======================= + +Here we will create and run a FDEM 1D inversion. + + + +.. plot:: + + from SimPEG import Examples + Examples.EM_FDEM_1D_Inversion.run() + +.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_1D_Inversion.py + :language: python + :linenos: diff --git a/docs/flow/index.rst b/docs/flow/index.rst index 3c8083fc..d4299b0d 100644 --- a/docs/flow/index.rst +++ b/docs/flow/index.rst @@ -41,7 +41,7 @@ Here we reproduce the results from Celia et al. (1990): Richards ======== -.. automodule:: simpegFLOW.Richards.Empirical +.. automodule:: SimPEG.FLOW.Richards.Empirical :show-inheritance: :members: :undoc-members: diff --git a/docs/images/simpegEM_noMath.png b/docs/images/simpegEM_noMath.png new file mode 100644 index 00000000..958f7003 Binary files /dev/null and b/docs/images/simpegEM_noMath.png differ diff --git a/docs/images/simpegEM_sensitivity_J_JTvec.png b/docs/images/simpegEM_sensitivity_J_JTvec.png new file mode 100644 index 00000000..f2e2d0e4 Binary files /dev/null and b/docs/images/simpegEM_sensitivity_J_JTvec.png differ diff --git a/docs/images/simpegEM_withMath.png b/docs/images/simpegEM_withMath.png new file mode 100644 index 00000000..4571c058 Binary files /dev/null and b/docs/images/simpegEM_withMath.png differ diff --git a/tests/base/test_regularization.py b/tests/base/test_regularization.py index af7da692..050c46ac 100644 --- a/tests/base/test_regularization.py +++ b/tests/base/test_regularization.py @@ -4,11 +4,17 @@ from SimPEG import * from scipy.sparse.linalg import dsolve import inspect +TOL = 1e-20 class RegularizationTests(unittest.TestCase): def setUp(self): - self.mesh2 = Mesh.TensorMesh([3, 2]) + hx, hy, hz = np.random.rand(10), np.random.rand(9), np.random.rand(8) + hx, hy, hz = hx/hx.sum(), hy/hy.sum(), hz/hz.sum() + mesh1 = Mesh.TensorMesh([hx]) + mesh2 = Mesh.TensorMesh([hx, hy]) + mesh3 = Mesh.TensorMesh([hx, hy, hz]) + self.meshlist = [mesh1,mesh2, mesh3] def test_regularization(self): for R in dir(Regularization): @@ -16,18 +22,63 @@ class RegularizationTests(unittest.TestCase): if not inspect.isclass(r): continue if not issubclass(r, Regularization.BaseRegularization): continue - # if 'Regularization' not in R: continue - mapping = r.mapPair(self.mesh2) - reg = r(self.mesh2, mapping=mapping) - m = np.random.rand(mapping.nP) - reg.mref = m[:]*np.mean(m) - print 'Check:', R - passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) - self.assertTrue(passed) - print 'Check 2 Deriv:', R - passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) - self.assertTrue(passed) + for i, mesh in enumerate(self.meshlist): + + print 'Testing %iD'%mesh.dim + + mapping = r.mapPair(mesh) + reg = r(mesh, mapping=mapping) + m = np.random.rand(mapping.nP) + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) + + def test_regularization_ActiveCells(self): + for R in dir(Regularization): + r = getattr(Regularization, R) + if not inspect.isclass(r): continue + if not issubclass(r, Regularization.BaseRegularization): + continue + + for i, mesh in enumerate(self.meshlist): + + print 'Testing Active Cells %iD'%(mesh.dim) + + if mesh.dim == 1: + indAct = Utils.mkvc(mesh.gridCC <= 0.8) + elif mesh.dim == 2: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5) + elif mesh.dim == 3: + indAct = Utils.mkvc(mesh.gridCC[:,-1] <= 2*np.sin(2*np.pi*mesh.gridCC[:,0])+0.5 * 2*np.sin(2*np.pi*mesh.gridCC[:,1])+0.5) + + mapping = Maps.IdentityMap(nP=indAct.nonzero()[0].size) + + reg = r(mesh, mapping=mapping, indActive=indAct) + m = np.random.rand(mesh.nC)[indAct] + reg.mref = np.ones_like(m)*np.mean(m) + + print 'Check: phi_m (mref) = %f' %reg.eval(reg.mref) + passed = reg.eval(reg.mref) < TOL + self.assertTrue(passed) + + print 'Check:', R + passed = Tests.checkDerivative(lambda m : [reg.eval(m), reg.evalDeriv(m)], m, plotIt=False) + self.assertTrue(passed) + + print 'Check 2 Deriv:', R + passed = Tests.checkDerivative(lambda m : [reg.evalDeriv(m), reg.eval2Deriv(m)], m, plotIt=False) + self.assertTrue(passed) if __name__ == '__main__': diff --git a/tests/examples/test_examples.py b/tests/examples/test_examples.py index 2e4803b1..edb5600c 100644 --- a/tests/examples/test_examples.py +++ b/tests/examples/test_examples.py @@ -1,8 +1,28 @@ import unittest import sys +import os from SimPEG import Examples import numpy as np +class compareInitFiles(unittest.TestCase): + def test_compareInitFiles(self): + print 'Checking that __init__.py up-to-date in SimPEG/Examples' + fName = os.path.abspath(__file__) + ExamplesDir = os.path.sep.join(fName.split(os.path.sep)[:-3] + ['SimPEG', 'Examples']) + + files = os.listdir(ExamplesDir) + + pyfiles = [] + [pyfiles.append(py.rstrip('.py')) for py in files if py.endswith('.py') and py != '__init__.py'] + + setdiff = set(pyfiles) - set(Examples.__examples__) + + print ' Any missing files? ', setdiff + + didpass = (setdiff == set()) + + self.assertTrue(didpass, "Examples not up to date, run 'python __init__.py' from SimPEG/Examples to update") + def get(test): def test_func(self): print '\nTesting %s.run(plotIt=False)\n'%test @@ -10,11 +30,11 @@ def get(test): self.assertTrue(True) return test_func attrs = dict() + for test in Examples.__examples__: attrs['test_'+test] = get(test) TestExamples = type('TestExamples', (unittest.TestCase,), attrs) - if __name__ == '__main__': unittest.main() diff --git a/tests/mesh/test_MeshIO.py b/tests/mesh/test_MeshIO.py new file mode 100644 index 00000000..e8dc0748 --- /dev/null +++ b/tests/mesh/test_MeshIO.py @@ -0,0 +1,100 @@ +import numpy as np +import unittest, os +import SimPEG as simpeg +from SimPEG.Mesh import TensorMesh, TreeMesh + + +class TestTensorMeshIO(unittest.TestCase): + + def setUp(self): + h = np.ones(16) + mesh = TensorMesh([h,2*h,3*h]) + self.mesh = mesh + + def test_UBCfiles(self): + + mesh = self.mesh + # Make a vector + vec = np.arange(mesh.nC) + # Write and read + mesh.writeUBC('temp.msh', {'arange.txt':vec}) + meshUBC = TensorMesh.readUBC('temp.msh') + vecUBC = meshUBC.readModelUBC('arange.txt') + + # The mesh + assert mesh.__str__() == meshUBC.__str__() + assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0 + assert np.sum(vec - vecUBC) == 0 + assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0) + + + vecUBC = mesh.readModelUBC('arange.txt') + assert np.sum(vec - vecUBC) == 0 + + mesh.writeModelUBC('arange2.txt', vec + 1) + vec2UBC = mesh.readModelUBC('arange2.txt') + assert np.sum(vec + 1 - vec2UBC) == 0 + + print 'IO of UBC tensor mesh files is working' + os.remove('temp.msh') + os.remove('arange.txt') + os.remove('arange2.txt') + + def test_VTKfiles(self): + mesh = self.mesh + vec = np.arange(mesh.nC) + + mesh.writeVTK('temp.vtr', {'arange.txt':vec}) + meshVTR, models = TensorMesh.readVTK('temp.vtr') + + assert mesh.__str__() == meshVTR.__str__() + assert np.all(np.array(mesh.h) - np.array(meshVTR.h) == 0) + + assert 'arange.txt' in models + vecVTK = models['arange.txt'] + assert np.sum(vec - vecVTK) == 0 + + print 'IO of VTR tensor mesh files is working' + os.remove('temp.vtr') + + +class TestOcTreeMeshIO(unittest.TestCase): + + def setUp(self): + h = np.ones(16) + mesh = TreeMesh([h,2*h,3*h]) + mesh.refine(3) + mesh._refineCell([0,0,0,3]) + mesh._refineCell([0,2,0,3]) + self.mesh = mesh + + def test_UBCfiles(self): + + mesh = self.mesh + # Make a vector + vec = np.arange(mesh.nC) + # Write and read + mesh.writeUBC('temp.msh', {'arange.txt':vec}) + meshUBC = TreeMesh.readUBC('temp.msh') + vecUBC = meshUBC.readModelUBC('arange.txt') + + # The mesh + assert mesh.__str__() == meshUBC.__str__() + assert np.sum(mesh.gridCC - meshUBC.gridCC) == 0 + assert np.sum(vec - vecUBC) == 0 + assert np.all(np.array(mesh.h) - np.array(meshUBC.h) == 0) + print 'IO of UBC octree files is working' + os.remove('temp.msh') + os.remove('arange.txt') + + def test_VTUfiles(self): + mesh = self.mesh + vec = np.arange(mesh.nC) + mesh.writeVTK('temp.vtu',{'arange':vec}) + print 'Writing of VTU files is working' + os.remove('temp.vtu') + + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/mesh/test_TreeMesh.py b/tests/mesh/test_TreeMesh.py index e624ce87..afad27d1 100644 --- a/tests/mesh/test_TreeMesh.py +++ b/tests/mesh/test_TreeMesh.py @@ -26,6 +26,27 @@ class TestSimpleQuadTree(unittest.TestCase): assert np.allclose(np.r_[M._areaFxFull, M._areaFyFull], M._deflationMatrix('F') * M.area) + def test_getitem(self): + M = Mesh.TreeMesh([4,4]) + M.refine(1) + assert M.nC == 4 + assert len(M) == M.nC + assert np.allclose(M[0].center, [0.25,0.25]) + actual = [[0,0],[0.5,0],[0,0.5],[0.5,0.5]] + for i, n in enumerate(M[0].nodes): + assert np.allclose(M._gridN[n,:], actual[i]) + + def test_getitem3D(self): + M = Mesh.TreeMesh([4,4,4]) + M.refine(1) + assert M.nC == 8 + assert len(M) == M.nC + assert np.allclose(M[0].center, [0.25,0.25,0.25]) + actual = [[0,0,0],[0.5,0,0],[0,0.5,0],[0.5,0.5,0], + [0,0,0.5],[0.5,0,0.5],[0,0.5,0.5],[0.5,0.5,0.5]] + for i, n in enumerate(M[0].nodes): + assert np.allclose(M._gridN[n,:], actual[i]) + def test_refine(self): M = Mesh.TreeMesh([4,4,4]) M.refine(1) diff --git a/tests/utils/test_Zero.py b/tests/utils/test_Zero.py index 594de6a6..7b3c6e5d 100644 --- a/tests/utils/test_Zero.py +++ b/tests/utils/test_Zero.py @@ -1,5 +1,5 @@ import unittest -from SimPEG.Utils import Zero, Identity, sdiag +from SimPEG.Utils import Zero, Identity, sdiag, mkvc from SimPEG import np, sp class Tests(unittest.TestCase): @@ -29,6 +29,11 @@ class Tests(unittest.TestCase): assert a == 1 self.assertRaises(ZeroDivisionError, lambda:3/z) + assert mkvc(z) == 0 + assert sdiag(z)*a == 0 + assert z.T == 0 + assert z.transpose == 0 + def test_mat_zero(self): z = Zero() S = sdiag(np.r_[2,3])