mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-18 12:40:30 +08:00
changed tx -> src
This commit is contained in:
Lindsey
+13
-13
@@ -1,13 +1,13 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class DipoleTx(Survey.BaseTx):
|
||||
"""A dipole transmitter, locA and locB are moved to the closest cell-centers"""
|
||||
class DipoleSrc(Survey.BaseSrc):
|
||||
"""A dipole source, locA and locB are moved to the closest cell-centers"""
|
||||
|
||||
current = 1
|
||||
|
||||
def __init__(self, locA, locB, rxList, **kwargs):
|
||||
super(DipoleTx, self).__init__((locA, locB), 'dipole', rxList, **kwargs)
|
||||
super(DipoleSrc, self).__init__((locA, locB), 'dipole', rxList, **kwargs)
|
||||
self._rhsDict = {}
|
||||
|
||||
def getRhs(self, mesh):
|
||||
@@ -21,7 +21,7 @@ class DipoleTx(Survey.BaseTx):
|
||||
|
||||
|
||||
class DipoleRx(Survey.BaseRx):
|
||||
"""A dipole transmitter, locA and locB are moved to the closest cell-centers"""
|
||||
"""A dipole source, locA and locB are moved to the closest cell-centers"""
|
||||
def __init__(self, locsM, locsN, **kwargs):
|
||||
locs = (locsM, locsN)
|
||||
assert locsM.shape == locsN.shape, 'locs must be the same shape.'
|
||||
@@ -46,8 +46,8 @@ class SurveyDC(Survey.BaseSurvey):
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
self.txList = txList
|
||||
def __init__(self, srcList, **kwargs):
|
||||
self.srcList = srcList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
self._rhsDict = {}
|
||||
self._Ps = {}
|
||||
@@ -64,7 +64,7 @@ class SurveyDC(Survey.BaseSurvey):
|
||||
|
||||
def getRhs(self, mesh):
|
||||
if mesh not in self._rhsDict:
|
||||
RHS = np.array([tx.getRhs(mesh) for tx in self.txList]).T
|
||||
RHS = np.array([src.getRhs(mesh) for src in self.srcList]).T
|
||||
self._rhsDict[mesh] = RHS
|
||||
return self._rhsDict[mesh]
|
||||
|
||||
@@ -72,9 +72,9 @@ class SurveyDC(Survey.BaseSurvey):
|
||||
if mesh in self._Ps:
|
||||
return self._Ps[mesh]
|
||||
|
||||
P_tx = [sp.vstack([rx.getP(mesh) for rx in tx.rxList]) for tx in self.txList]
|
||||
P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList]
|
||||
|
||||
self._Ps[mesh] = sp.block_diag(P_tx)
|
||||
self._Ps[mesh] = sp.block_diag(P_src)
|
||||
return self._Ps[mesh]
|
||||
|
||||
|
||||
@@ -176,7 +176,7 @@ class ProblemDC(Problem.BaseProblem):
|
||||
# Run forward simulation if $u$ not provided
|
||||
u = self.fields(self.curModel)
|
||||
else:
|
||||
shp = (self.mesh.nC, self.survey.nTx)
|
||||
shp = (self.mesh.nC, self.survey.nSrc)
|
||||
u = u.reshape(shp, order='F')
|
||||
|
||||
D = self.mesh.faceDiv
|
||||
@@ -186,8 +186,8 @@ class ProblemDC(Problem.BaseProblem):
|
||||
|
||||
# Take derivative of $C(m,u)$ w.r.t. $m$
|
||||
dCdm_x_v = np.empty_like(u)
|
||||
# loop over fields for each transmitter
|
||||
for i in range(self.survey.nTx):
|
||||
# loop over fields for each source
|
||||
for i in range(self.survey.nSrc):
|
||||
# Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
|
||||
dAdsig = D * self.dMdsig( G * u[:,i] )
|
||||
dCdm_x_v[:, i] = dAdsig * dsigdm_x_v
|
||||
@@ -207,7 +207,7 @@ class ProblemDC(Problem.BaseProblem):
|
||||
if u is None:
|
||||
u = self.fields(self.curModel)
|
||||
|
||||
shp = (self.mesh.nC, self.survey.nTx)
|
||||
shp = (self.mesh.nC, self.survey.nSrc)
|
||||
u = u.reshape(shp, order='F')
|
||||
P = self.survey.getP(self.mesh)
|
||||
PT_x_v = (P.T*v).reshape(shp, order='F')
|
||||
|
||||
@@ -23,40 +23,40 @@ def run(plotIt=False):
|
||||
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
|
||||
|
||||
rx = DC.DipoleRx(xyz_rxP, xyz_rxN)
|
||||
tx = DC.DipoleTx([-200, 0, -12.5],[+200, 0, -12.5], [rx])
|
||||
survey = DC.SurveyDC([tx])
|
||||
src = DC.DipoleSrc([-200, 0, -12.5],[+200, 0, -12.5], [rx])
|
||||
survey = DC.SurveyDC([src])
|
||||
problem = DC.ProblemDC(mesh)
|
||||
problem.pair(survey)
|
||||
data = survey.dpred(sigma)
|
||||
|
||||
def DChalf(txlocP, txlocN, rxloc, sigma, I=1.):
|
||||
rp = (txlocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
|
||||
rn = (txlocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
|
||||
def DChalf(srclocP, srclocN, rxloc, sigma, I=1.):
|
||||
rp = (srclocP.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
|
||||
rn = (srclocN.reshape([1,-1])).repeat(rxloc.shape[0], axis = 0)
|
||||
rP = np.sqrt(((rxloc-rp)**2).sum(axis=1))
|
||||
rN = np.sqrt(((rxloc-rn)**2).sum(axis=1))
|
||||
return I/(sigma*2.*np.pi)*(1/rP-1/rN)
|
||||
|
||||
data_analP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
|
||||
data_analN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
|
||||
data_anal = data_analP-data_analN
|
||||
Data_anal = data_anal.reshape((21, 21), order = 'F')
|
||||
data_anaP = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxP, sighalf)
|
||||
data_anaN = DChalf(np.r_[-200, 0, 0.],np.r_[+200, 0, 0.], xyz_rxN, sighalf)
|
||||
data_ana = data_anaP-data_anaN
|
||||
Data_ana = data_ana.reshape((21, 21), order = 'F')
|
||||
Data = data.reshape((21, 21), order = 'F')
|
||||
X = xyz_rxM[:,0].reshape((21, 21), order = 'F')
|
||||
Y = xyz_rxM[:,1].reshape((21, 21), order = 'F')
|
||||
|
||||
if plotIt:
|
||||
fig, ax = plt.subplots(1,2, figsize = (12, 5))
|
||||
vmin = np.r_[data, data_anal].min()
|
||||
vmax = np.r_[data, data_anal].max()
|
||||
vmin = np.r_[data, data_ana].min()
|
||||
vmax = np.r_[data, data_ana].max()
|
||||
dat1 = ax[1].contourf(X, Y, Data, 60, vmin = vmin, vmax = vmax)
|
||||
dat0 = ax[0].contourf(X, Y, Data_anal, 60, vmin = vmin, vmax = vmax)
|
||||
dat0 = ax[0].contourf(X, Y, Data_ana, 60, vmin = vmin, vmax = vmax)
|
||||
cb0 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[0])
|
||||
cb1 = plt.colorbar(dat1, orientation = 'horizontal', ax = ax[1])
|
||||
ax[1].set_title('Analytic')
|
||||
ax[0].set_title('Computed')
|
||||
plt.show()
|
||||
|
||||
return np.linalg.norm(data-data_anal)/np.linalg.norm(data_anal)
|
||||
return np.linalg.norm(data-data_ana)/np.linalg.norm(data_ana)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -3,7 +3,7 @@ import simpegDC as DC
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def getTxList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
def getSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
|
||||
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
|
||||
elocs -= (nElecs*aSpacing - aSpacing)/2
|
||||
@@ -24,19 +24,19 @@ def getTxList(nElecs, aSpacing, in2D=False, plotIt=False):
|
||||
plt.plot(elocs[WENNER[:,i]],s+'.')
|
||||
plt.show()
|
||||
|
||||
# Create transmitters and receivers
|
||||
# Create sources and receivers
|
||||
i = 0
|
||||
if in2D:
|
||||
getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
|
||||
else:
|
||||
getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0, 0]
|
||||
txList = []
|
||||
srcList = []
|
||||
for i in range(WENNER.shape[0]):
|
||||
rx = DC.DipoleRx(getLoc(i,1),getLoc(i,2))
|
||||
tx = DC.DipoleTx(getLoc(i,0),getLoc(i,3), [rx])
|
||||
txList += [tx]
|
||||
src = DC.DipoleSrc(getLoc(i,0),getLoc(i,3), [rx])
|
||||
srcList += [src]
|
||||
|
||||
return txList
|
||||
return srcList
|
||||
|
||||
|
||||
|
||||
@@ -53,8 +53,8 @@ def example(aSpacing=2.5, nElecs=10, plotIt=False):
|
||||
if plotIt:
|
||||
mesh.plotGrid(showIt=True)
|
||||
|
||||
txList = getTxList(nElecs, aSpacing, in2D=True)
|
||||
survey = DC.SurveyDC(txList)
|
||||
srcList = getSrcList(nElecs, aSpacing, in2D=True)
|
||||
survey = DC.SurveyDC(srcList)
|
||||
problem = DC.ProblemDC(mesh)
|
||||
problem.pair(survey)
|
||||
|
||||
|
||||
@@ -35,7 +35,7 @@ class DCProblemTests(unittest.TestCase):
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nTx)
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v, u=u))
|
||||
|
||||
Reference in New Issue
Block a user