wassname/simpeg
1
0
Fork 0
mirror of https://github.com/wassname/simpeg.git synced 2026-06-28 20:24:27 +08:00
Code Issues Packages Projects Releases Wiki Activity

merged in em/dev

Browse Source
This commit is contained in:
Lindsey Heagy
2016-02-22 11:26:30 -08:00
parent cd51ab8be7 52cc92f699
commit 341e902469
20 changed files with 136 additions and 96 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/EM/FDEM/FDEM.py
+2 -3
View File
@@ -92,8 +92,7 @@ class BaseFDEMProblem(BaseEMProblem):
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
df_dm_v = np.array(df_dm_v,dtype=complex)
Jv[src, rx] = rx.projectFieldsDeriv(src, self.mesh, u, df_dm_v)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
@@ -128,7 +127,7 @@ class BaseFDEMProblem(BaseEMProblem):
u_src = u[src, ftype]
for rx in src.rxList:
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
SimPEG/EM/FDEM/SurveyFDEM.py
+5 -5
View File
@@ -66,7 +66,7 @@ class Rx(SimPEG.Survey.BaseRx):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][2]
def projectFields(self, src, mesh, f):
def eval(self, src, mesh, f):
"""
Project fields to recievers to get data.
@@ -82,7 +82,7 @@ class Rx(SimPEG.Survey.BaseRx):
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False):
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
Derivative of projected fields with respect to the inversion model times a vector.
@@ -170,7 +170,7 @@ class Survey(SimPEG.Survey.BaseSurvey):
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
def eval(self, u):
"""
Project fields to receiver locations
:param Fields u: fields object
@@ -180,8 +180,8 @@ class Survey(SimPEG.Survey.BaseSurvey):
data = SimPEG.Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, u)
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def projectFieldsDeriv(self, u):
def evalDeriv(self, u):
raise Exception('Use Receivers to project fields deriv.')
SimPEG/EM/TDEM/SurveyTDEM.py
+6 -3
View File
@@ -48,15 +48,18 @@ class Rx(Survey.BaseTimeRx):
else:
return timeMesh.getInterpolationMat(self.times, self.projTLoc)
def projectFields(self, src, mesh, timeMesh, u):
def eval(self, src, mesh, timeMesh, u):
P = self.getP(mesh, timeMesh)
u_part = Utils.mkvc(u[src, self.projField, :])
return P*u_part
def projectFieldsDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
def evalDeriv(self, src, mesh, timeMesh, u, v, adjoint=False):
P = self.getP(mesh, timeMesh)
if not adjoint:
return P * Utils.mkvc(v[src, self.projField, :])
elif adjoint:
return P.T * v[src, self]
return P.T * v[src, self]
SimPEG/EM/TDEM_old/BaseTDEM.py
+2 -2
View File
@@ -128,7 +128,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
u = self.fields(m)
p = self.Gvec(m, v, u)
y = self.solveAh(m, p)
Jv = self.survey.projectFieldsDeriv(u, v=y)
Jv = self.survey.evalDeriv(u, v=y)
if self.verbose: print '%s\nDone calculating J(v)\n%s'%('*'*50,'*'*50)
return - mkvc(Jv)
@@ -155,7 +155,7 @@ class BaseTDEMProblem(BaseTimeProblem, BaseEMProblem):
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
p = self.survey.projectFieldsDeriv(u, v=v, adjoint=True)
p = self.survey.evalDeriv(u, v=v, adjoint=True)
y = self.solveAht(m, p)
w = self.Gtvec(m, y, u)
if self.verbose: print '%s\nDone calculating J^T(v)\n%s'%('*'*50,'*'*50)
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+2 -2
View File
@@ -21,8 +21,8 @@ def run(plotIt=True):
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
sig_half = 2e-2
sig_air = 1e-8
sig_layer = 1e-2
SimPEG/Examples/EM_TDEM_1D_Inversion.py
+2 -2
View File
@@ -19,8 +19,8 @@ def run(plotIt=True):
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-100.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
sig_half = 2e-3
sig_air = 1e-8
sig_layer = 1e-3
SimPEG/Examples/Inversion_Linear.py
+1 -1
View File
@@ -11,7 +11,7 @@ def run(N=100, plotIt=True):
"""
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
def eval(self, u):
return u
class LinearProblem(Problem.BaseProblem):
SimPEG/Examples/MT_3D_Foward.py
+1 -1
View File
@@ -52,7 +52,7 @@ def run(plotIt=True, nFreq=1):
# Calculate the data
fields = problem.fields(sig)
dataVec = survey.projectFields(fields)
dataVec = survey.eval(fields)
# Make the data
mtData = MT.Data(survey,dataVec)
SimPEG/FLOW/Richards/RichardsProblem.py
+10 -10
View File
@@ -8,7 +8,7 @@ class RichardsRx(Survey.BaseTimeRx):
knownRxTypes = ['saturation','pressureHead']
def projectFields(self, U, m, mapping, mesh, timeMesh):
def eval(self, U, m, mapping, mesh, timeMesh):
if self.rxType == 'pressureHead':
u = np.concatenate(U)
@@ -17,7 +17,7 @@ class RichardsRx(Survey.BaseTimeRx):
return self.getP(mesh, timeMesh) * u
def projectFieldsDeriv(self, U, m, mapping, mesh, timeMesh):
def evalDeriv(self, U, m, mapping, mesh, timeMesh):
P = self.getP(mesh, timeMesh)
if self.rxType == 'pressureHead':
@@ -57,13 +57,13 @@ class RichardsSurvey(Survey.BaseSurvey):
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.projectFields(u, m))
return Utils.mkvc(self.eval(u, m))
@Utils.requires('prob')
def projectFields(self, U, m):
def eval(self, U, m):
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFields(U, m,
Ds[ii] = rx.eval(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -71,11 +71,11 @@ class RichardsSurvey(Survey.BaseSurvey):
return np.concatenate(Ds)
@Utils.requires('prob')
def projectFieldsDeriv(self, U, m):
def evalDeriv(self, U, m):
"""The Derivative with respect to the fields."""
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFieldsDeriv(U, m,
Ds[ii] = rx.evalDeriv(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
B = np.array(sp.vstack(Bs).todense())
Ainv = self.Solver(A, **self.solverOpts)
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
AinvB = Ainv * B
z = np.zeros((self.mesh.nC, B.shape[1]))
zAinvB = np.vstack((z, AinvB))
@@ -277,7 +277,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
@Utils.timeIt
@@ -285,7 +285,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
if u is None:
u = self.field(m)
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
PTv = P.T*v
# This is done via backward substitution.
SimPEG/MT/BaseMT.py
+3 -3
View File
@@ -68,7 +68,7 @@ class BaseMTProblem(BaseFDEMProblem):
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
@@ -106,9 +106,9 @@ class BaseMTProblem(BaseFDEMProblem):
u_src = u[src, :]
for rx in src.rxList:
# Get the adjoint projectFieldsDeriv
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.projectFieldsDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
SimPEG/MT/SurveyMT.py
+5 -5
View File
@@ -59,7 +59,7 @@ class Rx(SimPEGsurvey.BaseRx):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][1]
def projectFields(self, src, mesh, f):
def eval(self, src, mesh, f):
'''
Project the fields to natural source data.
@@ -139,7 +139,7 @@ class Rx(SimPEGsurvey.BaseRx):
# print f_part
return f_part
def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False):
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
The derivative of the projection wrt u
@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
def eval(self, u):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, u)
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def projectFieldsDeriv(self, u):
def evalDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
#################
SimPEG/Maps.py
+41 -8
View File
@@ -4,6 +4,7 @@ from Tests import checkDerivative
from PropMaps import PropMap, Property
from numpy.polynomial import polynomial
from scipy.interpolate import UnivariateSpline
import warnings
class IdentityMap(object):
"""
@@ -296,11 +297,11 @@ class LogMap(IdentityMap):
def inverse(self, m):
return np.exp(Utils.mkvc(m))
class FullMap(IdentityMap):
class SurjectFull(IdentityMap):
"""
FullMap
SurjectFull
Given a scalar, the FullMap maps the value to the
Given a scalar, the SurjectFull maps the value to the
full model space.
"""
@@ -327,9 +328,15 @@ class FullMap(IdentityMap):
"""
return np.ones([self.mesh.nC,1])
class FullMap(SurjectFull):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`FullMap` is deprecated and will be removed in future versions. Use `SurjectFull` instead",
FutureWarning)
SurjectFull.__init__(self,mesh,**kwargs)
class Vertical1DMap(IdentityMap):
"""Vertical1DMap
class SurjectVertical1D(IdentityMap):
"""SurjectVertical1DMap
Given a 1D vector through the last dimension
of the mesh, this will extend to the full
@@ -369,8 +376,14 @@ class Vertical1DMap(IdentityMap):
), shape=(repNum, 1))
return sp.kron(sp.identity(self.nP), repVec)
class Vertical1DMap(SurjectVertical1D):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Vertical1DMap` is deprecated and will be removed in future versions. Use `SurjectVertical1D` instead",
FutureWarning)
SurjectVertical1D.__init__(self,mesh,**kwargs)
class Map2Dto3D(IdentityMap):
class Surject2Dto3D(IdentityMap):
"""Map2Dto3D
Given a 2D vector, this will extend to the full
@@ -425,6 +438,13 @@ class Map2Dto3D(IdentityMap):
), shape=(nC, nP))
return P
class Map2Dto3D(Surject2Dto3D):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Map2Dto3D` is deprecated and will be removed in future versions. Use `Surject2Dto3D` instead",
FutureWarning)
Surject2Dto3D.__init__(self,mesh,**kwargs)
class Mesh2Mesh(IdentityMap):
"""
Takes a model on one mesh are translates it to another mesh.
@@ -458,7 +478,7 @@ class Mesh2Mesh(IdentityMap):
return self.P
class ActiveCells(IdentityMap):
class InjectActiveCells(IdentityMap):
"""
Active model parameters.
@@ -506,7 +526,14 @@ class ActiveCells(IdentityMap):
def deriv(self, m):
return self.P
class ActiveCellsTopo(IdentityMap):
class ActiveCells(InjectActiveCells):
def __init__(self, mesh, indActive, valInactive, nC=None):
warnings.warn(
"`ActiveCells` is deprecated and will be removed in future versions. Use `InjectActiveCells` instead",
FutureWarning)
InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
class InjectActiveCellsTopo(IdentityMap):
"""
Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
@@ -577,6 +604,12 @@ class ActiveCellsTopo(IdentityMap):
def deriv(self, m):
return self.P
class ActiveCellsTopo(InjectActiveCellsTopo):
def __init__(self, mesh, indActive, valInactive, nC=None):
warnings.warn(
"`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead",
FutureWarning)
InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC)
class Weighting(IdentityMap):
"""
SimPEG/Survey.py
+7 -7
View File
@@ -307,12 +307,12 @@ class BaseSurvey(object):
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.projectFields(u))
return Utils.mkvc(self.eval(u))
@Utils.count
def projectFields(self, u):
"""projectFields(u)
def eval(self, u):
"""eval(u)
This function projects the fields onto the data space.
@@ -320,11 +320,11 @@ class BaseSurvey(object):
d_\\text{pred} = \mathbf{P} u(m)
"""
raise NotImplemented('projectFields is not yet implemented.')
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def projectFieldsDeriv(self, u):
"""projectFieldsDeriv(u)
def evalDeriv(self, u):
"""evalDeriv(u)
This function s the derivative of projects the fields onto the data space.
@@ -332,7 +332,7 @@ class BaseSurvey(object):
\\frac{\partial d_\\text{pred}}{\partial u} = \mathbf{P}
"""
raise NotImplemented('projectFields is not yet implemented.')
raise NotImplemented('eval is not yet implemented.')
@Utils.count
def residual(self, m, u=None):
tests/base/test_maps.py
+32 -27
View File
@@ -5,8 +5,8 @@ from scipy.sparse.linalg import dsolve
TOL = 1e-14
MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "Vertical1DMap", "Weighting", "FullMap"]
MAPS_TO_TEST_2D = ["CircleMap", "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
MAPS_TO_TEST_3D = [ "ComplexMap", "ExpMap", "IdentityMap", "SurjectVertical1D", "Weighting", "SurjectFull","FullMap","Vertical1DMap"]
class MapTests(unittest.TestCase):
@@ -52,7 +52,7 @@ class MapTests(unittest.TestCase):
def test_mapMultiplication(self):
M = Mesh.TensorMesh([2,3])
expMap = Maps.ExpMap(M)
vertMap = Maps.Vertical1DMap(M)
vertMap = Maps.SurjectVertical1D(M)
combo = expMap*vertMap
m = np.arange(3.0)
t_true = np.exp(np.r_[0,0,1,1,2,2.])
@@ -83,22 +83,23 @@ class MapTests(unittest.TestCase):
def test_activeCells(self):
M = Mesh.TensorMesh([2,4],'0C')
expMap = Maps.ExpMap(M)
actMap = Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
vertMap = Maps.Vertical1DMap(M)
combo = vertMap * actMap
m = np.r_[1,2.]
mod = Models.Model(m,combo)
# import matplotlib.pyplot as plt
# plt.colorbar(M.plotImage(mod.transform)[0])
# plt.show()
self.assertLess(np.linalg.norm(mod.transform - np.r_[1,1,2,2,10,10,10,10.]), TOL)
self.assertLess((mod.transformDeriv - combo.deriv(m)).toarray().sum(), TOL)
for actMap in [Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy), Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)]:
# actMap = Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
vertMap = Maps.SurjectVertical1D(M)
combo = vertMap * actMap
m = np.r_[1,2.]
mod = Models.Model(m,combo)
# import matplotlib.pyplot as plt
# plt.colorbar(M.plotImage(mod.transform)[0])
# plt.show()
self.assertLess(np.linalg.norm(mod.transform - np.r_[1,1,2,2,10,10,10,10.]), TOL)
self.assertLess((mod.transformDeriv - combo.deriv(m)).toarray().sum(), TOL)
def test_tripleMultiply(self):
M = Mesh.TensorMesh([2,4],'0C')
expMap = Maps.ExpMap(M)
vertMap = Maps.Vertical1DMap(M)
actMap = Maps.ActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
vertMap = Maps.SurjectVertical1D(M)
actMap = Maps.InjectActiveCells(M, M.vectorCCy <=0, 10, nC=M.nCy)
m = np.r_[1,2.]
t_true = np.exp(np.r_[1,1,2,2,10,10,10,10.])
self.assertLess(np.linalg.norm((expMap * vertMap * actMap * m)-t_true,np.inf),TOL)
@@ -115,29 +116,33 @@ class MapTests(unittest.TestCase):
M2 = Mesh.TensorMesh([2,4])
M3 = Mesh.TensorMesh([3,2,4])
m = np.random.rand(M2.nC)
m2to3 = Maps.Map2Dto3D(M3, normal='X')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[0,:,:] ) == m))
for m2to3 in [Maps.Surject2Dto3D(M3, normal='X'), Maps.Map2Dto3D(M3, normal='X')]:
# m2to3 = Maps.Surject2Dto3D(M3, normal='X')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[0,:,:] ) == m))
def test_map2Dto3D_y(self):
M2 = Mesh.TensorMesh([3,4])
M3 = Mesh.TensorMesh([3,2,4])
m = np.random.rand(M2.nC)
m2to3 = Maps.Map2Dto3D(M3, normal='Y')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,0,:] ) == m))
for m2to3 in [Maps.Surject2Dto3D(M3, normal='Y'),Maps.Map2Dto3D(M3, normal='Y')]:
# m2to3 = Maps.Surject2Dto3D(M3, normal='Y')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,0,:] ) == m))
def test_map2Dto3D_z(self):
M2 = Mesh.TensorMesh([3,2])
M3 = Mesh.TensorMesh([3,2,4])
m = np.random.rand(M2.nC)
m2to3 = Maps.Map2Dto3D(M3, normal='Z')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,:,0] ) == m))
for m2to3 in [Maps.Surject2Dto3D(M3, normal='Z'),Maps.Map2Dto3D(M3, normal='Z')]:
# m2to3 = Maps.Surject2Dto3D(M3, normal='Z')
m = np.arange(m2to3.nP)
self.assertTrue(m2to3.test())
self.assertTrue(np.all(Utils.mkvc( (m2to3 * m).reshape(M3.vnC,order='F')[:,:,0] ) == m))
if __name__ == '__main__':
tests/em/tdem/test_TDEM_b_DerivAdjoint.py
+4 -4
View File
@@ -18,8 +18,8 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = Mesh.CylMesh([hx,1,hy], '00C')
active = mesh.vectorCCz<0.
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
@@ -204,8 +204,8 @@ class TDEM_bDerivTests(unittest.TestCase):
d = Survey.Data(survey,v=d_vec)
# Check that d.T*Q*f = f.T*Q.T*d
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
V2 = f.tovec().dot(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec())
V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
V2 = f.tovec().dot(survey.evalDeriv(None, v=d, adjoint=True).tovec())
self.assertTrue((V1-V2)/np.abs(V1) < tol)
tests/em/tdem/test_TDEM_b_MultiSrc_DerivAdjoint.py
+4 -4
View File
@@ -17,8 +17,8 @@ class TDEM_bDerivTests(unittest.TestCase):
mesh = Mesh.CylMesh([hx,1,hy], '00C')
active = mesh.vectorCCz<0.
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
rxOffset = 40.
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-4,-3, 20), 'bz')
@@ -108,8 +108,8 @@ class TDEM_bDerivTests(unittest.TestCase):
d = Survey.Data(survey,v=d_vec)
# Check that d.T*Q*f = f.T*Q.T*d
V1 = d_vec.dot(survey.projectFieldsDeriv(None, v=f).tovec())
V2 = np.sum((f.tovec())*(survey.projectFieldsDeriv(None, v=d, adjoint=True).tovec()))
V1 = d_vec.dot(survey.evalDeriv(None, v=f).tovec())
V2 = np.sum((f.tovec())*(survey.evalDeriv(None, v=d, adjoint=True).tovec()))
self.assertTrue((V1-V2)/np.abs(V1) < 1e-6)
tests/em/tdem/test_TDEM_combos.py
+2 -2
View File
@@ -14,8 +14,8 @@ def getProb(meshType='CYL',rxTypes='bx,bz',nSrc=1):
mesh = Mesh.CylMesh([hx,1,hy], '00C')
active = mesh.vectorCCz<0.
activeMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * activeMap
activeMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * activeMap
rxOffset = 40.
tests/em/tdem/test_TDEM_forward_Analytic.py
+2 -2
View File
@@ -24,8 +24,8 @@ def halfSpaceProblemAnaDiff(meshType, sig_half=1e-2, rxOffset=50., bounds=[1e-5,
mesh = Mesh.TensorMesh([hx,hy,hz], 'CCC')
active = mesh.vectorCCz<0.
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
rx = EM.TDEM.RxTDEM(np.array([[rxOffset, 0., 0.]]), np.logspace(-5,-4, 21), 'bz')
src = EM.TDEM.SrcTDEM_VMD_MVP([rx], loc=np.array([0., 0., 0.]))
tests/mt/test_Problem1D_againstAnalyticHalfspace.py
+4 -4
View File
@@ -70,7 +70,7 @@ def appRes_TotalFieldNorm(sigmaHalf):
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
data = survey.eval(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
@@ -88,7 +88,7 @@ def appPhs_TotalFieldNorm(sigmaHalf):
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
data = survey.eval(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
@@ -106,7 +106,7 @@ def appRes_psFieldNorm(sigmaHalf):
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
data = survey.eval(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
@@ -124,7 +124,7 @@ def appPhs_psFieldNorm(sigmaHalf):
fields = problem.fields(sigma)
# Project the data
data = survey.projectFields(fields)
data = survey.eval(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
tests/mt/test_Problem3D_againstAnalytic.py
+1 -1
View File
@@ -210,7 +210,7 @@ def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
f = problem.fieldsPair(survey.mesh,survey)
f[src,'e_pxSolution'] = u[:len(u)/2]
f[src,'e_pySolution'] = u[len(u)/2::]
return rx.projectFields(src,survey.mesh,f), lambda t: rx.projectFieldsDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
return rx.eval(src,survey.mesh,f), lambda t: rx.evalDeriv(src,survey.mesh,f0,simpeg.mkvc(t,2))
return simpeg.Tests.checkDerivative(fun, u0, num=3, plotIt=False, eps=FLR)