diff --git a/.travis.yml b/.travis.yml index 39a3925b..c2c672bf 100644 --- a/.travis.yml +++ b/.travis.yml @@ -5,7 +5,11 @@ python: sudo: false env: - - TEST_DIR=tests/em + - TEST_DIR=tests/em/examples + - TEST_DIR=tests/em/fdem/forward + - TEST_DIR=tests/em/fdem/inverse/derivs + - TEST_DIR=tests/em/fdem/inverse/adjoint + - TEST_DIR=tests/em/tdem - TEST_DIR=tests/mesh - TEST_DIR=tests/flow - TEST_DIR=tests/utils diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 6d3df44f..cce32964 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -1,9 +1,9 @@ -from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver +from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 -from SurveyFDEM import SurveyFDEM -from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j +from SurveyFDEM import Survey as SurveyFDEM +from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j from SimPEG.EM.Base import BaseEMProblem -from SimPEG.EM.Utils.EMUtils import omega +from SimPEG.EM.Utils import omega class BaseFDEMProblem(BaseEMProblem): @@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem): \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}} - if using the E-B formulation (:code:`ProblemFDEM_e` - or :code:`ProblemFDEM_b`) or the magnetic field + if using the E-B formulation (:code:`Problem_e` + or :code:`Problem_b`) or the magnetic field \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: @@ -26,13 +26,13 @@ class BaseFDEMProblem(BaseEMProblem): \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} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`). + if using the H-J formulation (:code:`Problem_j` or :code:`Problem_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 + fieldsPair = Fields def fields(self, m=None): """ @@ -185,7 +185,7 @@ class BaseFDEMProblem(BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## -class ProblemFDEM_e(BaseFDEMProblem): +class Problem_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -205,7 +205,7 @@ class ProblemFDEM_e(BaseFDEMProblem): _fieldType = 'e' _eqLocs = 'FE' - fieldsPair = FieldsFDEM_e + fieldsPair = Fields_e def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -284,7 +284,7 @@ class ProblemFDEM_e(BaseFDEMProblem): return None -class ProblemFDEM_b(BaseFDEMProblem): +class Problem_b(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{e}\\\) using @@ -304,7 +304,7 @@ class ProblemFDEM_b(BaseFDEMProblem): _fieldType = 'b' _eqLocs = 'FE' - fieldsPair = FieldsFDEM_b + fieldsPair = Fields_b def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -425,7 +425,7 @@ class ProblemFDEM_b(BaseFDEMProblem): ########################################################################################## -class ProblemFDEM_j(BaseFDEMProblem): +class Problem_j(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{h}\\\) using @@ -446,7 +446,7 @@ class ProblemFDEM_j(BaseFDEMProblem): _fieldType = 'j' _eqLocs = 'EF' - fieldsPair = FieldsFDEM_j + fieldsPair = Fields_j def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -558,7 +558,7 @@ class ProblemFDEM_j(BaseFDEMProblem): -class ProblemFDEM_h(BaseFDEMProblem): +class Problem_h(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{j}\\\) using @@ -576,7 +576,7 @@ class ProblemFDEM_h(BaseFDEMProblem): _fieldType = 'h' _eqLocs = 'EF' - fieldsPair = FieldsFDEM_h + fieldsPair = Fields_h def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index bb786bd1..d83877c3 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -1,13 +1,16 @@ -from SimPEG import Survey, Problem, Utils, np, sp -from SimPEG.EM.Utils.EMUtils import omega +import numpy as np +import scipy.sparse as sp +import SimPEG +from SimPEG import Utils +from SimPEG.EM.Utils import omega -class FieldsFDEM(Problem.Fields): +class Fields(SimPEG.Problem.Fields): """Fancy Field Storage for a FDEM survey.""" knownFields = {} dtype = complex -class FieldsFDEM_e(FieldsFDEM): +class Fields_e(Fields): knownFields = {'eSolution':'E'} aliasFields = { 'e' : ['eSolution','E','_e'], @@ -19,7 +22,7 @@ class FieldsFDEM_e(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -89,7 +92,7 @@ class FieldsFDEM_e(FieldsFDEM): return self._bSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_b(FieldsFDEM): +class Fields_b(Fields): knownFields = {'bSolution':'F'} aliasFields = { 'b' : ['bSolution','F','_b'], @@ -101,7 +104,7 @@ class FieldsFDEM_b(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -190,7 +193,7 @@ class FieldsFDEM_b(FieldsFDEM): return self._eSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_j(FieldsFDEM): +class Fields_j(Fields): knownFields = {'jSolution':'F'} aliasFields = { 'j' : ['jSolution','F','_j'], @@ -202,7 +205,7 @@ class FieldsFDEM_j(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -293,7 +296,7 @@ class FieldsFDEM_j(FieldsFDEM): return self._hSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_h(FieldsFDEM): +class Fields_h(Fields): knownFields = {'hSolution':'E'} aliasFields = { 'h' : ['hSolution','E','_h'], @@ -305,7 +308,7 @@ class FieldsFDEM_h(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py new file mode 100644 index 00000000..5dcf9407 --- /dev/null +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -0,0 +1,347 @@ +from SimPEG import Survey, Problem, Utils, np, sp +from scipy.constants import mu_0 +from SimPEG.EM.Utils import * +# from SurveyFDEM import Rx + + +class BaseSrc(Survey.BaseSrc): + freq = None + # rxPair = Rx + integrate = True + + def eval(self, prob): + 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 bPrimary(self, prob): + return None + + def hPrimary(self, prob): + return None + + def ePrimary(self, prob): + return None + + def jPrimary(self, prob): + return None + + def S_m(self, prob): + return None + + def S_e(self, prob): + return None + + def S_mDeriv(self, prob, v, adjoint = False): + return None + + def S_eDeriv(self, prob, v, adjoint = False): + return None + + +class RawVec_e(BaseSrc): + """ + 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 + """ + + def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_e = np.array(S_e,dtype=complex) + self._ePrimary = ePrimary + self._bPrimary = bPrimary + self._hPrimary = hPrimary + self._jPrimary = jPrimary + self.freq = float(freq) + BaseSrc.__init__(self, rxList) + + def S_e(self, prob): + return self._S_e + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec_m(BaseSrc): + """ + 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 + """ + + def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_m = np.array(S_m,dtype=complex) + self.freq = float(freq) + self.integrate = integrate + self._ePrimary = np.array(ePrimary,dtype=complex) + self._bPrimary = np.array(bPrimary,dtype=complex) + self._hPrimary = np.array(hPrimary,dtype=complex) + self._jPrimary = np.array(jPrimary,dtype=complex) + + BaseSrc.__init__(self, rxList) + + def S_m(self, prob): + return self._S_m + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec(BaseSrc): + """ + 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 + """ + def __init__(self, rxList, freq, S_m, S_e, integrate = True): + self._S_m = np.array(S_m,dtype=complex) + self._S_e = np.array(S_e,dtype=complex) + self.freq = float(freq) + self.integrate = integrate + BaseSrc.__init__(self, rxList) + + def S_m(self, prob): + if prob._eqLocs is 'EF' and self.integrate is True: + return prob.Me * self._S_m + return self._S_m + + def S_e(self, prob): + if prob._eqLocs is 'FE' and self.integrate is True: + return prob.Me * self._S_e + return self._S_e + + +class MagDipole(BaseSrc): + + #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): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + self.integrate = False + BaseSrc.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob,b) + + def S_m(self, prob): + b_p = self.bPrimary(prob) + return -1j*omega(self.freq)*b_p + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class MagDipole_Bfield(BaseSrc): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + #TODO: neither does moment + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + BaseSrc.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl.T + + srcfct = MagneticDipoleFields + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,bz)) + else: + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,by,bz)) + + return b + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob, b) + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class CircularLoop(BaseSrc): + + #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): + self.freq = float(freq) + self.orientation = orientation + self.radius = radius + self.mu = mu + self.loc = loc + self.integrate = False + BaseSrc.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) + ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) + az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return 1./self.mu*b + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index dbda9f80..150a6c00 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,13 +1,13 @@ -from SimPEG import Survey, Problem, Utils, np, sp -from SimPEG.EM.Utils import SrcUtils -from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +import SimPEG +from SimPEG.EM.Utils import * from scipy.constants import mu_0 +import SrcFDEM as Src #################################################### # Receivers #################################################### -class RxFDEM(Survey.BaseRx): +class Rx(SimPEG.Survey.BaseRx): knownRxTypes = { 'exr':['e', 'Ex', 'real'], @@ -41,7 +41,7 @@ class RxFDEM(Survey.BaseRx): radius = None def __init__(self, locs, rxType): - Survey.BaseRx.__init__(self, locs, rxType) + SimPEG.Survey.BaseRx.__init__(self, locs, rxType) @property def projField(self): @@ -87,366 +87,21 @@ class RxFDEM(Survey.BaseRx): return Pv -#################################################### -# Sources -#################################################### - -class SrcFDEM(Survey.BaseSrc): - freq = None - rxPair = RxFDEM - integrate = True - - def eval(self, prob): - 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 bPrimary(self, prob): - return None - - def hPrimary(self, prob): - return None - - def ePrimary(self, prob): - return None - - def jPrimary(self, prob): - return None - - def S_m(self, prob): - return None - - def S_e(self, prob): - return None - - def S_mDeriv(self, prob, v, adjoint = False): - return None - - def S_eDeriv(self, prob, v, adjoint = False): - return None - - -class SrcFDEM_RawVec_e(SrcFDEM): - """ - 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 - """ - - def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_e = np.array(S_e,dtype=complex) - self._ePrimary = ePrimary - self._bPrimary = bPrimary - self._hPrimary = hPrimary - self._jPrimary = jPrimary - self.freq = float(freq) - SrcFDEM.__init__(self, rxList) - - def S_e(self, prob): - return self._S_e - - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary - - def jPrimary(self, prob): - return self._jPrimary - - -class SrcFDEM_RawVec_m(SrcFDEM): - """ - 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 - """ - - def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_m = np.array(S_m,dtype=complex) - self.freq = float(freq) - self.integrate = integrate - self._ePrimary = np.array(ePrimary,dtype=complex) - self._bPrimary = np.array(bPrimary,dtype=complex) - self._hPrimary = np.array(hPrimary,dtype=complex) - self._jPrimary = np.array(jPrimary,dtype=complex) - - SrcFDEM.__init__(self, rxList) - - def S_m(self, prob): - return self._S_m - - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary - - def jPrimary(self, prob): - return self._jPrimary - - -class SrcFDEM_RawVec(SrcFDEM): - """ - 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 - """ - def __init__(self, rxList, freq, S_m, S_e, integrate = True): - self._S_m = np.array(S_m,dtype=complex) - self._S_e = np.array(S_e,dtype=complex) - self.freq = float(freq) - self.integrate = integrate - SrcFDEM.__init__(self, rxList) - - def S_m(self, prob): - if prob._eqLocs is 'EF' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m - - def S_e(self, prob): - if prob._eqLocs is 'FE' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e - - -class SrcFDEM_MagDipole(SrcFDEM): - - #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): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - - else: - srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob,b) - - def S_m(self, prob): - b_p = self.bPrimary(prob) - return -1j*omega(self.freq)*b_p - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_MagDipole_Bfield(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - #TODO: neither does moment - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl.T - - srcfct = SrcUtils.MagneticDipoleFields - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,bz)) - else: - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,by,bz)) - - return b - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob, b) - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_CircularLoop(SrcFDEM): - - #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): - self.freq = float(freq) - self.orientation = orientation - self.radius = radius - self.mu = mu - self.loc = loc - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) - - else: - srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) - ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) - az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return 1./self.mu*b - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - #################################################### # Survey #################################################### -class SurveyFDEM(Survey.BaseSurvey): +class Survey(SimPEG.Survey.BaseSurvey): """ docstring for SurveyFDEM """ - srcPair = SrcFDEM + srcPair = Src.BaseSrc def __init__(self, srcList, **kwargs): # Sort these by frequency self.srcList = srcList - Survey.BaseSurvey.__init__(self, **kwargs) + SimPEG.Survey.BaseSurvey.__init__(self, **kwargs) _freqDict = {} for src in srcList: @@ -481,7 +136,7 @@ class SurveyFDEM(Survey.BaseSurvey): return self._freqDict[freq] def projectFields(self, u): - data = Survey.Data(self) + data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: data[src, rx] = rx.projectFields(src, self.mesh, u) diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 110b4d1e..978972f5 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyFDEM import * -from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h +from SurveyFDEM import Rx, Src, Survey +from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h from FieldsFDEM import * \ No newline at end of file diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 5e82b7b2..e36d76b8 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Solver, Problem from SimPEG.Problem import BaseTimeProblem -from SimPEG.EM.Utils import SrcUtils +from SimPEG.EM.Utils import * from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 9878438d..7f6e5c04 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey -from SimPEG.EM.Utils import SrcUtils +from SimPEG.EM.Utils import * from BaseTDEM import FieldsTDEM @@ -87,11 +87,11 @@ class SrcTDEM_VMD_MVP(SrcTDEM): """Vertical magnetic dipole, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') + MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) + MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: raise Exception('Unknown mesh for VMD') @@ -109,11 +109,11 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): """Circular Loop, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) + MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) + MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: raise Exception('Unknown mesh for CircularLoop') diff --git a/SimPEG/EM/Utils/SrcUtils.py b/SimPEG/EM/Utils/AnalyticUtils.py similarity index 100% rename from SimPEG/EM/Utils/SrcUtils.py rename to SimPEG/EM/Utils/AnalyticUtils.py diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index 6e430cf9..18dddde9 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,5 @@ # import Sources # import Ana # import Solver -import EMUtils -import SrcUtils \ No newline at end of file +from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py new file mode 100644 index 00000000..8c703083 --- /dev/null +++ b/SimPEG/EM/Utils/testingUtils.py @@ -0,0 +1,75 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 + +def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): + cs = 5. + ncx, ncy, ncz = 6, 6, 6 + npad = 3 + hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] + hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] + hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] + mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) + + mapping = Maps.ExpMap(mesh) + + x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source + XYZ = Utils.ndgrid(x,x,np.r_[0.]) + Rx0 = EM.FDEM.Rx(XYZ, comp) + + Src = [] + + for SrcType in SrcList: + if SrcType is 'MagDipole': + Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + elif SrcType is 'MagDipole_Bfield': + Src.append(EM.FDEM.Src.MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + elif SrcType is 'CircularLoop': + Src.append(EM.FDEM.Src.CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + elif SrcType is 'RawVec': + if fdemType is 'e' or fdemType is 'b': + S_m = np.zeros(mesh.nF) + S_e = np.zeros(mesh.nE) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + + elif fdemType is 'h' or fdemType is 'j': + S_m = np.zeros(mesh.nE) + S_e = np.zeros(mesh.nF) + S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. + S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + + if verbose: + print ' Fetching %s problem' % (fdemType) + + if fdemType == 'e': + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_e(mesh, mapping=mapping) + + elif fdemType == 'b': + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + + elif fdemType == 'j': + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_j(mesh, mapping=mapping) + + elif fdemType == 'h': + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_h(mesh, mapping=mapping) + + else: + raise NotImplementedError() + prb.pair(survey) + + try: + from pymatsolver import MumpsSolver + prb.Solver = MumpsSolver + except ImportError, e: + pass + + return prb \ No newline at end of file diff --git a/SimPEG/Utils/__init__.py b/SimPEG/Utils/__init__.py index 6637a138..5280ae79 100644 --- a/SimPEG/Utils/__init__.py +++ b/SimPEG/Utils/__init__.py @@ -7,4 +7,4 @@ from ipythonutils import easyAnimate as animate from CounterUtils import * import ModelBuilder import SolverUtils - +from coordutils import * diff --git a/SimPEG/Utils/codeutils.py b/SimPEG/Utils/codeutils.py index 0ba57b2e..4a9a28a7 100644 --- a/SimPEG/Utils/codeutils.py +++ b/SimPEG/Utils/codeutils.py @@ -3,10 +3,7 @@ import time import numpy as np from functools import wraps - -class SimPEGMetaClass(type): - def __new__(cls, name, bases, attrs): - return super(SimPEGMetaClass, cls).__new__(cls, name, bases, attrs) +SimPEGMetaClass = type def memProfileWrapper(towrap, *funNames): """ diff --git a/SimPEG/Utils/coordutils.py b/SimPEG/Utils/coordutils.py new file mode 100644 index 00000000..260e1a3b --- /dev/null +++ b/SimPEG/Utils/coordutils.py @@ -0,0 +1,62 @@ +import numpy as np +from SimPEG.Utils import mkvc + +def rotationMatrixFromNormals(v0,v1,tol=1e-20): + """ + Performs the minimum number of rotations to define a rotation from the direction indicated by the vector n0 to the direction indicated by n1. + The axis of rotation is n0 x n1 + https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula + + :param numpy.array v0: vector of length 3 + :param numpy.array v1: vector of length 3 + :param tol = 1e-20: tolerance. If the norm of the cross product between the two vectors is below this, no rotation is performed + :rtype: numpy.array, 3x3 + :return: rotation matrix which rotates the frame so that n0 is aligned with n1 + + """ + + # ensure both n0, n1 are vectors of length 1 + assert len(v0) == 3, "Length of n0 should be 3" + assert len(v1) == 3, "Length of n1 should be 3" + + # ensure both are true normals + n0 = v0*1./np.linalg.norm(v0) + n1 = v1*1./np.linalg.norm(v1) + + n0dotn1 = n0.dot(n1) + + # define the rotation axis, which is the cross product of the two vectors + rotAx = np.cross(n0,n1) + + if np.linalg.norm(rotAx) < tol: + return np.eye(3,dtype=float) + + rotAx *= 1./np.linalg.norm(rotAx) + + cosT = n0dotn1/(np.linalg.norm(n0)*np.linalg.norm(n1)) + sinT = np.sqrt(1.-n0dotn1**2) + + ux = np.array([[0., -rotAx[2], rotAx[1]], [rotAx[2], 0., -rotAx[0]], [-rotAx[1], rotAx[0], 0.]],dtype=float) + + return np.eye(3,dtype=float) + sinT*ux + (1.-cosT)*(ux.dot(ux)) + + +def rotatePointsFromNormals(XYZ,n0,n1,x0=np.r_[0.,0.,0.]): + """ + rotates a grid so that the vector n0 is aligned with the vector n1 + + :param numpy.array n0: vector of length 3, should have norm 1 + :param numpy.array n1: vector of length 3, should have norm 1 + :param numpy.array x0: vector of length 3, point about which we perform the rotation + :rtype: numpy.array, 3x3 + :return: rotation matrix which rotates the frame so that n0 is aligned with n1 + """ + + R = rotationMatrixFromNormals(n0, n1) + + assert XYZ.shape[1] == 3, "Grid XYZ should be 3 wide" + assert len(x0) == 3, "x0 should have length 3" + + X0 = np.ones([XYZ.shape[0],1])*mkvc(x0) + + return (XYZ - X0).dot(R.T) + X0 # equivalent to (R*(XYZ - X0)).T + X0 \ No newline at end of file diff --git a/tests/em/__init__.py b/tests/em/examples/__init__.py similarity index 100% rename from tests/em/__init__.py rename to tests/em/examples/__init__.py diff --git a/tests/em/test_Examples.py b/tests/em/examples/test_Examples.py similarity index 100% rename from tests/em/test_Examples.py rename to tests/em/examples/test_Examples.py diff --git a/tests/em/fdem/forward/__init__.py b/tests/em/fdem/forward/__init__.py new file mode 100644 index 00000000..38d84328 --- /dev/null +++ b/tests/em/fdem/forward/__init__.py @@ -0,0 +1,11 @@ +if __name__ == '__main__': + import os + import glob + import unittest + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py similarity index 93% rename from tests/em/test_FDEM_analytics.py rename to tests/em/fdem/forward/test_FDEM_analytics.py index cb906262..9786e7c8 100644 --- a/tests/em/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -28,12 +28,12 @@ class FDEM_analyticTests(unittest.TestCase): x = np.linspace(-10,10,5) XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0]) - rxList = EM.FDEM.RxFDEM(XYZ, 'exi') - Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) + rxList = EM.FDEM.Rx(XYZ, 'exi') + Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) - survey = EM.FDEM.SurveyFDEM([Src0]) + survey = EM.FDEM.Survey([Src0]) - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) prb.pair(survey) try: @@ -114,19 +114,19 @@ class FDEM_analyticTests(unittest.TestCase): de = np.zeros(mesh.nF,dtype=complex) de[s_ind] = 1./csz - de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)] + de_p = [EM.FDEM.Src.RawVec_e([],freq,de/mesh.area)] - dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)] + dm_p = [EM.FDEM.Src.MagDipole([],freq,src_loc)] # Pair the problem and survey - surveye = EM.FDEM.SurveyFDEM(de_p) - surveym = EM.FDEM.SurveyFDEM(dm_p) + surveye = EM.FDEM.Survey(de_p) + surveym = EM.FDEM.Survey(dm_p) mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))] - prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) - prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) + prbe = EM.FDEM.Problem_h(mesh, mapping=mapping) + prbm = EM.FDEM.Problem_e(mesh, mapping=mapping) prbe.pair(surveye) # pair problem and survey prbm.pair(surveym) diff --git a/tests/em/test_FDEMCasing.py b/tests/em/fdem/forward/test_FDEM_casing.py similarity index 100% rename from tests/em/test_FDEMCasing.py rename to tests/em/fdem/forward/test_FDEM_casing.py diff --git a/tests/em/fdem/forward/test_FDEM_forward.py b/tests/em/fdem/forward/test_FDEM_forward.py new file mode 100644 index 00000000..437f3708 --- /dev/null +++ b/tests/em/fdem/forward/test_FDEM_forward.py @@ -0,0 +1,127 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem + +testEB = True +testHJ = True + +verbose = False + +TOL = 1e-5 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 1e-1 +addrandoms = True + +SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] + + +def crossCheckTest(fdemType, comp): + + l2norm = lambda r: np.sqrt(r.dot(r)) + + prb1 = getFDEMProblem(fdemType, comp, SrcList, freq, verbose) + mesh = prb1.mesh + print 'Cross Checking Forward: %s formulation - %s' % (fdemType, comp) + m = np.log(np.ones(mesh.nC)*CONDUCTIVITY) + mu = np.log(np.ones(mesh.nC)*MU) + + if addrandoms is True: + m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(mesh.nC)*MU*1e-1 + + # prb1.PropMap.PropModel.mu = mu + # prb1.PropMap.PropModel.mui = 1./mu + survey1 = prb1.survey + d1 = survey1.dpred(m) + + if verbose: + print ' Problem 1 solved' + + if fdemType == 'e': + prb2 = getFDEMProblem('b', comp, SrcList, freq, verbose) + elif fdemType == 'b': + prb2 = getFDEMProblem('e', comp, SrcList, freq, verbose) + elif fdemType == 'j': + prb2 = getFDEMProblem('h', comp, SrcList, freq, verbose) + elif fdemType == 'h': + prb2 = getFDEMProblem('j', comp, SrcList, freq, verbose) + else: + raise NotImplementedError() + + # prb2.mu = mu + survey2 = prb2.survey + d2 = survey2.dpred(m) + + if verbose: + print ' Problem 2 solved' + + r = d2-d1 + l2r = l2norm(r) + + tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) + print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol + return l2r < tol + + +class FDEM_CrossCheck(unittest.TestCase): + if testEB: + def test_EB_CrossCheck_exr_Eform(self): + self.assertTrue(crossCheckTest('e', 'exr')) + 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')) + 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')) + + if testHJ: + 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')) + 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() \ No newline at end of file diff --git a/tests/em/fdem/inverse/adjoint/__init__.py b/tests/em/fdem/inverse/adjoint/__init__.py new file mode 100644 index 00000000..38d84328 --- /dev/null +++ b/tests/em/fdem/inverse/adjoint/__init__.py @@ -0,0 +1,11 @@ +if __name__ == '__main__': + import os + import glob + import unittest + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py b/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py new file mode 100644 index 00000000..f77f2131 --- /dev/null +++ b/tests/em/fdem/inverse/adjoint/test_FDEM_adjoint.py @@ -0,0 +1,156 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem + +testEB = True +testHJ = True + +verbose = False + +TOL = 1e-5 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 1e-1 +addrandoms = True + +SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + +def adjointTest(fdemType, comp): + prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + print 'Adjoint %s formulation - %s' % (fdemType, comp) + + m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) + mu = np.ones(prb.mesh.nC)*MU + + if addrandoms is True: + m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 + + survey = prb.survey + # prb.PropMap.PropModel.mu = mu + # prb.PropMap.PropModel.mui = 1./mu + u = prb.fields(m) + + v = np.random.rand(survey.nD) + w = np.random.rand(prb.mesh.nC) + + vJw = v.dot(prb.Jvec(m, w, u)) + wJtv = w.dot(prb.Jtvec(m, v, u)) + tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) + print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol + return np.abs(vJw - wJtv) < tol + +class FDEM_AdjointTests(unittest.TestCase): + if testEB: + def test_Jtvec_adjointTest_exr_Eform(self): + self.assertTrue(adjointTest('e', 'exr')) + def test_Jtvec_adjointTest_eyr_Eform(self): + self.assertTrue(adjointTest('e', 'eyr')) + def test_Jtvec_adjointTest_ezr_Eform(self): + self.assertTrue(adjointTest('e', 'ezr')) + def test_Jtvec_adjointTest_exi_Eform(self): + self.assertTrue(adjointTest('e', 'exi')) + def test_Jtvec_adjointTest_eyi_Eform(self): + self.assertTrue(adjointTest('e', 'eyi')) + def test_Jtvec_adjointTest_ezi_Eform(self): + self.assertTrue(adjointTest('e', 'ezi')) + + def test_Jtvec_adjointTest_bxr_Eform(self): + self.assertTrue(adjointTest('e', 'bxr')) + def test_Jtvec_adjointTest_byr_Eform(self): + self.assertTrue(adjointTest('e', 'byr')) + def test_Jtvec_adjointTest_bzr_Eform(self): + self.assertTrue(adjointTest('e', 'bzr')) + def test_Jtvec_adjointTest_bxi_Eform(self): + self.assertTrue(adjointTest('e', 'bxi')) + def test_Jtvec_adjointTest_byi_Eform(self): + self.assertTrue(adjointTest('e', 'byi')) + def test_Jtvec_adjointTest_bzi_Eform(self): + self.assertTrue(adjointTest('e', 'bzi')) + + def test_Jtvec_adjointTest_exr_Bform(self): + self.assertTrue(adjointTest('b', 'exr')) + def test_Jtvec_adjointTest_eyr_Bform(self): + self.assertTrue(adjointTest('b', 'eyr')) + def test_Jtvec_adjointTest_ezr_Bform(self): + self.assertTrue(adjointTest('b', 'ezr')) + def test_Jtvec_adjointTest_exi_Bform(self): + self.assertTrue(adjointTest('b', 'exi')) + def test_Jtvec_adjointTest_eyi_Bform(self): + self.assertTrue(adjointTest('b', 'eyi')) + def test_Jtvec_adjointTest_ezi_Bform(self): + self.assertTrue(adjointTest('b', 'ezi')) + def test_Jtvec_adjointTest_bxr_Bform(self): + self.assertTrue(adjointTest('b', 'bxr')) + def test_Jtvec_adjointTest_byr_Bform(self): + self.assertTrue(adjointTest('b', 'byr')) + def test_Jtvec_adjointTest_bzr_Bform(self): + self.assertTrue(adjointTest('b', 'bzr')) + def test_Jtvec_adjointTest_bxi_Bform(self): + self.assertTrue(adjointTest('b', 'bxi')) + def test_Jtvec_adjointTest_byi_Bform(self): + self.assertTrue(adjointTest('b', 'byi')) + def test_Jtvec_adjointTest_bzi_Bform(self): + self.assertTrue(adjointTest('b', 'bzi')) + + + if testHJ: + def test_Jtvec_adjointTest_jxr_Jform(self): + self.assertTrue(adjointTest('j', 'jxr')) + def test_Jtvec_adjointTest_jyr_Jform(self): + self.assertTrue(adjointTest('j', 'jyr')) + def test_Jtvec_adjointTest_jzr_Jform(self): + self.assertTrue(adjointTest('j', 'jzr')) + def test_Jtvec_adjointTest_jxi_Jform(self): + self.assertTrue(adjointTest('j', 'jxi')) + def test_Jtvec_adjointTest_jyi_Jform(self): + self.assertTrue(adjointTest('j', 'jyi')) + def test_Jtvec_adjointTest_jzi_Jform(self): + self.assertTrue(adjointTest('j', 'jzi')) + + def test_Jtvec_adjointTest_hxr_Jform(self): + self.assertTrue(adjointTest('j', 'hxr')) + def test_Jtvec_adjointTest_hyr_Jform(self): + self.assertTrue(adjointTest('j', 'hyr')) + def test_Jtvec_adjointTest_hzr_Jform(self): + self.assertTrue(adjointTest('j', 'hzr')) + def test_Jtvec_adjointTest_hxi_Jform(self): + self.assertTrue(adjointTest('j', 'hxi')) + def test_Jtvec_adjointTest_hyi_Jform(self): + self.assertTrue(adjointTest('j', 'hyi')) + def test_Jtvec_adjointTest_hzi_Jform(self): + self.assertTrue(adjointTest('j', 'hzi')) + + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'hxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'hyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'hzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'hxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'hyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'hzi')) + + def test_Jtvec_adjointTest_hxr_Hform(self): + self.assertTrue(adjointTest('h', 'jxr')) + def test_Jtvec_adjointTest_hyr_Hform(self): + self.assertTrue(adjointTest('h', 'jyr')) + def test_Jtvec_adjointTest_hzr_Hform(self): + self.assertTrue(adjointTest('h', 'jzr')) + def test_Jtvec_adjointTest_hxi_Hform(self): + self.assertTrue(adjointTest('h', 'jxi')) + def test_Jtvec_adjointTest_hyi_Hform(self): + self.assertTrue(adjointTest('h', 'jyi')) + def test_Jtvec_adjointTest_hzi_Hform(self): + self.assertTrue(adjointTest('h', 'jzi')) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/fdem/inverse/derivs/__init__.py b/tests/em/fdem/inverse/derivs/__init__.py new file mode 100644 index 00000000..38d84328 --- /dev/null +++ b/tests/em/fdem/inverse/derivs/__init__.py @@ -0,0 +1,11 @@ +if __name__ == '__main__': + import os + import glob + import unittest + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py new file mode 100644 index 00000000..52108c4e --- /dev/null +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -0,0 +1,154 @@ +import unittest +from SimPEG import * +from SimPEG import EM +import sys +from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem + +testDerivs = True +testEB = True +testHJ = True + +verbose = False + +TOL = 1e-5 +FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order +CONDUCTIVITY = 1e1 +MU = mu_0 +freq = 1e-1 +addrandoms = True + +SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' + + +def derivTest(fdemType, comp): + + prb = getFDEMProblem(fdemType, comp, [SrcType], freq) + print '%s formulation - %s' % (fdemType, comp) + x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) + mu = np.log(np.ones(prb.mesh.nC)*MU) + + if addrandoms is True: + x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 + mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1 + + # prb.PropMap.PropModel.mu = mu + # prb.PropMap.PropModel.mui = 1./mu + + survey = prb.survey + def fun(x): + return survey.dpred(x), lambda x: prb.Jvec(x0, x) + return Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR) + + +class FDEM_DerivTests(unittest.TestCase): + + if testEB: + def test_Jvec_exr_Eform(self): + self.assertTrue(derivTest('e', 'exr')) + def test_Jvec_eyr_Eform(self): + self.assertTrue(derivTest('e', 'eyr')) + def test_Jvec_ezr_Eform(self): + self.assertTrue(derivTest('e', 'ezr')) + def test_Jvec_exi_Eform(self): + self.assertTrue(derivTest('e', 'exi')) + def test_Jvec_eyi_Eform(self): + self.assertTrue(derivTest('e', 'eyi')) + def test_Jvec_ezi_Eform(self): + self.assertTrue(derivTest('e', 'ezi')) + + def test_Jvec_bxr_Eform(self): + self.assertTrue(derivTest('e', 'bxr')) + def test_Jvec_byr_Eform(self): + self.assertTrue(derivTest('e', 'byr')) + def test_Jvec_bzr_Eform(self): + self.assertTrue(derivTest('e', 'bzr')) + def test_Jvec_bxi_Eform(self): + self.assertTrue(derivTest('e', 'bxi')) + def test_Jvec_byi_Eform(self): + self.assertTrue(derivTest('e', 'byi')) + def test_Jvec_bzi_Eform(self): + self.assertTrue(derivTest('e', 'bzi')) + + def test_Jvec_exr_Bform(self): + self.assertTrue(derivTest('b', 'exr')) + def test_Jvec_eyr_Bform(self): + self.assertTrue(derivTest('b', 'eyr')) + def test_Jvec_ezr_Bform(self): + self.assertTrue(derivTest('b', 'ezr')) + def test_Jvec_exi_Bform(self): + self.assertTrue(derivTest('b', 'exi')) + def test_Jvec_eyi_Bform(self): + self.assertTrue(derivTest('b', 'eyi')) + def test_Jvec_ezi_Bform(self): + self.assertTrue(derivTest('b', 'ezi')) + + def test_Jvec_bxr_Bform(self): + self.assertTrue(derivTest('b', 'bxr')) + def test_Jvec_byr_Bform(self): + self.assertTrue(derivTest('b', 'byr')) + def test_Jvec_bzr_Bform(self): + self.assertTrue(derivTest('b', 'bzr')) + def test_Jvec_bxi_Bform(self): + self.assertTrue(derivTest('b', 'bxi')) + def test_Jvec_byi_Bform(self): + self.assertTrue(derivTest('b', 'byi')) + def test_Jvec_bzi_Bform(self): + self.assertTrue(derivTest('b', 'bzi')) + + if testHJ: + 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')) + + +if __name__ == '__main__': + unittest.main() diff --git a/tests/em/tdem/__init__.py b/tests/em/tdem/__init__.py new file mode 100644 index 00000000..38d84328 --- /dev/null +++ b/tests/em/tdem/__init__.py @@ -0,0 +1,11 @@ +if __name__ == '__main__': + import os + import glob + import unittest + test_file_strings = glob.glob('test_*.py') + module_strings = [str[0:len(str)-3] for str in test_file_strings] + suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str + in module_strings] + testSuite = unittest.TestSuite(suites) + + unittest.TextTestRunner(verbosity=2).run(testSuite) diff --git a/tests/em/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py similarity index 100% rename from tests/em/test_TDEM_b_DerivAdjoint.py rename to tests/em/tdem/test_TDEM_b_DerivAdjoint.py diff --git a/tests/em/test_TDEM_b_MultiSrc_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py similarity index 100% rename from tests/em/test_TDEM_b_MultiSrc_DerivAdjoint.py rename to tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py diff --git a/tests/em/test_TDEM_combos.py b/tests/em/tdem/test_TDEM_combos.py similarity index 100% rename from tests/em/test_TDEM_combos.py rename to tests/em/tdem/test_TDEM_combos.py diff --git a/tests/em/test_TDEM_forward_Analytic.py b/tests/em/tdem/test_TDEM_forward_Analytic.py similarity index 100% rename from tests/em/test_TDEM_forward_Analytic.py rename to tests/em/tdem/test_TDEM_forward_Analytic.py diff --git a/tests/em/test_FDEM.py b/tests/em/test_FDEM.py deleted file mode 100644 index 23bfa0d2..00000000 --- a/tests/em/test_FDEM.py +++ /dev/null @@ -1,478 +0,0 @@ -import unittest -from SimPEG import * -from SimPEG import EM -import sys -from scipy.constants import mu_0 - -testDerivs = True -testCrossCheck = True -testAdjoint = True -testEB = True -testHJ = True - -verbose = False - -TOL = 1e-5 -FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order -CONDUCTIVITY = 1e1 -MU = mu_0 -freq = 1e-1 -addrandoms = True - -SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' - - -def getProblem(fdemType, comp): - cs = 5. - ncx, ncy, ncz = 6, 6, 6 - npad = 3 - hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] - hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] - hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] - mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) - - mapping = Maps.ExpMap(mesh) - - x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source - XYZ = Utils.ndgrid(x,x,np.r_[0.]) - Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - - if SrcType is 'MagDipole': - Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - elif SrcType is 'MagDipole_Bfield': - Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - elif SrcType is 'CircularLoop': - Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - - if verbose: - print ' Fetching %s problem' % (fdemType) - - if fdemType == 'e': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nF) - S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) - - elif fdemType == 'b': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nF) - S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) - - elif fdemType == 'j': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nE) - S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) - - elif fdemType == 'h': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nE) - S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) - - else: - raise NotImplementedError() - prb.pair(survey) - - try: - from pymatsolver import MumpsSolver - prb.Solver = MumpsSolver - except ImportError, e: - pass - - return prb - -def adjointTest(fdemType, comp): - prb = getProblem(fdemType, comp) - print 'Adjoint %s formulation - %s' % (fdemType, comp) - - m = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) - mu = np.ones(prb.mesh.nC)*MU - - if addrandoms is True: - m = m + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1 - - survey = prb.survey - # prb.PropMap.PropModel.mu = mu - # prb.PropMap.PropModel.mui = 1./mu - u = prb.fields(m) - - v = np.random.rand(survey.nD) - w = np.random.rand(prb.mesh.nC) - - vJw = v.dot(prb.Jvec(m, w, u)) - wJtv = w.dot(prb.Jtvec(m, v, u)) - tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR]) - print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol - return np.abs(vJw - wJtv) < tol - - -def derivTest(fdemType, comp): - - prb = getProblem(fdemType, comp) - print '%s formulation - %s' % (fdemType, comp) - x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) - mu = np.log(np.ones(prb.mesh.nC)*MU) - - if addrandoms is True: - x0 = x0 + np.random.randn(prb.mapping.nP)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(prb.mapping.nP)*MU*1e-1 - - # prb.PropMap.PropModel.mu = mu - # prb.PropMap.PropModel.mui = 1./mu - - survey = prb.survey - def fun(x): - return survey.dpred(x), lambda x: prb.Jvec(x0, x) - 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)*np.log(CONDUCTIVITY)*1e-1 - mu = mu + np.random.randn(mesh.nC)*MU*1e-1 - - # prb1.PropMap.PropModel.mu = mu - # prb1.PropMap.PropModel.mui = 1./mu - survey1 = prb1.survey - d1 = survey1.dpred(m) - - if verbose: - print ' Problem 1 solved' - - if fdemType == 'e': - prb2 = 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 - d2 = survey2.dpred(m) - - if verbose: - print ' Problem 2 solved' - - r = d2-d1 - l2r = l2norm(r) - - tol = np.max([TOL*(10**int(np.log10(l2norm(d1)))),FLR]) - print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol - return l2r < tol - - - - - -class FDEM_DerivTests(unittest.TestCase): - - if testDerivs: - if testEB: - def test_Jvec_exr_Eform(self): - self.assertTrue(derivTest('e', 'exr')) - def test_Jvec_eyr_Eform(self): - self.assertTrue(derivTest('e', 'eyr')) - def test_Jvec_ezr_Eform(self): - self.assertTrue(derivTest('e', 'ezr')) - def test_Jvec_exi_Eform(self): - self.assertTrue(derivTest('e', 'exi')) - def test_Jvec_eyi_Eform(self): - self.assertTrue(derivTest('e', 'eyi')) - def test_Jvec_ezi_Eform(self): - self.assertTrue(derivTest('e', 'ezi')) - - def test_Jvec_bxr_Eform(self): - self.assertTrue(derivTest('e', 'bxr')) - def test_Jvec_byr_Eform(self): - self.assertTrue(derivTest('e', 'byr')) - def test_Jvec_bzr_Eform(self): - self.assertTrue(derivTest('e', 'bzr')) - def test_Jvec_bxi_Eform(self): - self.assertTrue(derivTest('e', 'bxi')) - def test_Jvec_byi_Eform(self): - self.assertTrue(derivTest('e', 'byi')) - def test_Jvec_bzi_Eform(self): - self.assertTrue(derivTest('e', 'bzi')) - - def test_Jvec_exr_Bform(self): - self.assertTrue(derivTest('b', 'exr')) - def test_Jvec_eyr_Bform(self): - self.assertTrue(derivTest('b', 'eyr')) - def test_Jvec_ezr_Bform(self): - self.assertTrue(derivTest('b', 'ezr')) - def test_Jvec_exi_Bform(self): - self.assertTrue(derivTest('b', 'exi')) - def test_Jvec_eyi_Bform(self): - self.assertTrue(derivTest('b', 'eyi')) - def test_Jvec_ezi_Bform(self): - self.assertTrue(derivTest('b', 'ezi')) - - def test_Jvec_bxr_Bform(self): - self.assertTrue(derivTest('b', 'bxr')) - def test_Jvec_byr_Bform(self): - self.assertTrue(derivTest('b', 'byr')) - def test_Jvec_bzr_Bform(self): - self.assertTrue(derivTest('b', 'bzr')) - def test_Jvec_bxi_Bform(self): - self.assertTrue(derivTest('b', 'bxi')) - def test_Jvec_byi_Bform(self): - self.assertTrue(derivTest('b', 'byi')) - def test_Jvec_bzi_Bform(self): - self.assertTrue(derivTest('b', 'bzi')) - - if testHJ: - 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')) - - - if testAdjoint: - if testEB: - def test_Jtvec_adjointTest_exr_Eform(self): - self.assertTrue(adjointTest('e', 'exr')) - def test_Jtvec_adjointTest_eyr_Eform(self): - self.assertTrue(adjointTest('e', 'eyr')) - def test_Jtvec_adjointTest_ezr_Eform(self): - self.assertTrue(adjointTest('e', 'ezr')) - def test_Jtvec_adjointTest_exi_Eform(self): - self.assertTrue(adjointTest('e', 'exi')) - def test_Jtvec_adjointTest_eyi_Eform(self): - self.assertTrue(adjointTest('e', 'eyi')) - def test_Jtvec_adjointTest_ezi_Eform(self): - self.assertTrue(adjointTest('e', 'ezi')) - - def test_Jtvec_adjointTest_bxr_Eform(self): - self.assertTrue(adjointTest('e', 'bxr')) - def test_Jtvec_adjointTest_byr_Eform(self): - self.assertTrue(adjointTest('e', 'byr')) - def test_Jtvec_adjointTest_bzr_Eform(self): - self.assertTrue(adjointTest('e', 'bzr')) - def test_Jtvec_adjointTest_bxi_Eform(self): - self.assertTrue(adjointTest('e', 'bxi')) - def test_Jtvec_adjointTest_byi_Eform(self): - self.assertTrue(adjointTest('e', 'byi')) - def test_Jtvec_adjointTest_bzi_Eform(self): - self.assertTrue(adjointTest('e', 'bzi')) - - def test_Jtvec_adjointTest_exr_Bform(self): - self.assertTrue(adjointTest('b', 'exr')) - def test_Jtvec_adjointTest_eyr_Bform(self): - self.assertTrue(adjointTest('b', 'eyr')) - def test_Jtvec_adjointTest_ezr_Bform(self): - self.assertTrue(adjointTest('b', 'ezr')) - def test_Jtvec_adjointTest_exi_Bform(self): - self.assertTrue(adjointTest('b', 'exi')) - def test_Jtvec_adjointTest_eyi_Bform(self): - self.assertTrue(adjointTest('b', 'eyi')) - def test_Jtvec_adjointTest_ezi_Bform(self): - self.assertTrue(adjointTest('b', 'ezi')) - def test_Jtvec_adjointTest_bxr_Bform(self): - self.assertTrue(adjointTest('b', 'bxr')) - def test_Jtvec_adjointTest_byr_Bform(self): - self.assertTrue(adjointTest('b', 'byr')) - def test_Jtvec_adjointTest_bzr_Bform(self): - self.assertTrue(adjointTest('b', 'bzr')) - def test_Jtvec_adjointTest_bxi_Bform(self): - self.assertTrue(adjointTest('b', 'bxi')) - def test_Jtvec_adjointTest_byi_Bform(self): - self.assertTrue(adjointTest('b', 'byi')) - def test_Jtvec_adjointTest_bzi_Bform(self): - self.assertTrue(adjointTest('b', 'bzi')) - - - if testHJ: - def test_Jtvec_adjointTest_jxr_Jform(self): - self.assertTrue(adjointTest('j', 'jxr')) - def test_Jtvec_adjointTest_jyr_Jform(self): - self.assertTrue(adjointTest('j', 'jyr')) - def test_Jtvec_adjointTest_jzr_Jform(self): - self.assertTrue(adjointTest('j', 'jzr')) - def test_Jtvec_adjointTest_jxi_Jform(self): - self.assertTrue(adjointTest('j', 'jxi')) - def test_Jtvec_adjointTest_jyi_Jform(self): - self.assertTrue(adjointTest('j', 'jyi')) - def test_Jtvec_adjointTest_jzi_Jform(self): - self.assertTrue(adjointTest('j', 'jzi')) - - def test_Jtvec_adjointTest_hxr_Jform(self): - self.assertTrue(adjointTest('j', 'hxr')) - def test_Jtvec_adjointTest_hyr_Jform(self): - self.assertTrue(adjointTest('j', 'hyr')) - def test_Jtvec_adjointTest_hzr_Jform(self): - self.assertTrue(adjointTest('j', 'hzr')) - def test_Jtvec_adjointTest_hxi_Jform(self): - self.assertTrue(adjointTest('j', 'hxi')) - def test_Jtvec_adjointTest_hyi_Jform(self): - self.assertTrue(adjointTest('j', 'hyi')) - def test_Jtvec_adjointTest_hzi_Jform(self): - self.assertTrue(adjointTest('j', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'hxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'hyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'hzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'hxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'hyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'hzi')) - - def test_Jtvec_adjointTest_hxr_Hform(self): - self.assertTrue(adjointTest('h', 'jxr')) - def test_Jtvec_adjointTest_hyr_Hform(self): - self.assertTrue(adjointTest('h', 'jyr')) - def test_Jtvec_adjointTest_hzr_Hform(self): - self.assertTrue(adjointTest('h', 'jzr')) - def test_Jtvec_adjointTest_hxi_Hform(self): - self.assertTrue(adjointTest('h', 'jxi')) - def test_Jtvec_adjointTest_hyi_Hform(self): - self.assertTrue(adjointTest('h', 'jyi')) - def test_Jtvec_adjointTest_hzi_Hform(self): - self.assertTrue(adjointTest('h', 'jzi')) - - - if testCrossCheck: - if testEB: - 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')) - 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')) - - if testHJ: - 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')) - 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() diff --git a/tests/utils/test_coordutils.py b/tests/utils/test_coordutils.py new file mode 100644 index 00000000..24f8451c --- /dev/null +++ b/tests/utils/test_coordutils.py @@ -0,0 +1,56 @@ +import unittest, os +import numpy as np +from SimPEG import Utils + +tol = 1e-15 + +class coorUtilsTest(unittest.TestCase): + + def test_rotationMatrixFromNormals(self): + np.random.seed(0) + v0 = np.random.rand(3) + v0 *= 1./np.linalg.norm(v0) + + np.random.seed(5) + v1 = np.random.rand(3) + v1 *= 1./np.linalg.norm(v1) + + Rf = Utils.coordutils.rotationMatrixFromNormals(v0,v1) + Ri = Utils.coordutils.rotationMatrixFromNormals(v1,v0) + + self.assertTrue(np.linalg.norm(Utils.mkvc(Rf.dot(v0) - v1)) < tol) + self.assertTrue(np.linalg.norm(Utils.mkvc(Ri.dot(v1) - v0)) < tol) + + def test_rotatePointsFromNormals(self): + np.random.seed(10) + v0 = np.random.rand(3) + v0*= 1./np.linalg.norm(v0) + + np.random.seed(15) + v1 = np.random.rand(3) + v1*= 1./np.linalg.norm(v1) + + v2 = Utils.mkvc(Utils.coordutils.rotatePointsFromNormals(Utils.mkvc(v0,2).T,v0,v1)) + + self.assertTrue(np.linalg.norm(v2-v1) < tol) + + def test_rotateMatrixFromNormals(self): + np.random.seed(20) + n0 = np.random.rand(3) + n0 *= 1./np.linalg.norm(n0) + + np.random.seed(25) + n1 = np.random.rand(3) + n1 *= 1./np.linalg.norm(n1) + + np.random.seed(30) + scale = np.random.rand(100,1) + XYZ0 = scale * n0 + XYZ1 = scale * n1 + + XYZ2 = Utils.coordutils.rotatePointsFromNormals(XYZ0,n0,n1) + self.assertTrue(np.linalg.norm(Utils.mkvc(XYZ1) - Utils.mkvc(XYZ2))/np.linalg.norm(Utils.mkvc(XYZ1)) < tol) + +if __name__ == '__main__': + unittest.main() +