From 51342e7cc81bfb12297d8d5e4f6806855bebf271 Mon Sep 17 00:00:00 2001 From: Lindsey Heagy Date: Wed, 29 Apr 2015 16:32:50 -0700 Subject: [PATCH] Fairly major source refactor: - removed the folder Sources and put those routines inside of SrcUtils. - Broke apart calls for MagDipole, MagDipole_B, CircularLoop --- simpegEM/FDEM/SurveyFDEM.py | 269 +++++++++++------- simpegEM/Sources/CircularLoop.py | 101 ------- simpegEM/Sources/__init__.py | 3 - simpegEM/TDEM/BaseTDEM.py | 2 +- simpegEM/TDEM/SurveyTDEM.py | 2 +- simpegEM/Tests/test_FDEM.py | 2 +- .../magneticDipole.py => Utils/SrcUtils.py} | 106 ++++++- simpegEM/Utils/__init__.py | 3 +- simpegEM/__init__.py | 1 - 9 files changed, 280 insertions(+), 209 deletions(-) delete mode 100644 simpegEM/Sources/CircularLoop.py delete mode 100644 simpegEM/Sources/__init__.py rename simpegEM/{Sources/magneticDipole.py => Utils/SrcUtils.py} (51%) diff --git a/simpegEM/FDEM/SurveyFDEM.py b/simpegEM/FDEM/SurveyFDEM.py index 1106267e..f22f38cd 100644 --- a/simpegEM/FDEM/SurveyFDEM.py +++ b/simpegEM/FDEM/SurveyFDEM.py @@ -1,5 +1,5 @@ from SimPEG import Survey, Problem, Utils, np, sp -from simpegEM import Sources +from simpegEM.Utils import SrcUtils from simpegEM.Utils.EMUtils import omega @@ -91,134 +91,207 @@ class RxFDEM(Survey.BaseRx): # Sources #################################################### -# class SrcFDEM(Survey.BaseSrc): -# freq = None -# rxPair = RxFDEM -# knownSrcTypes = {} - - class SrcFDEM(Survey.BaseSrc): - #TODO: Break these out into Classes of Sources. - - freq = None #: Frequency (float) - + freq = None rxPair = RxFDEM + knownSrcTypes = ['Simple', 'MagDipole'] #TODO: Do we want to just classify them by Magnetic, Electric, Both? - knownSrcTypes = ['VMD', 'VMD_B', 'CircularLoop', 'Simple'] - radius = None +class SrcFDEM_Simple_e(SrcFDEM): + """ + Simple electric source. It is defined by the user provided vector S_e - def __init__(self, loc, srcType, freq, rxList): + :param numpy.array S_e: electric source term + :param float freq: frequency + :param rxList: receiver list + """ + + def __init__(self, S_e, freq, rxList): + self.S_e = S_e self.freq = float(freq) - Survey.BaseSrc.__init__(self, loc, srcType, rxList) + SrcFDEM.__init__(self, None, 'Simple', rxList) def getSource(self, prob): + return None, self.S_e - src = self - freq = src.freq - solType = prob._fieldType # Hack, should just ask whether j_m, j_g are defined on edges or faces - - if solType == 'e' or solType == 'b': - gridEJx = prob.mesh.gridEx - gridEJy = prob.mesh.gridEy - gridEJz = prob.mesh.gridEz - nEJ = prob.mesh.nE - - gridBHx = prob.mesh.gridFx - gridBHy = prob.mesh.gridFy - gridBHz = prob.mesh.gridFz - nBH = prob.mesh.nF + def getSourceDeriv(self, prob, v, adjoint = False): + return None, None +class SrcFDEM_Simple_m(SrcFDEM): + """ + Simple 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, S_m, freq, rxList): + self.S_m = S_m + self.freq = float(freq) + SrcFDEM.__init__(self, None, 'Simple', rxList) + + def getSource(self, prob): + return self.S_m, None + + def getSourceDeriv(self, prob, v, adjoint = False): + return None, None + + +class SrcFDEM_Simple(SrcFDEM): + """ + Simple 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, S_m, S_e, freq, rxList): + self.S_m = S_m + self.S_e = S_e + SrcFDEM.__init__(self, None, 'Simple', rxList) + + def getSource(self, prob): + return self.S_m, self.S_e + + def getSourceDeriv(self, prob, v, adjoint=None): + return None, None + + +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, loc, freq, rxList, orientation='Z'): + self.freq = float(freq) + self.orientation = orientation + SrcFDEM.__init__(self, loc, 'MagDipole', rxList) + + def getSource(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 - mui = prob.MfMui - - elif solType == 'h' or solType == 'j': - gridEJx = prob.mesh.gridFx - gridEJy = prob.mesh.gridFy - gridEJz = prob.mesh.gridFz - nEJ = prob.mesh.nF - - gridBHx = prob.mesh.gridEx - gridBHy = prob.mesh.gridEy - gridBHz = prob.mesh.gridEz - nBH = prob.mesh.nE + elif eqLocs is 'EF': + gridX = prob.mesh.gridFx + gridY = prob.mesh.gridFy + gridZ = prob.mesh.gridFz C = prob.mesh.edgeCurl.T - mui = prob.MeMuI - - else: - NotImplementedError('Only E or F sources') if prob.mesh._meshType is 'CYL': if not prob.mesh.isSymmetric: + # TODO ? raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!') - - if src.srcType == 'VMD': - SRC = Sources.MagneticDipoleVectorPotential(src.loc, gridEJy, 'y') - elif src.srcType == 'CircularLoop': - SRC = Sources.MagneticLoopVectorPotential(src.loc, gridEJy, 'y', src.radius) - else: - raise NotImplementedError('Only VMD and CircularLoop') - - elif prob.mesh._meshType is 'TENSOR': - - if src.srcType == 'VMD': - srcfct = Sources.MagneticDipoleVectorPotential - SRCx = srcfct(src.loc, gridEJx, 'x') - SRCy = srcfct(src.loc, gridEJy, 'y') - SRCz = srcfct(src.loc, gridEJz, 'z') - - elif src.srcType == 'VMD_B': - srcfct = Sources.MagneticDipoleFields - SRCx = srcfct(src.loc, gridBHx, 'x') - SRCy = srcfct(src.loc, gridBHy, 'y') - SRCz = srcfct(src.loc, gridBHz, 'z') - - elif src.srcType == 'CircularLoop': - srcfct = Sources.MagneticLoopVectorPotential - SRCx = srcfct(src.loc, gridEJx, 'x', src.radius) - SRCy = srcfct(src.loc, gridEJy, 'y', src.radius) - SRCz = srcfct(src.loc, gridEJz, 'z', src.radius) - else: - - raise NotImplemented('%s srcType is not implemented' % src.srcType) - SRC = np.concatenate((SRCx, SRCy, SRCz)) + a = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y') else: - raise Exception('Unknown mesh for VMD') + srcfct = SrcUtils.MagneticDipoleVectorPotential + ax = srcfct(self.loc, gridX, 'x') + ay = srcfct(self.loc, gridY, 'y') + az = srcfct(self.loc, gridZ, 'z') + a = np.concatenate((ax, ay, az)) - # b-forumlation - if src.srcType == 'VMD_B': - b_0 = SRC + S_m = -1j*omega(self.freq)*C*a + + return S_m, None + + + def getSourceDeriv(self, prob, v, adjoint=None): + return None, None + + +class MagDipole_Bfield(SrcFDEM): + + #TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that + def __init__(self, loc, freq, rxList, orientation='Z'): + self.freq = float(freq) + self.orientation = orientation + SrcFDEM.__init__(self, loc, 'MagDipole', rxList) + + def getSource(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 = SrcUtils.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') + bz = srcfct(self.loc, gridZ, 'z') + b = np.concatenate((bx,bz)) else: - a = SRC - b_0 = C*a + bx = srcfct(self.loc, gridX, 'x') + by = srcfct(self.loc, gridY, 'y') + bz = srcfct(self.loc, gridZ, 'z') + b = np.concatenate((bx,by,bz)) - return -1j*omega(freq)*b_0, None + return -1j*omega(self.freq)*b, None -class SimpleSrcFDEM_e(SrcFDEM): + def getSourceDeriv(self, prob, v, adjoint=None): + return None, None - def __init__(self, vec, freq, rxList): - self.vec = vec + +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, loc, freq, rxList, orientation='Z', radius = 1.): self.freq = float(freq) - SrcFDEM.__init__(self, None, 'Simple', freq, rxList) + self.orientation = orientation + self.radius = radius + SrcFDEM.__init__(self, loc, 'MagDipole', rxList) def getSource(self, prob): - return None, self.vec + 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 = SrcUtils.MagneticDipoleVectorPotential(src.loc, gridY, 'y', self.radius) + + 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) + a = np.concatenate((ax, ay, az)) + + return -1j*omega(self.freq)*C*a -class SimpleSrcFDEM_m(SrcFDEM): - - def __init__(self, vec, freq, rxList): - self.vec = vec - self.freq = float(freq) - SrcFDEM.__init__(self, None, 'Simple', freq, rxList) - - def getSource(self, prob): - return self.vec, None - +#################################################### +# Survey +#################################################### class SurveyFDEM(Survey.BaseSurvey): """ diff --git a/simpegEM/Sources/CircularLoop.py b/simpegEM/Sources/CircularLoop.py deleted file mode 100644 index feb55bc8..00000000 --- a/simpegEM/Sources/CircularLoop.py +++ /dev/null @@ -1,101 +0,0 @@ -from SimPEG import * -from scipy.special import ellipk, ellipe - -def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): - """ - Calculate the vector potential of horizontal circular loop - at given locations - - :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) - :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh - :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list - :param numpy.ndarray I: Input current of the loop - :param numpy.ndarray radius: radius of the loop - :rtype: numpy.ndarray - :return: The vector potential each dipole at each observation location - """ - - if type(component) in [list, tuple]: - out = range(len(component)) - for i, comp in enumerate(component): - out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius) - 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) - - srcLoc = np.atleast_2d(srcLoc) - obsLoc = np.atleast_2d(obsLoc) - - n = obsLoc.shape[0] - nSrc = srcLoc.shape[0] - - if component=='z': - A = np.zeros((n, nSrc)) - if nSrc ==1: - return A.flatten() - return A - - else: - - A = np.zeros((n, nSrc)) - for i in range (nSrc): - x = obsLoc[:, 0] - srcLoc[i, 0] - y = obsLoc[:, 1] - srcLoc[i, 1] - z = obsLoc[:, 2] - srcLoc[i, 2] - r = np.sqrt(x**2 + y**2) - m = (4 * radius * r) / ((radius + r)**2 + z**2) - m[m > 1.] = 1. - # m might be slightly larger than 1 due to rounding errors - # but ellipke requires 0 <= m <= 1 - K = ellipk(m) - E = ellipe(m) - ind = (r > 0) & (m < 1) - # % 1/r singular at r = 0 and K(m) singular at m = 1 - Aphi = np.zeros(n) - # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi. - Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind]) - if component == 'x': - A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] ) - elif component == 'y': - A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] ) - else: - raise ValueError('Invalid component') - - if nSrc == 1: - return A.flatten() - return A - -if __name__ == '__main__': - from SimPEG import Mesh - import matplotlib.pyplot as plt - cs = 20 - ncx, ncy, ncz = 41, 41, 40 - hx = np.ones(ncx)*cs - hy = np.ones(ncy)*cs - hz = np.ones(ncz)*cs - mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC') - srcLoc = np.r_[0., 0., 0.] - Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200) - Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200) - Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200) - A = np.r_[Ax, Ay, Az] - B0 = mesh.edgeCurl*A - J0 = mesh.edgeCurl.T*B0 - - # mesh.plotImage(A, vType = 'Ex') - # mesh.plotImage(A, vType = 'Ey') - - mesh.plotImage(B0, vType = 'Fx') - mesh.plotImage(B0, vType = 'Fy') - mesh.plotImage(B0, vType = 'Fz') - - # # mesh.plotImage(J0, vType = 'Ex') - # mesh.plotImage(J0, vType = 'Ey') - # mesh.plotImage(J0, vType = 'Ez') - - plt.show() - - diff --git a/simpegEM/Sources/__init__.py b/simpegEM/Sources/__init__.py deleted file mode 100644 index 1f04379e..00000000 --- a/simpegEM/Sources/__init__.py +++ /dev/null @@ -1,3 +0,0 @@ -from magneticDipole import MagneticDipoleVectorPotential -from CircularLoop import MagneticLoopVectorPotential -from magneticDipole import MagneticDipoleFields diff --git a/simpegEM/TDEM/BaseTDEM.py b/simpegEM/TDEM/BaseTDEM.py index 3f695416..a4955945 100644 --- a/simpegEM/TDEM/BaseTDEM.py +++ b/simpegEM/TDEM/BaseTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Solver, Problem from SimPEG.Problem import BaseTimeProblem -from simpegEM import Sources +from simpegEM.Utils import SrcUtils from scipy.constants import mu_0 from SimPEG.Utils import sdiag, mkvc from SimPEG import Utils, Mesh diff --git a/simpegEM/TDEM/SurveyTDEM.py b/simpegEM/TDEM/SurveyTDEM.py index 02be5d02..668b43a0 100644 --- a/simpegEM/TDEM/SurveyTDEM.py +++ b/simpegEM/TDEM/SurveyTDEM.py @@ -1,6 +1,6 @@ from SimPEG import Utils, Survey, np from SimPEG.Survey import BaseSurvey -from simpegEM import Sources +from simpegEM.Utils import SrcUtils from BaseTDEM import FieldsTDEM diff --git a/simpegEM/Tests/test_FDEM.py b/simpegEM/Tests/test_FDEM.py index 8f4a8365..1c5f4230 100644 --- a/simpegEM/Tests/test_FDEM.py +++ b/simpegEM/Tests/test_FDEM.py @@ -35,7 +35,7 @@ def getProblem(fdemType, comp): 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) - Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0]) + Src0 = EM.FDEM.SrcFDEM_MagDipole(np.r_[0.,0.,0.], freq, [Rx0]) survey = EM.FDEM.SurveyFDEM([Src0]) diff --git a/simpegEM/Sources/magneticDipole.py b/simpegEM/Utils/SrcUtils.py similarity index 51% rename from simpegEM/Sources/magneticDipole.py rename to simpegEM/Utils/SrcUtils.py index 526fc1ac..5ff26467 100644 --- a/simpegEM/Sources/magneticDipole.py +++ b/simpegEM/Utils/SrcUtils.py @@ -1,6 +1,6 @@ -import numpy as np +from SimPEG import * +from scipy.special import ellipk, ellipe from scipy.constants import mu_0, pi -from SimPEG import Mesh def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)): """ @@ -53,6 +53,7 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0 return A.flatten() return A + def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): """ Calculate the vector potential of a set of magnetic dipoles @@ -98,3 +99,104 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.): if nSrc == 1: return B.flatten() return B + + + +def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius): + """ + Calculate the vector potential of horizontal circular loop + at given locations + + :param numpy.ndarray srcLoc: Location of the source(s) (x, y, z) + :param numpy.ndarray,SimPEG.Mesh obsLoc: Where the potentials will be calculated (x, y, z) or a SimPEG Mesh + :param str,list component: The component to calculate - 'x', 'y', or 'z' if an array, or grid type if mesh, can be a list + :param numpy.ndarray I: Input current of the loop + :param numpy.ndarray radius: radius of the loop + :rtype: numpy.ndarray + :return: The vector potential each dipole at each observation location + """ + + if type(component) in [list, tuple]: + out = range(len(component)) + for i, comp in enumerate(component): + out[i] = MagneticLoopVectorPotential(srcLoc, obsLoc, comp, radius) + 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) + + srcLoc = np.atleast_2d(srcLoc) + obsLoc = np.atleast_2d(obsLoc) + + n = obsLoc.shape[0] + nSrc = srcLoc.shape[0] + + if component=='z': + A = np.zeros((n, nSrc)) + if nSrc ==1: + return A.flatten() + return A + + else: + + A = np.zeros((n, nSrc)) + for i in range (nSrc): + x = obsLoc[:, 0] - srcLoc[i, 0] + y = obsLoc[:, 1] - srcLoc[i, 1] + z = obsLoc[:, 2] - srcLoc[i, 2] + r = np.sqrt(x**2 + y**2) + m = (4 * radius * r) / ((radius + r)**2 + z**2) + m[m > 1.] = 1. + # m might be slightly larger than 1 due to rounding errors + # but ellipke requires 0 <= m <= 1 + K = ellipk(m) + E = ellipe(m) + ind = (r > 0) & (m < 1) + # % 1/r singular at r = 0 and K(m) singular at m = 1 + Aphi = np.zeros(n) + # % Common factor is (mu * I) / pi with I = 1 and mu = 4e-7 * pi. + Aphi[ind] = 4e-7 / np.sqrt(m[ind]) * np.sqrt(radius / r[ind]) *((1. - m[ind] / 2.) * K[ind] - E[ind]) + if component == 'x': + A[ind, i] = Aphi[ind] * (-y[ind] / r[ind] ) + elif component == 'y': + A[ind, i] = Aphi[ind] * ( x[ind] / r[ind] ) + else: + raise ValueError('Invalid component') + + if nSrc == 1: + return A.flatten() + return A + +if __name__ == '__main__': + from SimPEG import Mesh + import matplotlib.pyplot as plt + cs = 20 + ncx, ncy, ncz = 41, 41, 40 + hx = np.ones(ncx)*cs + hy = np.ones(ncy)*cs + hz = np.ones(ncz)*cs + mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC') + srcLoc = np.r_[0., 0., 0.] + Ax = MagneticLoopVectorPotential(srcLoc, mesh.gridEx, 'x', 200) + Ay = MagneticLoopVectorPotential(srcLoc, mesh.gridEy, 'y', 200) + Az = MagneticLoopVectorPotential(srcLoc, mesh.gridEz, 'z', 200) + A = np.r_[Ax, Ay, Az] + B0 = mesh.edgeCurl*A + J0 = mesh.edgeCurl.T*B0 + + # mesh.plotImage(A, vType = 'Ex') + # mesh.plotImage(A, vType = 'Ey') + + mesh.plotImage(B0, vType = 'Fx') + mesh.plotImage(B0, vType = 'Fy') + mesh.plotImage(B0, vType = 'Fz') + + # # mesh.plotImage(J0, vType = 'Ex') + # mesh.plotImage(J0, vType = 'Ey') + # mesh.plotImage(J0, vType = 'Ez') + + plt.show() + + diff --git a/simpegEM/Utils/__init__.py b/simpegEM/Utils/__init__.py index 608ff235..6e430cf9 100644 --- a/simpegEM/Utils/__init__.py +++ b/simpegEM/Utils/__init__.py @@ -1,4 +1,5 @@ # import Sources # import Ana # import Solver -import EMUtils \ No newline at end of file +import EMUtils +import SrcUtils \ No newline at end of file diff --git a/simpegEM/__init__.py b/simpegEM/__init__.py index 381c18c8..6a1ca774 100644 --- a/simpegEM/__init__.py +++ b/simpegEM/__init__.py @@ -2,6 +2,5 @@ import TDEM import FDEM import Base -import Sources import Analytics import Utils