wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-08-04 13:13:40 +08:00
Code Issues Packages Projects Releases Wiki Activity

tx -> src

Browse Source
This commit is contained in:
Lindsey
2015-04-17 16:41:54 -07:00
parent 08d90bbb67
commit dc7cc1c716
18 changed files with 248 additions and 238 deletions
Download Patch File Download Diff File Expand all files Collapse all files
docs/examples/FDEM_b.py
+7 -7
View File
@@ -19,9 +19,9 @@ model = Model.LogModel(mesh)
x = np.linspace(-10,10,5)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
rxList = EM.FDEM.RxListFDEM(XYZ, 'Ex')
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, rxList)
survey = EM.FDEM.SurveyFDEM([Tx0])
survey = EM.FDEM.SurveyFDEM([Src0])
prb = EM.FDEM.ProblemFDEM_b(model)
prb.pair(survey)
@@ -31,26 +31,26 @@ sigma = np.ones(mesh.nC)*sig
sigma[mesh.gridCC[:,2] > 0] = 1e-8
m = np.log(sigma)
skin = 500*np.sqrt(1/(sig*Tx0.freq))
skin = 500*np.sqrt(1/(sig*Src0.freq))
print 'The skin depth is: %4.2f m' % skin
prb.Solver = Utils.SolverUtils.DSolverWrap(sp.linalg.spsolve, factorize=False, checkAccuracy=True)
u = prb.fields(m)
plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Tx0, 'b'].real)), 'Fz'))
plt.colorbar(mesh.plotImage(np.log10(np.abs(u[Src0, 'b'].real)), 'Fz'))
bfz = mesh.r(u[Tx0, 'b'],'F','Fz','M')
bfz = mesh.r(u[Src0, 'b'],'F','Fz','M')
x = np.linspace(-55,55,12)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
P = mesh.getInterpolationMat(XYZ, 'Fz')
an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Tx0.freq, sig)
an = EM.Utils.Ana.FEM.hzAnalyticDipoleF(x, Src0.freq, sig)
plt.figure(2)
plt.plot(x,np.log10(np.abs(P*np.imag(u[Tx0, 'b']))))
plt.plot(x,np.log10(np.abs(P*np.imag(u[Src0, 'b']))))
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
plt.xlabel('Distance, m')
plt.ylabel('Log10 Response imag($B_z$)')
simpegEM/Analytics/FDEM.py
+4 -4
View File
@@ -38,7 +38,7 @@ def hzAnalyticDipoleF(r, freq, sigma, secondary=True):
return hz
def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'):
def AnalyticMagDipoleWholeSpace(XYZ, srcLoc, sig, f, m=1., orientation='X'):
"""
Analytical solution for a dipole in a whole-space.
@@ -67,9 +67,9 @@ def AnalyticMagDipoleWholeSpace(XYZ, txLoc, sig, f, m=1., orientation='X'):
XYZ = Utils.asArray_N_x_Dim(XYZ, 3)
dx = XYZ[:,0]-txLoc[0]
dy = XYZ[:,1]-txLoc[1]
dz = XYZ[:,2]-txLoc[2]
dx = XYZ[:,0]-srcLoc[0]
dy = XYZ[:,1]-srcLoc[1]
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 )
simpegEM/Base.py
-4
View File
@@ -113,10 +113,6 @@ class BaseEMProblem(Problem.BaseProblem):
sigma = self.curModel.transform
self._MfSigmai = self.mesh.getFaceInnerProduct(1/sigma)
return self._MfSigmai
@property
def dMfSigmai_dsig(self):
return self._dMfSigmai_dsig
deleteTheseOnModelUpdate = ['_MeSigma', '_MeSigmaI','_MfSigmai']
simpegEM/Examples/CylInversion.py
+2 -2
View File
@@ -36,8 +36,8 @@ if plotIt:
rxOffset=1e-3
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 30]]), np.logspace(-5,-3, 31), 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 80]), 'VMD_MVP', [rx])
survey = EM.TDEM.SurveyTDEM([tx])
src = EM.TDEM.SrcTDEM(np.array([0., 0., 80]), 'VMD_MVP', [rx])
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
prb.Solver = SolverLU
simpegEM/FDEM/FDEM.py
+28 -28
View File
@@ -28,8 +28,8 @@ class BaseFDEMProblem(BaseEMProblem):
rhs = RHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
sol = Ainv * rhs
Txs = self.survey.getTransmitters(freq)
F[Txs, self._fieldType] = sol
Srcs = self.survey.getSources(freq)
F[Srcs, self._fieldType] = sol
return F
@@ -45,19 +45,19 @@ class BaseFDEMProblem(BaseEMProblem):
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts)
for tx in self.survey.getTransmitters(freq):
u_tx = u[tx, self.solType]
w = self.getADeriv(freq, u_tx, v)
for src in self.survey.getSources(freq):
u_src = u[src, self.solType]
w = self.getADeriv(freq, u_src, v)
Ainvw = Ainv * w
for rx in tx.rxList:
for rx in src.rxList:
fAinvw = self.calcFields(Ainvw, freq, rx.projField)
P = lambda v: rx.projectFieldsDeriv(tx, self.mesh, u, v)
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v)
Jv[tx, rx] = - P(fAinvw)
Jv[src, rx] = - P(fAinvw)
df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, v)
df_dm = self.calcFieldsDeriv(u_src, freq, rx.projField, v)
if df_dm is not None:
Jv[tx, rx] += P(df_dm)
Jv[src, rx] += P(df_dm)
return Utils.mkvc(Jv)
@@ -77,17 +77,17 @@ class BaseFDEMProblem(BaseEMProblem):
AT = self.getA(freq).T
ATinv = self.Solver(AT, **self.solverOpts)
for tx in self.survey.getTransmitters(freq):
u_tx = u[tx, self.solType]
for src in self.survey.getSources(freq):
u_src = u[src, self.solType]
for rx in tx.rxList:
PTv = rx.projectFieldsDeriv(tx, self.mesh, u, v[tx, rx], adjoint=True)
for rx in src.rxList:
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True)
fPTv = self.calcFields(PTv, freq, rx.projField, adjoint=True)
w = ATinv * fPTv
Jtv_rx = - self.getADeriv(freq, u_tx, w, adjoint=True)
Jtv_rx = - self.getADeriv(freq, u_src, w, adjoint=True)
df_dm = self.calcFieldsDeriv(u_tx, freq, rx.projField, PTv, adjoint=True)
df_dm = self.calcFieldsDeriv(u_src, freq, rx.projField, PTv, adjoint=True)
if df_dm is not None:
Jtv_rx += df_dm
@@ -105,19 +105,19 @@ class BaseFDEMProblem(BaseEMProblem):
def getSource(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE or nF, nTx)
:rtype: numpy.ndarray (nE or nF, nSrc)
:return: RHS
"""
Txs = self.survey.getTransmitters(freq)
Srcs = self.survey.getSources(freq)
if self._eqLocs is 'FE':
S_m = 1j*np.zeros((self.mesh.nF,len(Txs)))
S_e = 1j*np.zeros((self.mesh.nE,len(Txs)))
S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
elif self._eqLocs is 'EF':
S_m = 1j*np.zeros((self.mesh.nE,len(Txs)))
S_e = 1j*np.zeros((self.mesh.nF,len(Txs)))
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
for i, tx in enumerate(Txs):
smi, sei = tx.getSource(self)
for i, src in enumerate(Srcs):
smi, sei = src.getSource(self)
if smi is not None:
S_m[:,i] = smi
if sei is not None:
@@ -185,7 +185,7 @@ class ProblemFDEM_e(BaseFDEMProblem):
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
@@ -250,7 +250,7 @@ class ProblemFDEM_b(BaseFDEMProblem):
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
@@ -358,7 +358,7 @@ class ProblemFDEM_j(BaseFDEMProblem):
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
@@ -441,7 +441,7 @@ class ProblemFDEM_h(BaseFDEMProblem):
def getRHS(self, freq):
"""
:param float freq: Frequency
:rtype: numpy.ndarray (nE, nTx)
:rtype: numpy.ndarray (nE, nSrc)
:return: RHS
"""
simpegEM/FDEM/SurveyFDEM.py
+67 -53
View File
@@ -3,6 +3,10 @@ from simpegEM import Sources
from simpegEM.Utils.EMUtils import omega
####################################################
# Receivers
####################################################
class RxFDEM(Survey.BaseRx):
knownRxTypes = {
@@ -54,15 +58,15 @@ class RxFDEM(Survey.BaseRx):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][2]
def projectFields(self, tx, mesh, u):
def projectFields(self, src, mesh, u):
P = self.getP(mesh)
u_part_complex = u[tx, self.projField]
u_part_complex = u[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
def projectFieldsDeriv(self, tx, mesh, u, v, adjoint=False):
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
P = self.getP(mesh)
if not adjoint:
@@ -82,26 +86,36 @@ class RxFDEM(Survey.BaseRx):
return Pv
# SrcFDEM
class TxFDEM(Survey.BaseTx):
####################################################
# 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)
rxPair = RxFDEM
knownTxTypes = ['VMD', 'VMD_B', 'CircularLoop', 'Simple']
knownSrcTypes = ['VMD', 'VMD_B', 'CircularLoop', 'Simple']
radius = None
def __init__(self, loc, txType, freq, rxList):
def __init__(self, loc, srcType, freq, rxList):
self.freq = float(freq)
Survey.BaseTx.__init__(self, loc, txType, rxList)
Survey.BaseSrc.__init__(self, loc, srcType, rxList)
def getSource(self, prob):
tx = self
freq = tx.freq
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':
@@ -141,42 +155,42 @@ class TxFDEM(Survey.BaseTx):
if not prob.mesh.isSymmetric:
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
if tx.txType == 'VMD':
SRC = Sources.MagneticDipoleVectorPotential(tx.loc, gridEJy, 'y')
elif tx.txType == 'CircularLoop':
SRC = Sources.MagneticLoopVectorPotential(tx.loc, gridEJy, 'y', tx.radius)
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 tx.txType == 'VMD':
src = Sources.MagneticDipoleVectorPotential
SRCx = src(tx.loc, gridEJx, 'x')
SRCy = src(tx.loc, gridEJy, 'y')
SRCz = src(tx.loc, gridEJz, 'z')
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 tx.txType == 'VMD_B':
src = Sources.MagneticDipoleFields
SRCx = src(tx.loc, gridBHx, 'x')
SRCy = src(tx.loc, gridBHy, 'y')
SRCz = src(tx.loc, gridBHz, '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 tx.txType == 'CircularLoop':
src = Sources.MagneticLoopVectorPotential
SRCx = src(tx.loc, gridEJx, 'x', tx.radius)
SRCy = src(tx.loc, gridEJy, 'y', tx.radius)
SRCz = src(tx.loc, gridEJz, 'z', tx.radius)
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 txType is not implemented' % tx.txType)
raise NotImplemented('%s srcType is not implemented' % src.srcType)
SRC = np.concatenate((SRCx, SRCy, SRCz))
else:
raise Exception('Unknown mesh for VMD')
# b-forumlation
if tx.txType == 'VMD_B':
if src.srcType == 'VMD_B':
b_0 = SRC
else:
a = SRC
@@ -184,23 +198,23 @@ class TxFDEM(Survey.BaseTx):
return -1j*omega(freq)*b_0, None
class SimpleTxFDEM_g(TxFDEM):
class SimpleSrcFDEM_e(SrcFDEM):
def __init__(self, vec, freq, rxList):
self.vec = vec
self.freq = float(freq)
TxFDEM.__init__(self, None, 'Simple', freq, rxList)
SrcFDEM.__init__(self, None, 'Simple', freq, rxList)
def getSource(self, prob):
return None, self.vec
class SimpleTxFDEM_m(TxFDEM):
class SimpleSrcFDEM_m(SrcFDEM):
def __init__(self, vec, freq, rxList):
self.vec = vec
self.freq = float(freq)
TxFDEM.__init__(self, None, 'Simple', freq, rxList)
SrcFDEM.__init__(self, None, 'Simple', freq, rxList)
def getSource(self, prob):
return self.vec, None
@@ -211,18 +225,18 @@ class SurveyFDEM(Survey.BaseSurvey):
docstring for SurveyFDEM
"""
txPair = TxFDEM
srcPair = SrcFDEM
def __init__(self, txList, **kwargs):
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.txList = txList
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
_freqDict = {}
for tx in txList:
if tx.freq not in _freqDict:
_freqDict[tx.freq] = []
_freqDict[tx.freq] += [tx]
for src in srcList:
if src.freq not in _freqDict:
_freqDict[src.freq] = []
_freqDict[src.freq] += [src]
self._freqDict = _freqDict
self._freqs = sorted([f for f in self._freqDict])
@@ -238,24 +252,24 @@ class SurveyFDEM(Survey.BaseSurvey):
return len(self._freqDict)
@property
def nTxByFreq(self):
if getattr(self, '_nTxByFreq', None) is None:
self._nTxByFreq = {}
def nSrcByFreq(self):
if getattr(self, '_nSrcByFreq', None) is None:
self._nSrcByFreq = {}
for freq in self.freqs:
self._nTxByFreq[freq] = len(self.getTransmitters(freq))
return self._nTxByFreq
self._nSrcByFreq[freq] = len(self.getSources(freq))
return self._nSrcByFreq
def getTransmitters(self, freq):
"""Returns the transmitters associated with a specific frequency."""
def getSources(self, freq):
"""Returns the sources associated with a specific frequency."""
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
data = Survey.Data(self)
for tx in self.txList:
for rx in tx.rxList:
data[tx, rx] = rx.projectFields(tx, self.mesh, u)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, u)
return data
def projectFieldsDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
raise Exception('Use Sources to project fields deriv.')
simpegEM/Sources/CircularLoop.py
+18 -18
View File
@@ -1,12 +1,12 @@
from SimPEG import *
from scipy.special import ellipk, ellipe
def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius):
def MagneticLoopVectorPotential(srcLoc, obsLoc, component, radius):
"""
Calculate the vector potential of horizontal circular loop
at given locations
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
: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
@@ -18,33 +18,33 @@ def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius):
if type(component) in [list, tuple]:
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticLoopVectorPotential(txLoc, obsLoc, comp, radius)
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(txLoc, getattr(mesh,'grid'+component), component[1], radius)
return MagneticLoopVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], radius)
txLoc = np.atleast_2d(txLoc)
srcLoc = np.atleast_2d(srcLoc)
obsLoc = np.atleast_2d(obsLoc)
n = obsLoc.shape[0]
nTx = txLoc.shape[0]
nSrc = srcLoc.shape[0]
if component=='z':
A = np.zeros((n, nTx))
if nTx ==1:
A = np.zeros((n, nSrc))
if nSrc ==1:
return A.flatten()
return A
else:
A = np.zeros((n, nTx))
for i in range (nTx):
x = obsLoc[:, 0] - txLoc[i, 0]
y = obsLoc[:, 1] - txLoc[i, 1]
z = obsLoc[:, 2] - txLoc[i, 2]
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.
@@ -64,7 +64,7 @@ def MagneticLoopVectorPotential(txLoc, obsLoc, component, radius):
else:
raise ValueError('Invalid component')
if nTx == 1:
if nSrc == 1:
return A.flatten()
return A
@@ -77,10 +77,10 @@ if __name__ == '__main__':
hy = np.ones(ncy)*cs
hz = np.ones(ncz)*cs
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCC')
txLoc = np.r_[0., 0., 0.]
Ax = MagneticLoopVectorPotential(txLoc, mesh.gridEx, 'x', 200)
Ay = MagneticLoopVectorPotential(txLoc, mesh.gridEy, 'y', 200)
Az = MagneticLoopVectorPotential(txLoc, mesh.gridEz, 'z', 200)
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
simpegEM/Sources/magneticDipole.py
+18 -18
View File
@@ -2,12 +2,12 @@ import numpy as np
from scipy.constants import mu_0, pi
from SimPEG import Mesh
def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
def MagneticDipoleVectorPotential(srcLoc, obsLoc, component, dipoleMoment=(0., 0., 1.)):
"""
Calculate the vector potential of a set of magnetic dipoles
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
: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 dipoleMoment: The vector dipole moment
@@ -18,13 +18,13 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0.
if type(component) in [list, tuple]:
out = range(len(component))
for i, comp in enumerate(component):
out[i] = MagneticDipoleVectorPotential(txLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
out[i] = MagneticDipoleVectorPotential(srcLoc, obsLoc, comp, dipoleMoment=dipoleMoment)
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 MagneticDipoleVectorPotential(txLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment)
return MagneticDipoleVectorPotential(srcLoc, getattr(mesh,'grid'+component), component[1], dipoleMoment=dipoleMoment)
if component == 'x':
dimInd = 0
@@ -35,30 +35,30 @@ def MagneticDipoleVectorPotential(txLoc, obsLoc, component, dipoleMoment=(0., 0.
else:
raise ValueError('Invalid component')
txLoc = np.atleast_2d(txLoc)
srcLoc = np.atleast_2d(srcLoc)
obsLoc = np.atleast_2d(obsLoc)
dipoleMoment = np.atleast_2d(dipoleMoment)
nEdges = obsLoc.shape[0]
nTx = txLoc.shape[0]
nSrc = srcLoc.shape[0]
m = np.array(dipoleMoment).repeat(nEdges, axis=0)
A = np.empty((nEdges, nTx))
for i in range(nTx):
dR = obsLoc - txLoc[i, np.newaxis].repeat(nEdges, axis=0)
A = np.empty((nEdges, nSrc))
for i in range(nSrc):
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)
if nTx == 1:
if nSrc == 1:
return A.flatten()
return A
def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.):
def MagneticDipoleFields(srcLoc, obsLoc, component, dipoleMoment=1.):
"""
Calculate the vector potential of a set of magnetic dipoles
at given locations 'ref. <http://en.wikipedia.org/wiki/Dipole#Magnetic_vector_potential>'
:param numpy.ndarray txLoc: Location of the transmitter(s) (x, y, z)
: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)
@@ -75,17 +75,17 @@ def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.):
else:
raise ValueError('Invalid component')
txLoc = np.atleast_2d(txLoc)
srcLoc = np.atleast_2d(srcLoc)
obsLoc = np.atleast_2d(obsLoc)
dipoleMoment = np.atleast_2d(dipoleMoment)
nFaces = obsLoc.shape[0]
nTx = txLoc.shape[0]
nSrc = srcLoc.shape[0]
m = np.array(dipoleMoment).repeat(nFaces, axis=0)
B = np.empty((nFaces, nTx))
for i in range(nTx):
dR = obsLoc - txLoc[i, np.newaxis].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)
@@ -95,6 +95,6 @@ def MagneticDipoleFields(txLoc, obsLoc, component, dipoleMoment=1.):
B[:, i] = +(mu_0/(4*pi)) /(r**3) * (3*dR[:,2]**2/r**2-1)
else:
raise Exception("Not Implemented")
if nTx == 1:
if nSrc == 1:
return B.flatten()
return B
simpegEM/TDEM/BaseTDEM.py
+6 -6
View File
@@ -13,19 +13,19 @@ class FieldsTDEM(Problem.TimeFields):
knownFields = {'b': 'F', 'e': 'E'}
def tovec(self):
nTx, nF, nE = self.survey.nTx, self.mesh.nF, self.mesh.nE
u = np.empty(0 if nTx == 1 else (0, nTx))
nSrc, nF, nE = self.survey.nSrc, self.mesh.nF, self.mesh.nE
u = np.empty(0 if nSrc == 1 else (0, nSrc))
for i in range(self.survey.prob.nT):
if 'b' in self:
b = self[:,'b',i+1]
else:
b = np.zeros(nF if nTx == 1 else (nF, nTx))
b = np.zeros(nF if nSrc == 1 else (nF, nSrc))
if 'e' in self:
e = self[:,'e',i+1]
else:
e = np.zeros(nE if nTx == 1 else (nE, nTx))
e = np.zeros(nE if nSrc == 1 else (nE, nSrc))
u = np.concatenate((u, b, e))
return Utils.mkvc(u)
@@ -42,9 +42,9 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
self.curModel = m
# Create a fields storage object
F = self._FieldsForward_pair(self.mesh, self.survey)
for tx in self.survey.txList:
for src in self.survey.srcList:
# Set the initial conditions
F[tx,:,0] = tx.getInitialFields(self.mesh)
F[src,:,0] = src.getInitialFields(self.mesh)
F = self.forward(m, self.getRHS, F=F)
if self.verbose: print '%s\nDone calculating fields(m)\n%s'%('*'*50,'*'*50)
return F
simpegEM/TDEM/SurveyTDEM.py
+22 -22
View File
@@ -51,27 +51,27 @@ class RxTDEM(Survey.BaseTimeRx):
else:
return timeMesh.getInterpolationMat(self.times, self.projTLoc)
def projectFields(self, tx, mesh, timeMesh, u):
def projectFields(self, src, mesh, timeMesh, u):
P = self.getP(mesh, timeMesh)
u_part = Utils.mkvc(u[tx, self.projField, :])
u_part = Utils.mkvc(u[src, self.projField, :])
return P*u_part
def projectFieldsDeriv(self, tx, mesh, timeMesh, u, v, adjoint=False):
def projectFieldsDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
P = self.getP(mesh, timeMesh)
if not adjoint:
return P * Utils.mkvc(v[tx, self.projField, :])
return P * Utils.mkvc(v[src, self.projField, :])
elif adjoint:
return P.T * v[tx, self]
return P.T * v[src, self]
class TxTDEM(Survey.BaseTx):
class SrcTDEM(Survey.BaseSrc):
rxPair = RxTDEM
radius = None
knownTxTypes = ['VMD_MVP', 'CircularLoop_MVP']
knownSrcTypes = ['VMD_MVP', 'CircularLoop_MVP']
def getInitialFields(self, mesh):
F0 = getattr(self, '_getInitialFields_' + self.txType)(mesh)
F0 = getattr(self, '_getInitialFields_' + self.srcType)(mesh)
return F0
def _getInitialFields_VMD_MVP(self, mesh):
@@ -109,18 +109,18 @@ class SurveyTDEM(Survey.BaseSurvey):
"""
docstring for SurveyTDEM
"""
txPair = TxTDEM
srcPair = SrcTDEM
def __init__(self, txList, **kwargs):
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.txList = txList
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
def projectFields(self, u):
data = Survey.Data(self)
for tx in self.txList:
for rx in tx.rxList:
data[tx, rx] = rx.projectFields(tx, self.mesh, self.prob.timeMesh, u)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
return data
def projectFieldsDeriv(self, u, v=None, adjoint=False):
@@ -128,20 +128,20 @@ class SurveyTDEM(Survey.BaseSurvey):
if not adjoint:
data = Survey.Data(self)
for tx in self.txList:
for rx in tx.rxList:
data[tx, rx] = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
return data
else:
f = FieldsTDEM(self.mesh, self)
for tx in self.txList:
for rx in tx.rxList:
Ptv = rx.projectFieldsDeriv(tx, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
for src in self.srcList:
for rx in src.rxList:
Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
Ptv = Ptv.reshape((-1, self.prob.timeMesh.nN), order='F')
if rx.projField not in f: # first time we are projecting
f[tx, rx.projField, :] = Ptv
f[src, rx.projField, :] = Ptv
else: # there are already fields, so let's add to them!
f[tx, rx.projField, :] += Ptv
f[src, rx.projField, :] += Ptv
return f
simpegEM/TDEM/TDEM_b.py
+14 -14
View File
@@ -14,7 +14,7 @@ class FieldsTDEM_e_from_b(FieldsTDEM):
self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
self.MfMui = self.survey.prob.MfMui
def e_from_b(self, b, txInd, timeInd):
def e_from_b(self, b, srcInd, timeInd):
# TODO: implement non-zero js
return self.MeSigmaI*(self.edgeCurlT*(self.MfMui*b))
@@ -32,10 +32,10 @@ class FieldsTDEM_e_from_b_Ah(FieldsTDEM):
self.edgeCurlT = self.survey.prob.mesh.edgeCurl.T
self.MfMui = self.survey.prob.MfMui
def e_from_b(self, y_b, txInd, tInd):
def e_from_b(self, y_b, srcInd, tInd):
y_e = self.MeSigmaI*(self.edgeCurlT*(self.MfMui*y_b))
if 'e' in self.p:
y_e = y_e - self.MeSigmaI*self.p[txInd,'e',tInd]
y_e = y_e - self.MeSigmaI*self.p[srcInd,'e',tInd]
return y_e
class ProblemTDEM_b(BaseTDEMProblem):
@@ -73,7 +73,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
def getRHS(self, tInd, F):
dt = self.timeSteps[tInd]
B_n = np.c_[[F[tx,'b',tInd] for tx in self.survey.txList]].T
B_n = np.c_[[F[src,'b',tInd] for src in self.survey.srcList]].T
RHS = (1.0/dt)*self.MfMui*B_n
return RHS
@@ -95,7 +95,7 @@ class ProblemTDEM_b(BaseTDEMProblem):
u = self.fields(m)
self.curModel = m
# Note: Fields has shape (nF/E, nTx, nT+1)
# Note: Fields has shape (nF/E, nSrc, nT+1)
# However, p will only really fill (:,:,1:nT+1)
# meaning the 'initial fields' are zero (:,:,0)
p = FieldsTDEM(self.mesh, self.survey)
@@ -112,9 +112,9 @@ class ProblemTDEM_b(BaseTDEMProblem):
# TODO: G[1] may be dependent on the model
# for a galvanic source (deriv of the dc problem)
#
# Do multiplication for all tx in self.survey.txList
for tx in self.survey.txList:
p[tx, 'e', i] = - dMdsig(u[tx,'e',i]) * dsigdm_x_v
# Do multiplication for all src in self.survey.srcList
for src in self.survey.srcList:
p[src, 'e', i] = - dMdsig(u[src,'e',i]) * dsigdm_x_v
return p
def Gtvec(self, m, vec, u=None):
@@ -133,14 +133,14 @@ class ProblemTDEM_b(BaseTDEMProblem):
dMdsig = self.mesh.getEdgeInnerProductDeriv(self.curModel.transform)
dsigdm = self.curModel.transformDeriv
nTx = self.survey.nTx
nSrc = self.survey.nSrc
VUs = None
# Here we can do internal multiplications of Gt*v and then multiply by MsigDeriv.T in one go.
for i in range(1,self.nT+1):
vu = None
for tx in self.survey.txList:
vutx = dMdsig(u[tx,'e',i]).T * vec[tx,'e',i]
vu = vutx if vu is None else vu + vutx
for src in self.survey.srcList:
vusrc = dMdsig(u[src,'e',i]).T * vec[src,'e',i]
vu = vusrc if vu is None else vu + vusrc
VUs = vu if VUs is None else VUs + vu
p = -dsigdm.T*VUs
return p
@@ -241,8 +241,8 @@ class ProblemTDEM_b(BaseTDEMProblem):
# 1 (tInd=1 uses fields 2 and 3)
def AhtRHS(tInd, y):
nTx, nF = self.survey.nTx, self.mesh.nF
rhs = np.zeros(nF if nTx == 1 else (nF, nTx))
nSrc, nF = self.survey.nSrc, self.mesh.nF
rhs = np.zeros(nF if nSrc == 1 else (nF, nSrc))
if 'e' in p:
rhs += self.MfMui*(self.mesh.edgeCurl*(self.MeSigmaI*p[:,'e',tInd+1]))
simpegEM/TDEM/__init__.py
+1 -1
View File
@@ -1,3 +1,3 @@
from SurveyTDEM import SurveyTDEM, RxTDEM, TxTDEM
from SurveyTDEM import SurveyTDEM, RxTDEM, SrcTDEM
from BaseTDEM import BaseTDEMProblem, FieldsTDEM
from TDEM_b import ProblemTDEM_b
simpegEM/Tests/test_FDEM.py
+3 -3
View File
@@ -32,12 +32,12 @@ def getProblem(fdemType, comp):
mapping = Maps.ExpMap(mesh)
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the transmitter
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)
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', freq, [Rx0])
survey = EM.FDEM.SurveyFDEM([Tx0])
survey = EM.FDEM.SurveyFDEM([Src0])
if verbose:
simpegEM/Tests/test_FDEM_analytics.py
+7 -7
View File
@@ -22,9 +22,9 @@ 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')
Tx0 = EM.FDEM.TxFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, [rxList])
Src0 = EM.FDEM.SrcFDEM(np.r_[0.,0.,0.], 'VMD', 1e2, [rxList])
survey = EM.FDEM.SurveyFDEM([Tx0])
survey = EM.FDEM.SurveyFDEM([Src0])
prb = EM.FDEM.ProblemFDEM_b(mesh, mapping=mapping)
prb.pair(survey)
@@ -43,7 +43,7 @@ class FDEM_analyticTests(unittest.TestCase):
self.prb = prb
self.mesh = mesh
self.m = m
self.Tx0 = Tx0
self.Src0 = Src0
self.sig = sig
def test_Transect(self):
@@ -51,20 +51,20 @@ class FDEM_analyticTests(unittest.TestCase):
u = self.prb.fields(self.m)
bfz = self.mesh.r(u[self.Tx0, 'b'],'F','Fz','M')
bfz = self.mesh.r(u[self.Src0, 'b'],'F','Fz','M')
x = np.linspace(-55,55,12)
XYZ = Utils.ndgrid(x,np.r_[0],np.r_[0])
P = self.mesh.getInterpolationMat(XYZ, 'Fz')
an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Tx0.freq, self.sig)
an = EM.Analytics.FDEM.hzAnalyticDipoleF(x, self.Src0.freq, self.sig)
diff = np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']) - mu_0*np.imag(an)))
diff = np.log10(np.abs(P*np.imag(u[self.Src0, 'b']) - mu_0*np.imag(an)))
if plotIt:
import matplotlib.pyplot as plt
plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Tx0, 'b']))))
plt.plot(x,np.log10(np.abs(P*np.imag(u[self.Src0, 'b']))))
plt.plot(x,np.log10(np.abs(mu_0*np.imag(an))), 'r')
plt.plot(x,diff,'g')
plt.show()
simpegEM/Tests/test_FieldsObject.py
+40 -40
View File
@@ -8,86 +8,86 @@ class FieldsTest(unittest.TestCase):
mesh = Mesh.TensorMesh([np.ones(n)*5 for n in [10,11,12]],[0,0,-30])
x = np.linspace(5,10,3)
XYZ = Utils.ndgrid(x,x,np.r_[0.])
txLoc = np.r_[0,0,0.]
srcLoc = np.r_[0,0,0.]
rxList0 = EM.FDEM.RxFDEM(XYZ, 'exi')
Tx0 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList0])
Src0 = EM.FDEM.SrcFDEM(srcLoc, 'VMD', 3., [rxList0])
rxList1 = EM.FDEM.RxFDEM(XYZ, 'bxi')
Tx1 = EM.FDEM.TxFDEM(txLoc, 'VMD', 3., [rxList1])
Src1 = EM.FDEM.SrcFDEM(srcLoc, 'VMD', 3., [rxList1])
rxList2 = EM.FDEM.RxFDEM(XYZ, 'bxi')
Tx2 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList2])
Src2 = EM.FDEM.SrcFDEM(srcLoc, 'VMD', 2., [rxList2])
rxList3 = EM.FDEM.RxFDEM(XYZ, 'bxi')
Tx3 = EM.FDEM.TxFDEM(txLoc, 'VMD', 2., [rxList3])
Tx4 = EM.FDEM.TxFDEM(txLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3])
txList = [Tx0,Tx1,Tx2,Tx3,Tx4]
survey = EM.FDEM.SurveyFDEM(txList)
Src3 = EM.FDEM.SrcFDEM(srcLoc, 'VMD', 2., [rxList3])
Src4 = EM.FDEM.SrcFDEM(srcLoc, 'VMD', 1., [rxList0, rxList1, rxList2, rxList3])
srcList = [Src0,Src1,Src2,Src3,Src4]
survey = EM.FDEM.SurveyFDEM(srcList)
self.F = EM.FDEM.FieldsFDEM(mesh, survey)
self.Tx0 = Tx0
self.Tx1 = Tx1
self.Src0 = Src0
self.Src1 = Src1
self.mesh = mesh
self.XYZ = XYZ
def test_SetGet(self):
F = self.F
for freq in F.survey.freqs:
nFreq = F.survey.nTxByFreq[freq]
Txs = F.survey.getTransmitters(freq)
nFreq = F.survey.nSrcByFreq[freq]
Srcs = F.survey.getSources(freq)
e = np.random.rand(F.mesh.nE, nFreq)
F[Txs, 'e'] = e
F[Srcs, 'e'] = e
b = np.random.rand(F.mesh.nF, nFreq)
F[Txs, 'b'] = b
F[Srcs, 'b'] = b
if nFreq == 1:
F[Txs, 'b'] = Utils.mkvc(b)
F[Srcs, 'b'] = Utils.mkvc(b)
if e.shape[1] == 1:
e, b = Utils.mkvc(e), Utils.mkvc(b)
self.assertTrue(np.all(F[Txs, 'e'] == e))
self.assertTrue(np.all(F[Txs, 'b'] == b))
F[Txs] = {'b':b,'e':e}
self.assertTrue(np.all(F[Txs, 'e'] == e))
self.assertTrue(np.all(F[Txs, 'b'] == b))
self.assertTrue(np.all(F[Srcs, 'e'] == e))
self.assertTrue(np.all(F[Srcs, 'b'] == b))
F[Srcs] = {'b':b,'e':e}
self.assertTrue(np.all(F[Srcs, 'e'] == e))
self.assertTrue(np.all(F[Srcs, 'b'] == b))
lastFreq = F[Txs]
lastFreq = F[Srcs]
self.assertTrue(type(lastFreq) is dict)
self.assertTrue(sorted([k for k in lastFreq]) == ['b','e'])
self.assertTrue(np.all(lastFreq['b'] == b))
self.assertTrue(np.all(lastFreq['e'] == e))
Tx_f3 = F.survey.getTransmitters(3.)
self.assertTrue(F[Tx_f3,'b'].shape == (F.mesh.nF, 2))
Src_f3 = F.survey.getSources(3.)
self.assertTrue(F[Src_f3,'b'].shape == (F.mesh.nF, 2))
b = np.random.rand(F.mesh.nF, 2)
Tx_f0 = F.survey.getTransmitters(self.Tx0.freq)
F[Tx_f0,'b'] = b
self.assertTrue(F[self.Tx0]['b'].shape == (F.mesh.nF,))
self.assertTrue(F[self.Tx0,'b'].shape == (F.mesh.nF,))
self.assertTrue(np.all(F[self.Tx0,'b'] == b[:,0]))
self.assertTrue(np.all(F[self.Tx1,'b'] == b[:,1]))
Src_f0 = F.survey.getSources(self.Src0.freq)
F[Src_f0,'b'] = b
self.assertTrue(F[self.Src0]['b'].shape == (F.mesh.nF,))
self.assertTrue(F[self.Src0,'b'].shape == (F.mesh.nF,))
self.assertTrue(np.all(F[self.Src0,'b'] == b[:,0]))
self.assertTrue(np.all(F[self.Src1,'b'] == b[:,1]))
def test_assertions(self):
freq = self.F.survey.freqs[0]
Txs = self.F.survey.getTransmitters(freq)
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nTxByFreq[freq])
def fun(): self.F[Txs, 'b'] = bWrongSize
Srcs = self.F.survey.getSources(freq)
bWrongSize = np.random.rand(self.F.mesh.nE, self.F.survey.nSrcByFreq[freq])
def fun(): self.F[Srcs, 'b'] = bWrongSize
self.assertRaises(ValueError, fun)
def fun(): self.F[-999.]
self.assertRaises(KeyError, fun)
def fun(): self.F['notRight']
self.assertRaises(KeyError, fun)
def fun(): self.F[Txs,'notThere']
def fun(): self.F[Srcs,'notThere']
self.assertRaises(KeyError, fun)
def test_FieldProjections(self):
F = self.F
for freq in F.survey.freqs:
nFreq = F.survey.nTxByFreq[freq]
Txs = F.survey.getTransmitters(freq)
nFreq = F.survey.nSrcByFreq[freq]
Srcs = F.survey.getSources(freq)
e = np.random.rand(F.mesh.nE, nFreq)
b = np.random.rand(F.mesh.nF, nFreq)
F[Txs] = {'b':b,'e':e}
F[Srcs] = {'b':b,'e':e}
Txs = F.survey.getTransmitters(freq)
for ii, tx in enumerate(Txs):
for jj, rx in enumerate(tx.rxList):
dat = rx.projectFields(tx, self.mesh, F)
Srcs = F.survey.getSources(freq)
for ii, src in enumerate(Srcs):
for jj, rx in enumerate(src.rxList):
dat = rx.projectFields(src, self.mesh, F)
self.assertTrue(dat.dtype == float)
fieldType = rx.projField
u = {'b':b[:,ii], 'e': e[:,ii]}[fieldType]
simpegEM/Tests/test_TDEM_b_DerivAdjoint.py
+2 -2
View File
@@ -22,9 +22,9 @@ class TDEM_bDerivTests(unittest.TestCase):
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
survey = EM.TDEM.SurveyTDEM([tx])
survey = EM.TDEM.SurveyTDEM([src])
self.prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
# self.prb.timeSteps = [1e-5]
simpegEM/Tests/test_TDEM_combos.py
+7 -7
View File
@@ -4,7 +4,7 @@ import simpegEM as EM
plotIt = False
def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
cs = 5.
ncx = 20
ncy = 6
@@ -19,12 +19,12 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nTx=1):
rxOffset = 40.
txs = []
for ii in range(nTx):
srcs = []
for ii in range(nSrc):
rxs = [EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20 + ii), rxType) for rxType in rxTypes.split(',')]
txs += [EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)]
srcs += [EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', rxs)]
survey = EM.TDEM.SurveyTDEM(txs)
survey = EM.TDEM.SurveyTDEM(srcs)
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
# prb.timeSteps = [1e-5]
@@ -68,8 +68,8 @@ class TDEM_bDerivTests(unittest.TestCase):
def test_Jvec_bxbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz')))
def test_Adjoint_bxbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz')))
def test_Jvec_bxbz_2tx(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nTx=2)))
def test_Adjoint_bxbz_2tx(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nTx=2)))
def test_Jvec_bxbz_2src(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz',nSrc=2)))
def test_Adjoint_bxbz_2src(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz',nSrc=2)))
def test_Jvec_bxbzbz(self): self.assertTrue(dotestJvec(*getProb(rxTypes='bx,bz,bz')))
def test_Adjoint_bxbzbz(self): self.assertLess(*dotestAdjoint(*getProb(rxTypes='bx,bz,bz')))
simpegEM/Tests/test_TDEM_forward_Analytic.py
+2 -2
View File
@@ -28,9 +28,9 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
tx = EM.TDEM.TxTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
src = EM.TDEM.SrcTDEM(np.array([0., 0., 0.]), 'VMD_MVP', [rx])
survey = EM.TDEM.SurveyTDEM([tx])
survey = EM.TDEM.SurveyTDEM([src])
prb = EM.TDEM.ProblemTDEM_b(mesh, mapping=mapping)
prb.Solver = MumpsSolver