diff --git a/simpegEM/Analytics/FDEM.py b/simpegEM/Analytics/FDEM.py index a6c051cd..9abb0a15 100644 --- a/simpegEM/Analytics/FDEM.py +++ b/simpegEM/Analytics/FDEM.py @@ -42,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0): return hz -def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X', mu = mu_0): +def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0): """ Analytical solution for a dipole in a whole-space. @@ -79,7 +79,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X', mu = k = np.sqrt( -1j*2.*np.pi*f*mu*sig ) kr = k*r - front = m / (4.*pi * r**3.) * np.exp(-1j*kr) + front = moment / (4.*pi * r**3.) * np.exp(-1j*kr) mid = -kr**2. + 3.*1j*kr + 3. if orientation.upper() == 'X': diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index ff60cb95..8a6937e1 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -217,13 +217,13 @@ 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) + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment) else: srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', mu=self.mu) - ay = srcfct(self.loc, gridY, 'y', mu=self.mu) - az = srcfct(self.loc, gridZ, 'z', mu=self.mu) + 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 @@ -286,13 +286,13 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - bx = srcfct(self.loc, gridX, 'x', mu=self.mu) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu) + 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) - by = srcfct(self.loc, gridY, 'y', mu=self.mu) - bz = srcfct(self.loc, gridZ, 'z', mu=self.mu) + 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 @@ -305,14 +305,33 @@ class SrcFDEM_MagDipole_Bfield(SrcFDEM): 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.): + 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 SrcFDEM.__init__(self, rxList) def bPrimary(self,prob): @@ -334,25 +353,42 @@ class SrcFDEM_CircularLoop(SrcFDEM): if not prob.mesh.isSymmetric: # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y', self.radius) + a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', moment=self.radius, mu=self.mu) else: srcfct = SrcUtils.MagneticDipoleVectorPotential - ax = srcfct(self.loc, gridX, 'x', self.radius) - ay = srcfct(self.loc, gridY, 'y', self.radius) - az = srcfct(self.loc, gridZ, 'z', self.radius) + 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 h_from_b + 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/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index bcddd91d..e3cbfc7b 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -36,7 +36,8 @@ def getProblem(fdemType, comp): XYZ = Utils.ndgrid(x,x,np.r_[0.]) Rx0 = EM.FDEM.RxFDEM(XYZ, comp) Src0 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) - + Src1 = EM.FDEM.SrcFDEM_MagDipole_Bfield([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) + Src2 = EM.FDEM.SrcFDEM_CircularLoop([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.]) if verbose: print ' Fetching %s problem' % (fdemType) @@ -46,32 +47,32 @@ def getProblem(fdemType, comp): 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. - Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) prb = EM.FDEM.ProblemFDEM_e(mesh, mapping=mapping) elif fdemType == '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. - Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping) elif fdemType == '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. - Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) prb = EM.FDEM.ProblemFDEM_j(mesh, mapping=mapping) elif fdemType == 'h': 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. - Src1 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) - survey = EM.FDEM.SurveyFDEM([Src0,Src1]) + Src3 = EM.FDEM.SrcFDEM_RawVec([Rx0], freq[0], S_m, S_e) + survey = EM.FDEM.SurveyFDEM([Src0,Src1,Src2,Src3]) prb = EM.FDEM.ProblemFDEM_h(mesh, mapping=mapping) else: raise NotImplementedError() diff --git a/simpegEM/Utils/SrcUtils.py b/simpegEM/Utils/SrcUtils.py index 7c672bf4..1827f6b2 100644 --- a/simpegEM/Utils/SrcUtils.py +++ b/simpegEM/Utils/SrcUtils.py @@ -2,7 +2,7 @@ from SimPEG import * from scipy.special import ellipk, ellipe from scipy.constants import mu_0, pi -def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.), mu = mu_0): +def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, moment=1., dipoleMoment=(0., 0., 1.), mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -15,6 +15,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 :return: The vector potential each dipole at each observation location """ #TODO: break this out! + if type(component) in [list, tuple]: out = range(len(component)) for i, comp in enumerate(component): @@ -54,7 +55,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 return A -def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): +def MagneticDipoleFields(srcLoc, obsLoc, component, moment=1., mu = mu_0): """ Calculate the vector potential of a set of magnetic dipoles at given locations 'ref. ' @@ -62,7 +63,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) :param numpy.ndarray obsLoc: Where the potentials will be calculated (x, y, z) :param str component: The component to calculate - 'x', 'y', or 'z' - :param numpy.ndarray dipoleMoment: The vector dipole moment (vertical) + :param numpy.ndarray moment: The vector dipole moment (vertical) :rtype: numpy.ndarray :return: The vector potential each dipole at each observation location """ @@ -78,12 +79,12 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): srcLoc = np.atleast_2d(srcLoc) obsLoc = np.atleast_2d(obsLoc) - dipoleMoment = np.atleast_2d(dipoleMoment) + moment = np.atleast_2d(moment) nFaces = obsLoc.shape[0] nSrc = srcLoc.shape[0] - m = np.array(dipoleMoment).repeat(nFaces, axis=0) + m = np.array(moment).repeat(nFaces, axis=0) B = np.empty((nFaces, nSrc)) for i in range(nSrc): dR = obsLoc - srcLoc[i, np.newaxis].repeat(nFaces, axis=0) @@ -102,7 +103,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1., mu = mu_0): -def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): +def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius, mu=mu_0): """ Calculate the vector potential of horizontal circular loop at given locations @@ -119,13 +120,13 @@ def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): if type(component) in [list, tuple]: out = range(len(component)) for i, comp in enumerate(component): - out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius) + out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius, mu) return np.concatenate(out) if isinstance(obsLoc, Mesh.BaseMesh): mesh = obsLoc assert component in ['Ex','Ey','Ez','Fx','Fy','Fz'], "Components must be in: ['Ex','Ey','Ez','Fx','Fy','Fz']" - return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius) + return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius, mu) srcLoc = np.atleast_2d(srcLoc) obsLoc = np.atleast_2d(obsLoc)