wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-07-10 12:39:23 +08:00
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #24 from simpeg/feat/sourceRefactor

Browse Source 
Feat/source refactor
This commit is contained in:
Lindsey
2015-06-30 14:10:05 -07:00
parent 527133c709 5ec86300f6
commit f96cf0418a
7 changed files with 198 additions and 88 deletions
Download Patch File Download Diff File Expand all files Collapse all files
simpegEM/Analytics/FDEM.py
+29 -16
View File
@@ -5,7 +5,8 @@ from scipy.special import erf
import matplotlib.pyplot as plt
from SimPEG import Utils
def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
def hzAnalyticDipoleF(r, freq, sigma, secondary=True, mu=mu_0):
"""
4.56 in Ward and Hohmann
@@ -25,7 +26,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
"""
r = np.abs(r)
k = np.sqrt(-1j*2.*np.pi*freq*mu_0*sigma)
k = np.sqrt(-1j*2.*np.pi*freq*mu*sigma)
m = 1
front = m / (2. * np.pi * (k**2) * (r**5) )
@@ -41,7 +42,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
return hz
def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, moment=1., orientation='X', mu = mu_0):
"""
Analytical solution for a dipole in a whole-space.
@@ -75,10 +76,10 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig )
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':
@@ -96,9 +97,9 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
Hy = front*( (dy*dz/r**2.) * mid )
Hz = front*( (dz/r)**2. * mid + (kr**2. - 1j*kr - 1.) )
Bx = mu_0*Hx
By = mu_0*Hy
Bz = mu_0*Hz
Bx = mu*Hx
By = mu*Hy
Bz = mu*Hz
if Bx.ndim is 1:
Bx = Utils.mkvc(Bx,2)
@@ -112,7 +113,7 @@ def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
return Bx, By, Bz
def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, current=1., length=1., orientation='X', mu=mu_0):
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
dx = XYZ[:,0]-srcLoc[0]
@@ -120,21 +121,33 @@ def ElectricDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
dz = XYZ[:,2]-srcLoc[2]
r = np.sqrt( dx**2. + dy**2. + dz**2.)
k = np.sqrt( -1j*2.*np.pi*f*mu_0*sig )
k = np.sqrt( -1j*2.*np.pi*f*mu*sig )
kr = k*r
front = moment / (4. * np.pi * sig * r**3) * exp(-1j*k*r)
front = current * length / (4. * np.pi * sig * r**3) * np.exp(-1j*k*r)
mid = -k**2 * r**2 + 3*1j*k*r + 3
Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r))
Ey = front*(dx*dy / r**2)*mid
Ez = front*(dx*dz / r**2)*mid
# Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r))
# Ey = front*(dx*dy / r**2)*mid
# Ez = front*(dx*dz / r**2)*mid
if orientation.upper() == 'X':
Ex = front*((dx**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ey = front*(dx*dy / r**2)*mid
Ez = front*(dx*dz / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Y':
return Ez, Ex, Ey
# x--> y, y--> z, z-->x
Ey = front*((dy**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ez = front*(dy*dz / r**2)*mid
Ex = front*(dy*dx / r**2)*mid
return Ex, Ey, Ez
elif orientation.upper() == 'Z':
return Ey, Ez, Ex
# x --> z, y --> x, z --> y
Ez = front*((dz**2 / r**2)*mid + (k**2 * r**2 -1j*k*r-1.))
Ex = front*(dz*dx / r**2)*mid
Ey = front*(dz*dy / r**2)*mid
return Ex, Ey, Ez
# return Ey, Ez, Ex
simpegEM/Base.py
+1 -1
View File
@@ -141,5 +141,5 @@ class BaseEMProblem(Problem.BaseProblem):
# TODO: This isn't going to work yet
# TODO: This should take a vector
def dMfRhoIDeriv(self,u):
def MfRhoIDeriv(self,u):
return self.mesh.getFaceInnerProductDeriv(self.curModel.rho, invMat=True)(u) * self.curModel.rhoDeriv
simpegEM/FDEM/FDEM.py
+4 -8
View File
@@ -19,7 +19,7 @@ class BaseFDEMProblem(BaseEMProblem):
surveyPair = SurveyFDEM
fieldsPair = FieldsFDEM
def fields(self, m):
def fields(self, m=None):
self.curModel = m
F = self.fieldsPair(self.mesh, self.survey)
@@ -151,10 +151,6 @@ class BaseFDEMProblem(BaseEMProblem):
return S_m, S_e
def getSourceTermDeriv(self,freq,m,v,u=None,adjoint=False):
raise NotImplementedError('getSourceTermDeriv not implemented yet')
return None, None
##########################################################################################
################################ E-B Formulation #########################################
@@ -321,14 +317,14 @@ class ProblemFDEM_b(BaseFDEMProblem):
def getRHSDeriv_m(self, src, v, adjoint=False):
C = self.mesh.edgeCurl
S_m, S_e = self.getSourceTerm(src.freq)
S_m, S_e = src.eval(self)
MfMui = self.MfMui
if self._makeASymmetric and adjoint:
v = self.MfMui * v
if S_e is not None:
MeSigmaIDeriv = self.MeSigmaIDeriv(S_e)
MeSigmaIDeriv = self.MeSigmaIDeriv(Utils.mkvc(S_e))
if not adjoint:
RHSderiv = C * (MeSigmaIDeriv * v)
elif adjoint:
@@ -577,7 +573,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, src, v, adjoint=False):
_, S_e = self.getSourceTerm(src.freq)
_, S_e = src.eval(self)
C = self.mesh.edgeCurl
MfRho = self.MfRho
MfRhoDeriv = self.MfRhoDeriv(S_e)
simpegEM/FDEM/FieldsFDEM.py
+23 -20
View File
@@ -7,7 +7,6 @@ class FieldsFDEM(Problem.Fields):
knownFields = {}
dtype = complex
class FieldsFDEM_e(FieldsFDEM):
knownFields = {'eSolution':'E'}
aliasFields = {
@@ -23,6 +22,7 @@ class FieldsFDEM_e(FieldsFDEM):
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
# def getDeriv_u(self, fieldsList, src, v, adjoint=False):
@@ -32,7 +32,7 @@ class FieldsFDEM_e(FieldsFDEM):
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
ep = src.ePrimary(self.prob)
if ep is not None:
ePrimary[:,i] = ep
return ePrimary
@@ -53,7 +53,7 @@ class FieldsFDEM_e(FieldsFDEM):
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self._edgeCurl.shape[0],eSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
bp = src.bPrimary(self.prob)
if bp is not None:
bPrimary[:,i] += bp
return bPrimary
@@ -63,7 +63,7 @@ class FieldsFDEM_e(FieldsFDEM):
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
S_m, _ = src.eval(self.survey.prob)
S_m, _ = src.eval(self.prob)
if S_m is not None:
b[:,i] += 1./(1j*omega(src.freq)) * S_m
return b
@@ -75,7 +75,7 @@ class FieldsFDEM_e(FieldsFDEM):
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
S_mDeriv, _ = src.evalDeriv(self.survey.prob, adjoint)
S_mDeriv, _ = src.evalDeriv(self.prob, adjoint)
S_mDeriv = S_mDeriv(v)
if S_mDeriv is not None:
return 1./(1j * omega(src.freq)) * S_mDeriv
@@ -108,6 +108,7 @@ class FieldsFDEM_b(FieldsFDEM):
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._MeSigmaI = self.survey.prob.MeSigmaI
self._MfMui = self.survey.prob.MfMui
@@ -116,7 +117,7 @@ class FieldsFDEM_b(FieldsFDEM):
def _bPrimary(self, bSolution, srcList):
bPrimary = np.zeros_like(bSolution)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
bp = src.bPrimary(self.prob)
if bp is not None:
bPrimary[:,i] = bp
return bPrimary
@@ -137,7 +138,7 @@ class FieldsFDEM_b(FieldsFDEM):
def _ePrimary(self, bSolution, srcList):
ePrimary = np.zeros([self._edgeCurl.shape[1],bSolution.shape[1]],dtype = complex)
for i,src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
ep = src.ePrimary(self.prob)
if ep is not None:
ePrimary[:,i] = ep
return ePrimary
@@ -145,9 +146,9 @@ class FieldsFDEM_b(FieldsFDEM):
def _eSecondary(self, bSolution, srcList):
e = self._MeSigmaI * ( self._edgeCurl.T * ( self._MfMui * bSolution))
for i,src in enumerate(srcList):
_,S_e = src.eval(self.survey.prob)
_,S_e = src.eval(self.prob)
if S_e is not None:
e += -self._MeSigmaI*S_e
e[:,i] += -self._MeSigmaI*S_e
return e
def _eSecondaryDeriv_u(self, src, v, adjoint=False):
@@ -158,18 +159,18 @@ class FieldsFDEM_b(FieldsFDEM):
def _eSecondaryDeriv_m(self, src, v, adjoint=False):
bSolution = self[[src],'bSolution']
_,S_e = src.eval(self.survey.prob)
_,S_e = src.eval(self.prob)
w = self._edgeCurl.T * (self._MfMui * bSolution)
if S_e is not None:
w += -S_e
w += -Utils.mkvc(S_e,2)
if not adjoint:
de_dm = self._MeSigmaIDeriv(w) * v
elif adjoint:
de_dm = self._MeSigmaIDeriv(w).T * v
_, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
_, S_eDeriv = src.evalDeriv(self.prob, adjoint)
Se_Deriv = S_eDeriv(v)
if Se_Deriv is not None:
@@ -203,6 +204,7 @@ class FieldsFDEM_j(FieldsFDEM):
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._MeMuI = self.survey.prob.MeMuI
self._MfRho = self.survey.prob.MfRho
@@ -211,7 +213,7 @@ class FieldsFDEM_j(FieldsFDEM):
def _jPrimary(self, jSolution, srcList):
jPrimary = np.zeros_like(jSolution,dtype = complex)
for i, src in enumerate(srcList):
jp = src.jPrimary(self.survey.prob)
jp = src.jPrimary(self.prob)
if jp is not None:
jPrimary[:,i] += jp
return jPrimary
@@ -232,7 +234,7 @@ class FieldsFDEM_j(FieldsFDEM):
def _hPrimary(self, jSolution, srcList):
hPrimary = np.zeros([self._edgeCurl.shape[1],jSolution.shape[1]],dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.survey.prob)
hp = src.hPrimary(self.prob)
if hp is not None:
hPrimary[:,i] = hp
return hPrimary
@@ -244,7 +246,7 @@ class FieldsFDEM_j(FieldsFDEM):
h = MeMuI * (C.T * (MfRho * jSolution) )
for i, src in enumerate(srcList):
h[:,i] *= -1./(1j*omega(src.freq))
S_m,_ = src.eval(self.survey.prob)
S_m,_ = src.eval(self.prob)
if S_m is not None:
h[:,i] += 1./(1j*omega(src.freq)) * MeMuI * S_m
return h
@@ -270,7 +272,7 @@ class FieldsFDEM_j(FieldsFDEM):
elif adjoint:
hDeriv_m = -1./(1j*omega(src.freq)) * MfRhoDeriv(jSolution).T * ( C * (MeMuI.T * v ) )
S_mDeriv,_ = src.evalDeriv(self.survey.prob, adjoint)
S_mDeriv,_ = src.evalDeriv(self.prob, adjoint)
if not adjoint:
S_mDeriv = S_mDeriv(v)
@@ -309,6 +311,7 @@ class FieldsFDEM_h(FieldsFDEM):
FieldsFDEM.__init__(self,mesh,survey,**kwargs)
def startup(self):
self.prob = self.survey.prob
self._edgeCurl = self.survey.prob.mesh.edgeCurl
self._MeMuI = self.survey.prob.MeMuI
self._MfRho = self.survey.prob.MfRho
@@ -316,7 +319,7 @@ class FieldsFDEM_h(FieldsFDEM):
def _hPrimary(self, hSolution, srcList):
hPrimary = np.zeros_like(hSolution,dtype = complex)
for i, src in enumerate(srcList):
hp = src.hPrimary(self.survey.prob)
hp = src.hPrimary(self.prob)
if hp is not None:
hPrimary[:,i] += hp
return hPrimary
@@ -337,7 +340,7 @@ class FieldsFDEM_h(FieldsFDEM):
def _jPrimary(self, hSolution, srcList):
jPrimary = np.zeros([self._edgeCurl.shape[0], hSolution.shape[1]])
for i, src in enumerate(srcList):
jp = src.jPrimary(self.survey.prob)
jp = src.jPrimary(self.prob)
if jp is not None:
jPrimary[:,i] = jp
return jPrimary
@@ -345,7 +348,7 @@ class FieldsFDEM_h(FieldsFDEM):
def _jSecondary(self, hSolution, srcList):
j = self._edgeCurl*hSolution
for i, src in enumerate(srcList):
_,S_e = src.eval(self.survey.prob)
_,S_e = src.eval(self.prob)
if S_e is not None:
j[:,i] += -S_e
return j
@@ -357,7 +360,7 @@ class FieldsFDEM_h(FieldsFDEM):
return self._edgeCurl.T*v
def _jSecondaryDeriv_m(self, src, v, adjoint=False):
_,S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
_,S_eDeriv = src.evalDeriv(self.prob, adjoint)
S_eDeriv = S_eDeriv(v)
if S_eDeriv is not None:
return -S_eDeriv
simpegEM/FDEM/SurveyFDEM.py
+72 -18
View File
@@ -1,7 +1,7 @@
from SimPEG import Survey, Problem, Utils, np, sp
from simpegEM.Utils import SrcUtils
from simpegEM.Utils.EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
from scipy.constants import mu_0
####################################################
# Receivers
@@ -189,11 +189,12 @@ class SrcFDEM_RawVec(SrcFDEM):
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.):
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):
@@ -216,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')
a = SrcUtils.MagneticDipoleVectorPotential(self.loc, gridY, 'y', mu=self.mu, moment=self.moment)
else:
srcfct = SrcUtils.MagneticDipoleVectorPotential
ax = srcfct(self.loc, gridX, 'x')
ay = srcfct(self.loc, gridY, 'y')
az = srcfct(self.loc, gridZ, 'z')
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
@@ -235,17 +236,34 @@ class SrcFDEM_MagDipole(SrcFDEM):
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.):
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):
@@ -268,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')
bz = srcfct(self.loc, gridZ, 'z')
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')
by = srcfct(self.loc, gridY, 'y')
bz = srcfct(self.loc, gridZ, 'z')
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
@@ -287,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):
@@ -316,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
+56 -13
View File
@@ -17,9 +17,11 @@ TOL = 1e-4
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = [1e-1, 2e-1]
freq = 1e-1
addrandoms = True
SrcType = 'MagDipole' #or 'MAgDipole_Bfield', 'CircularLoop', 'RawVec'
def getProblem(fdemType, comp):
cs = 5.
@@ -35,23 +37,61 @@ 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_MagDipole([Rx0], freq=freq[0], loc=np.r_[0.,0.,0.])
Src1 = EM.FDEM.SrcFDEM_MagDipole([Rx0], freq=freq[1], loc=np.r_[0.,0.,0.])
survey = EM.FDEM.SurveyFDEM([Src0, Src1])
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)
@@ -69,15 +109,15 @@ def adjointTest(fdemType, comp):
print 'Adjoint %s formulation - %s' % (fdemType, comp)
m = np.log(np.ones(prb.mesh.nC)*CONDUCTIVITY)
mu = np.log(np.ones(prb.mesh.nC))*MU
mu = np.ones(prb.mesh.nC)*MU
if addrandoms is True:
m = m + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
m = m + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
prb.mu = mu
survey = prb.survey
prb.PropMap.PropModel.mu = mu
prb.PropMap.PropModel.mui = 1./mu
u = prb.fields(m)
v = np.random.rand(survey.nD)
@@ -99,10 +139,12 @@ def derivTest(fdemType, comp):
mu = np.log(np.ones(prb.mesh.nC)*MU)
if addrandoms is True:
x0 = x0 + np.random.randn(prb.mesh.nC)*CONDUCTIVITY*1e-1
x0 = x0 + np.random.randn(prb.mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(prb.mesh.nC)*MU*1e-1
prb.mu = mu
prb.PropMap.PropModel.mu = mu
prb.PropMap.PropModel.mui = 1./mu
survey = prb.survey
def fun(x):
return survey.dpred(x), lambda x: prb.Jvec(x0, x)
@@ -120,10 +162,11 @@ def crossCheckTest(fdemType, comp):
mu = np.log(np.ones(mesh.nC)*MU)
if addrandoms is True:
m = m + np.random.randn(mesh.nC)*CONDUCTIVITY*1e-1
m = m + np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu = mu + np.random.randn(mesh.nC)*MU*1e-1
prb1.mu = mu
prb1.PropMap.PropModel.mu = mu
prb1.PropMap.PropModel.mui = 1./mu
survey1 = prb1.survey
d1 = survey1.dpred(m)
simpegEM/Utils/SrcUtils.py
+13 -12
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.)):
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):
@@ -48,13 +49,13 @@ def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0
dR = obsLoc - srcLoc[i, np.newaxis].repeat(nEdges, axis=0)
mCr = np.cross(m, dR)
r = np.sqrt((dR**2).sum(axis=1))
A[:, i] = +(mu_0/(4*pi)) * mCr[:,dimInd]/(r**3)
A[:, i] = +(mu/(4*pi)) * mCr[:,dimInd]/(r**3)
if nSrc == 1:
return A.flatten()
return A
def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
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.):
: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,22 +79,22 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
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)
r = np.sqrt((dR**2).sum(axis=1))
if dimInd == 0:
B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,0]/r**2)
elif dimInd == 1:
B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]*dR[:,1]/r**2)
elif dimInd == 2:
B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
B[:, i] = +(mu/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
else:
raise Exception("Not Implemented")
if nSrc == 1:
@@ -102,7 +103,7 @@ def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
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)