From 25e21608df0cdf7fc11e974e44ffdd20bb858e08 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 08:34:46 -0800 Subject: [PATCH 01/14] import all EM utils under Utils namespace --- SimPEG/EM/Utils/__init__.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index 6e430cf9..487bd909 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,5 @@ # import Sources # import Ana # import Solver -import EMUtils -import SrcUtils \ No newline at end of file +from EMUtils import * +from SrcUtils import * \ No newline at end of file From 3b62957d08b08f9e43f75ff557536f9d417cd986 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 08:44:34 -0800 Subject: [PATCH 02/14] re-named SrcUtils to AnalyticUtils, cleaned up namespace for imports to close #159 --- SimPEG/EM/FDEM/FDEM.py | 2 +- SimPEG/EM/FDEM/FieldsFDEM.py | 2 +- SimPEG/EM/FDEM/SurveyFDEM.py | 13 ++++++------- SimPEG/EM/Utils/{SrcUtils.py => AnalyticUtils.py} | 0 SimPEG/EM/Utils/__init__.py | 4 ++-- 5 files changed, 10 insertions(+), 11 deletions(-) rename SimPEG/EM/Utils/{SrcUtils.py => AnalyticUtils.py} (100%) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 6d3df44f..9d3be6b2 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -3,7 +3,7 @@ from scipy.constants import mu_0 from SurveyFDEM import SurveyFDEM from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j from SimPEG.EM.Base import BaseEMProblem -from SimPEG.EM.Utils.EMUtils import omega +from SimPEG.EM.Utils import omega class BaseFDEMProblem(BaseEMProblem): diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index bb786bd1..53f90c4b 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -1,5 +1,5 @@ from SimPEG import Survey, Problem, Utils, np, sp -from SimPEG.EM.Utils.EMUtils import omega +from SimPEG.EM.Utils import omega class FieldsFDEM(Problem.Fields): diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index dbda9f80..94cac712 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,6 +1,5 @@ from SimPEG import Survey, Problem, Utils, np, sp -from SimPEG.EM.Utils import SrcUtils -from SimPEG.EM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from SimPEG.EM.Utils import * from scipy.constants import mu_0 #################################################### @@ -258,10 +257,10 @@ class SrcFDEM_MagDipole(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) else: - srcfct = SrcUtils.MagneticDipoleVectorPotential + 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) @@ -322,7 +321,7 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): gridZ = prob.mesh.gridEz C = prob.mesh.edgeCurl.T - srcfct = SrcUtils.MagneticDipoleFields + srcfct = MagneticDipoleFields if prob.mesh._meshType is 'CYL': if not prob.mesh.isSymmetric: # TODO ? @@ -395,10 +394,10 @@ class SrcFDEM_CircularLoop(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) else: - srcfct = SrcUtils.MagneticDipoleVectorPotential + srcfct = MagneticDipoleVectorPotential ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) diff --git a/SimPEG/EM/Utils/SrcUtils.py b/SimPEG/EM/Utils/AnalyticUtils.py similarity index 100% rename from SimPEG/EM/Utils/SrcUtils.py rename to SimPEG/EM/Utils/AnalyticUtils.py diff --git a/SimPEG/EM/Utils/__init__.py b/SimPEG/EM/Utils/__init__.py index 487bd909..18dddde9 100644 --- a/SimPEG/EM/Utils/__init__.py +++ b/SimPEG/EM/Utils/__init__.py @@ -1,5 +1,5 @@ # import Sources # import Ana # import Solver -from EMUtils import * -from SrcUtils import * \ No newline at end of file +from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b +from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential \ No newline at end of file From 19bcdbfbf5361f6d9cb48199a8b086349b230e0e Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 10:39:08 -0800 Subject: [PATCH 03/14] EM.FDEM.SrcFDEM_XXX --> EM.FDEM.Src.XXX --- SimPEG/EM/FDEM/SrcFDEM.py | 348 ++++++++++++++++++ SimPEG/EM/FDEM/SurveyFDEM.py | 621 ++++++++++++++++---------------- SimPEG/EM/FDEM/__init__.py | 2 +- SimPEG/EM/Utils/testingUtils.py | 10 +- 4 files changed, 665 insertions(+), 316 deletions(-) create mode 100644 SimPEG/EM/FDEM/SrcFDEM.py diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py new file mode 100644 index 00000000..f7b95925 --- /dev/null +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -0,0 +1,348 @@ +from SimPEG import Survey, Problem, Utils, np, sp +# import SimPEG.EM as EM +from SimPEG.EM.Utils import * +from scipy.constants import mu_0 +# from SurveyFDEM import RxFDEM + + +class BaseSrcFDEM(Survey.BaseSrc): + freq = None + # rxPair = EM.FDEM.RxFDEM + integrate = True + + def eval(self, prob): + S_m = self.S_m(prob) + S_e = self.S_e(prob) + return S_m, S_e + + def evalDeriv(self, prob, v, adjoint=False): + return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) + + def bPrimary(self, prob): + return None + + def hPrimary(self, prob): + return None + + def ePrimary(self, prob): + return None + + def jPrimary(self, prob): + return None + + def S_m(self, prob): + return None + + def S_e(self, prob): + return None + + def S_mDeriv(self, prob, v, adjoint = False): + return None + + def S_eDeriv(self, prob, v, adjoint = False): + return None + + +class RawVec_e(BaseSrcFDEM): + """ + RawVec electric source. It is defined by the user provided vector S_e + + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_e = np.array(S_e,dtype=complex) + self._ePrimary = ePrimary + self._bPrimary = bPrimary + self._hPrimary = hPrimary + self._jPrimary = jPrimary + self.freq = float(freq) + BaseSrcFDEM.__init__(self, rxList) + + def S_e(self, prob): + return self._S_e + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec_m(BaseSrcFDEM): + """ + RawVec magnetic source. It is defined by the user provided vector S_m + + :param numpy.array S_m: magnetic source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): + self._S_m = np.array(S_m,dtype=complex) + self.freq = float(freq) + self.integrate = integrate + self._ePrimary = np.array(ePrimary,dtype=complex) + self._bPrimary = np.array(bPrimary,dtype=complex) + self._hPrimary = np.array(hPrimary,dtype=complex) + self._jPrimary = np.array(jPrimary,dtype=complex) + + BaseSrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + return self._S_m + + def ePrimary(self, prob): + return self._ePrimary + + def bPrimary(self, prob): + return self._bPrimary + + def hPrimary(self, prob): + return self._hPrimary + + def jPrimary(self, prob): + return self._jPrimary + + +class RawVec(BaseSrcFDEM): + """ + RawVec source. It is defined by the user provided vectors S_m, S_e + + :param numpy.array S_m: magnetic source term + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + def __init__(self, rxList, freq, S_m, S_e, integrate = True): + self._S_m = np.array(S_m,dtype=complex) + self._S_e = np.array(S_e,dtype=complex) + self.freq = float(freq) + self.integrate = integrate + BaseSrcFDEM.__init__(self, rxList) + + def S_m(self, prob): + if prob._eqLocs is 'EF' and self.integrate is True: + return prob.Me * self._S_m + return self._S_m + + def S_e(self, prob): + if prob._eqLocs is 'FE' and self.integrate is True: + return prob.Me * self._S_e + return self._S_e + + +class MagDipole(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + self.integrate = False + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob,b) + + def S_m(self, prob): + b_p = self.bPrimary(prob) + return -1j*omega(self.freq)*b_p + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class MagDipole_Bfield(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + #TODO: neither does moment + def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): + self.freq = float(freq) + self.loc = loc + self.orientation = orientation + self.moment = moment + self.mu = mu + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl.T + + srcfct = MagneticDipoleFields + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,bz)) + else: + bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) + by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) + b = np.concatenate((bx,by,bz)) + + return b + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return h_from_b(prob, b) + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + +class CircularLoop(BaseSrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): + self.freq = float(freq) + self.orientation = orientation + self.radius = radius + self.mu = mu + self.loc = loc + self.integrate = False + BaseSrcFDEM.__init__(self, rxList) + + def bPrimary(self, prob): + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + gridX = prob.mesh.gridEx + gridY = prob.mesh.gridEy + gridZ = prob.mesh.gridEz + C = prob.mesh.edgeCurl + + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz + C = prob.mesh.edgeCurl.T + + if prob.mesh._meshType is 'CYL': + if not prob.mesh.isSymmetric: + # TODO ? + raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') + a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + + else: + srcfct = MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) + ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) + az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) + a = np.concatenate((ax, ay, az)) + + return C*a + + def hPrimary(self, prob): + b = self.bPrimary(prob) + return 1./self.mu*b + + def S_m(self, prob): + b = self.bPrimary(prob) + return -1j*omega(self.freq)*b + + def S_e(self, prob): + if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): + return None + else: + eqLocs = prob._eqLocs + + if eqLocs is 'FE': + mui_s = prob.curModel.mui - 1./self.mu + MMui_s = prob.mesh.getFaceInnerProduct(mui_s) + C = prob.mesh.edgeCurl + elif eqLocs is 'EF': + mu_s = prob.curModel.mu - self.mu + MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) + C = prob.mesh.edgeCurl.T + + return -C.T * (MMui_s * self.bPrimary(prob)) + + diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 94cac712..1cdc104c 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,6 +1,7 @@ from SimPEG import Survey, Problem, Utils, np, sp from SimPEG.EM.Utils import * from scipy.constants import mu_0 +import SrcFDEM as Src #################################################### # Receivers @@ -90,345 +91,345 @@ class RxFDEM(Survey.BaseRx): # Sources #################################################### -class SrcFDEM(Survey.BaseSrc): - freq = None - rxPair = RxFDEM - integrate = True +# class SrcFDEM(Survey.BaseSrc): +# freq = None +# rxPair = RxFDEM +# integrate = True - def eval(self, prob): - S_m = self.S_m(prob) - S_e = self.S_e(prob) - return S_m, S_e +# def eval(self, prob): +# S_m = self.S_m(prob) +# S_e = self.S_e(prob) +# return S_m, S_e - def evalDeriv(self, prob, v, adjoint=False): - return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) +# def evalDeriv(self, prob, v, adjoint=False): +# return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) - def bPrimary(self, prob): - return None +# def bPrimary(self, prob): +# return None - def hPrimary(self, prob): - return None +# def hPrimary(self, prob): +# return None - def ePrimary(self, prob): - return None +# def ePrimary(self, prob): +# return None - def jPrimary(self, prob): - return None +# def jPrimary(self, prob): +# return None - def S_m(self, prob): - return None +# def S_m(self, prob): +# return None - def S_e(self, prob): - return None +# def S_e(self, prob): +# return None - def S_mDeriv(self, prob, v, adjoint = False): - return None +# def S_mDeriv(self, prob, v, adjoint = False): +# return None - def S_eDeriv(self, prob, v, adjoint = False): - return None +# def S_eDeriv(self, prob, v, adjoint = False): +# return None -class SrcFDEM_RawVec_e(SrcFDEM): - """ - RawVec electric source. It is defined by the user provided vector S_e +# class SrcFDEM_RawVec_e(SrcFDEM): +# """ +# RawVec electric source. It is defined by the user provided vector S_e - :param numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list - """ +# :param numpy.array S_e: electric source term +# :param float freq: frequency +# :param rxList: receiver list +# """ - def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_e = np.array(S_e,dtype=complex) - self._ePrimary = ePrimary - self._bPrimary = bPrimary - self._hPrimary = hPrimary - self._jPrimary = jPrimary - self.freq = float(freq) - SrcFDEM.__init__(self, rxList) +# def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): +# self._S_e = np.array(S_e,dtype=complex) +# self._ePrimary = ePrimary +# self._bPrimary = bPrimary +# self._hPrimary = hPrimary +# self._jPrimary = jPrimary +# self.freq = float(freq) +# SrcFDEM.__init__(self, rxList) - def S_e(self, prob): - return self._S_e +# def S_e(self, prob): +# return self._S_e - def ePrimary(self, prob): - return self._ePrimary +# def ePrimary(self, prob): +# return self._ePrimary - def bPrimary(self, prob): - return self._bPrimary +# def bPrimary(self, prob): +# return self._bPrimary - def hPrimary(self, prob): - return self._hPrimary +# def hPrimary(self, prob): +# return self._hPrimary - def jPrimary(self, prob): - return self._jPrimary +# def jPrimary(self, prob): +# return self._jPrimary -class SrcFDEM_RawVec_m(SrcFDEM): - """ - RawVec magnetic source. It is defined by the user provided vector S_m +# class SrcFDEM_RawVec_m(SrcFDEM): +# """ +# RawVec magnetic source. It is defined by the user provided vector S_m - :param numpy.array S_m: magnetic source term - :param float freq: frequency - :param rxList: receiver list - """ +# :param numpy.array S_m: magnetic source term +# :param float freq: frequency +# :param rxList: receiver list +# """ - def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): - self._S_m = np.array(S_m,dtype=complex) - self.freq = float(freq) - self.integrate = integrate - self._ePrimary = np.array(ePrimary,dtype=complex) - self._bPrimary = np.array(bPrimary,dtype=complex) - self._hPrimary = np.array(hPrimary,dtype=complex) - self._jPrimary = np.array(jPrimary,dtype=complex) +# def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): +# self._S_m = np.array(S_m,dtype=complex) +# self.freq = float(freq) +# self.integrate = integrate +# self._ePrimary = np.array(ePrimary,dtype=complex) +# self._bPrimary = np.array(bPrimary,dtype=complex) +# self._hPrimary = np.array(hPrimary,dtype=complex) +# self._jPrimary = np.array(jPrimary,dtype=complex) - SrcFDEM.__init__(self, rxList) +# SrcFDEM.__init__(self, rxList) - def S_m(self, prob): - return self._S_m +# def S_m(self, prob): +# return self._S_m - def ePrimary(self, prob): - return self._ePrimary - - def bPrimary(self, prob): - return self._bPrimary - - def hPrimary(self, prob): - return self._hPrimary +# def ePrimary(self, prob): +# return self._ePrimary + +# def bPrimary(self, prob): +# return self._bPrimary + +# def hPrimary(self, prob): +# return self._hPrimary - def jPrimary(self, prob): - return self._jPrimary +# def jPrimary(self, prob): +# return self._jPrimary -class SrcFDEM_RawVec(SrcFDEM): - """ - RawVec source. It is defined by the user provided vectors S_m, S_e +# class SrcFDEM_RawVec(SrcFDEM): +# """ +# RawVec source. It is defined by the user provided vectors S_m, S_e - :param numpy.array S_m: magnetic source term - :param numpy.array S_e: electric source term - :param float freq: frequency - :param rxList: receiver list - """ - def __init__(self, rxList, freq, S_m, S_e, integrate = True): - self._S_m = np.array(S_m,dtype=complex) - self._S_e = np.array(S_e,dtype=complex) - self.freq = float(freq) - self.integrate = integrate - SrcFDEM.__init__(self, rxList) - - def S_m(self, prob): - if prob._eqLocs is 'EF' and self.integrate is True: - return prob.Me * self._S_m - return self._S_m - - def S_e(self, prob): - if prob._eqLocs is 'FE' and self.integrate is True: - return prob.Me * self._S_e - return self._S_e - - -class SrcFDEM_MagDipole(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - - else: - srcfct = MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob,b) - - def S_m(self, prob): - b_p = self.bPrimary(prob) - return -1j*omega(self.freq)*b_p - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_MagDipole_Bfield(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - #TODO: neither does moment - def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): - self.freq = float(freq) - self.loc = loc - self.orientation = orientation - self.moment = moment - self.mu = mu - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl.T - - srcfct = MagneticDipoleFields - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,bz)) - else: - bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) - by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) - b = np.concatenate((bx,by,bz)) - - return b - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return h_from_b(prob, b) - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) - - -class SrcFDEM_CircularLoop(SrcFDEM): - - #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that - def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): - self.freq = float(freq) - self.orientation = orientation - self.radius = radius - self.mu = mu - self.loc = loc - self.integrate = False - SrcFDEM.__init__(self, rxList) - - def bPrimary(self, prob): - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - gridX = prob.mesh.gridEx - gridY = prob.mesh.gridEy - gridZ = prob.mesh.gridEz - C = prob.mesh.edgeCurl - - elif eqLocs is 'EF': - gridX = prob.mesh.gridFx - gridY = prob.mesh.gridFy - gridZ = prob.mesh.gridFz - C = prob.mesh.edgeCurl.T - - if prob.mesh._meshType is 'CYL': - if not prob.mesh.isSymmetric: - # TODO ? - raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) - - else: - srcfct = MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) - ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) - az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) - a = np.concatenate((ax, ay, az)) - - return C*a - - def hPrimary(self, prob): - b = self.bPrimary(prob) - return 1./self.mu*b - - def S_m(self, prob): - b = self.bPrimary(prob) - return -1j*omega(self.freq)*b - - def S_e(self, prob): - if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): - return None - else: - eqLocs = prob._eqLocs - - if eqLocs is 'FE': - mui_s = prob.curModel.mui - 1./self.mu - MMui_s = prob.mesh.getFaceInnerProduct(mui_s) - C = prob.mesh.edgeCurl - elif eqLocs is 'EF': - mu_s = prob.curModel.mu - self.mu - MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) - C = prob.mesh.edgeCurl.T - - return -C.T * (MMui_s * self.bPrimary(prob)) +# :param numpy.array S_m: magnetic source term +# :param numpy.array S_e: electric source term +# :param float freq: frequency +# :param rxList: receiver list +# """ +# def __init__(self, rxList, freq, S_m, S_e, integrate = True): +# self._S_m = np.array(S_m,dtype=complex) +# self._S_e = np.array(S_e,dtype=complex) +# self.freq = float(freq) +# self.integrate = integrate +# SrcFDEM.__init__(self, rxList) + +# def S_m(self, prob): +# if prob._eqLocs is 'EF' and self.integrate is True: +# return prob.Me * self._S_m +# return self._S_m + +# def S_e(self, prob): +# if prob._eqLocs is 'FE' and self.integrate is True: +# return prob.Me * self._S_e +# return self._S_e + + +# class SrcFDEM_MagDipole(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): +# self.freq = float(freq) +# self.loc = loc +# self.orientation = orientation +# self.moment = moment +# self.mu = mu +# self.integrate = False +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl.T + + +# if prob.mesh._meshType is 'CYL': +# if not prob.mesh.isSymmetric: +# # TODO ? +# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') +# a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) + +# else: +# srcfct = MagneticDipoleVectorPotential +# ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) +# az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# a = np.concatenate((ax, ay, az)) + +# return C*a + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return h_from_b(prob,b) + +# def S_m(self, prob): +# b_p = self.bPrimary(prob) +# return -1j*omega(self.freq)*b_p + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) + + +# class SrcFDEM_MagDipole_Bfield(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# #TODO: neither does moment +# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): +# self.freq = float(freq) +# self.loc = loc +# self.orientation = orientation +# self.moment = moment +# self.mu = mu +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl.T + +# srcfct = MagneticDipoleFields +# if prob.mesh._meshType is 'CYL': +# if not prob.mesh.isSymmetric: +# # TODO ? +# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') +# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# b = np.concatenate((bx,bz)) +# else: +# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) +# by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) +# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) +# b = np.concatenate((bx,by,bz)) + +# return b + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return h_from_b(prob, b) + +# def S_m(self, prob): +# b = self.bPrimary(prob) +# return -1j*omega(self.freq)*b + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) + + +# class SrcFDEM_CircularLoop(SrcFDEM): + +# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that +# def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): +# self.freq = float(freq) +# self.orientation = orientation +# self.radius = radius +# self.mu = mu +# self.loc = loc +# self.integrate = False +# SrcFDEM.__init__(self, rxList) + +# def bPrimary(self, prob): +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# gridX = prob.mesh.gridEx +# gridY = prob.mesh.gridEy +# gridZ = prob.mesh.gridEz +# C = prob.mesh.edgeCurl + +# elif eqLocs is 'EF': +# gridX = prob.mesh.gridFx +# gridY = prob.mesh.gridFy +# gridZ = prob.mesh.gridFz +# C = prob.mesh.edgeCurl.T + +# if prob.mesh._meshType is 'CYL': +# if not prob.mesh.isSymmetric: +# # TODO ? +# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') +# a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) + +# else: +# srcfct = MagneticDipoleVectorPotential +# ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) +# ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) +# az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) +# a = np.concatenate((ax, ay, az)) + +# return C*a + +# def hPrimary(self, prob): +# b = self.bPrimary(prob) +# return 1./self.mu*b + +# def S_m(self, prob): +# b = self.bPrimary(prob) +# return -1j*omega(self.freq)*b + +# def S_e(self, prob): +# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): +# return None +# else: +# eqLocs = prob._eqLocs + +# if eqLocs is 'FE': +# mui_s = prob.curModel.mui - 1./self.mu +# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) +# C = prob.mesh.edgeCurl +# elif eqLocs is 'EF': +# mu_s = prob.curModel.mu - self.mu +# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) +# C = prob.mesh.edgeCurl.T + +# return -C.T * (MMui_s * self.bPrimary(prob)) #################################################### @@ -440,7 +441,7 @@ class SurveyFDEM(Survey.BaseSurvey): docstring for SurveyFDEM """ - srcPair = SrcFDEM + srcPair = Src.BaseSrcFDEM def __init__(self, srcList, **kwargs): # Sort these by frequency diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 110b4d1e..9a262ffb 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyFDEM import * +from SurveyFDEM import RxFDEM, Src, SurveyFDEM from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h from FieldsFDEM import * \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index b4570c0d..cb1e1d26 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -23,25 +23,25 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): for SrcType in SrcList: if SrcType is 'MagDipole': - Src.append(EM.FDEM.SrcFDEM_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.SrcFDEM_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.SrcFDEM_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])] = 1. S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) elif fdemType is 'h' or fdemType is 'j': S_m = np.zeros(mesh.nE) S_e = np.zeros(mesh.nF) S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src.append(EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e)) + Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e)) if verbose: print ' Fetching %s problem' % (fdemType) From aed5f5ad52be82943848dd965a63597844e682ea Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 12:15:32 -0800 Subject: [PATCH 04/14] fixed SrcUtils import in BaseTDEM --- SimPEG/EM/TDEM/BaseTDEM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/TDEM/BaseTDEM.py b/SimPEG/EM/TDEM/BaseTDEM.py index 5e82b7b2..e36d76b8 100644 --- a/SimPEG/EM/TDEM/BaseTDEM.py +++ b/SimPEG/EM/TDEM/BaseTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Solver, Problem from SimPEG.Problem import BaseTimeProblem -from SimPEG.EM.Utils import SrcUtils +from SimPEG.EM.Utils import * from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh From 356a5b103dbe311428ec2c44651133787099accf Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 12:20:07 -0800 Subject: [PATCH 05/14] and in SurveyTDEM --- SimPEG/EM/TDEM/SurveyTDEM.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 9878438d..0c87334f 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey -from SimPEG.EM.Utils import SrcUtils +from SimPEG.EM.Utils import * from BaseTDEM import FieldsTDEM From 2da82b7d190ecd500124283fe43b20537963e49a Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 12:24:12 -0800 Subject: [PATCH 06/14] fixed SrcUtils call in SurveyTDEM --- SimPEG/EM/TDEM/SurveyTDEM.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/TDEM/SurveyTDEM.py b/SimPEG/EM/TDEM/SurveyTDEM.py index 0c87334f..7f6e5c04 100644 --- a/SimPEG/EM/TDEM/SurveyTDEM.py +++ b/SimPEG/EM/TDEM/SurveyTDEM.py @@ -87,11 +87,11 @@ class SrcTDEM_VMD_MVP(SrcTDEM): """Vertical magnetic dipole, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') + MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey') else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = SrcUtils.MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) + MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez']) else: raise Exception('Unknown mesh for VMD') @@ -109,11 +109,11 @@ class SrcTDEM_CircularLoop_MVP(SrcTDEM): """Circular Loop, magnetic vector potential""" if mesh._meshType is 'CYL': if mesh.isSymmetric: - MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) + MVP = MagneticLoopVectorPotential(self.loc, mesh, 'Ey', self.radius) else: raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') elif mesh._meshType is 'TENSOR': - MVP = SrcUtils.MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) + MVP = MagneticLoopVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'], self.radius) else: raise Exception('Unknown mesh for CircularLoop') From 52375ecd0619e9417726a334d716ece8dd95f110 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 12:32:09 -0800 Subject: [PATCH 07/14] fixed fdem deriv test --- .../fdem/inverse/derivs/test_FDEM_derivs.py | 85 +------------------ 1 file changed, 2 insertions(+), 83 deletions(-) diff --git a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py index 32b171c9..52108c4e 100644 --- a/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py +++ b/tests/em/fdem/inverse/derivs/test_FDEM_derivs.py @@ -3,6 +3,7 @@ from SimPEG import * from SimPEG import EM import sys from scipy.constants import mu_0 +from SimPEG.EM.Utils.testingUtils import getFDEMProblem testDerivs = True testEB = True @@ -20,91 +21,9 @@ addrandoms = True SrcType = 'RawVec' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec' -def getProblem(fdemType, comp): - cs = 5. - ncx, ncy, ncz = 6, 6, 6 - npad = 3 - hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)] - hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)] - hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)] - mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C']) - - mapping = Maps.ExpMap(mesh) - - x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source - XYZ = Utils.ndgrid(x,x,np.r_[0.]) - Rx0 = EM.FDEM.RxFDEM(XYZ, comp) - - if SrcType is 'MagDipole': - Src = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - elif SrcType is 'MagDipole_Bfield': - Src = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - elif SrcType is 'CircularLoop': - Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq, loc=np.r_[0.,0.,0.]) - - if verbose: - print ' Fetching %s problem' % (fdemType) - - if fdemType == 'e': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nF) - S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) - - elif fdemType == 'b': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nF) - S_e = np.zeros(mesh.nE) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) - - elif fdemType == 'j': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nE) - S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) - - elif fdemType == 'h': - if SrcType is 'RawVec': - S_m = np.zeros(mesh.nE) - S_e = np.zeros(mesh.nF) - S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1. - S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1. - Src = EM.FDEM.SrcFDEM_RawVec([Rx0], freq, S_m, S_e) - - survey = EM.FDEM.SurveyFDEM([Src]) - prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) - - else: - raise NotImplementedError() - prb.pair(survey) - - try: - from pymatsolver import MumpsSolver - prb.Solver = MumpsSolver - except ImportError, e: - pass - - return prb - - def derivTest(fdemType, comp): - prb = getProblem(fdemType, comp) + prb = getFDEMProblem(fdemType, comp, [SrcType], freq) print '%s formulation - %s' % (fdemType, comp) x0 = np.log(np.ones(prb.mapping.nP)*CONDUCTIVITY) mu = np.log(np.ones(prb.mesh.nC)*MU) From 9d92f62cfaf28fa4844517954606533f60a8f882 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Fri, 13 Nov 2015 13:50:10 -0800 Subject: [PATCH 08/14] analytic test update --- tests/em/fdem/forward/test_FDEM_analytics.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index cb906262..0f9d64f8 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -29,7 +29,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.RxFDEM(XYZ, 'exi') - Src0 = EM.FDEM.SrcFDEM_MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) + Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.SurveyFDEM([Src0]) @@ -114,9 +114,9 @@ class FDEM_analyticTests(unittest.TestCase): de = np.zeros(mesh.nF,dtype=complex) de[s_ind] = 1./csz - de_p = [EM.FDEM.SrcFDEM_RawVec_e([],freq,de/mesh.area)] + de_p = [EM.FDEM.Src.RawVec_e([],freq,de/mesh.area)] - dm_p = [EM.FDEM.SrcFDEM_MagDipole([],freq,src_loc)] + dm_p = [EM.FDEM.Src.MagDipole([],freq,src_loc)] # Pair the problem and survey From c041ad3de6dea293f624951ef02283238fb577fd Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 08:29:27 -0800 Subject: [PATCH 09/14] EM.FDEM.RxFDEM --> EM.FDEM.Rx --- SimPEG/EM/FDEM/SurveyFDEM.py | 347 +------------------------------- SimPEG/EM/FDEM/__init__.py | 2 +- SimPEG/EM/Utils/testingUtils.py | 2 +- 3 files changed, 3 insertions(+), 348 deletions(-) diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 1cdc104c..2af72d81 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -7,7 +7,7 @@ import SrcFDEM as Src # Receivers #################################################### -class RxFDEM(Survey.BaseRx): +class Rx(Survey.BaseRx): knownRxTypes = { 'exr':['e', 'Ex', 'real'], @@ -87,351 +87,6 @@ class RxFDEM(Survey.BaseRx): return Pv -#################################################### -# Sources -#################################################### - -# class SrcFDEM(Survey.BaseSrc): -# freq = None -# rxPair = RxFDEM -# integrate = True - -# def eval(self, prob): -# S_m = self.S_m(prob) -# S_e = self.S_e(prob) -# return S_m, S_e - -# def evalDeriv(self, prob, v, adjoint=False): -# return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint) - -# def bPrimary(self, prob): -# return None - -# def hPrimary(self, prob): -# return None - -# def ePrimary(self, prob): -# return None - -# def jPrimary(self, prob): -# return None - -# def S_m(self, prob): -# return None - -# def S_e(self, prob): -# return None - -# def S_mDeriv(self, prob, v, adjoint = False): -# return None - -# def S_eDeriv(self, prob, v, adjoint = False): -# return None - - -# class SrcFDEM_RawVec_e(SrcFDEM): -# """ -# RawVec electric source. It is defined by the user provided vector S_e - -# :param numpy.array S_e: electric source term -# :param float freq: frequency -# :param rxList: receiver list -# """ - -# def __init__(self, rxList, freq, S_e, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): -# self._S_e = np.array(S_e,dtype=complex) -# self._ePrimary = ePrimary -# self._bPrimary = bPrimary -# self._hPrimary = hPrimary -# self._jPrimary = jPrimary -# self.freq = float(freq) -# SrcFDEM.__init__(self, rxList) - -# def S_e(self, prob): -# return self._S_e - -# def ePrimary(self, prob): -# return self._ePrimary - -# def bPrimary(self, prob): -# return self._bPrimary - -# def hPrimary(self, prob): -# return self._hPrimary - -# def jPrimary(self, prob): -# return self._jPrimary - - -# class SrcFDEM_RawVec_m(SrcFDEM): -# """ -# RawVec magnetic source. It is defined by the user provided vector S_m - -# :param numpy.array S_m: magnetic source term -# :param float freq: frequency -# :param rxList: receiver list -# """ - -# def __init__(self, rxList, freq, S_m, integrate = True, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None): -# self._S_m = np.array(S_m,dtype=complex) -# self.freq = float(freq) -# self.integrate = integrate -# self._ePrimary = np.array(ePrimary,dtype=complex) -# self._bPrimary = np.array(bPrimary,dtype=complex) -# self._hPrimary = np.array(hPrimary,dtype=complex) -# self._jPrimary = np.array(jPrimary,dtype=complex) - -# SrcFDEM.__init__(self, rxList) - -# def S_m(self, prob): -# return self._S_m - -# def ePrimary(self, prob): -# return self._ePrimary - -# def bPrimary(self, prob): -# return self._bPrimary - -# def hPrimary(self, prob): -# return self._hPrimary - -# def jPrimary(self, prob): -# return self._jPrimary - - -# class SrcFDEM_RawVec(SrcFDEM): -# """ -# RawVec source. It is defined by the user provided vectors S_m, S_e - -# :param numpy.array S_m: magnetic source term -# :param numpy.array S_e: electric source term -# :param float freq: frequency -# :param rxList: receiver list -# """ -# def __init__(self, rxList, freq, S_m, S_e, integrate = True): -# self._S_m = np.array(S_m,dtype=complex) -# self._S_e = np.array(S_e,dtype=complex) -# self.freq = float(freq) -# self.integrate = integrate -# SrcFDEM.__init__(self, rxList) - -# def S_m(self, prob): -# if prob._eqLocs is 'EF' and self.integrate is True: -# return prob.Me * self._S_m -# return self._S_m - -# def S_e(self, prob): -# if prob._eqLocs is 'FE' and self.integrate is True: -# return prob.Me * self._S_e -# return self._S_e - - -# class SrcFDEM_MagDipole(SrcFDEM): - -# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that -# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): -# self.freq = float(freq) -# self.loc = loc -# self.orientation = orientation -# self.moment = moment -# self.mu = mu -# self.integrate = False -# SrcFDEM.__init__(self, rxList) - -# def bPrimary(self, prob): -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# gridX = prob.mesh.gridEx -# gridY = prob.mesh.gridEy -# gridZ = prob.mesh.gridEz -# C = prob.mesh.edgeCurl - -# elif eqLocs is 'EF': -# gridX = prob.mesh.gridFx -# gridY = prob.mesh.gridFy -# gridZ = prob.mesh.gridFz -# C = prob.mesh.edgeCurl.T - - -# if prob.mesh._meshType is 'CYL': -# if not prob.mesh.isSymmetric: -# # TODO ? -# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') -# a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) - -# else: -# srcfct = MagneticDipoleVectorPotential -# ax = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) -# ay = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) -# az = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) -# a = np.concatenate((ax, ay, az)) - -# return C*a - -# def hPrimary(self, prob): -# b = self.bPrimary(prob) -# return h_from_b(prob,b) - -# def S_m(self, prob): -# b_p = self.bPrimary(prob) -# return -1j*omega(self.freq)*b_p - -# def S_e(self, prob): -# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): -# return None -# else: -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# mui_s = prob.curModel.mui - 1./self.mu -# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) -# C = prob.mesh.edgeCurl -# elif eqLocs is 'EF': -# mu_s = prob.curModel.mu - self.mu -# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) -# C = prob.mesh.edgeCurl.T - -# return -C.T * (MMui_s * self.bPrimary(prob)) - - -# class SrcFDEM_MagDipole_Bfield(SrcFDEM): - -# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that -# #TODO: neither does moment -# def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): -# self.freq = float(freq) -# self.loc = loc -# self.orientation = orientation -# self.moment = moment -# self.mu = mu -# SrcFDEM.__init__(self, rxList) - -# def bPrimary(self, prob): -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# gridX = prob.mesh.gridFx -# gridY = prob.mesh.gridFy -# gridZ = prob.mesh.gridFz -# C = prob.mesh.edgeCurl - -# elif eqLocs is 'EF': -# gridX = prob.mesh.gridEx -# gridY = prob.mesh.gridEy -# gridZ = prob.mesh.gridEz -# C = prob.mesh.edgeCurl.T - -# srcfct = MagneticDipoleFields -# if prob.mesh._meshType is 'CYL': -# if not prob.mesh.isSymmetric: -# # TODO ? -# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') -# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) -# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) -# b = np.concatenate((bx,bz)) -# else: -# bx = srcfct(self.loc, gridX, 'x', mu=self.mu, moment=self.moment) -# by = srcfct(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) -# bz = srcfct(self.loc, gridZ, 'z', mu=self.mu, moment=self.moment) -# b = np.concatenate((bx,by,bz)) - -# return b - -# def hPrimary(self, prob): -# b = self.bPrimary(prob) -# return h_from_b(prob, b) - -# def S_m(self, prob): -# b = self.bPrimary(prob) -# return -1j*omega(self.freq)*b - -# def S_e(self, prob): -# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): -# return None -# else: -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# mui_s = prob.curModel.mui - 1./self.mu -# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) -# C = prob.mesh.edgeCurl -# elif eqLocs is 'EF': -# mu_s = prob.curModel.mu - self.mu -# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) -# C = prob.mesh.edgeCurl.T - -# return -C.T * (MMui_s * self.bPrimary(prob)) - - -# class SrcFDEM_CircularLoop(SrcFDEM): - -# #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that -# def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): -# self.freq = float(freq) -# self.orientation = orientation -# self.radius = radius -# self.mu = mu -# self.loc = loc -# self.integrate = False -# SrcFDEM.__init__(self, rxList) - -# def bPrimary(self, prob): -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# gridX = prob.mesh.gridEx -# gridY = prob.mesh.gridEy -# gridZ = prob.mesh.gridEz -# C = prob.mesh.edgeCurl - -# elif eqLocs is 'EF': -# gridX = prob.mesh.gridFx -# gridY = prob.mesh.gridFy -# gridZ = prob.mesh.gridFz -# C = prob.mesh.edgeCurl.T - -# if prob.mesh._meshType is 'CYL': -# if not prob.mesh.isSymmetric: -# # TODO ? -# raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') -# a = MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) - -# else: -# srcfct = MagneticDipoleVectorPotential -# ax = srcfct(self.loc, gridX, 'x', self.radius, mu=self.mu) -# ay = srcfct(self.loc, gridY, 'y', self.radius, mu=self.mu) -# az = srcfct(self.loc, gridZ, 'z', self.radius, mu=self.mu) -# a = np.concatenate((ax, ay, az)) - -# return C*a - -# def hPrimary(self, prob): -# b = self.bPrimary(prob) -# return 1./self.mu*b - -# def S_m(self, prob): -# b = self.bPrimary(prob) -# return -1j*omega(self.freq)*b - -# def S_e(self, prob): -# if all(np.r_[self.mu] == np.r_[prob.curModel.mu]): -# return None -# else: -# eqLocs = prob._eqLocs - -# if eqLocs is 'FE': -# mui_s = prob.curModel.mui - 1./self.mu -# MMui_s = prob.mesh.getFaceInnerProduct(mui_s) -# C = prob.mesh.edgeCurl -# elif eqLocs is 'EF': -# mu_s = prob.curModel.mu - self.mu -# MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True) -# C = prob.mesh.edgeCurl.T - -# return -C.T * (MMui_s * self.bPrimary(prob)) - - #################################################### # Survey #################################################### diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 9a262ffb..1beef79c 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyFDEM import RxFDEM, Src, SurveyFDEM +from SurveyFDEM import Rx, Src, SurveyFDEM from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h from FieldsFDEM import * \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index cb1e1d26..315dc3f1 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -17,7 +17,7 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source XYZ = Utils.ndgrid(x,x,np.r_[0.]) - Rx0 = EM.FDEM.RxFDEM(XYZ, comp) + Rx0 = EM.FDEM.Rx(XYZ, comp) Src = [] From fe2e29031cb221f5664fdb4b0935d2e5e112a6c6 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 21:30:09 -0800 Subject: [PATCH 10/14] updated Rx naming in FDEM analytics test --- tests/em/fdem/forward/test_FDEM_analytics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/em/fdem/forward/test_FDEM_analytics.py b/tests/em/fdem/forward/test_FDEM_analytics.py index 0f9d64f8..7b4d4439 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.RxFDEM(XYZ, 'exi') + rxList = EM.FDEM.Rx(XYZ, 'exi') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) survey = EM.FDEM.SurveyFDEM([Src0]) From 2fbe0af1086b8ce4ab1545dffc3b606e1ebb4175 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 21:47:38 -0800 Subject: [PATCH 11/14] BaseSrcFDEM --> BaseSrc --- SimPEG/EM/FDEM/SrcFDEM.py | 26 +++++++++++++------------- SimPEG/EM/FDEM/SurveyFDEM.py | 2 +- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index f7b95925..93b94c3a 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -5,7 +5,7 @@ from scipy.constants import mu_0 # from SurveyFDEM import RxFDEM -class BaseSrcFDEM(Survey.BaseSrc): +class BaseSrc(Survey.BaseSrc): freq = None # rxPair = EM.FDEM.RxFDEM integrate = True @@ -43,7 +43,7 @@ class BaseSrcFDEM(Survey.BaseSrc): return None -class RawVec_e(BaseSrcFDEM): +class RawVec_e(BaseSrc): """ RawVec electric source. It is defined by the user provided vector S_e @@ -59,7 +59,7 @@ class RawVec_e(BaseSrcFDEM): self._hPrimary = hPrimary self._jPrimary = jPrimary self.freq = float(freq) - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def S_e(self, prob): return self._S_e @@ -77,7 +77,7 @@ class RawVec_e(BaseSrcFDEM): return self._jPrimary -class RawVec_m(BaseSrcFDEM): +class RawVec_m(BaseSrc): """ RawVec magnetic source. It is defined by the user provided vector S_m @@ -95,7 +95,7 @@ class RawVec_m(BaseSrcFDEM): self._hPrimary = np.array(hPrimary,dtype=complex) self._jPrimary = np.array(jPrimary,dtype=complex) - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def S_m(self, prob): return self._S_m @@ -113,7 +113,7 @@ class RawVec_m(BaseSrcFDEM): return self._jPrimary -class RawVec(BaseSrcFDEM): +class RawVec(BaseSrc): """ RawVec source. It is defined by the user provided vectors S_m, S_e @@ -127,7 +127,7 @@ class RawVec(BaseSrcFDEM): self._S_e = np.array(S_e,dtype=complex) self.freq = float(freq) self.integrate = integrate - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def S_m(self, prob): if prob._eqLocs is 'EF' and self.integrate is True: @@ -140,7 +140,7 @@ class RawVec(BaseSrcFDEM): return self._S_e -class MagDipole(BaseSrcFDEM): +class MagDipole(BaseSrc): #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0): @@ -150,7 +150,7 @@ class MagDipole(BaseSrcFDEM): self.moment = moment self.mu = mu self.integrate = False - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def bPrimary(self, prob): eqLocs = prob._eqLocs @@ -209,7 +209,7 @@ class MagDipole(BaseSrcFDEM): return -C.T * (MMui_s * self.bPrimary(prob)) -class MagDipole_Bfield(BaseSrcFDEM): +class MagDipole_Bfield(BaseSrc): #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that #TODO: neither does moment @@ -219,7 +219,7 @@ class MagDipole_Bfield(BaseSrcFDEM): self.orientation = orientation self.moment = moment self.mu = mu - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def bPrimary(self, prob): eqLocs = prob._eqLocs @@ -278,7 +278,7 @@ class MagDipole_Bfield(BaseSrcFDEM): return -C.T * (MMui_s * self.bPrimary(prob)) -class CircularLoop(BaseSrcFDEM): +class CircularLoop(BaseSrc): #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0): @@ -288,7 +288,7 @@ class CircularLoop(BaseSrcFDEM): self.mu = mu self.loc = loc self.integrate = False - BaseSrcFDEM.__init__(self, rxList) + BaseSrc.__init__(self, rxList) def bPrimary(self, prob): eqLocs = prob._eqLocs diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index 2af72d81..a475b550 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -96,7 +96,7 @@ class SurveyFDEM(Survey.BaseSurvey): docstring for SurveyFDEM """ - srcPair = Src.BaseSrcFDEM + srcPair = Src.BaseSrc def __init__(self, srcList, **kwargs): # Sort these by frequency From 7023e999bcc669f8ba15f3610bbfbf3cdb24c675 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 22:17:38 -0800 Subject: [PATCH 12/14] SurveyFDEM --> Survey, FieldsFDEM --> Fields --- SimPEG/EM/FDEM/FDEM.py | 30 ++++++++++---------- SimPEG/EM/FDEM/FieldsFDEM.py | 23 ++++++++------- SimPEG/EM/FDEM/SurveyFDEM.py | 12 ++++---- SimPEG/EM/FDEM/__init__.py | 4 +-- SimPEG/EM/Utils/testingUtils.py | 16 +++++------ tests/em/fdem/forward/test_FDEM_analytics.py | 12 ++++---- 6 files changed, 50 insertions(+), 47 deletions(-) diff --git a/SimPEG/EM/FDEM/FDEM.py b/SimPEG/EM/FDEM/FDEM.py index 9d3be6b2..cce32964 100644 --- a/SimPEG/EM/FDEM/FDEM.py +++ b/SimPEG/EM/FDEM/FDEM.py @@ -1,7 +1,7 @@ -from SimPEG import Survey, Problem, Utils, np, sp, Solver as SimpegSolver +from SimPEG import Problem, Utils, np, sp, Solver as SimpegSolver from scipy.constants import mu_0 -from SurveyFDEM import SurveyFDEM -from FieldsFDEM import FieldsFDEM, FieldsFDEM_e, FieldsFDEM_b, FieldsFDEM_h, FieldsFDEM_j +from SurveyFDEM import Survey as SurveyFDEM +from FieldsFDEM import Fields, Fields_e, Fields_b, Fields_h, Fields_j from SimPEG.EM.Base import BaseEMProblem from SimPEG.EM.Utils import omega @@ -17,8 +17,8 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\ {\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{M^e} \mathbf{s_e}} - if using the E-B formulation (:code:`ProblemFDEM_e` - or :code:`ProblemFDEM_b`) or the magnetic field + if using the E-B formulation (:code:`Problem_e` + or :code:`Problem_b`) or the magnetic field \\\(\\\mathbf{h}\\\) and current density \\\(\\\mathbf{j}\\\) .. math :: @@ -26,13 +26,13 @@ class BaseFDEMProblem(BaseEMProblem): \mathbf{C}^T \mathbf{M_{\\rho}^f} \mathbf{j} + i \omega \mathbf{M_{\mu}^e} \mathbf{h} = \mathbf{M^e} \mathbf{s_m} \\\\ \mathbf{C} \mathbf{h} - \mathbf{j} = \mathbf{s_e} - if using the H-J formulation (:code:`ProblemFDEM_j` or :code:`ProblemFDEM_h`). + if using the H-J formulation (:code:`Problem_j` or :code:`Problem_h`). The problem performs the elimination so that we are solving the system for \\\(\\\mathbf{e},\\\mathbf{b},\\\mathbf{j} \\\) or \\\(\\\mathbf{h}\\\) """ surveyPair = SurveyFDEM - fieldsPair = FieldsFDEM + fieldsPair = Fields def fields(self, m=None): """ @@ -185,7 +185,7 @@ class BaseFDEMProblem(BaseEMProblem): ################################ E-B Formulation ######################################### ########################################################################################## -class ProblemFDEM_e(BaseFDEMProblem): +class Problem_e(BaseFDEMProblem): """ By eliminating the magnetic flux density using @@ -205,7 +205,7 @@ class ProblemFDEM_e(BaseFDEMProblem): _fieldType = 'e' _eqLocs = 'FE' - fieldsPair = FieldsFDEM_e + fieldsPair = Fields_e def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -284,7 +284,7 @@ class ProblemFDEM_e(BaseFDEMProblem): return None -class ProblemFDEM_b(BaseFDEMProblem): +class Problem_b(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{e}\\\) using @@ -304,7 +304,7 @@ class ProblemFDEM_b(BaseFDEMProblem): _fieldType = 'b' _eqLocs = 'FE' - fieldsPair = FieldsFDEM_b + fieldsPair = Fields_b def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -425,7 +425,7 @@ class ProblemFDEM_b(BaseFDEMProblem): ########################################################################################## -class ProblemFDEM_j(BaseFDEMProblem): +class Problem_j(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{h}\\\) using @@ -446,7 +446,7 @@ class ProblemFDEM_j(BaseFDEMProblem): _fieldType = 'j' _eqLocs = 'EF' - fieldsPair = FieldsFDEM_j + fieldsPair = Fields_j def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) @@ -558,7 +558,7 @@ class ProblemFDEM_j(BaseFDEMProblem): -class ProblemFDEM_h(BaseFDEMProblem): +class Problem_h(BaseFDEMProblem): """ We eliminate \\\(\\\mathbf{j}\\\) using @@ -576,7 +576,7 @@ class ProblemFDEM_h(BaseFDEMProblem): _fieldType = 'h' _eqLocs = 'EF' - fieldsPair = FieldsFDEM_h + fieldsPair = Fields_h def __init__(self, mesh, **kwargs): BaseFDEMProblem.__init__(self, mesh, **kwargs) diff --git a/SimPEG/EM/FDEM/FieldsFDEM.py b/SimPEG/EM/FDEM/FieldsFDEM.py index 53f90c4b..d83877c3 100644 --- a/SimPEG/EM/FDEM/FieldsFDEM.py +++ b/SimPEG/EM/FDEM/FieldsFDEM.py @@ -1,13 +1,16 @@ -from SimPEG import Survey, Problem, Utils, np, sp +import numpy as np +import scipy.sparse as sp +import SimPEG +from SimPEG import Utils from SimPEG.EM.Utils import omega -class FieldsFDEM(Problem.Fields): +class Fields(SimPEG.Problem.Fields): """Fancy Field Storage for a FDEM survey.""" knownFields = {} dtype = complex -class FieldsFDEM_e(FieldsFDEM): +class Fields_e(Fields): knownFields = {'eSolution':'E'} aliasFields = { 'e' : ['eSolution','E','_e'], @@ -19,7 +22,7 @@ class FieldsFDEM_e(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -89,7 +92,7 @@ class FieldsFDEM_e(FieldsFDEM): return self._bSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_b(FieldsFDEM): +class Fields_b(Fields): knownFields = {'bSolution':'F'} aliasFields = { 'b' : ['bSolution','F','_b'], @@ -101,7 +104,7 @@ class FieldsFDEM_b(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -190,7 +193,7 @@ class FieldsFDEM_b(FieldsFDEM): return self._eSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_j(FieldsFDEM): +class Fields_j(Fields): knownFields = {'jSolution':'F'} aliasFields = { 'j' : ['jSolution','F','_j'], @@ -202,7 +205,7 @@ class FieldsFDEM_j(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob @@ -293,7 +296,7 @@ class FieldsFDEM_j(FieldsFDEM): return self._hSecondaryDeriv_m(src, v, adjoint) -class FieldsFDEM_h(FieldsFDEM): +class Fields_h(Fields): knownFields = {'hSolution':'E'} aliasFields = { 'h' : ['hSolution','E','_h'], @@ -305,7 +308,7 @@ class FieldsFDEM_h(FieldsFDEM): } def __init__(self,mesh,survey,**kwargs): - FieldsFDEM.__init__(self,mesh,survey,**kwargs) + Fields.__init__(self,mesh,survey,**kwargs) def startup(self): self.prob = self.survey.prob diff --git a/SimPEG/EM/FDEM/SurveyFDEM.py b/SimPEG/EM/FDEM/SurveyFDEM.py index a475b550..150a6c00 100644 --- a/SimPEG/EM/FDEM/SurveyFDEM.py +++ b/SimPEG/EM/FDEM/SurveyFDEM.py @@ -1,4 +1,4 @@ -from SimPEG import Survey, Problem, Utils, np, sp +import SimPEG from SimPEG.EM.Utils import * from scipy.constants import mu_0 import SrcFDEM as Src @@ -7,7 +7,7 @@ import SrcFDEM as Src # Receivers #################################################### -class Rx(Survey.BaseRx): +class Rx(SimPEG.Survey.BaseRx): knownRxTypes = { 'exr':['e', 'Ex', 'real'], @@ -41,7 +41,7 @@ class Rx(Survey.BaseRx): radius = None def __init__(self, locs, rxType): - Survey.BaseRx.__init__(self, locs, rxType) + SimPEG.Survey.BaseRx.__init__(self, locs, rxType) @property def projField(self): @@ -91,7 +91,7 @@ class Rx(Survey.BaseRx): # Survey #################################################### -class SurveyFDEM(Survey.BaseSurvey): +class Survey(SimPEG.Survey.BaseSurvey): """ docstring for SurveyFDEM """ @@ -101,7 +101,7 @@ class SurveyFDEM(Survey.BaseSurvey): def __init__(self, srcList, **kwargs): # Sort these by frequency self.srcList = srcList - Survey.BaseSurvey.__init__(self, **kwargs) + SimPEG.Survey.BaseSurvey.__init__(self, **kwargs) _freqDict = {} for src in srcList: @@ -136,7 +136,7 @@ class SurveyFDEM(Survey.BaseSurvey): return self._freqDict[freq] def projectFields(self, u): - data = Survey.Data(self) + data = SimPEG.Survey.Data(self) for src in self.srcList: for rx in src.rxList: data[src, rx] = rx.projectFields(src, self.mesh, u) diff --git a/SimPEG/EM/FDEM/__init__.py b/SimPEG/EM/FDEM/__init__.py index 1beef79c..978972f5 100644 --- a/SimPEG/EM/FDEM/__init__.py +++ b/SimPEG/EM/FDEM/__init__.py @@ -1,3 +1,3 @@ -from SurveyFDEM import Rx, Src, SurveyFDEM -from FDEM import BaseFDEMProblem, ProblemFDEM_e, ProblemFDEM_b, ProblemFDEM_j, ProblemFDEM_h +from SurveyFDEM import Rx, Src, Survey +from FDEM import BaseFDEMProblem, Problem_e, Problem_b, Problem_j, Problem_h from FieldsFDEM import * \ No newline at end of file diff --git a/SimPEG/EM/Utils/testingUtils.py b/SimPEG/EM/Utils/testingUtils.py index 315dc3f1..8c703083 100644 --- a/SimPEG/EM/Utils/testingUtils.py +++ b/SimPEG/EM/Utils/testingUtils.py @@ -47,20 +47,20 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False): print ' Fetching %s problem' % (fdemType) if fdemType == 'e': - survey = EM.FDEM.SurveyFDEM(Src) - prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_e(mesh, mapping=mapping) elif fdemType == 'b': - survey = EM.FDEM.SurveyFDEM(Src) - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) elif fdemType == 'j': - survey = EM.FDEM.SurveyFDEM(Src) - prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_j(mesh, mapping=mapping) elif fdemType == 'h': - survey = EM.FDEM.SurveyFDEM(Src) - prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) + survey = EM.FDEM.Survey(Src) + prb = EM.FDEM.Problem_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 7b4d4439..9786e7c8 100644 --- a/tests/em/fdem/forward/test_FDEM_analytics.py +++ b/tests/em/fdem/forward/test_FDEM_analytics.py @@ -31,9 +31,9 @@ class FDEM_analyticTests(unittest.TestCase): rxList = EM.FDEM.Rx(XYZ, 'exi') Src0 = EM.FDEM.Src.MagDipole([rxList],loc=np.r_[0.,0.,0.], freq=freq) - survey = EM.FDEM.SurveyFDEM([Src0]) + survey = EM.FDEM.Survey([Src0]) - prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) + prb = EM.FDEM.Problem_b(mesh, mapping=mapping) prb.pair(survey) try: @@ -120,13 +120,13 @@ class FDEM_analyticTests(unittest.TestCase): # Pair the problem and survey - surveye = EM.FDEM.SurveyFDEM(de_p) - surveym = EM.FDEM.SurveyFDEM(dm_p) + surveye = EM.FDEM.Survey(de_p) + surveym = EM.FDEM.Survey(dm_p) mapping = [('sigma', Maps.IdentityMap(mesh)),('mu', Maps.IdentityMap(mesh))] - prbe = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) - prbm = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) + prbe = EM.FDEM.Problem_h(mesh, mapping=mapping) + prbm = EM.FDEM.Problem_e(mesh, mapping=mapping) prbe.pair(surveye) # pair problem and survey prbm.pair(surveym) From 699f6452b2b54ed644a295faeee7bd8950ac6232 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 22:29:44 -0800 Subject: [PATCH 13/14] cleaned up imports for source --- SimPEG/EM/FDEM/SrcFDEM.py | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 93b94c3a..7c305090 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -1,13 +1,11 @@ from SimPEG import Survey, Problem, Utils, np, sp -# import SimPEG.EM as EM -from SimPEG.EM.Utils import * from scipy.constants import mu_0 -# from SurveyFDEM import RxFDEM +# from SurveyFDEM import Rx class BaseSrc(Survey.BaseSrc): freq = None - # rxPair = EM.FDEM.RxFDEM + # rxPair = Rx integrate = True def eval(self, prob): From 622152c0c8e081fe8ef1c749221854f35e7fb9b0 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Tue, 17 Nov 2015 22:33:59 -0800 Subject: [PATCH 14/14] import simpegEM utils --- SimPEG/EM/FDEM/SrcFDEM.py | 1 + 1 file changed, 1 insertion(+) diff --git a/SimPEG/EM/FDEM/SrcFDEM.py b/SimPEG/EM/FDEM/SrcFDEM.py index 7c305090..5dcf9407 100644 --- a/SimPEG/EM/FDEM/SrcFDEM.py +++ b/SimPEG/EM/FDEM/SrcFDEM.py @@ -1,5 +1,6 @@ from SimPEG import Survey, Problem, Utils, np, sp from scipy.constants import mu_0 +from SimPEG.EM.Utils import * # from SurveyFDEM import Rx