From 936a7aaadcf5292a2fa99f18d18f57f78fb4a318 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 29 Mar 2016 21:32:08 -0700 Subject: [PATCH 01/16] make integrate = False default for all sources --- SimPEG/EM/FDEM/SrcFDEM.py | 26 ++++++++++++-------------- SimPEG/EM/Utils/testingUtils.py | 8 ++++---- 2 files changed, 16 insertions(+), 18 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 87967dd5..6da924be 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -9,8 +9,10 @@ class BaseSrc(Survey.BaseSrc): """ freq = None - # rxPair = RxFDEM - integrate = True + integrate = False + + def __init__(self, rxList, **kwargs): + Survey.BaseSrc.__init__(self, rxList, **kwargs) def eval(self, prob): """ @@ -136,13 +138,12 @@ 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): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - self.integrate = integrate BaseSrc.__init__(self, rxList) @@ -166,13 +167,12 @@ 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()): self._s_m = np.array(s_m, dtype=complex) self.freq = float(freq) - self.integrate = integrate BaseSrc.__init__(self, rxList) @@ -197,14 +197,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 +277,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 +541,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/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 569f8e6d..e5e75c06 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -20,7 +20,7 @@ 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) @@ -43,14 +43,14 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): 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) @@ -90,7 +90,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 From d8eeb7cd0511d9a30341b6f58d40da6905d42750 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 May 2016 11:18:36 -0700 Subject: [PATCH 02/16] use Problem3D_assumption, Fields3D_assumption --- SimPEG/EM/Base.py | 18 +++++++++++++++ SimPEG/EM/FDEM/FDEM.py | 24 ++++++++++---------- SimPEG/EM/FDEM/FieldsFDEM.py | 16 ++++++------- SimPEG/EM/FDEM/__init__.py | 4 ++-- SimPEG/EM/Utils/testingUtils.py | 8 +++---- tests/em/fdem/forward/test_FDEM_analytics.py | 6 ++--- tests/em/fdem/forward/test_FDEM_forwardHB.py | 4 ++-- 7 files changed, 49 insertions(+), 31 deletions(-) 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/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index caca7602..78e7fa61 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.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}\\\) """ @@ -177,7 +177,7 @@ class BaseFDEMProblem(BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## -class Problem_e(BaseFDEMProblem): +class Problem3D_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -199,7 +199,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 +288,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 +310,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 +436,7 @@ class Problem_b(BaseFDEMProblem): ########################################################################################## -class Problem_j(BaseFDEMProblem): +class Problem3D_j(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{h}\\\) using @@ -458,7 +458,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 +577,7 @@ class Problem_j(BaseFDEMProblem): -class Problem_h(BaseFDEMProblem): +class Problem3D_h(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{j}\\\) using @@ -596,7 +596,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/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/__init__.py b/SimPEG/EM/FDEM/__init__.py index 978972f5..1565d062 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ 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 FDEM 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 e5e75c06..00461e2d 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -57,19 +57,19 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): 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() diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index 9786e7c8..f596d638 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -33,7 +33,7 @@ class FDEM_analyticTests(unittest.TestCase): 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() From 52747c09262d5fa835d15f445ef9a3fccd4dba4d Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 May 2016 11:35:28 -0700 Subject: [PATCH 03/16] update example --- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index e76b2439..824c4f95 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -42,7 +42,7 @@ def run(plotIt=True): ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5) - rxOffset=10. + rxOffset=10. bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') freqs = np.logspace(1,3,10) @@ -52,7 +52,7 @@ def run(plotIt=True): [srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs] survey = EM.FDEM.Survey(srcList) - prb = EM.FDEM.Problem_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem3D_b(mesh, mapping=mapping) try: from pymatsolver import MumpsSolver From f7c46ed83be33915ca91ad64265a90758941f151 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 May 2016 12:41:06 -0700 Subject: [PATCH 04/16] Rx classes for FDEM --- SimPEG/EM/FDEM/RxFDEM.py | 105 ++++++++++++++++ SimPEG/EM/FDEM/SurveyFDEM.py | 121 +------------------ SimPEG/EM/FDEM/__init__.py | 4 +- SimPEG/EM/Utils/testingUtils.py | 17 ++- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 2 +- tests/em/fdem/forward/test_FDEM_analytics.py | 2 +- 6 files changed, 123 insertions(+), 128 deletions(-) create mode 100644 SimPEG/EM/FDEM/RxFDEM.py diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py new file mode 100644 index 00000000..d0445ff8 --- /dev/null +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -0,0 +1,105 @@ +import SimPEG +from SimPEG.EM.Utils import * +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 + +class BaseRx(SimPEG.Survey.BaseRx): + """ + Frequency domain receivers + + :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) + """ + + def __init__(self, locs, orientation=None, real_or_imag=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(real_or_imag in ['real', 'imag']), "'real_or_imag' must be 'real' or 'imag', not %s"%real_or_imag + + self.projComp = orientation + self.real_or_imag = real_or_imag + + 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] + # get the real or imag component + f_part = getattr(f_part_complex, self.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, self.real_or_imag) + elif adjoint: + Pv_real = P.T * v + + # real_or_imag = self.projComp + if self.real_or_imag == 'imag': + Pv = 1j*Pv_real + elif self.real_or_imag == 'real': + Pv = Pv_real.astype(complex) + else: + raise NotImplementedError('must be real or imag') + + return Pv + + +class eField(BaseRx): + + def __init__(self, locs, orientation=None, real_or_imag=None): + self.projField = 'e' + BaseRx.__init__(self, locs, orientation, real_or_imag) + + +class bField(BaseRx): + + def __init__(self, locs, orientation=None, real_or_imag=None): + self.projField = 'b' + BaseRx.__init__(self, locs, orientation, real_or_imag) + + +class hField(BaseRx): + + def __init__(self, locs, orientation=None, real_or_imag=None): + self.projField = 'h' + BaseRx.__init__(self, locs, orientation, real_or_imag) + + +class jField(BaseRx): + + def __init__(self, locs, orientation=None, real_or_imag=None): + self.projField = 'j' + BaseRx.__init__(self, locs, orientation, real_or_imag) 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 1565d062..1701fe3e 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,5 @@ -from SurveyFDEM import Rx, Src, Survey +from SurveyFDEM import Survey +import SrcFDEM as Src +import RxFDEM as Rx from FDEM 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 00461e2d..c3fc50d2 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -26,31 +26,36 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): 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, comp[0] + 'Field') + 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, mesh.getEdgeInnerProduct()*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, mesh.getEdgeInnerProduct()*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) diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 824c4f95..1764f28c 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -43,7 +43,7 @@ def run(plotIt=True): rxOffset=10. - bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi') + bzi = EM.FDEM.Rx.bField(np.array([[rxOffset, 0., 1e-3]]), orientation='z', real_or_imag='imag') freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index f596d638..0ea43ca7 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -28,7 +28,7 @@ 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.eField(XYZ, orientation='x', real_or_imag='imag') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.Survey([Src0]) From cb042ac9380a20eb431d6b39ef6c8eff554ac94c Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 May 2016 13:00:29 -0700 Subject: [PATCH 05/16] cleanup imports, docstrings --- SimPEG/EM/FDEM/RxFDEM.py | 43 ++++++++++++++++++++++++++++++---------- 1 file changed, 32 insertions(+), 11 deletions(-) diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py index d0445ff8..ef58d807 100644 --- a/SimPEG/EM/FDEM/RxFDEM.py +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -1,17 +1,13 @@ import SimPEG -from SimPEG.EM.Utils import * -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 class BaseRx(SimPEG.Survey.BaseRx): """ - Frequency domain receivers + 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 real_or_imag: real or imaginary component 'real' or 'imag' """ def __init__(self, locs, orientation=None, real_or_imag=None): @@ -40,8 +36,7 @@ class BaseRx(SimPEG.Survey.BaseRx): P = self.getP(mesh, self.projGLoc(f)) f_part_complex = f[src, self.projField] - # get the real or imag component - f_part = getattr(f_part_complex, self.real_or_imag) + f_part = getattr(f_part_complex, self.real_or_imag) # get the real or imag component return P*f_part @@ -61,12 +56,10 @@ class BaseRx(SimPEG.Survey.BaseRx): if not adjoint: Pv_complex = P * v - # real_or_imag = self.projComp Pv = getattr(Pv_complex, self.real_or_imag) elif adjoint: Pv_real = P.T * v - # real_or_imag = self.projComp if self.real_or_imag == 'imag': Pv = 1j*Pv_real elif self.real_or_imag == 'real': @@ -78,6 +71,13 @@ class BaseRx(SimPEG.Survey.BaseRx): class eField(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 real_or_imag: real or imaginary component 'real' or 'imag' + """ def __init__(self, locs, orientation=None, real_or_imag=None): self.projField = 'e' @@ -85,6 +85,13 @@ class eField(BaseRx): class bField(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 real_or_imag: real or imaginary component 'real' or 'imag' + """ def __init__(self, locs, orientation=None, real_or_imag=None): self.projField = 'b' @@ -92,6 +99,13 @@ class bField(BaseRx): class hField(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 real_or_imag: real or imaginary component 'real' or 'imag' + """ def __init__(self, locs, orientation=None, real_or_imag=None): self.projField = 'h' @@ -99,6 +113,13 @@ class hField(BaseRx): class jField(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 real_or_imag: real or imaginary component 'real' or 'imag' + """ def __init__(self, locs, orientation=None, real_or_imag=None): self.projField = 'j' From 0a714663d38d255595346dd656a2bdf36cf81c3a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Sun, 8 May 2016 13:12:37 -0700 Subject: [PATCH 06/16] update Jtvec to work with Rx classes --- SimPEG/EM/FDEM/FDEM.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 78e7fa61..026230c6 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -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.real_or_imag is 'real': Jtv += np.array(df_dmT, dtype=complex).real - elif real_or_imag is 'imag': + elif rx.real_or_imag is 'imag': Jtv += - np.array(df_dmT, dtype=complex).real else: raise Exception('Must be real or imag') From 73c219ff5c3ce9af371a479a050783223cfbe89a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 9 May 2016 12:29:52 -0700 Subject: [PATCH 07/16] updated Problem naming in casing example --- SimPEG/Examples/EM_Schenkel_Morrison_Casing.py | 2 +- docs/examples/DC_Forward_PseudoSection.rst | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/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 From 11e6b452c96cf419e0f36224975c7c8fdbba8a8b Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 17:26:16 -0700 Subject: [PATCH 08/16] renamed FDEM.py to ProblemFDEM.py, changed real_or_imag to component --- SimPEG/EM/FDEM/{FDEM.py => ProblemFDEM.py} | 4 +-- SimPEG/EM/FDEM/RxFDEM.py | 40 +++++++++++----------- SimPEG/EM/FDEM/__init__.py | 2 +- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 2 +- 4 files changed, 24 insertions(+), 24 deletions(-) rename SimPEG/EM/FDEM/{FDEM.py => ProblemFDEM.py} (99%) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/ProblemFDEM.py similarity index 99% rename from SimPEG/EM/FDEM/FDEM.py rename to SimPEG/EM/FDEM/ProblemFDEM.py index 026230c6..75b05bf9 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/ProblemFDEM.py @@ -137,9 +137,9 @@ class BaseFDEMProblem(BaseEMProblem): df_dmT = df_dmT + du_dmT # TODO: this should be taken care of by the reciever? - if rx.real_or_imag is 'real': + if rx.component is 'real': Jtv += np.array(df_dmT, dtype=complex).real - elif rx.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') diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py index ef58d807..f15040fb 100644 --- a/SimPEG/EM/FDEM/RxFDEM.py +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -7,15 +7,15 @@ class BaseRx(SimPEG.Survey.BaseRx): :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) :param string orientation: receiver orientation 'x', 'y' or 'z' - :param string real_or_imag: real or imaginary component 'real' or 'imag' + :param string component: real or imaginary component 'real' or 'imag' """ - def __init__(self, locs, orientation=None, real_or_imag=None): + 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(real_or_imag in ['real', 'imag']), "'real_or_imag' must be 'real' or 'imag', not %s"%real_or_imag + assert(component in ['real', 'imag']), "'component' must be 'real' or 'imag', not %s"%component self.projComp = orientation - self.real_or_imag = real_or_imag + self.component = component SimPEG.Survey.BaseRx.__init__(self, locs, rxType=None) #TODO: remove rxType from baseRx @@ -36,7 +36,7 @@ class BaseRx(SimPEG.Survey.BaseRx): P = self.getP(mesh, self.projGLoc(f)) f_part_complex = f[src, self.projField] - f_part = getattr(f_part_complex, self.real_or_imag) # get the real or imag component + f_part = getattr(f_part_complex, self.component) # get the real or imag component return P*f_part @@ -56,13 +56,13 @@ class BaseRx(SimPEG.Survey.BaseRx): if not adjoint: Pv_complex = P * v - Pv = getattr(Pv_complex, self.real_or_imag) + Pv = getattr(Pv_complex, self.component) elif adjoint: Pv_real = P.T * v - if self.real_or_imag == 'imag': + if self.component == 'imag': Pv = 1j*Pv_real - elif self.real_or_imag == 'real': + elif self.component == 'real': Pv = Pv_real.astype(complex) else: raise NotImplementedError('must be real or imag') @@ -76,12 +76,12 @@ class eField(BaseRx): :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) :param string orientation: receiver orientation 'x', 'y' or 'z' - :param string real_or_imag: real or imaginary component 'real' or 'imag' + :param string component: real or imaginary component 'real' or 'imag' """ - def __init__(self, locs, orientation=None, real_or_imag=None): + def __init__(self, locs, orientation=None, component=None): self.projField = 'e' - BaseRx.__init__(self, locs, orientation, real_or_imag) + BaseRx.__init__(self, locs, orientation, component) class bField(BaseRx): @@ -90,12 +90,12 @@ class bField(BaseRx): :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) :param string orientation: receiver orientation 'x', 'y' or 'z' - :param string real_or_imag: real or imaginary component 'real' or 'imag' + :param string component: real or imaginary component 'real' or 'imag' """ - def __init__(self, locs, orientation=None, real_or_imag=None): + def __init__(self, locs, orientation=None, component=None): self.projField = 'b' - BaseRx.__init__(self, locs, orientation, real_or_imag) + BaseRx.__init__(self, locs, orientation, component) class hField(BaseRx): @@ -104,12 +104,12 @@ class hField(BaseRx): :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) :param string orientation: receiver orientation 'x', 'y' or 'z' - :param string real_or_imag: real or imaginary component 'real' or 'imag' + :param string component: real or imaginary component 'real' or 'imag' """ - def __init__(self, locs, orientation=None, real_or_imag=None): + def __init__(self, locs, orientation=None, component=None): self.projField = 'h' - BaseRx.__init__(self, locs, orientation, real_or_imag) + BaseRx.__init__(self, locs, orientation, component) class jField(BaseRx): @@ -118,9 +118,9 @@ class jField(BaseRx): :param numpy.ndarray locs: receiver locations (ie. :code:`np.r_[x,y,z]`) :param string orientation: receiver orientation 'x', 'y' or 'z' - :param string real_or_imag: real or imaginary component 'real' or 'imag' + :param string component: real or imaginary component 'real' or 'imag' """ - def __init__(self, locs, orientation=None, real_or_imag=None): + def __init__(self, locs, orientation=None, component=None): self.projField = 'j' - BaseRx.__init__(self, locs, orientation, real_or_imag) + BaseRx.__init__(self, locs, orientation, component) diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 1701fe3e..c4ff6451 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,5 +1,5 @@ from SurveyFDEM import Survey import SrcFDEM as Src import RxFDEM as Rx -from FDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h +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/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 0a5c59a0..95e3f09f 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -43,7 +43,7 @@ def run(plotIt=True): rxOffset=10. - bzi = EM.FDEM.Rx.bField(np.array([[rxOffset, 0., 1e-3]]), orientation='z', real_or_imag='imag') + bzi = EM.FDEM.Rx.bField(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag') freqs = np.logspace(1,3,10) srcLoc = np.array([0., 0., 10.]) From c1b1c2467f737aa6462070826b4ee5277f343dc4 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 10 May 2016 19:57:16 -0700 Subject: [PATCH 09/16] import from ProblemFDEM in baseMT, fixed a missed real_or_imag --> component --- SimPEG/MT/BaseMT.py | 2 +- tests/em/fdem/forward/test_FDEM_analytics.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) 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/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index 0ea43ca7..b9de4a26 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -28,7 +28,7 @@ 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.eField(XYZ, orientation='x', real_or_imag='imag') + rxList = EM.FDEM.Rx.eField(XYZ, orientation='x', component='imag') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.Survey([Src0]) From a690cab13132681f7b9972f3c72547883af5c422 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 11 May 2016 09:05:13 -0700 Subject: [PATCH 10/16] simple field receivers are `Point` receivers --- SimPEG/EM/FDEM/RxFDEM.py | 16 ++++++++-------- SimPEG/EM/Utils/testingUtils.py | 2 +- SimPEG/Examples/EM_FDEM_1D_Inversion.py | 2 +- tests/em/fdem/forward/test_FDEM_analytics.py | 2 +- 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/SimPEG/EM/FDEM/RxFDEM.py b/SimPEG/EM/FDEM/RxFDEM.py index f15040fb..53d6c722 100644 --- a/SimPEG/EM/FDEM/RxFDEM.py +++ b/SimPEG/EM/FDEM/RxFDEM.py @@ -70,7 +70,7 @@ class BaseRx(SimPEG.Survey.BaseRx): return Pv -class eField(BaseRx): +class Point_e(BaseRx): """ Electric field FDEM receiver @@ -81,10 +81,10 @@ class eField(BaseRx): def __init__(self, locs, orientation=None, component=None): self.projField = 'e' - BaseRx.__init__(self, locs, orientation, component) + super(Point_e, self).__init__(locs, orientation, component) -class bField(BaseRx): +class Point_b(BaseRx): """ Magnetic flux FDEM receiver @@ -95,10 +95,10 @@ class bField(BaseRx): def __init__(self, locs, orientation=None, component=None): self.projField = 'b' - BaseRx.__init__(self, locs, orientation, component) + super(Point_b, self).__init__(locs, orientation, component) -class hField(BaseRx): +class Point_h(BaseRx): """ Magnetic field FDEM receiver @@ -109,10 +109,10 @@ class hField(BaseRx): def __init__(self, locs, orientation=None, component=None): self.projField = 'h' - BaseRx.__init__(self, locs, orientation, component) + super(Point_h, self).__init__(locs, orientation, component) -class jField(BaseRx): +class Point_j(BaseRx): """ Current density FDEM receiver @@ -123,4 +123,4 @@ class jField(BaseRx): def __init__(self, locs, orientation=None, component=None): self.projField = 'j' - BaseRx.__init__(self, locs, orientation, component) + super(Point_j, self).__init__(locs, orientation, component) diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index c3fc50d2..22c925b6 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -26,7 +26,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False): 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 = getattr(EM.FDEM.Rx, comp[0] + 'Field') + Rx0 = getattr(EM.FDEM.Rx, 'Point_' + comp[0]) if comp[2] == 'r': real_or_imag = 'real' elif comp[2] == 'i': diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index 95e3f09f..4275903d 100644 --- a/SimPEG/Examples/EM_FDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_FDEM_1D_Inversion.py @@ -43,7 +43,7 @@ def run(plotIt=True): rxOffset=10. - bzi = EM.FDEM.Rx.bField(np.array([[rxOffset, 0., 1e-3]]), orientation='z', component='imag') + 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.]) diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index b9de4a26..6f283666 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -28,7 +28,7 @@ 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.eField(XYZ, orientation='x', component='imag') + 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]) From 029171fb1d7652207f7ba10423d0f3a8ceef89d1 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 11 May 2016 09:09:26 -0700 Subject: [PATCH 11/16] use .format for strings --- SimPEG/EM/FDEM/ProblemFDEM.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/FDEM/ProblemFDEM.py b/SimPEG/EM/FDEM/ProblemFDEM.py index 75b05bf9..2aed1d91 100644 --- a/SimPEG/EM/FDEM/ProblemFDEM.py +++ b/SimPEG/EM/FDEM/ProblemFDEM.py @@ -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 From c88263234bd0d5cac71d6ddb4bf9a8d50dbcbd45 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 May 2016 23:56:06 -0700 Subject: [PATCH 12/16] rename FDEM --> ProblemFDEM --- SimPEG/EM/FDEM/{FDEM.py => ProblemFDEM.py} | 0 SimPEG/EM/FDEM/__init__.py | 2 +- 2 files changed, 1 insertion(+), 1 deletion(-) rename SimPEG/EM/FDEM/{FDEM.py => ProblemFDEM.py} (100%) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/ProblemFDEM.py similarity index 100% rename from SimPEG/EM/FDEM/FDEM.py rename to SimPEG/EM/FDEM/ProblemFDEM.py diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 1701fe3e..c4ff6451 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,5 +1,5 @@ from SurveyFDEM import Survey import SrcFDEM as Src import RxFDEM as Rx -from FDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h +from ProblemFDEM import Problem3D_e, Problem3D_b, Problem3D_j, Problem3D_h from FieldsFDEM import Fields3D_e, Fields3D_b, Fields3D_j, Fields3D_h From 10c8791514e9bf636d71331adbda2f63327fa8a5 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 18 May 2016 00:33:30 -0700 Subject: [PATCH 13/16] update base MT to import ProblemFDEM --- SimPEG/MT/BaseMT.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From e25b496ab0117d2a6b4f8e73ba3f6ca96a09c7c0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 23 May 2016 12:07:14 -0700 Subject: [PATCH 14/16] allow kwarg input of primary fields --- SimPEG/EM/FDEM/SrcFDEM.py | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 6da924be..2c97adf1 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -10,6 +10,10 @@ class BaseSrc(Survey.BaseSrc): freq = None integrate = False + _ePrimary = None + _bPrimary = None + _hPrimary = None + _jPrimary = None def __init__(self, rxList, **kwargs): Survey.BaseSrc.__init__(self, rxList, **kwargs) @@ -53,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): """ @@ -63,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): """ @@ -73,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): """ @@ -83,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): """ @@ -144,6 +156,7 @@ class RawVec_e(BaseSrc): def __init__(self, rxList, freq, s_e): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) + self._ePrimary = ePrimary BaseSrc.__init__(self, rxList) From beca0203dfcd0cfc7197512a90004f4cab2a5491 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 23 May 2016 12:11:03 -0700 Subject: [PATCH 15/16] typo fix --- SimPEG/EM/FDEM/SrcFDEM.py | 1 - 1 file changed, 1 deletion(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 2c97adf1..e7650164 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -156,7 +156,6 @@ class RawVec_e(BaseSrc): def __init__(self, rxList, freq, s_e): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - self._ePrimary = ePrimary BaseSrc.__init__(self, rxList) From 2c87a50d29b9888cba3f2701d19c9a44924432dd Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Mon, 23 May 2016 12:21:29 -0700 Subject: [PATCH 16/16] add kwargs to raw vec e,m --- SimPEG/EM/FDEM/SrcFDEM.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index e7650164..9bb6ffce 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -153,11 +153,11 @@ class RawVec_e(BaseSrc): :param bool integrate: Integrate the source term (multiply by Me) [False] """ - def __init__(self, rxList, freq, s_e): + def __init__(self, rxList, freq, s_e, **kwargs): self._s_e = np.array(s_e, dtype=complex) self.freq = float(freq) - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_e(self, prob): """ @@ -182,11 +182,11 @@ class RawVec_m(BaseSrc): :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) - BaseSrc.__init__(self, rxList) + BaseSrc.__init__(self, rxList, **kwargs) def s_m(self, prob): """