diff --git a/simpegEM/Base.py b/simpegEM/Base.py index c22252fa..d1e6bb08 100644 --- a/simpegEM/Base.py +++ b/simpegEM/Base.py @@ -17,6 +17,14 @@ class BaseEMProblem(Problem.BaseProblem): verbose = False + #################################################### + # Make A Symmetric + #################################################### + @property + def _makeASymmetric(self): + if getattr(self, '__makeASymmetric', None) is None: + self.__makeASymmetric = True + return self.__makeASymmetric #################################################### # Mu Model @@ -30,6 +38,10 @@ class BaseEMProblem(Problem.BaseProblem): def mu(self, value): if getattr(self, '_MfMui', None) is not None: del self._MfMui + if getattr(self, '_MeMu', None) is not None: + del delf._MeMu + if getattr(self, '_MeMuI', None) is not None: + del self._MeMuI self._mu = value @@ -44,6 +56,20 @@ class BaseEMProblem(Problem.BaseProblem): self._MfMui = self.mesh.getFaceInnerProduct(1/self.mu) return self._MfMui + @property + def MeMuI(self): + # TODO: Assuming isotropic mu + if getattr(self, '_MeMuI', None) is None: + self._MeMuI = self.mesh.getEdgeInnerProduct(self.mu, invMat=True) + return self._MeMuI + + @property + def MeMu(self): + #TODO: Assuming isotropic mu + if getattr(self, '_MeMu', None) is None: + self._MeMu = self.mesh.getEdgeInnerProduct(self.mu) + return self._MeMu + @property def Me(self): if getattr(self, '_Me', None) is None: @@ -66,7 +92,16 @@ class BaseEMProblem(Problem.BaseProblem): self._MeSigmaI = self.mesh.getEdgeInnerProduct(sigma, invMat=True) return self._MeSigmaI - deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI'] + @property + def MfSigmai(self): + #TODO: hardcoded to sigma as the model + #TODO: hardcoded to sigma diagonal + if getattr(self, '_MfSigmai', None) is None: + sigma = self.curModel.transform + self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma) + return self._MfSigmai + + deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai'] def fields(self, m): self.curModel = m diff --git a/simpegEM/FDEM/FDEM.py b/simpegEM/FDEM/FDEM.py index c4137486..1a45d139 100644 --- a/simpegEM/FDEM/FDEM.py +++ b/simpegEM/FDEM/FDEM.py @@ -8,6 +8,102 @@ def omega(freq): """Change frequency to angular frequency, omega""" return 2.*np.pi*freq +def getSource(self,freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + Txs = self.survey.getTransmitters(freq) + rhs = range(len(Txs)) + + solType = self.solType + + if solType == 'e' or solType == 'b': + gridEJx = self.mesh.gridEx + gridEJy = self.mesh.gridEy + gridEJz = self.mesh.gridEz + nEJ = self.mesh.nE + + gridBHx = self.mesh.gridFx + gridBHy = self.mesh.gridFy + gridBHz = self.mesh.gridFz + nBH = self.mesh.nF + + + C = self.mesh.edgeCurl + mui = self.MfMui + + elif solType == 'h' or solType == 'j': + gridEJx = self.mesh.gridFx + gridEJy = self.mesh.gridFy + gridEJz = self.mesh.gridFz + nEJ = self.mesh.nF + + gridBHx = self.mesh.gridEx + gridBHy = self.mesh.gridEy + gridBHz = self.mesh.gridEz + nBH = self.mesh.nE + + C = self.mesh.edgeCurl.T + mui = self.MeMuI + + else: + NotImplementedError('Only E or F sources') + + for i, tx in enumerate(Txs): + if self.mesh._meshType is 'CYL': + if self.mesh.isSymmetric: + if tx.txType == 'VMD': + SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y') + elif tx.txType =='CircularLoop': + SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius) + else: + raise NotImplementedError('Only VMD and CircularLoop') + else: + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + + elif self.mesh._meshType is 'TENSOR': + + if tx.txType == 'VMD': + src = Sources.MagneticDipoleVectorPotential + SRCx = src(tx.loc, gridEJx, 'x') + SRCy = src(tx.loc, gridEJy, 'y') + SRCz = src(tx.loc, gridEJz, 'z') + + elif tx.txType == 'VMD_B': + src = Sources.MagneticDipoleFields + SRCx = src(tx.loc, gridBHx, 'x') + SRCy = src(tx.loc, gridBHy, 'y') + SRCz = src(tx.loc, gridBHz, 'z') + + elif tx.txType == 'CircularLoop': + src = Sources.MagneticLoopVectorPotential + SRCx = src(tx.loc, gridEJx, 'x', tx.radius) + SRCy = src(tx.loc, gridEJy, 'y', tx.radius) + SRCz = src(tx.loc, gridEJz, 'z', tx.radius) + else: + + raise NotImplemented('%s txType is not implemented' % tx.txType) + SRC = np.concatenate((SRCx, SRCy, SRCz)) + + else: + raise Exception('Unknown mesh for VMD') + + rhs[i] = SRC + + # b-forumlation + if tx.txType == 'VMD_B': + b_0 = np.concatenate(rhs).reshape((nBH, len(Txs)), order='E') + else: + a = np.concatenate(rhs).reshape((nEJ, len(Txs)), order='F') + b_0 = C*a + + if solType == 'b' or solType == 'h': + return b_0 + elif solType == 'e' or solType == 'j': + return C.T*mui*b_0 + class BaseFDEMProblem(BaseEMProblem): """ We start by looking at Maxwell's equations in the electric field \\(\\vec{E}\\) and the magnetic flux density \\(\\vec{B}\\): @@ -18,7 +114,6 @@ class BaseFDEMProblem(BaseEMProblem): \\nabla \\times \\mu^{-1} \\vec{B} - \\sigma \\vec{E} = \\vec{J_s} """ - surveyPair = SurveyFDEM def forward(self, m, RHS, CalcFields): @@ -105,6 +200,12 @@ class BaseFDEMProblem(BaseEMProblem): return Jtv + def getSource(self,freq): + return self.getSource(freq) + +########################################################################################## +################################ E-B Formulation ######################################### +########################################################################################## class ProblemFDEM_e(BaseFDEMProblem): """ @@ -141,6 +242,7 @@ class ProblemFDEM_e(BaseFDEMProblem): return C.T*mui*C + 1j*omega(freq)*sig + def getADeriv(self, freq, u, v, adjoint=False): sig = self.curModel.transform dsig_dm = self.curModel.transformDeriv @@ -157,29 +259,8 @@ class ProblemFDEM_e(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x') - SRCy = src(tx.loc, self.mesh.gridEy, 'y') - SRCz = src(tx.loc, self.mesh.gridEz, 'z') - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - else: - raise NotImplemented('%s txType is not implemented' % tx.txType) - rhs[i] = np.concatenate((SRCx, SRCy, SRCz)) - - a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - mui = self.MfMui - C = self.mesh.edgeCurl - - j_s = C.T*mui*C*a + j_s = getSource(self,freq) return -1j*omega(freq)*j_s def calcFields(self, sol, freq, fieldType, adjoint=False): @@ -221,8 +302,13 @@ class ProblemFDEM_b(BaseFDEMProblem): mui = self.MfMui sigI = self.MeSigmaI C = self.mesh.edgeCurl + iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) - return mui*C*sigI*C.T*mui + 1j*omega(freq)*mui + A = C*sigI*C.T*mui + iomega + + if self._makeASymmetric is True: + return mui.T*A + return A def getADeriv(self, freq, u, v, adjoint=False): @@ -237,9 +323,14 @@ class ProblemFDEM_b(BaseFDEMProblem): dMe_dsig = self.mesh.getEdgeInnerProductDeriv(sig)(vec) if adjoint: - return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * ( mui.T * v ) ) ) ) + if self._makeASymmetric is True: + v = mui * v + return dsig_dm.T * ( dMe_dsig.T * ( dMeSigmaI_dI.T * ( C.T * v ) ) ) + + if self._makeASymmetric is True: + return mui.T * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) + return C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) - return mui * ( C * ( dMeSigmaI_dI * ( dMe_dsig * ( dsig_dm * v ) ) ) ) def getRHS(self, freq): """ @@ -247,59 +338,14 @@ class ProblemFDEM_b(BaseFDEMProblem): :rtype: numpy.ndarray (nE, nTx) :return: RHS """ - Txs = self.survey.getTransmitters(freq) - rhs = range(len(Txs)) - for i, tx in enumerate(Txs): - - if self.mesh._meshType is 'CYL': - if self.mesh.isSymmetric: - if tx.txType == 'VMD': - SRC = Sources.MagneticDipoleVectorPotential(tx.loc, self.mesh.gridEy, 'y') - elif tx.txType =='CircularLoop': - SRC = Sources.MagneticLoopVectorPotential(tx.loc, self.mesh.gridEy, 'y', tx.radius) - else: - raise NotImplementedError('Only VMD and CircularLoop') - else: - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - elif self.mesh._meshType is 'TENSOR': + b_0 = getSource(self,freq) - if tx.txType == 'VMD': - src = Sources.MagneticDipoleVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x') - SRCy = src(tx.loc, self.mesh.gridEy, 'y') - SRCz = src(tx.loc, self.mesh.gridEz, 'z') - - elif tx.txType == 'VMD_B': - src = Sources.MagneticDipoleFields - SRCx = src(tx.loc, self.mesh.gridFx, 'x') - SRCy = src(tx.loc, self.mesh.gridFy, 'y') - SRCz = src(tx.loc, self.mesh.gridFz, 'z') - - elif tx.txType == 'CircularLoop': - src = Sources.MagneticLoopVectorPotential - SRCx = src(tx.loc, self.mesh.gridEx, 'x', tx.radius) - SRCy = src(tx.loc, self.mesh.gridEy, 'y', tx.radius) - SRCz = src(tx.loc, self.mesh.gridEz, 'z', tx.radius) - else: - - raise NotImplemented('%s txType is not implemented' % tx.txType) - SRC = np.concatenate((SRCx, SRCy, SRCz)) - - else: - raise Exception('Unknown mesh for VMD') - - rhs[i] = SRC - - mui = self.MfMui - if tx.txType == 'VMD_B': - b_0 = np.concatenate(rhs).reshape((self.mesh.nF, len(Txs)), order='F') - else: - a = np.concatenate(rhs).reshape((self.mesh.nE, len(Txs)), order='F') - C = self.mesh.edgeCurl - b_0 = C*a - - return -1j*omega(freq)*mui*b_0 + rhs = -1j*omega(freq)*b_0 + if self._makeASymmetric is True: + mui = self.MfMui + return mui.T*rhs + return rhs def calcFields(self, sol, freq, fieldType, adjoint=False): b = sol @@ -335,3 +381,225 @@ class ProblemFDEM_b(BaseFDEMProblem): return None raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + +########################################################################################## +################################ H-J Formulation ######################################### +########################################################################################## + + +class ProblemFDEM_j(BaseFDEMProblem): + """ + Using the H-J formulation of Maxwell's equations + + .. math:: + \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 + \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + + Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + + For this implementation, we solve for J using \( \\vec{H} = - (i\\omega\\mu)^{-1} \\nabla \\times \\sigma^{-1} \\vec{J} \) : + + .. 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} + + .. note:: + This implementation does not yet work with full anisotropy!! + + """ + + solType = 'j' + storeTheseFields = ['j','h'] + + def __init__(self, model, **kwargs): + BaseFDEMProblem.__init__(self, model, **kwargs) + + 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 + + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + MeMuI = self.MeMuI + MfSigi = self.MfSigmai + C = self.mesh.edgeCurl + iomega = 1j * omega(freq) * sp.eye(self.mesh.nF) + + A = C * MeMuI * C.T * MfSigi + iomega + + if self._makeASymmetric is True: + return MfSigi.T*A + return A + + + def getADeriv(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:: + \\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}} + """ + + MeMuI = self.MeMuI + MfSigi = self.MfSigmai + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(u) + + if adjoint: + if self._makeASymmetric is True: + v = MfSigi * v + return dsig_dm.T * ( dsigi_dsig.T *( dMf_dsigi.T * ( C * ( MeMuI.T * ( C.T * v ) ) ) ) ) + + if self._makeASymmetric is True: + return MfSigi.T * ( C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) ) + return C * ( MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) ) + + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + j_s = getSource(self,freq) + rhs = -1j*omega(freq)*j_s + + if self._makeASymmetric is True: + MfSigi = self.MfSigmai + return MfSigi.T*rhs + return rhs + + def calcFields(self, sol, freq, fieldType, adjoint=False): + j = sol + if fieldType == 'j': + return j + elif fieldType == 'h': + MeMuI = self.MeMuI + C = self.mesh.edgeCurl + MfSigi = self.MfSigmai + if not adjoint: + h = -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( MfSigi * j ) ) + else: + h = -(1./(1j*omega(freq))) * MfSigi.T * ( C * ( MeMuI.T * j ) ) + return h + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + j = sol + if fieldType == 'j': + return None + elif fieldType == 'h': + MeMuI = self.MeMuI + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(j) + sigi = self.MfSigmai + if not adjoint: + return -(1./(1j*omega(freq))) * MeMuI * ( C.T * ( dMf_dsigi * ( dsigi_dsig * ( dsig_dm * v ) ) ) ) + else: + return -(1./(1j*omega(freq))) * dsig_dm.T * ( dsigi_dsig.T * ( dMf_dsigi.T * ( C * ( MeMuI.T * v ) ) ) ) + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + + + +# Solving for h! - using primary- secondary approach +class ProblemFDEM_h(BaseFDEMProblem): + """ + Using the H-J formulation of Maxwell's equations + + .. math:: + \\nabla \\times \\sigma^{-1} \\vec{J} + i\\omega\\mu\\vec{H} = 0 + \\nabla \\times \\vec{H} - \\vec{J} = \\vec{J_s} + + Since \(\\vec{J}\) is a flux and \(\\vec{H}\) is a field, we discretize \(\\vec{J}\) on faces and \(\\vec{H}\) on edges. + + 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!! + + """ + + solType = 'h' + storeTheseFields = ['j','h'] + + def __init__(self, model, **kwargs): + BaseFDEMProblem.__init__(self, model, **kwargs) + + def getA(self, freq): + """ + :param float freq: Frequency + :rtype: scipy.sparse.csr_matrix + :return: A + """ + + MeMu = self.MeMu + MfSigi = self.MfSigmai + C = self.mesh.edgeCurl + + return C.T * MfSigi * C + 1j*omega(freq)*MeMu + + def getADeriv(self, freq, u, v, adjoint=False): + + MeMu = self.MeMu + C = self.mesh.edgeCurl + sig = self.curModel.transform + sigi = 1/sig + dsig_dm = self.curModel.transformDeriv + dsigi_dsig = -Utils.sdiag(sigi)**2 + + dMf_dsigi = self.mesh.getFaceInnerProductDeriv(sigi)(C*u) + + if adjoint: + return (dsig_dm.T * (dsigi_dsig.T * (dMf_dsigi.T * (C * v)))) + return (C.T * (dMf_dsigi * (dsigi_dsig * (dsig_dm * v)))) + + + + def getRHS(self, freq): + """ + :param float freq: Frequency + :rtype: numpy.ndarray (nE, nTx) + :return: RHS + """ + b_0 = getSource(self,freq) + return -1j*omega(freq)*b_0 + + def calcFields(self, sol, freq, fieldType, adjoint=False): + h = sol + if fieldType == 'j': + C = self.mesh.edgeCurl + if adjoint: + return C.T*h + return C*h + elif fieldType == 'h': + return h + raise NotImplementedError('fieldType "%s" is not implemented.' % fieldType) + + def calcFieldsDeriv(self, sol, freq, fieldType, v, adjoint=False): + return None \ No newline at end of file diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 9e87465f..e9000b7e 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -16,6 +16,20 @@ class RxFDEM(Survey.BaseRx): 'bxi':['b', 'Fx', 'imag'], 'byi':['b', 'Fy', 'imag'], 'bzi':['b', 'Fz', 'imag'], + + 'jxr':['j', 'Fx', 'real'], + 'jyr':['j', 'Fy', 'real'], + 'jzr':['j', 'Fz', 'real'], + 'jxi':['j', 'Fx', 'imag'], + 'jyi':['j', 'Fy', 'imag'], + 'jzi':['j', 'Fz', 'imag'], + + 'hxr':['h', 'Ex', 'real'], + 'hyr':['h', 'Ey', 'real'], + 'hzr':['h', 'Ez', 'real'], + 'hxi':['h', 'Ex', 'imag'], + 'hyi':['h', 'Ey', 'imag'], + 'hzi':['h', 'Ez', 'imag'], } radius = None @@ -84,7 +98,7 @@ class TxFDEM(Survey.BaseTx): class FieldsFDEM(Problem.Fields): """Fancy Field Storage for a FDEM survey.""" - knownFields = {'b': 'F', 'e': 'E'} + knownFields = {'b': 'F', 'e': 'E', 'j': 'F', 'h': 'E'} dtype = complex diff --git a/simpegEM/FDEM/__init__.py b/simpegEM/FDEM/__init__.py index 70124686..562b9218 100644 --- a/simpegEM/FDEM/__init__.py +++ b/simpegEM/FDEM/__init__.py @@ -1,2 +1,2 @@ from SurveyFDEM import * -from FDEM import ProblemFDEM_e, ProblemFDEM_b +from FDEM import ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index e1d8f1c5..3a1c4c92 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -5,9 +5,11 @@ import sys from scipy.constants import mu_0 TOL = 1e-4 -CONDUCTIVITY = 1e3 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 MU = mu_0 -addrandoms = True +freq = 1e-1 +addrandoms = True def getProblem(fdemType, comp): cs = 5. @@ -20,17 +22,22 @@ def getProblem(fdemType, comp): mapping = Maps.ExpMap(mesh) - x = np.linspace(-30,30,6) - XYZ = Utils.ndgrid(x,x,np.r_[0]) + x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter + XYZ = Utils.ndgrid(x,x,np.r_[0.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - Tx0 = EM.FDEM.TxFDEM(np.r_[4.,2.,2.], 'VMD', 1e-2, [Rx0]) + Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0]) survey = EM.FDEM.SurveyFDEM([Tx0]) + if fdemType == 'e': prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) elif fdemType == 'b': prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + elif fdemType == 'j': + prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + elif fdemType == 'h': + prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) else: raise NotImplementedError() prb.pair(survey) @@ -51,19 +58,20 @@ def adjointTest(fdemType, comp): mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3 + m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 prb.mu = mu survey = prb.survey + u = prb.fields(m) + v = np.random.rand(survey.nD) w = np.random.rand(prb.mapping.nP) - u = prb.fields(m) vJw = v.dot(prb.Jvec(m, w, u=u)) wJtv = w.dot(prb.Jtvec(m, v, u=u)) - tol = TOL*(10**int(np.log10(np.abs(vJw)))) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol return np.abs(vJw - wJtv) < tol @@ -74,14 +82,64 @@ def derivTest(fdemType, comp): mu = np.log(np.ones(prb.mesh.nC)*MU) if addrandoms is True: - x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-3 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-3 + x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - prb.mu = mu + prb.mu = mu survey = prb.survey def fun(x): return survey.dpred(x), lambda x: prb.Jvec(x0, x) - return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=1e-25) + return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) + + +def crossCheckTest(fdemType, comp): + + l2norm = lambda r: np.sqrt(r.dot(r)) + + prb1 = getProblem(fdemType, comp) + mesh = prb1.mesh + print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp) + m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) + mu = np.log(np.ones(mesh.nC)*MU) + + if addrandoms is True: + m = m + np.random.randn(mesh.nC)*CONDUCTIVITY*1e-1 + mu = mu + np.random.randn(mesh.nC)*MU*1e-1 + + prb1.mu = mu + survey1 = prb1.survey + + u1 = prb1.fields(m) + d1 = Utils.mkvc(survey1.projectFields(u1)) + + prb1.unpair + + if fdemType == 'e': + prb2 = getProblem('b', comp) + elif fdemType == 'b': + prb2 = getProblem('e', comp) + elif fdemType == 'j': + prb2 = getProblem('h', comp) + elif fdemType == 'h': + prb2 = getProblem('j', comp) + else: + raise NotImplementedError() + + prb2.mu = mu + survey2 = prb2.survey + + u2 = prb2.fields(m) + d2 = Utils.mkvc(survey2.projectFields(u2)) + + r = d2-d1 + l2r = l2norm(r) + + tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) + print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol + return l2r < tol + + + class FDEM_DerivTests(unittest.TestCase): @@ -189,6 +247,162 @@ class FDEM_DerivTests(unittest.TestCase): self.assertTrue(adjointTest('b', 'bzi')) + def test_Jvec_jxr_Jform(self): + self.assertTrue(derivTest('j', 'jxr')) + def test_Jvec_jyr_Jform(self): + self.assertTrue(derivTest('j', 'jyr')) + def test_Jvec_jzr_Jform(self): + self.assertTrue(derivTest('j', 'jzr')) + def test_Jvec_jxi_Jform(self): + self.assertTrue(derivTest('j', 'jxi')) + def test_Jvec_jyi_Jform(self): + self.assertTrue(derivTest('j', 'jyi')) + def test_Jvec_jzi_Jform(self): + self.assertTrue(derivTest('j', 'jzi')) + + def test_Jvec_hxr_Jform(self): + self.assertTrue(derivTest('j', 'hxr')) + def test_Jvec_hyr_Jform(self): + self.assertTrue(derivTest('j', 'hyr')) + def test_Jvec_hzr_Jform(self): + self.assertTrue(derivTest('j', 'hzr')) + def test_Jvec_hxi_Jform(self): + self.assertTrue(derivTest('j', 'hxi')) + def test_Jvec_hyi_Jform(self): + self.assertTrue(derivTest('j', 'hyi')) + def test_Jvec_hzi_Jform(self): + self.assertTrue(derivTest('j', 'hzi')) + + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'hxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'hyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'hzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'hxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'hyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'hzi')) + + def test_Jvec_hxr_Hform(self): + self.assertTrue(derivTest('h', 'jxr')) + def test_Jvec_hyr_Hform(self): + self.assertTrue(derivTest('h', 'jyr')) + def test_Jvec_hzr_Hform(self): + self.assertTrue(derivTest('h', 'jzr')) + def test_Jvec_hxi_Hform(self): + self.assertTrue(derivTest('h', 'jxi')) + def test_Jvec_hyi_Hform(self): + self.assertTrue(derivTest('h', 'jyi')) + def test_Jvec_hzi_Hform(self): + self.assertTrue(derivTest('h', 'jzi')) + + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) + + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + self.assertTrue(adjointTest('j', 'hzi')) + + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) + + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'jxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'jyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'jzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'jxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'jyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'jzi')) + + + def test_EB_CrossCheck_exr_Eform(self): + self.assertTrue(crossCheckTest('e', 'exr')) + def test_EB_CrossCheck_eyr_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyr')) + def test_EB_CrossCheck_ezr_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezr')) + def test_EB_CrossCheck_exi_Eform(self): + self.assertTrue(crossCheckTest('e', 'exi')) + def test_EB_CrossCheck_eyi_Eform(self): + self.assertTrue(crossCheckTest('e', 'eyi')) + def test_EB_CrossCheck_ezi_Eform(self): + self.assertTrue(crossCheckTest('e', 'ezi')) + + def test_EB_CrossCheck_bxr_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxr')) + def test_EB_CrossCheck_byr_Eform(self): + self.assertTrue(crossCheckTest('e', 'byr')) + # def test_EB_CrossCheck_bzr_Eform(self): + # self.assertTrue(crossCheckTest('e', 'bzr')) # Doesn't make sense to test this for p-s approach + def test_EB_CrossCheck_bxi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bxi')) + def test_EB_CrossCheck_byi_Eform(self): + self.assertTrue(crossCheckTest('e', 'byi')) + def test_EB_CrossCheck_bzi_Eform(self): + self.assertTrue(crossCheckTest('e', 'bzi')) + + def test_HJ_CrossCheck_jxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxr')) + def test_HJ_CrossCheck_jyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyr')) + def test_HJ_CrossCheck_jzr_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzr')) + def test_HJ_CrossCheck_jxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jxi')) + def test_HJ_CrossCheck_jyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jyi')) + def test_HJ_CrossCheck_jzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'jzi')) + + def test_HJ_CrossCheck_hxr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxr')) + def test_HJ_CrossCheck_hyr_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyr')) + # def test_HJ_CrossCheck_hzr_Jform(self): + # self.assertTrue(crossCheckTest('j', 'hzr')) # Doesn't make sense to test this for p-s approach + def test_HJ_CrossCheck_hxi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hxi')) + def test_HJ_CrossCheck_hyi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hyi')) + def test_HJ_CrossCheck_hzi_Jform(self): + self.assertTrue(crossCheckTest('j', 'hzi')) if __name__ == '__main__': unittest.main()