diff --git a/SimPEG/EM/TDEM/SrcTDEM.py b/SimPEG/EM/TDEM/SrcTDEM.py new file mode 100644 index 00000000..ff4f9c66 --- /dev/null +++ b/SimPEG/EM/TDEM/SrcTDEM.py @@ -0,0 +1,204 @@ +import SimPEG +from SimPEG import np, Utils +from SimPEG.Utils import Zero, Identity +from scipy.constants import mu_0 +from SimPEG.EM.Utils import * + +#################################################### +# Sources +#################################################### + +class BaseWaveform(object): + + def __init__(self, offTime=0., hasInitialFields=False): + self.offTime = offTime + self.hasInitialFields = hasInitialFields + + def _assertMatchesPair(self, pair): + assert (isinstance(self, pair) + ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__) + + def eval(self, time): + raise NotImplementedError + + def evalDeriv(self, time): + raise NotImplementedError # needed for E-formulation + + +class StepOffWaveform(BaseWaveform): + + def __init__(self, offTime=0.): + BaseWaveform.__init__(self, offTime, hasInitialFields=True) + + def eval(self, time): + return 0. + + +class RawWaveform(BaseWaveform): + + def __init__(self, offTime=0.): + BaseWaveform.__init__(self, offTime, hasInitialFields=True) + + def eval(self, time): + raise NotImplementedError('RawWaveform has not been implemented, you should write it!') + + +class TriangularWaveform(BaseWaveform): + + def __init__(self, offTime=0.): + BaseWaveform.__init__(self, offTime, hasInitialFields=True) + + def eval(self, time): + raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!') + + + +class BaseSrc(SimPEG.Survey.BaseSrc): + + # rxPair = Rx + integrate = True + waveformPair = BaseWaveform + + @property + def waveform(self): + "A waveform instance is not None" + return getattr(self, '_waveform', None) + @waveform.setter + def waveform(self, val): + if self.waveform is None: + val._assertMatchesPair(self.waveformPair) + self._mapping = val + else: + self._mapping = self.PropMap(val) + + + def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs): + self.waveform = waveform + SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) + + + def bInitial(self, prob): + return Zero() + + def bInitialDeriv(self, prob, v=None, adjoint=False): + return Zero() + + def eInitial(self, prob): + return Zero() + + def eInitialDeriv(self, prob, v=None, adjoint=False): + return Zero() + + def eval(self, prob, time): + S_m = self.S_m(prob, time) + S_e = self.S_e(prob, time) + return S_m, S_e + + def evalDeriv(self, prob, time, v=None, adjoint=False): + if v is not None: + return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint) + else: + return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint) + + def S_m(self, prob, time): + return Zero() + + def S_e(self, prob, time): + return Zero() + + def S_mDeriv(self, prob, time, v=None, adjoint=False): + return Zero() + + def S_eDeriv(self, prob, time, v=None, adjoint=False): + return Zero() + + +class MagDipole(BaseSrc): + + waveform = None + loc = None + orientation = 'Z' + moment = 1. + mu = mu_0 + + def __init__(self, rxList, **kwargs): + assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." + self.integrate = False + BaseSrc.__init__(self, rxList, **kwargs) + + def _bfromVectorPotential(self, prob): + if prob._eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif prob._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: + 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 bInitial(self, prob): + + if self.waveform.hasInitialFields is False: + return Zero() + + return self._bfromVectorPotential(prob) + + def eInitial(self, prob): + + if self.waveform.hasInitialFields is False: + return Zero() + + b = self.bInitial(prob) + MeSigmaI = prob.MeSigmaI + MfMui = prob.MfMui + C = prob.mesh.edgeCurl + + return MeSigmaI * (C.T * (MfMui * b)) + + def eInitialDeriv(self, prob, v=None, adjoint=False): + + if self.waveform.hasInitialFields is False: + return Zero() + + b = self.bInitial(prob) + MeSigmaIDeriv = prob.MeSigmaIDeriv + MfMui = prob.MfMui + C = prob.mesh.edgeCurl + S_e = self.S_e(prob, prob.t0) + + # S_e doesn't depend on the model + + if adjoint: + return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v + + return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v + + + def S_m(self, prob, time): + if self.waveform.hasInitialFields is False: + raise NotImplementedError + return Zero() + + def S_e(self, prob, time): + if self.waveform.hasInitialFields is False: + raise NotImplementedError + return Zero() diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index eb534978..6acf3c04 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -3,6 +3,8 @@ from SimPEG import np, Utils from SimPEG.Utils import Zero, Identity from scipy.constants import mu_0 from SimPEG.EM.Utils import * +import SrcTDEM as Src + #################################################### # Receivers @@ -70,204 +72,6 @@ class Rx(SimPEG.Survey.BaseTimeRx): # newshape = (len(dP_dF_T)/timeMesh.nN, timeMesh.nN ) return P.T * v #np.reshape(dP_dF_T, newshape, order='F') -#################################################### -# Sources -#################################################### - -class BaseWaveform(object): - - def __init__(self, offTime=0., hasInitialFields=False): - self.offTime = offTime - self.hasInitialFields = hasInitialFields - - def _assertMatchesPair(self, pair): - assert (isinstance(self, pair) - ), "Waveform object must be an instance of a %s BaseWaveform class."%(pair.__name__) - - def eval(self, time): - raise NotImplementedError - - def evalDeriv(self, time): - raise NotImplementedError # needed for E-formulation - - -class StepOffWaveform(BaseWaveform): - - def __init__(self, offTime=0.): - BaseWaveform.__init__(self, offTime, hasInitialFields=True) - - def eval(self, time): - return 0. - - -class RawWaveform(BaseWaveform): - - def __init__(self, offTime=0.): - BaseWaveform.__init__(self, offTime, hasInitialFields=True) - - def eval(self, time): - raise NotImplementedError('RawWaveform has not been implemented, you should write it!') - - -class TriangularWaveform(BaseWaveform): - - def __init__(self, offTime=0.): - BaseWaveform.__init__(self, offTime, hasInitialFields=True) - - def eval(self, time): - raise NotImplementedError('TriangularWaveform has not been implemented, you should write it!') - - - -class BaseSrc(SimPEG.Survey.BaseSrc): - - rxPair = Rx - integrate = True - waveformPair = BaseWaveform - - @property - def waveform(self): - "A waveform instance is not None" - return getattr(self, '_waveform', None) - @waveform.setter - def waveform(self, val): - if self.waveform is None: - val._assertMatchesPair(self.waveformPair) - self._mapping = val - else: - self._mapping = self.PropMap(val) - - - def __init__(self, rxList, waveform = StepOffWaveform(), **kwargs): - self.waveform = waveform - SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs) - - - def bInitial(self, prob): - return Zero() - - def bInitialDeriv(self, prob, v=None, adjoint=False): - return Zero() - - def eInitial(self, prob): - return Zero() - - def eInitialDeriv(self, prob, v=None, adjoint=False): - return Zero() - - def eval(self, prob, time): - S_m = self.S_m(prob, time) - S_e = self.S_e(prob, time) - return S_m, S_e - - def evalDeriv(self, prob, time, v=None, adjoint=False): - if v is not None: - return self.S_mDeriv(prob, time, v, adjoint), self.S_eDeriv(prob, time, v, adjoint) - else: - return lambda v: self.S_mDeriv(prob, time, v, adjoint), lambda v: self.S_eDeriv(prob, time, v, adjoint) - - def S_m(self, prob, time): - return Zero() - - def S_e(self, prob, time): - return Zero() - - def S_mDeriv(self, prob, time, v=None, adjoint=False): - return Zero() - - def S_eDeriv(self, prob, time, v=None, adjoint=False): - return Zero() - - -class MagDipole(BaseSrc): - - waveform = None - loc = None - orientation = 'Z' - moment = 1. - mu = mu_0 - - def __init__(self, rxList, **kwargs): - assert self.orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..." - self.integrate = False - BaseSrc.__init__(self, rxList, **kwargs) - - def _bfromVectorPotential(self, prob): - if prob._eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif prob._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: - 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 bInitial(self, prob): - - if self.waveform.hasInitialFields is False: - return Zero() - - return self._bfromVectorPotential(prob) - - def eInitial(self, prob): - - if self.waveform.hasInitialFields is False: - return Zero() - - b = self.bInitial(prob) - MeSigmaI = prob.MeSigmaI - MfMui = prob.MfMui - C = prob.mesh.edgeCurl - - return MeSigmaI * (C.T * (MfMui * b)) - - def eInitialDeriv(self, prob, v=None, adjoint=False): - - if self.waveform.hasInitialFields is False: - return Zero() - - b = self.bInitial(prob) - MeSigmaIDeriv = prob.MeSigmaIDeriv - MfMui = prob.MfMui - C = prob.mesh.edgeCurl - S_e = self.S_e(prob, prob.t0) - - # S_e doesn't depend on the model - - if adjoint: - return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ).T * v - - return MeSigmaIDeriv( -S_e + C.T * ( MfMui * b ) ) * v - - - def S_m(self, prob, time): - if self.waveform.hasInitialFields is False: - raise NotImplementedError - return Zero() - - def S_e(self, prob, time): - if self.waveform.hasInitialFields is False: - raise NotImplementedError - return Zero() #################################################### # Survey @@ -278,7 +82,7 @@ class Survey(SimPEG.Survey.BaseSurvey): Time domain electromagnetic survey """ - srcPair = BaseSrc + srcPair = Src.BaseSrc rxPair = Rx def __init__(self, srcList, **kwargs): diff --git a/SimPEG/EM/TDEM/__init__.py b/SimPEG/EM/TDEM/__init__.py index 13010109..3c22d42d 100644 --- a/SimPEG/EM/TDEM/__init__.py +++ b/SimPEG/EM/TDEM/__init__.py @@ -1,3 +1,3 @@ from TDEM import BaseTDEMProblem, Problem_b, Problem_e from FieldsTDEM import Fields, Fields_b -from SurveyTDEM import Survey, BaseSrc, Rx +from SurveyTDEM import Survey, Src, Rx diff --git a/SimPEG/EM/TDEM_old/SurveyTDEM.py b/SimPEG/EM/TDEM_old/SurveyTDEM.py index 338fd484..6b5aec2b 100644 --- a/SimPEG/EM/TDEM_old/SurveyTDEM.py +++ b/SimPEG/EM/TDEM_old/SurveyTDEM.py @@ -2,7 +2,7 @@ from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey from SimPEG.EM.Utils import * from BaseTDEM import FieldsTDEM - +import SrcTDEM as Src class RxTDEM(Survey.BaseTimeRx): @@ -87,7 +87,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): def getInitialFields(self, mesh): """Vertical magnetic dipole, magnetic vector potential""" if self.waveformType == "STEPOFF": - print ">> Step waveform: Non-zero initial condition" + print ">> Step waveform: Non-zero initial condition" if mesh._meshType is 'CYL': if mesh.isSymmetric: MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') @@ -96,8 +96,8 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') - return {"b": mesh.edgeCurl*MVP} + raise Exception('Unknown mesh for VMD') + return {"b": mesh.edgeCurl*MVP} elif self.waveformType == "GENERAL": print ">> General waveform: Zero initial condition" return {"b": np.zeros(mesh.nF)} @@ -113,7 +113,7 @@ class SrcTDEM_VMD_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: - raise Exception('Unknown mesh for VMD') + raise Exception('Unknown mesh for VMD') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP @@ -122,7 +122,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): self.loc = loc self.radius = radius self.waveformType = waveformType - SrcTDEM.__init__(self,rxList) + SrcTDEM.__init__(self,rxList) def getInitialFields(self, mesh): """Circular Loop, magnetic vector potential""" @@ -153,7 +153,7 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): elif mesh._meshType is 'TENSOR': MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: - raise Exception('Unknown mesh for CircularLoop') + raise Exception('Unknown mesh for CircularLoop') return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP diff --git a/SimPEG/Examples/EM_FDEM_1D_Inversion.py b/SimPEG/Examples/EM_FDEM_1D_Inversion.py index e76b2439..1556bf63 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) diff --git a/SimPEG/Examples/EM_TDEM_1D_Inversion.py b/SimPEG/Examples/EM_TDEM_1D_Inversion.py index c92b972d..71ea0619 100644 --- a/SimPEG/Examples/EM_TDEM_1D_Inversion.py +++ b/SimPEG/Examples/EM_TDEM_1D_Inversion.py @@ -43,7 +43,7 @@ def run(plotIt=True): rxOffset=1e-3 rx = EM.TDEM.Rx(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz') - src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 80])) + src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 80])) survey = EM.TDEM.Survey([src]) prb = EM.TDEM.Problem_b(mesh, mapping=mapping) diff --git a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py index c3ac40ba..085b4f1a 100644 --- a/tests/em/tdem/test_TDEM_b_DerivAdjoint.py +++ b/tests/em/tdem/test_TDEM_b_DerivAdjoint.py @@ -24,7 +24,7 @@ def setUp(prbtype='b', rxcomp='bz'): rxOffset = 10. rx = EM.TDEM.Rx(np.array([[rxOffset, 0., -1e-2]]), np.logspace(-4,-3, 20), rxcomp) #,] - src = EM.TDEM.SurveyTDEM.MagDipole([rx], loc=np.array([0., 0., 0.])) + src = EM.TDEM.Src.MagDipole([rx], loc=np.array([0., 0., 0.])) survey = EM.TDEM.Survey([src])