diff --git a/SimPEG/EM/Base.py b/SimPEG/EM/Base.py index a16cdb91..431086a4 100644 --- a/SimPEG/EM/Base.py +++ b/SimPEG/EM/Base.py @@ -61,6 +61,15 @@ class BaseEMProblem(Problem.BaseProblem): self._Me = self.mesh.getEdgeInnerProduct() return self._Me + @property + def MeI(self): + """ + Edge inner product matrix + """ + if getattr(self, '_MeI', None) is None: + self._MeI = self.mesh.getEdgeInnerProduct(invMat=True) + return self._MeI + @property def Mf(self): """ @@ -70,6 +79,15 @@ class BaseEMProblem(Problem.BaseProblem): self._Mf = self.mesh.getFaceInnerProduct() return self._Mf + @property + def MfI(self): + """ + Face inner product matrix + """ + if getattr(self, '_MfI', None) is None: + self._MfI = self.mesh.getFaceInnerProduct(invMat=True) + return self._MfI + # ----- Magnetic Permeability ----- # @property diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index e2193973..47f37616 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -160,9 +160,9 @@ class Fields(SimPEG.Problem.Fields): return self._jDeriv_u(src, v, adjoint), self._jDeriv_m(src, v, adjoint) return np.array(self._jDeriv_u(src, du_dm_v, adjoint) + self._jDeriv_m(src, v, adjoint), dtype = complex) -class Fields_e(Fields): +class Fields3D_e(Fields): """ - Fields object for Problem_e. + Fields object for Problem3D_e. :param Mesh mesh: mesh :param Survey survey: survey @@ -426,9 +426,9 @@ class Fields_e(Fields): -class Fields_b(Fields): +class Fields3D_b(Fields): """ - Fields object for Problem_b. + Fields object for Problem3D_b. :param Mesh mesh: mesh :param Survey survey: survey @@ -693,9 +693,9 @@ class Fields_b(Fields): return Zero() -class Fields_j(Fields): +class Fields3D_j(Fields): """ - Fields object for Problem_j. + Fields object for Problem3D_j. :param Mesh mesh: mesh :param Survey survey: survey @@ -988,9 +988,9 @@ class Fields_j(Fields): return 1./(1j * omega(src.freq)) * VI * (self._aveE2CCV * ( s_mDeriv(v) - self._edgeCurl.T * ( self._MfRhoDeriv(jSolution) * v ) ) ) -class Fields_h(Fields): +class Fields3D_h(Fields): """ - Fields object for Problem_h. + Fields object for Problem3D_h. :param Mesh mesh: mesh :param Survey survey: survey diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/ProblemFDEM.py similarity index 95% rename from SimPEG/EM/FDEM/FDEM.py rename to SimPEG/EM/FDEM/ProblemFDEM.py index caca7602..2aed1d91 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/ProblemFDEM.py @@ -1,7 +1,7 @@ from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 from SurveyFDEM import Survey as SurveyFDEM -from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j +from FieldsFDEM import Fields, Fields3D_e, Fields3D_b, Fields3D_h, Fields3D_j from SimPEG.EM.Base import BaseEMProblem from SimPEG.EM.Utils import omega @@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ {\mathbf{C}^{\\top} \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}} - if using the E-B formulation (:code:`Problem_e` - or :code:`Problem_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. + if using the E-B formulation (:code:`Problem3D_e` + or :code:`Problem3D_b`). Note that in this case, :math:`\mathbf{s_e}` is an integrated quantity. If we write Maxwell's equations in terms of \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) @@ -28,7 +28,7 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C}^{\\top} \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. + if using the H-J formulation (:code:`Problem3D_j` or :code:`Problem3D_h`). Note that here, :math:`\mathbf{s_m}` is an integrated quantity. The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ @@ -87,7 +87,7 @@ class BaseFDEMProblem(BaseEMProblem): du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v ) for rx in src.rxList: - df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_dmFun = getattr(f, '_{0}Deriv'.format(rx.projField), None) df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False) Jv[src, rx] = rx.evalDeriv(src, self.mesh, f, df_dm_v) Ainv.clean() @@ -125,7 +125,7 @@ class BaseFDEMProblem(BaseEMProblem): for rx in src.rxList: PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx], adjoint=True) # wrt f, need possibility wrt m - df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None) + df_duTFun = getattr(f, '_{0}Deriv'.format(rx.projField), None) df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True) ATinvdf_duT = ATinv * df_duT @@ -137,10 +137,9 @@ class BaseFDEMProblem(BaseEMProblem): df_dmT = df_dmT + du_dmT # TODO: this should be taken care of by the reciever? - real_or_imag = rx.projComp - if real_or_imag is 'real': + if rx.component is 'real': Jtv += np.array(df_dmT, dtype=complex).real - elif real_or_imag is 'imag': + elif rx.component is 'imag': Jtv += - np.array(df_dmT, dtype=complex).real else: raise Exception('Must be real or imag') @@ -177,7 +176,7 @@ class BaseFDEMProblem(BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## -class Problem_e(BaseFDEMProblem): +class Problem3D_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -199,7 +198,7 @@ class Problem_e(BaseFDEMProblem): _solutionType = 'eSolution' _formulation = 'EB' - fieldsPair = Fields_e + fieldsPair = Fields3D_e def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -288,7 +287,7 @@ class Problem_e(BaseFDEMProblem): return C.T * (MfMui * s_mDeriv(v)) -1j * omega(freq) * s_eDeriv(v) -class Problem_b(BaseFDEMProblem): +class Problem3D_b(BaseFDEMProblem): """ We eliminate :math:`\mathbf{e}` using @@ -310,7 +309,7 @@ class Problem_b(BaseFDEMProblem): _solutionType = 'bSolution' _formulation = 'EB' - fieldsPair = Fields_b + fieldsPair = Fields3D_b def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -436,7 +435,7 @@ class Problem_b(BaseFDEMProblem): ########################################################################################## -class Problem_j(BaseFDEMProblem): +class Problem3D_j(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{h}\\\) using @@ -458,7 +457,7 @@ class Problem_j(BaseFDEMProblem): _solutionType = 'jSolution' _formulation = 'HJ' - fieldsPair = Fields_j + fieldsPair = Fields3D_j def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -577,7 +576,7 @@ class Problem_j(BaseFDEMProblem): -class Problem_h(BaseFDEMProblem): +class Problem3D_h(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{j}\\\) using @@ -596,7 +595,7 @@ class Problem_h(BaseFDEMProblem): _solutionType = 'hSolution' _formulation = 'HJ' - fieldsPair = Fields_h + fieldsPair = Fields3D_h def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py new file mode 100644 index 00000000..53d6c722 --- /dev/null +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -0,0 +1,126 @@ +import SimPEG +from SimPEG import sp + +class BaseRx(SimPEG.Survey.BaseRx): + """ + Frequency domain receiver base class + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + assert(orientation in ['x','y','z']), "Orientation %s not known. Orientation must be in 'x', 'y', 'z'. Arbitrary orientations have not yet been implemented."%orientation + assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component + + self.projComp = orientation + self.component = component + + SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx + + def projGLoc(self, u): + """Grid Location projection (e.g. Ex Fy ...)""" + return u._GLoc(self.projField) + self.projComp + + def eval(self, src, mesh, f): + """ + Project fields to recievers to get data. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields f: fields object + :rtype: numpy.ndarray + :return: fields projected to recievers + """ + + P = self.getP(mesh, self.projGLoc(f)) + f_part_complex = f[src, self.projField] + f_part = getattr(f_part_complex, self.component) # get the real or imag component + + return P*f_part + + def evalDeriv(self, src, mesh, f, v, adjoint=False): + """ + Derivative of projected fields with respect to the inversion model times a vector. + + :param Source src: FDEM source + :param Mesh mesh: mesh used + :param Fields f: fields object + :param numpy.ndarray v: vector to multiply + :rtype: numpy.ndarray + :return: fields projected to recievers + """ + + P = self.getP(mesh, self.projGLoc(f)) + + if not adjoint: + Pv_complex = P * v + Pv = getattr(Pv_complex, self.component) + elif adjoint: + Pv_real = P.T * v + + if self.component == 'imag': + Pv = 1j*Pv_real + elif self.component == 'real': + Pv = Pv_real.astype(complex) + else: + raise NotImplementedError('must be real or imag') + + return Pv + + +class Point_e(BaseRx): + """ + Electric field FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'e' + super(Point_e, self).__init__(locs, orientation, component) + + +class Point_b(BaseRx): + """ + Magnetic flux FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'b' + super(Point_b, self).__init__(locs, orientation, component) + + +class Point_h(BaseRx): + """ + Magnetic field FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'h' + super(Point_h, self).__init__(locs, orientation, component) + + +class Point_j(BaseRx): + """ + Current density FDEM receiver + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + :param string orientation: receiver orientation 'x', 'y' or 'z' + :param string component: real or imaginary component 'real' or 'imag' + """ + + def __init__(self, locs, orientation=None, component=None): + self.projField = 'j' + super(Point_j, self).__init__(locs, orientation, component) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 87967dd5..9bb6ffce 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -9,8 +9,14 @@ class BaseSrc(Survey.BaseSrc): """ freq = None - # rxPair = RxFDEM - integrate = True + integrate = False + _ePrimary = None + _bPrimary = None + _hPrimary = None + _jPrimary = None + + def __init__(self, rxList, **kwargs): + Survey.BaseSrc.__init__(self, rxList, **kwargs) def eval(self, prob): """ @@ -51,7 +57,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary magnetic flux density """ - return Zero() + if self._bPrimary is None: + return Zero() + return self._bPrimary def hPrimary(self, prob): """ @@ -61,7 +69,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary magnetic field """ - return Zero() + if self._hPrimary is None: + return Zero() + return self._hPrimary def ePrimary(self, prob): """ @@ -71,7 +81,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary electric field """ - return Zero() + if self._ePrimary is None: + return Zero() + return self._ePrimary def jPrimary(self, prob): """ @@ -81,7 +93,9 @@ class BaseSrc(Survey.BaseSrc): :rtype: numpy.ndarray :return: primary current density """ - return Zero() + if self._jPrimary is None: + return Zero() + return self._jPrimary def s_m(self, prob): """ @@ -136,15 +150,14 @@ class RawVec_e(BaseSrc): :param list rxList: receiver list :param float freq: frequency :param numpy.array s_e: electric source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + def __init__(self, rxList, freq, s_e, **kwargs): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - self.integrate = integrate - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_e(self, prob): """ @@ -166,15 +179,14 @@ class RawVec_m(BaseSrc): :param float freq: frequency :param rxList: receiver list :param numpy.array s_m: magnetic source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): + def __init__(self, rxList, freq, s_m, **kwargs): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()): self._s_m = np.array(s_m, dtype=complex) self.freq = float(freq) - self.integrate = integrate - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_m(self, prob): """ @@ -197,14 +209,13 @@ class RawVec(BaseSrc): :param float freq: frequency :param numpy.array s_m: magnetic source term :param numpy.array s_e: electric source term - :param bool integrate: Integrate the source term (multiply by Me) [True] + :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_m, s_e, integrate=True): + def __init__(self, rxList, freq, s_m, s_e, **kwargs): 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) + BaseSrc.__init__(self, rxList, **kwargs) def s_m(self, prob): """ @@ -278,14 +289,13 @@ class MagDipole(BaseSrc): :param float mu: background magnetic permeability """ - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0): + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0, **kwargs): self.freq = float(freq) self.loc = loc self.orientation = orientation assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." self.moment = moment self.mu = mu - self.integrate = False BaseSrc.__init__(self, rxList) def bPrimary(self, prob): @@ -543,7 +553,7 @@ class CircularLoop(BaseSrc): 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) + a = MagneticLoopVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) else: srcfct = MagneticDipoleVectorPotential diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 1552a12c..46ae2523 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -4,126 +4,9 @@ from SimPEG.EM.Base import BaseEMSurvey from scipy.constants import mu_0 from SimPEG.Utils import Zero, Identity import SrcFDEM as Src +import RxFDEM as Rx from SimPEG import sp - -#################################################### -# Receivers -#################################################### - -class Rx(SimPEG.Survey.BaseRx): - """ - Frequency domain receivers - - :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) - :param string rxType: reciever type from knownRxTypes - """ - - knownRxTypes = { - 'exr':['e', 'x', 'real'], - 'eyr':['e', 'y', 'real'], - 'ezr':['e', 'z', 'real'], - 'exi':['e', 'x', 'imag'], - 'eyi':['e', 'y', 'imag'], - 'ezi':['e', 'z', 'imag'], - - 'bxr':['b', 'x', 'real'], - 'byr':['b', 'y', 'real'], - 'bzr':['b', 'z', 'real'], - 'bxi':['b', 'x', 'imag'], - 'byi':['b', 'y', 'imag'], - 'bzi':['b', 'z', 'imag'], - - 'jxr':['j', 'x', 'real'], - 'jyr':['j', 'y', 'real'], - 'jzr':['j', 'z', 'real'], - 'jxi':['j', 'x', 'imag'], - 'jyi':['j', 'y', 'imag'], - 'jzi':['j', 'z', 'imag'], - - 'hxr':['h', 'x', 'real'], - 'hyr':['h', 'y', 'real'], - 'hzr':['h', 'z', 'real'], - 'hxi':['h', 'x', 'imag'], - 'hyi':['h', 'y', 'imag'], - 'hzi':['h', 'z', 'imag'], - } - radius = None - - def __init__(self, locs, rxType): - SimPEG.Survey.BaseRx.__init__(self, locs, rxType) - - @property - def projField(self): - """Field Type projection (e.g. e b ...)""" - return self.knownRxTypes[self.rxType][0] - - @property - def projComp(self): - """Component projection (real/imag)""" - return self.knownRxTypes[self.rxType][2] - - def projGLoc(self, u): - """Grid Location projection (e.g. Ex Fy ...)""" - return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1] - - def eval(self, src, mesh, f): - """ - Project fields to recievers to get data. - - :param Source src: FDEM source - :param Mesh mesh: mesh used - :param Fields f: fields object - :rtype: numpy.ndarray - :return: fields projected to recievers - """ - # projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0]) - # projGLoc += self.knownRxTypes[self.rxType][1] - - P = self.getP(mesh, self.projGLoc(f)) - f_part_complex = f[src, self.projField] - # get the real or imag component - real_or_imag = self.projComp - f_part = getattr(f_part_complex, real_or_imag) - - return P*f_part - - def evalDeriv(self, src, mesh, f, v, adjoint=False): - """ - Derivative of projected fields with respect to the inversion model times a vector. - - :param Source src: FDEM source - :param Mesh mesh: mesh used - :param Fields f: fields object - :param numpy.ndarray v: vector to multiply - :rtype: numpy.ndarray - :return: fields projected to recievers - """ - - P = self.getP(mesh, self.projGLoc(f)) - - if not adjoint: - Pv_complex = P * v - real_or_imag = self.projComp - Pv = getattr(Pv_complex, real_or_imag) - elif adjoint: - Pv_real = P.T * v - - real_or_imag = self.projComp - if real_or_imag == 'imag': - Pv = 1j*Pv_real - elif real_or_imag == 'real': - Pv = Pv_real.astype(complex) - else: - raise NotImplementedError('must be real or imag') - - return Pv - - -#################################################### -# Survey -#################################################### - class Survey(BaseEMSurvey): """ Frequency domain electromagnetic survey @@ -132,7 +15,7 @@ class Survey(BaseEMSurvey): """ srcPair = Src.BaseSrc - rxPair = Rx + rxPair = Rx.BaseRx def __init__(self, srcList, **kwargs): # Sort these by frequency diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 978972f5..c4ff6451 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,5 @@ -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 +from SurveyFDEM import Survey +import SrcFDEM as Src +import RxFDEM as Rx +from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h +from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 569f8e6d..22c925b6 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -20,56 +20,61 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) if useMu is True: - mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] + mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))] else: mapping = Maps.ExpMap(mesh) x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5)) - Rx0 = EM.FDEM.Rx(XYZ, comp) + Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0]) + if comp[2] == 'r': + real_or_imag = 'real' + elif comp[2] == 'i': + real_or_imag = 'imag' + rx0 = Rx0(XYZ, comp[1], 'imag') Src = [] for SrcType in SrcList: if SrcType is 'MagDipole': - Src.append(EM.FDEM.Src.MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.])) + 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.])) + 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.])) + 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])] = 1e-3 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3 - Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([rx0], freq, S_m, mesh.getEdgeInnerProduct()*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])] = 1e-3 S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3 - Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([rx0], freq, mesh.getEdgeInnerProduct()*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) + prb = EM.FDEM.Problem3D_e(mesh, mapping=mapping) elif fdemType == 'b': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) elif fdemType == 'j': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_j(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_j(mesh, mapping=mapping) elif fdemType == 'h': survey = EM.FDEM.Survey(Src) - prb = EM.FDEM.Problem_h(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_h(mesh, mapping=mapping) else: raise NotImplementedError() @@ -90,7 +95,7 @@ def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useM prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose) mesh = prb1.mesh print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp) - + logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY) mu = np.ones(mesh.nC)*MU diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 29f51ed4..4275903d 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -42,8 +42,8 @@ def run(plotIt=True): ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) - rxOffset=10. - bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') + rxOffset=10. + bzi = EM.FDEM.Rx.Point_b(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag') freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) @@ -51,7 +51,7 @@ def run(plotIt=True): srcList = [EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z') for freq in freqs] survey = EM.FDEM.Survey(srcList) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) try: from pymatsolver import MumpsSolver diff --git a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py index 76af4a3d..930b452a 100644 --- a/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py +++ b/SimPEG/Examples/EM_Schenkel_Morrison_Casing.py @@ -215,7 +215,7 @@ def run(plotIt=True): # ------------ Problem and Survey --------------- survey = FDEM.Survey(sg_p + dg_p) mapping = [('sigma', Maps.IdentityMap(mesh))] - problem = FDEM.Problem_h(mesh, mapping=mapping) + problem = FDEM.Problem3D_h(mesh, mapping=mapping) problem.pair(survey) # ------------- Solve --------------------------- diff --git a/SimPEG/MT/BaseMT.py b/SimPEG/MT/BaseMT.py index c201dfb0..579f590f 100644 --- a/SimPEG/MT/BaseMT.py +++ b/SimPEG/MT/BaseMT.py @@ -1,5 +1,5 @@ from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np -from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem +from SimPEG.EM.FDEM.ProblemFDEM import BaseFDEMProblem from SurveyMT import Survey, Data from FieldsMT import BaseMTFields diff --git a/docs/examples/DC_Forward_PseudoSection.rst b/docs/examples/DC_Forward_PseudoSection.rst index 80cf0307..29855d4a 100644 --- a/docs/examples/DC_Forward_PseudoSection.rst +++ b/docs/examples/DC_Forward_PseudoSection.rst @@ -22,7 +22,7 @@ radi = Radius of spheres [r1,r2] param = Conductivity of background and two spheres [m0,m1,m2] stype = survey type "pdp" (pole dipole) or "dpdp" (dipole dipole) dtype = Data type "appr" (app res) | "appc" (app cond) | "volt" (potential) -Created by @fourndo on Mon Feb 01 19:28:06 2016 +Created by @fourndo diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index 9786e7c8..6f283666 100644 --- a/tests/em/fdem/forward/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.Rx(XYZ, 'exi') + rxList = EM.FDEM.Rx.Point_e(XYZ, orientation='x', component='imag') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.Survey([Src0]) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) prb.pair(survey) try: @@ -125,8 +125,8 @@ class FDEM_analyticTests(unittest.TestCase): mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))] - prbe = EM.FDEM.Problem_h(mesh, mapping=mapping) - prbm = EM.FDEM.Problem_e(mesh, mapping=mapping) + prbe = EM.FDEM.Problem3D_h(mesh, mapping=mapping) + prbm = EM.FDEM.Problem3D_e(mesh, mapping=mapping) prbe.pair(surveye) # pair problem and survey prbm.pair(surveym) diff --git a/tests/em/fdem/forward/test_FDEM_forwardHB.py b/tests/em/fdem/forward/test_FDEM_forwardHB.py index 545a5014..8fce615b 100644 --- a/tests/em/fdem/forward/test_FDEM_forwardHB.py +++ b/tests/em/fdem/forward/test_FDEM_forwardHB.py @@ -12,7 +12,7 @@ testBH = True verbose = False TOLEJHB = 1 # averaging and more sensitive to boundary condition violations (ie. the impact of violating the boundary conditions in each case is different.) -#TODO: choose better testing parameters to lower this +#TODO: choose better testing parameters to lower this SrcList = ['RawVec', 'MagDipole_Bfield', 'MagDipole', 'CircularLoop'] @@ -125,4 +125,4 @@ class FDEM_CrossCheck(unittest.TestCase): self.assertTrue(crossCheckTest(SrcList, 'b', 'h', 'hzi', verbose=verbose, TOL=TOLEJHB)) if __name__ == '__main__': - unittest.main() \ No newline at end of file + unittest.main()