wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-02 18:48:59 +08:00
Code Issues Packages Projects Releases Wiki Activity

fixed loop source and mag dipole fields for variable mu, added both to testing

Browse Source
This commit is contained in:
Lindsey Heagy
2015-06-23 23:02:33 -07:00
parent ac3e7765fe
commit d9e8996336
4 changed files with 72 additions and 34 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Analytics/FDEM.py
+2 -2
View File
@@ -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':
simpegEM/FDEM/SurveyFDEM.py
+51 -15
View File
@@ -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
simpegEM/Tests/test_FDEM.py
+10 -9
View File
@@ -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()
simpegEM/Utils/SrcUtils.py
+9 -8
View File
@@ -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. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
@@ -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. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
@@ -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)