diff --git a/docs/api_FDEM.rst b/docs/api_FDEM.rst index 2c063216..1848c569 100644 --- a/docs/api_FDEM.rst +++ b/docs/api_FDEM.rst @@ -36,7 +36,7 @@ In the frequency domain, Maxwell's equations are given by .. math :: \curl \vec{E} = - i \omega \vec{B} \\ - \curl \vec{H} = \vec{J} + i \omega \vec{D} + \vec{J}_s \\ + \curl \vec{H} = \vec{J} + i \omega \vec{D} + \vec{S} \\ \div \vec{B} = 0 \\ @@ -50,7 +50,7 @@ where: - \\(\\vec{D}\\) : electric displacement / electric flux density (\\(C/m^2\\)) - \\(\\vec{J}\\) : electric current density (\\(A/m^2\\)) - \\(\\rho_f\\) : free charge density -The source term is \\(\\vec{J}_s\\) +The source term is \\(\\vec{S}\\) Constitutive Relations @@ -71,7 +71,6 @@ where: - \\(\\varepsilon\\) : dielectric permittivity \\(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\\) -For a more complete discussion of physical properties see `GPG `_ Quasi-static Approximation @@ -80,17 +79,65 @@ For the frequency range typical of most geophysical surveys, the contribution of .. math :: \nabla \times \vec{E} = -i \omega \vec{B} \\ - \nabla \times \vec{H} = \vec{J} + \vec{J}_s + \nabla \times \vec{H} = \vec{J} + \vec{S} -Fields from a Dipole --------------------- +Implementation in simpegEM +========================== +We consider two formulations in simpegEM, 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: -Forward Problem -=============== +.. math :: + \nabla \times \vec{E} + i \omega \vec{B} = \vec{S}_m \\ + \nabla \times \mu^{-1} \vec{B} - \sigma \vec{E} = \vec{S}_e -Inverse Problem -=============== +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 + + +Discretizing +------------ +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) + +.. figure:: ./images/finitevolrealestate.png + :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 % + +Note that resistivity is the inverse of conductivity, \\(\\rho = \\sigma^{-1}\\). + + +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{s_e} + +H-J Formulation: +**************** + +.. math :: + \mathbf{C^T} \mathbf{M^f_\rho} \mathbf{j} + i \omega \mathbf{M^e_\mu} \mathbf{h} = \mathbf{s_m} \\ + \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} + + +.. Forward Problem +.. =============== + +.. Inverse Problem +.. =============== API === diff --git a/docs/images/ebjhdiscretizations.png b/docs/images/ebjhdiscretizations.png new file mode 100644 index 00000000..85bfbd36 Binary files /dev/null and b/docs/images/ebjhdiscretizations.png differ diff --git a/docs/images/finitevolrealestate.png b/docs/images/finitevolrealestate.png new file mode 100644 index 00000000..38a53c71 Binary files /dev/null and b/docs/images/finitevolrealestate.png differ diff --git a/simpegEM/Base.py b/simpegEM/Base.py index 690b2463..32018f7e 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -2,6 +2,10 @@ from SimPEG import Survey, Problem, Utils, Models, Maps, PropMaps, np, sp, Solve from scipy.constants import mu_0 class EMPropMap(Maps.PropMap): + """ + Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m) + """ + sigma = Maps.Property("Electrical Conductivity", defaultInvProp = True, propertyLink=('rho',Maps.ReciprocalMap)) mu = Maps.Property("Inverse Magnetic Permeability", defaultVal = mu_0, propertyLink=('mui',Maps.ReciprocalMap)) @@ -50,12 +54,18 @@ class BaseEMProblem(Problem.BaseProblem): @property def Me(self): + """ + Edge inner product matrix + """ if getattr(self, '_Me', None) is None: self._Me = self.mesh.getEdgeInnerProduct() return self._Me @property def Mf(self): + """ + Face inner product matrix + """ if getattr(self, '_Mf', None) is None: self._Mf = self.mesh.getFaceInnerProduct() return self._Mf @@ -64,32 +74,48 @@ class BaseEMProblem(Problem.BaseProblem): # ----- Magnetic Permeability ----- # @property def MfMui(self): + """ + Face inner product matrix for \\(\\mu^{-1}\\). Used in the E-B formulation + """ if getattr(self, '_MfMui', None) is None: self._MfMui = self.mesh.getFaceInnerProduct(self.curModel.mui) return self._MfMui @property def MfMuiI(self): + """ + Inverse of :code:`MfMui`. + """ if getattr(self, '_MfMuiI', None) is None: self._MfMuiI = self.mesh.getFaceInnerProduct(self.curModel.mui, invMat=True) return self._MfMuiI @property def MeMu(self): + """ + Edge inner product matrix for \\(\\mu\\). Used in the H-J formulation + """ if getattr(self, '_MeMu', None) is None: self._MeMu = self.mesh.getEdgeInnerProduct(self.curModel.mu) return self._MeMu @property def MeMuI(self): + """ + Inverse of :code:`MeMu` + """ if getattr(self, '_MeMuI', None) is None: self._MeMuI = self.mesh.getEdgeInnerProduct(self.curModel.mu, invMat=True) return self._MeMuI + # ----- Electrical Conductivity ----- # #TODO: hardcoded to sigma as the model @property def MeSigma(self): + """ + Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation + """ if getattr(self, '_MeSigma', None) is None: self._MeSigma = self.mesh.getEdgeInnerProduct(self.curModel.sigma) return self._MeSigma @@ -97,13 +123,16 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MeSigmaDeriv(self, u): """ - Deriv of MeSigma wrt sigma + Derivative of MeSigma with respect to the model """ return self.mesh.getEdgeInnerProductDeriv(self.curModel.sigma)(u) * self.curModel.sigmaDeriv @property def MeSigmaI(self): + """ + Inverse of the edge inner product matrix for \\(\\sigma\\). + """ if getattr(self, '_MeSigmaI', None) is None: self._MeSigmaI = self.mesh.getEdgeInnerProduct(self.curModel.sigma, invMat=True) return self._MeSigmaI @@ -111,7 +140,7 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This should take a vector def MeSigmaIDeriv(self, u): """ - Deriv of MeSigma wrt sigma + Derivative of :code:`MeSigma` with respect to the model """ # TODO: only works for diagonal tensors. getEdgeInnerProductDeriv, invMat=True should be implemented in SimPEG @@ -124,17 +153,26 @@ class BaseEMProblem(Problem.BaseProblem): @property def MfRho(self): + """ + Face inner product matrix for \\(\\rho\\). Used in the H-J formulation + """ if getattr(self, '_MfRho', None) is None: self._MfRho = self.mesh.getFaceInnerProduct(self.curModel.rho) return self._MfRho # TODO: This should take a vector def MfRhoDeriv(self,u): + """ + Derivative of :code:`MfRho` with respect to the model. + """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho)(u) * (-Utils.sdiag(self.curModel.rho**2) * self.curModel.sigmaDeriv) # self.curModel.rhoDeriv @property def MfRhoI(self): + """ + Inverse of :code:`MfRho` + """ if getattr(self, '_MfRhoI', None) is None: self._MfRhoI = self.mesh.getFaceInnerProduct(self.curModel.rho, invMat=True) return self._MfRhoI @@ -142,4 +180,7 @@ class BaseEMProblem(Problem.BaseProblem): # TODO: This isn't going to work yet # TODO: This should take a vector def MfRhoIDeriv(self,u): + """ + Derivative of :code:`MfRhoI` with respect to the model. + """ return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index c54137fd..5baca01e 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -8,18 +8,33 @@ from simpegEM.Utils.EMUtils import omega class BaseFDEMProblem(BaseEMProblem): """ - We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\): + We start by looking at Maxwell's equations in the electric field \\(\\mathbf{e}\\) and the magnetic flux density \\(\\mathbf{b}\\): - .. math:: + .. math :: - \\nabla \\times \\vec{E} + i \\omega \\vec{B} = \\vec{S_m} \\\\ - \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{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{s_e} + if using the E-B formulation (:code:`ProblemFDEM_e` or :code:`ProblemFDEM_b`) or the magnetic field \\(\\mathbf{h}\\) and current density \\(\\mathbf{j}\\) + + .. math :: + + \\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{j} + i \\omega \\mathbf{M^e_{\\mu}} \\mathbf{h} = \\mathbf{s_m} \\\\ + \\mathbf{C} \\mathbf{h} - \\mathbf{j} = \\mathbf{s_e} + + if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`). + + The problem performs the elimination so that we are solving the system for \\(\\mathbf{e},\\mathbf{b},\\mathbf{j} or \\mathbf{h}\\) """ + surveyPair = SurveyFDEM fieldsPair = FieldsFDEM def fields(self, m=None): + """ + Solve the forward problem for the fields. + """ + self.curModel = m F = self.fieldsPair(self.mesh, self.survey) @@ -35,6 +50,10 @@ class BaseFDEMProblem(BaseEMProblem): return F def Jvec(self, m, v, f=None): + """ + Sensitivity times a vector + """ + if f is None: f = self.fields(m) @@ -75,6 +94,10 @@ class BaseFDEMProblem(BaseEMProblem): return Utils.mkvc(Jv) def Jtvec(self, m, v, f=None): + """ + Sensitivity transpose times a vector + """ + if f is None: f = self.fields(m) @@ -130,9 +153,11 @@ class BaseFDEMProblem(BaseEMProblem): def getSourceTerm(self, freq): """ + 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: RHS + :return: S_m, S_e """ Srcs = self.survey.getSrcByFreq(freq) if self._eqLocs is 'FE': @@ -160,19 +185,19 @@ class ProblemFDEM_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using - .. math:: + .. math :: - \\vec{B} = \\frac{-1}{i\\omega}\\nabla\\times\\vec{E}, - - we can write Maxwell's equations as a second order system in \\ \\vec{E} \\ only: - - .. math:: - - \\nabla \\times \\mu^{-1} \\nabla \\times \\vec{E} + i \\omega \\sigma \\vec{E} = \\vec{J_s} - - This is the definition of the Forward Problem using the E-formulation of Maxwell's equations. + \\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: + + .. math :: + + \\left(\\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}} \\mathbf{C} + i \\omega \\mathbf{M^e_{\\sigma}} \\right)\\mathbf{e} + = \\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}}\\mathbf{s_m} -i\\omega\\mathbf{s_e} \\\\ + + which we solve for \\(\\mathbf{e}\\). """ _fieldType = 'e' @@ -184,6 +209,9 @@ class ProblemFDEM_e(BaseFDEMProblem): def getA(self, freq): """ + .. math :: + \\mathbf{A} = \\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}} \\mathbf{C} + i \\omega \\mathbf{M^e_{\\sigma}} + :param float freq: Frequency :rtype: scipy.sparse.csr_matrix :return: A @@ -206,6 +234,9 @@ class ProblemFDEM_e(BaseFDEMProblem): def getRHS(self, freq): """ + .. math :: + \\mathbf{RHS} = \\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}}\\mathbf{s_m} -i\\omega\\mathbf{s_e} + :param float freq: Frequency :rtype: numpy.ndarray (nE, nSrc) :return: RHS @@ -251,8 +282,20 @@ class ProblemFDEM_e(BaseFDEMProblem): class ProblemFDEM_b(BaseFDEMProblem): """ - Solving for b! + We eliminate \\(\\mathbf{e}\\) using + + .. math :: + \\mathbf{e} = \\mathbf{M^e_{\\sigma}}^{-1} \\left(\\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}} \\mathbf{b} - \\mathbf{s_e}\\right) + + and solve for \\(\\mathbf{b}\\) using: + + .. math :: + \\left(\\mathbf{C} \\mathbf{M^e_{\\sigma}}^{-1} \\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}} + i \\omega \\right)\\mathbf{b} = \\mathbf{s_m} + \\mathbf{M^e_{\\sigma}}^{-1}\\mathbf{s_e} + + .. note :: + The inverse problem will not work with full anisotropy """ + _fieldType = 'b' _eqLocs = 'FE' fieldsPair = FieldsFDEM_b @@ -262,10 +305,14 @@ class ProblemFDEM_b(BaseFDEMProblem): def getA(self, freq): """ + .. math :: + \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{\\sigma}}^{-1} \\mathbf{C}^T \\mathbf{M^f_{\\mu^{-1}} + i \\omega + :param float freq: Frequency :rtype: scipy.sparse.csr_matrix :return: A """ + MfMui = self.MfMui MeSigmaI = self.MeSigmaI C = self.mesh.edgeCurl @@ -298,6 +345,9 @@ class ProblemFDEM_b(BaseFDEMProblem): def getRHS(self, freq): """ + .. math :: + \\mathbf{RHS} = \\mathbf{s_m} + \\mathbf{M^e_{\\sigma}}^{-1}\\mathbf{s_e} + :param float freq: Frequency :rtype: numpy.ndarray (nE, nSrc) :return: RHS @@ -369,23 +419,16 @@ class ProblemFDEM_b(BaseFDEMProblem): class ProblemFDEM_j(BaseFDEMProblem): """ - Using the H-J formulation of Maxwell's equations + We eliminate \\(\\mathbf{h}\\) using - .. math:: - \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 - \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + .. math :: - Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + \\mathbf{h} = \\frac{1}{i \\omega} \\mathbf{M^e_{\\mu}}^{-1} \\left(- \\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{j} + \\mathbf{s_m} \\right) - For this implementation, we solve for J using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) : + and solve for \\(\\mathbf{j}\\) using - .. math:: - \\nabla \\times ( \\mu^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} ) + i\\omega \\vec{J} = - i\\omega\\vec{J_s} - - We discretize this to: - - .. math:: - (\\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega ) \\mathbf{j} = - i\\omega \\mathbf{j_s} + .. math :: + \\left(\\mathbf{C} \\mathbf{M^e_{\\mu}}^{-1} \\mathbf{C}^T \\mathbf{M^f_{\\rho}} + i \\omega\\right))\\mathbf{j} = \\mathbf{C} \\mathbf{M^e_{\\mu}}^{-1}\\mathbf{s_m} -i\\omega\\mathbf{s_e} .. note:: This implementation does not yet work with full anisotropy!! @@ -401,9 +444,8 @@ class ProblemFDEM_j(BaseFDEMProblem): def getA(self, freq): """ - Here, we form the operator \(\\mathbf{A}\) to solce - .. math:: - \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega + .. math :: + \\mathbf{A} = \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C}^T \\mathbf{M^f_{\\sigma^{-1}}} + i\\omega :param float freq: Frequency :rtype: scipy.sparse.csr_matrix @@ -425,7 +467,7 @@ class ProblemFDEM_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 - .. math:: + .. math :: \\frac{\mathbf{A(\\sigma)} \mathbf{v}}{d \\mathbf{m}} &= \\mathbf{C} \\mathbf{M^e_{mu^{-1}}} \\mathbf{C^T} \\frac{d \\mathbf{M^f_{\\sigma^{-1}}}}{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}} """ @@ -447,6 +489,8 @@ class ProblemFDEM_j(BaseFDEMProblem): def getRHS(self, freq): """ + .. math :: + \\mathbf{RHS} = \\mathbf{C} \\mathbf{M^e_{\\mu}}^{-1}\\mathbf{s_m} -i\\omega\\mathbf{s_e} :param float freq: Frequency :rtype: numpy.ndarray (nE, nSrc) :return: RHS @@ -505,26 +549,15 @@ class ProblemFDEM_j(BaseFDEMProblem): class ProblemFDEM_h(BaseFDEMProblem): """ - Using the H-J formulation of Maxwell's equations + We eliminate \\(\\mathbf{j}\\) using - .. math:: - \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 - \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + .. math :: + \\mathbf{j} = \\mathbf{C} \\mathbf{h} - \\mathbf{s_e} - Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + and solve for \\(\\mathbf{h}\\) using - For this implementation, we solve for J using \( \\vec{J} = \\nabla \\times \\vec{H} - \\vec{J_s} \) - - .. math:: - \\nabla \\times \\sigma^{-1} \\nabla \\times \\vec{H} + i\\omega\\mu\\vec{H} = \\nabla \\times \\sigma^{-1} \\vec{J_s} - - We discretize and solve - - .. math:: - (\\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\mathbf{C} + i\\omega \\mathbf{M_{\mu}} ) \\mathbf{h} = \\mathbf{C^T} \\mathbf{M^f_{\\sigma^{-1}}} \\vec{J_s} - - .. note:: - This implementation does not yet work with full anisotropy!! + .. math :: + \\left(\\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{C} + i \\omega \\mathbf{M^e_{\\mu}}\\right) \\mathbf{h} = \\mathbf{s_m} + \\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{s_e} """ @@ -537,6 +570,8 @@ class ProblemFDEM_h(BaseFDEMProblem): def getA(self, freq): """ + .. math :: + \\mathbf{A} = \\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{C} + i \\omega \\mathbf{M^e_{\\mu}} :param float freq: Frequency :rtype: scipy.sparse.csr_matrix :return: A @@ -560,6 +595,9 @@ class ProblemFDEM_h(BaseFDEMProblem): def getRHS(self, freq): """ + .. math :: + \\mathbf{RHS} = \\mathbf{s_m} + \\mathbf{C}^T \\mathbf{M^f_{\\rho}} \\mathbf{s_e} + :param float freq: Frequency :rtype: numpy.ndarray (nE, nSrc) :return: RHS