Merge branch 'dev' into docs

This commit is contained in:
Lindsey Heagy
2016-03-06 21:50:27 -08:00
66 changed files with 6439 additions and 126 deletions
+4
View File
@@ -18,7 +18,9 @@ env:
- TEST_DIR="tests/mesh tests/base tests/utils"
- TEST_DIR=tests/em/fdem/inverse/derivs
- TEST_DIR=tests/em/tdem
- TEST_DIR=tests/dcip
- TEST_DIR=tests/flow
- TEST_DIR=tests/mt
- TEST_DIR=tests/examples
- TEST_DIR=tests/em/fdem/inverse/adjoint
- TEST_DIR=tests/em/fdem/forward
@@ -57,3 +59,5 @@ notifications:
email:
- rowanc1@gmail.com
- lindseyheagy@gmail.com
- gkrosen@gmail.com
- sgkang09@gmail.com
+294
View File
@@ -0,0 +1,294 @@
from SimPEG import *
class FieldsDC_CC(Problem.Fields):
knownFields = {'phi_sol':'CC'}
aliasFields = {
'phi' : ['phi_sol','CC','_phi'],
'e' : ['phi_sol','F','_e'],
'j' : ['phi_sol','F','_j']
}
def __init__(self,mesh,survey,**kwargs):
super(FieldsDC_CC, self).__init__(mesh, survey, **kwargs)
def startup(self):
self._cellGrad = self.survey.prob.mesh.cellGrad
self._Mfinv = self.survey.prob.mesh.getFaceInnerProduct(invMat=True)
def _phi(self, phi_sol, srcList):
phi = phi_sol
# for i, src in enumerate(srcList):
# phi_p = src.phi_p(self.survey.prob)
# if phi_p is not None:
# phi[:,i] += phi_p
return phi
def _e(self, phi_sol, srcList):
e = -self._cellGrad*phi_sol
# for i, src in enumerate(srcList):
# e_p = src.e_p(self.survey.prob)
# if e_p is not None:
# e[:,i] += e_p
return e
def _j(self, phi_sol, srcList):
j = -self._Mfinv*self.survey.prob.Msig*self._cellGrad*phi_sol
# for i, src in enumerate(srcList):
# j_p = src.j_p(self.survey.prob)
# if j_p is not None:
# j[:,i] += j_p
return j
class SrcDipole(Survey.BaseSrc):
"""A dipole source, locA and locB are moved to the closest cell-centers"""
current = 1
loc = None
# _rhsDict = None
def __init__(self, rxList, locA, locB, **kwargs):
self.loc = (locA, locB)
super(SrcDipole, self).__init__(rxList, **kwargs)
def eval(self, prob):
# Recompute rhs
# if getattr(self, '_rhsDict', None) is None:
# self._rhsDict = {}
# if mesh not in self._rhsDict:
pts = [self.loc[0], self.loc[1]]
inds = Utils.closestPoints(prob.mesh, pts)
q = np.zeros(prob.mesh.nC)
q[inds] = - self.current * ( np.r_[1., -1.] / prob.mesh.vol[inds] )
# self._rhsDict[mesh] = q
# return self._rhsDict[mesh]
return q
class RxDipole(Survey.BaseRx):
"""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.'
super(RxDipole, self).__init__(locs, 'dipole', storeProjections=False, **kwargs)
@property
def nD(self):
"""Number of data in the receiver."""
return self.locs[0].shape[0]
def getP(self, mesh):
P0 = mesh.getInterpolationMat(self.locs[0], self.projGLoc)
P1 = mesh.getInterpolationMat(self.locs[1], self.projGLoc)
return P0 - P1
class SurveyDC(Survey.BaseSurvey):
"""
**SurveyDC**
Geophysical DC resistivity data.
"""
uncert = None
def __init__(self, srcList, **kwargs):
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
# self._rhsDict = {}
self._Ps = {}
def eval(self, u):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
"""
P = self.getP(self.prob.mesh)
return P*mkvc(u[self.srcList, 'phi_sol'])
def getP(self, mesh):
if mesh in self._Ps:
return self._Ps[mesh]
P_src = [sp.vstack([rx.getP(mesh) for rx in src.rxList]) for src in self.srcList]
self._Ps[mesh] = sp.block_diag(P_src)
return self._Ps[mesh]
class ProblemDC_CC(Problem.BaseProblem):
"""
**ProblemDC**
Geophysical DC resistivity problem.
"""
surveyPair = SurveyDC
Solver = Solver
fieldsPair = FieldsDC_CC
Ainv = None
def __init__(self, mesh, **kwargs):
Problem.BaseProblem.__init__(self, mesh)
self.mesh.setCellGradBC('neumann')
Utils.setKwargs(self, **kwargs)
deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
@property
def Msig(self):
if getattr(self, '_Msig', None) is None:
sigma = self.curModel.transform
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
def dMdsig(self):
if getattr(self, '_dMdsig', None) is None:
sigma = self.curModel.transform
Av = self.mesh.aveF2CC
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
return self._dMdsig
@property
def A(self):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
Where M() is the mass matrix and mT is the model transform.
"""
if getattr(self, '_A', None) is None:
D = self.mesh.faceDiv
G = self.mesh.cellGrad
self._A = D*self.Msig*G
# Remove the null space from the matrix.
self._A[0,0] /= self.mesh.vol[0]
self._A = self._A.tocsc()
return self._A
def getRHS(self):
# if self.mesh not in self._rhsDict:
RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
# self._rhsDict[mesh] = RHS
# return self._rhsDict[mesh]
return RHS
def fields(self, m):
F = self.fieldsPair(self.mesh, self.survey)
self.curModel = m
A = self.A
self.Ainv = self.Solver(A, **self.solverOpts)
RHS = self.getRHS()
Phi = self.Ainv * RHS
Srcs = self.survey.srcList
F[Srcs, 'phi_sol'] = Phi
return F
def Jvec(self, m, v, u=None):
"""
:param numpy.array m: model
:param numpy.array v: vector to multiply
:param numpy.array u: fields
:rtype: numpy.array
:return: Jv
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
\\nabla_u (A(m)u - q) = A(m)
\\nabla_m (A(m)u - q) = G\\text{sdiag}(Du)\\nabla_m(M(mT(m)))
Where M() is the mass matrix and mT is the model transform.
.. math::
J = - P \left( \\nabla_u c(m, u) \\right)^{-1} \\nabla_m c(m, u)
J(v) = - P ( A(m)^{-1} ( G\\text{sdiag}(Du)\\nabla_m(M(mT(m))) v ) )
"""
# Set current model; clear dependent property $\mathbf{A(m)}$
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
D = self.mesh.faceDiv
G = self.mesh.cellGrad
# Derivative of model transform, $\deriv{\sigma}{\m}$
dsigdm_x_v = self.curModel.transformDeriv * v
# Take derivative of $C(m,u)$ w.r.t. $m$
dCdm_x_v = np.empty_like(u)
# 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
# Take derivative of $C(m,u)$ w.r.t. $u$
dA_du = self.A
# Solve for $\deriv{u}{m}$
# dCdu_inv = self.Solver(dCdu, **self.solverOpts)
if self.Ainv is None:
self.Ainv = self.Solver(dA_du, **self.solverOpts)
P = self.survey.getP(self.mesh)
Jv = - P * mkvc( self.Ainv * dCdm_x_v )
return Jv
def Jtvec(self, m, v, u=None):
self.curModel = m
sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
if u is None:
# Run forward simulation if $u$ not provided
u = self.fields(self.curModel)[self.survey.srcList, 'phi_sol']
else:
u = u[self.survey.srcList, 'phi_sol']
shp = u.shape
P = self.survey.getP(self.mesh)
PT_x_v = (P.T*v).reshape(shp, order='F')
D = self.mesh.faceDiv
G = self.mesh.cellGrad
dA_du = self.A
mT_dm = self.mapping.deriv(m)
# We probably always need this due to the linesearch .. (?)
self.Ainv = self.Solver(dA_du.T, **self.solverOpts)
# if self.Ainv is None:
# self.Ainv = self.Solver(dCdu, **self.solverOpts)
w = self.Ainv * PT_x_v
Jtv = 0
for i, ui in enumerate(u.T): # loop over each column
Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
Jtv = - mT_dm.T * ( Jtv )
return Jtv
+182
View File
@@ -0,0 +1,182 @@
from SimPEG import *
from BaseDC import SurveyDC, FieldsDC_CC
class SurveyIP(SurveyDC):
"""
**SurveyDC**
Geophysical DC resistivity data.
"""
def __init__(self, srcList, **kwargs):
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
self._Ps = {}
def dpred(self, m, u=None):
"""
Predicted data.
.. math::
d_\\text{pred} = Pu(m)
"""
return self.prob.forward(m)
class ProblemIP(Problem.BaseProblem):
"""
**ProblemIP**
Geophysical IP resistivity problem.
"""
surveyPair = SurveyDC
Solver = Solver
sigma = None
Ainv = None
u = None
def __init__(self, mesh, **kwargs):
Problem.BaseProblem.__init__(self, mesh)
self.mesh.setCellGradBC('neumann')
Utils.setKwargs(self, **kwargs)
# deleteTheseOnModelUpdate = ['_A', '_Msig', '_dMdsig']
@property
def Msig(self):
if getattr(self, '_Msig', None) is None:
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
self._Msig = Utils.sdiag(1/(self.mesh.dim * Av.T * (1/sigma)))
return self._Msig
@property
def dMdsig(self):
if getattr(self, '_dMdsig', None) is None:
# sigma = self.curModel.transform
sigma = self.sigma
Av = self.mesh.aveF2CC
dMdprop = self.mesh.dim * Utils.sdiag(self.Msig.diagonal()**2) * Av.T * Utils.sdiag(1./sigma**2)
self._dMdsig = lambda Gu: Utils.sdiag(Gu) * dMdprop
return self._dMdsig
@property
def A(self):
"""
Makes the matrix A(m) for the DC resistivity problem.
:param numpy.array m: model
:rtype: scipy.csc_matrix
:return: A(m)
.. math::
c(m,u) = A(m)u - q = G\\text{sdiag}(M(mT(m)))Du - q = 0
Where M() is the mass matrix and mT is the model transform.
"""
if getattr(self, '_A', None) is None:
D = self.mesh.faceDiv
G = self.mesh.cellGrad
self._A = D*self.Msig*G
# Remove the null space from the matrix.
self._A[-1,-1] /= self.mesh.vol[-1]
self._A = self._A.tocsc()
return self._A
def getRHS(self):
# if self.mesh not in self._rhsDict:
RHS = np.array([src.eval(self) for src in self.survey.srcList]).T
# self._rhsDict[mesh] = RHS
# return self._rhsDict[mesh]
return RHS
def fields(self, m):
if self.u is None:
A = self.A
if self.Ainv == None:
self.Ainv = self.Solver(A, **self.solverOpts)
Q = self.getRHS()
self.u = self.Ainv * Q
return self.u
def forward(self, m, u=None):
# Set current model; clear dependent property $\mathbf{A(m)}$
self.curModel = m
# sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
sigma = self.sigma
if self.u is None:
# Run forward simulation if $u$ not provided
u = self.fields(sigma)
shp = (self.mesh.nC, self.survey.nSrc)
u = self.u.reshape(shp, order='F')
D = self.mesh.faceDiv
G = self.mesh.cellGrad
# Derivative of model transform, $\deriv{\sigma}{\m}$
# dsigdm_x_v = self.curModel.transformDeriv * v
dsigdm_x_v = Utils.sdiag(sigma) * self.curModel.transformDeriv * m
# Take derivative of $C(m,u)$ w.r.t. $m$
dCdm_x_v = np.empty_like(u)
# 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
# Take derivative of $C(m,u)$ w.r.t. $u$
if self.Ainv == None:
self.Ainv = self.Solver(A, **self.solverOpts)
# dCdu = self.A
# Solve for $\deriv{u}{m}$
# dCdu_inv = self.Solver(dCdu, **self.solverOpts)
P = self.survey.getP(self.mesh)
J_x_v = - P * mkvc( self.Ainv * dCdm_x_v )
return -J_x_v
def Jvec(self, m, v, u=None):
return self.forward(v)
def Jtvec(self, m, v, u=None):
self.curModel = m
# sigma = self.curModel.transform # $\sigma = \mathcal{M}(\m)$
sigma = self.sigma
if self.u is None:
u = self.fields(sigma)
else:
u = self.u
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')
D = self.mesh.faceDiv
G = self.mesh.cellGrad
A = self.A
mT_dm = Utils.sdiag(sigma)*self.mapping.deriv(m)
# mT_dm = self.mapping.deriv(m)
# dCdu = A.T
# Ainv = self.Solver(dCdu, **self.solverOpts)
# if self.Ainv == None:
self.Ainv = self.Solver(A.T, **self.solverOpts)
w = self.Ainv * PT_x_v
Jtv = 0
for i, ui in enumerate(u.T): # loop over each column
Jtv += self.dMdsig( G * ui ).T * ( D.T * w[:,i] )
Jtv = - mT_dm.T * ( Jtv )
return -Jtv
+934
View File
@@ -0,0 +1,934 @@
from SimPEG import np
import BaseDC as DC
import BaseDC as IP
def getActiveindfromTopo(mesh, topo):
# def genActiveindfromTopo(mesh, topo):
"""
Get active indices from topography
"""
from scipy.interpolate import NearestNDInterpolator
if mesh.dim==3:
nCxy = mesh.nCx*mesh.nCy
Zcc = mesh.gridCC[:,2].reshape((nCxy, mesh.nCz), order='F')
Ftopo = NearestNDInterpolator(topo[:,:2], topo[:,2])
XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
XY.shape
topo = Ftopo(XY)
actind = []
for ixy in range(nCxy):
actind.append(topo[ixy] <= Zcc[ixy,:])
else:
raise NotImplementedError("Only 3D is working")
return Utils.mkvc(np.vstack(actind))
def gettopoCC(mesh, airind):
# def gettopoCC(mesh, airind):
"""
Get topography from active indices of mesh.
"""
mesh2D = Mesh.TensorMesh([mesh.hx, mesh.hy], mesh.x0[:2])
zc = mesh.gridCC[:,2]
AIRIND = airind.reshape((mesh.vnC[0]*mesh.vnC[1],mesh.vnC[2]), order='F')
ZC = zc.reshape((mesh.vnC[0]*mesh.vnC[1], mesh.vnC[2]), order='F')
topo = np.zeros(ZC.shape[0])
topoCC = np.zeros(ZC.shape[0])
for i in range(ZC.shape[0]):
ind = np.argmax(ZC[i,:][~AIRIND[i,:]])
topo[i] = ZC[i,:][~AIRIND[i,:]].max() + mesh.hz[~AIRIND[i,:]][ind]*0.5
topoCC[i] = ZC[i,:][~AIRIND[i,:]].max()
XY = Utils.ndgrid(mesh.vectorCCx, mesh.vectorCCy)
return mesh2D, topoCC
def readUBC_DC3Dobstopo(filename,mesh,topo,probType="CC"):
"""
Seogi's personal readObs function.
"""
text_file = open(filename, "r")
lines = text_file.readlines()
text_file.close()
SRC = []
DATA = []
srcLists = []
isrc = 0
# airind = getActiveindfromTopo(mesh, topo)
# mesh2D, topoCC = gettopoCC(mesh, airind)
for line in lines:
if "!" in line.split(): continue
elif line == '\n': continue
elif line == ' \n': continue
temp = map(float, line.split())
# Read a line for the current electrode
if len(temp) == 5: # SRC: Only X and Y are provided (assume no topography)
#TODO consider topography and assign the closest cell center in the earth
if isrc == 0:
DATA_temp = []
else:
DATA.append(np.asarray(DATA_temp))
DATA_temp = []
indM = Utils.closestPoints(mesh2D, DATA[isrc-1][:,1:3])
indN = Utils.closestPoints(mesh2D, DATA[isrc-1][:,3:5])
rx = DCIP.RxDipole(np.c_[DATA[isrc-1][:,1:3], topoCC[indM]], np.c_[DATA[isrc-1][:,3:5], topoCC[indN]])
temp = np.asarray(temp)
if [SRC[isrc-1][0], SRC[isrc-1][1]] == [SRC[isrc-1][2], SRC[isrc-1][3]]:
indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]])
tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[mesh.vectorCCx.max(), mesh.vectorCCy.max(), topoCC[-1]])
else:
indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]])
indB = Utils.closestPoints(mesh2D, [SRC[isrc-1][2], SRC[isrc-1][3]])
tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[SRC[isrc-1][2], SRC[isrc-1][3], topoCC[indB]])
srcLists.append(tx)
SRC.append(temp)
isrc += 1
elif len(temp) == 7: # SRC: X, Y and Z are provided
SRC.append(temp)
isrc += 1
elif len(temp) == 6: #
DATA_temp.append(np.r_[isrc, np.asarray(temp)])
elif len(temp) > 7:
DATA_temp.append(np.r_[isrc, np.asarray(temp)])
DATA.append(np.asarray(DATA_temp))
DATA_temp = []
indM = Utils.closestPoints(mesh2D, DATA[isrc-1][:,1:3])
indN = Utils.closestPoints(mesh2D, DATA[isrc-1][:,3:5])
rx = DCIP.RxDipole(np.c_[DATA[isrc-1][:,1:3], topoCC[indM]], np.c_[DATA[isrc-1][:,3:5], topoCC[indN]])
temp = np.asarray(temp)
if [SRC[isrc-1][0], SRC[isrc-1][1]] == [SRC[isrc-1][2], SRC[isrc-1][3]]:
indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]])
tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[mesh.vectorCCx.max(), mesh.vectorCCy.max(), topoCC[-1]])
else:
indA = Utils.closestPoints(mesh2D, [SRC[isrc-1][0], SRC[isrc-1][1]])
indB = Utils.closestPoints(mesh2D, [SRC[isrc-1][2], SRC[isrc-1][3]])
tx = DCIP.SrcDipole([rx], [SRC[isrc-1][0], SRC[isrc-1][1], topoCC[indA]],[SRC[isrc-1][2], SRC[isrc-1][3], topoCC[indB]])
srcLists.append(tx)
text_file.close()
survey = DCIP.SurveyDC(srcLists)
# Do we need this?
SRC = np.asarray(SRC)
DATA = np.vstack(DATA)
survey.dobs = np.vstack(DATA)[:,-2]
return {'DCsurvey':survey, 'airind':airind, 'topoCC':topoCC, 'SRC':SRC}
def readUBC_DC2DModel(fileName):
"""
Read UBC GIF 2DTensor model and generate 2D Tensor model in simpeg
Input:
:param fileName, path to the UBC GIF 2D model file
Output:
:param SimPEG TensorMesh 2D object
:return
Created on Thu Nov 12 13:14:10 2015
@author: dominiquef
"""
from SimPEG import np, mkvc
# Open fileand skip header... assume that we know the mesh already
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
dim = np.array(obsfile[0].split(),dtype=float)
temp = np.array(obsfile[1].split(),dtype=float)
if len(temp) > 1:
model = np.zeros(dim)
for ii in range(len(obsfile)-1):
mm = np.array(obsfile[ii+1].split(),dtype=float)
model[:,ii] = mm
model = model[:,::-1]
else:
if len(obsfile[1:])==1:
mm = np.array(obsfile[1:].split(),dtype=float)
else:
mm = np.array(obsfile[1:],dtype=float)
# Permute the second dimension to flip the order
model = mm.reshape(dim[1],dim[0])
model = model[::-1,:]
model = np.transpose(model, (1, 0))
model = mkvc(model)
return model
def plot_pseudoSection(DCsurvey, axs, stype):
"""
Read list of 2D tx-rx location and plot a speudo-section of apparent
resistivity.
Assumes flat topo for now...
Input:
:param d2D, z0
:switch stype -> Either 'pdp' (pole-dipole) | 'dpdp' (dipole-dipole)
Output:
:figure scatter plot overlayed on image
Edited Feb 17th, 2016
@author: dominiquef
"""
from SimPEG import np
from scipy.interpolate import griddata
import pylab as plt
# Set depth to 0 for now
z0 = 0.
# Pre-allocate
midx = []
midz = []
rho = []
count = 0 # Counter for data
for ii in range(DCsurvey.nSrc):
Tx = DCsurvey.srcList[ii].loc
Rx = DCsurvey.srcList[ii].rxList[0].locs
nD = DCsurvey.srcList[ii].rxList[0].nD
data = DCsurvey.dobs[count:count+nD]
count += nD
# Get distances between each poles A-B-M-N
MA = np.abs(Tx[0][0] - Rx[0][:,0])
MB = np.abs(Tx[1][0] - Rx[0][:,0])
NB = np.abs(Tx[1][0] - Rx[1][:,0])
NA = np.abs(Tx[0][0] - Rx[1][:,0])
MN = np.abs(Rx[1][:,0] - Rx[0][:,0])
# Create mid-point location
Cmid = (Tx[0][0] + Tx[1][0])/2
Pmid = (Rx[0][:,0] + Rx[1][:,0])/2
# Compute pant leg of apparent rho
if stype == 'pdp':
leg = data * 2*np.pi * MA * ( MA + MN ) / MN
leg = np.log10(abs(1/leg))
elif stype == 'dpdp':
leg = data * 2*np.pi / ( 1/MA - 1/MB - 1/NB + 1/NA )
midx = np.hstack([midx, ( Cmid + Pmid )/2 ])
midz = np.hstack([midz, -np.abs(Cmid-Pmid)/2 + z0 ])
rho = np.hstack([rho,leg])
ax = axs
# Grid points
grid_x, grid_z = np.mgrid[np.min(midx):np.max(midx), np.min(midz):np.max(midz)]
grid_rho = griddata(np.c_[midx,midz], rho.T, (grid_x, grid_z), method='linear')
plt.imshow(grid_rho.T, extent = (np.min(midx),np.max(midx),np.min(midz),np.max(midz)), origin='lower', alpha=0.8, vmin = np.min(rho), vmax = np.max(rho))
cbar = plt.colorbar(format = '%.2f',fraction=0.04,orientation="horizontal")
cmin,cmax = cbar.get_clim()
ticks = np.linspace(cmin,cmax,3)
cbar.set_ticks(ticks)
# Plot apparent resistivity
plt.scatter(midx,midz,s=50,c=rho.T)
ax.set_xticklabels([])
ax.set_ylabel('Z')
ax.yaxis.tick_right()
ax.yaxis.set_label_position('right')
plt.gca().set_aspect('equal', adjustable='box')
return ax
def gen_DCIPsurvey(endl, mesh, stype, a, b, n):
"""
Load in endpoints and survey specifications to generate Tx, Rx location
stations.
Assumes flat topo for now...
Input:
:param endl -> input endpoints [x1, y1, z1, x2, y2, z2]
:object mesh -> SimPEG mesh object
:switch stype -> "dpdp" (dipole-dipole) | "pdp" (pole-dipole) | 'gradient'
: param a, n -> pole seperation, number of rx dipoles per tx
Output:
:param Tx, Rx -> List objects for each tx location
Lines: P1x, P1y, P1z, P2x, P2y, P2z
Created on Wed December 9th, 2015
@author: dominiquef
!! Require clean up to deal with DCsurvey
"""
from SimPEG import np
def xy_2_r(x1,x2,y1,y2):
r = np.sqrt( np.sum((x2 - x1)**2 + (y2 - y1)**2) )
return r
## Evenly distribute electrodes and put on surface
# Mesure survey length and direction
dl_len = xy_2_r(endl[0,0],endl[1,0],endl[0,1],endl[1,1])
dl_x = ( endl[1,0] - endl[0,0] ) / dl_len
dl_y = ( endl[1,1] - endl[0,1] ) / dl_len
nstn = np.floor( dl_len / a )
# Compute discrete pole location along line
stn_x = endl[0,0] + np.array(range(int(nstn)))*dl_x*a
stn_y = endl[0,1] + np.array(range(int(nstn)))*dl_y*a
# Create line of P1 locations
M = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
# Create line of P2 locations
N = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
## Build list of Tx-Rx locations depending on survey type
# Dipole-dipole: Moving tx with [a] spacing -> [AB a MN1 a MN2 ... a MNn]
# Pole-dipole: Moving pole on one end -> [A a MN1 a MN2 ... MNn a B]
Tx = []
Rx = []
SrcList = []
if stype != 'gradient':
for ii in range(0, int(nstn)-1):
if stype == 'dpdp':
tx = np.c_[M[ii,:],N[ii,:]]
elif stype == 'pdp':
tx = np.c_[M[ii,:],M[ii,:]]
# Rx.append(np.c_[M[ii+1:indx,:],N[ii+1:indx,:]])
# Current elctrode seperation
AB = xy_2_r(tx[0,1],endl[1,0],tx[1,1],endl[1,1])
# Number of receivers to fit
nstn = np.min([np.floor( (AB - b) / a ) , n])
# Check if there is enough space, else break the loop
if nstn <= 0:
continue
# Compute discrete pole location along line
stn_x = N[ii,0] + dl_x*b + np.array(range(int(nstn)))*dl_x*a
stn_y = N[ii,1] + dl_y*b + np.array(range(int(nstn)))*dl_y*a
# Create receiver poles
# Create line of P1 locations
P1 = np.c_[stn_x, stn_y, np.ones(nstn).T*mesh.vectorNz[-1]]
# Create line of P2 locations
P2 = np.c_[stn_x+a*dl_x, stn_y+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
Rx.append(np.c_[P1,P2])
rxClass = DC.RxDipole(P1, P2)
Tx.append(tx)
if stype == 'dpdp':
srcClass = DC.SrcDipole([rxClass], M[ii,:],N[ii,:])
elif stype == 'pdp':
srcClass = DC.SrcDipole([rxClass], M[ii,:],M[ii,:])
SrcList.append(srcClass)
#==============================================================================
# elif re.match(stype,'dpdp'):
#
# for ii in range(0, int(nstn)-2):
#
# indx = np.min([ii+n+1,nstn])
# Tx.append(np.c_[M[ii,:],N[ii,:]])
# Rx.append(np.c_[M[ii+2:indx,:],N[ii+2:indx,:]])
#==============================================================================
elif stype == 'gradient':
# Gradient survey only requires Tx at end of line and creates a square
# grid of receivers at in the middle at a pre-set minimum distance
Tx.append(np.c_[M[0,:],N[-1,:]])
# Get the edge limit of survey area
min_x = endl[0,0] + dl_x * b
min_y = endl[0,1] + dl_y * b
max_x = endl[1,0] - dl_x * b
max_y = endl[1,1] - dl_y * b
box_l = np.sqrt( (min_x - max_x)**2 + (min_y - max_y)**2 )
box_w = box_l/2.
nstn = np.floor( box_l / a )
# Compute discrete pole location along line
stn_x = min_x + np.array(range(int(nstn)))*dl_x*a
stn_y = min_y + np.array(range(int(nstn)))*dl_y*a
# Define number of cross lines
nlin = int(np.floor( box_w / a ))
lind = range(-nlin,nlin+1)
ngrad = nstn * len(lind)
rx = np.zeros([ngrad,6])
for ii in range( len(lind) ):
# Move line in perpendicular direction by dipole spacing
lxx = stn_x - lind[ii]*a*dl_y
lyy = stn_y + lind[ii]*a*dl_x
M = np.c_[ lxx, lyy , np.ones(nstn).T*mesh.vectorNz[-1]]
N = np.c_[ lxx+a*dl_x, lyy+a*dl_y, np.ones(nstn).T*mesh.vectorNz[-1]]
rx[(ii*nstn):((ii+1)*nstn),:] = np.c_[M,N]
Rx.append(rx)
rxClass = DC.RxDipole(rx[:,:3], rx[:,3:])
srcClass = DC.SrcDipole([rxClass], M[0,:], N[-1,:])
SrcList.append(srcClass)
else:
print """stype must be either 'pdp', 'dpdp' or 'gradient'. """
survey = DC.SurveyDC(SrcList)
return survey, Tx, Rx
def writeUBC_DCobs(fileName, DCsurvey, dtype, stype):
"""
Write UBC GIF DCIP 2D or 3D observation file
Input:
:string fileName -> including path where the file is written out
:DCsurvey -> DC survey class object
:string dtype -> either '2D' | '3D'
:string stype -> either 'SURFACE' | 'GENERAL'
Output:
:param UBC2D-Data file
:return
Last edit: February 16th, 2016
@author: dominiquef
"""
from SimPEG import mkvc
assert (dtype=='2D') | (dtype=='3D'), "Data must be either '2D' | '3D'"
assert (stype=='SURFACE') | (stype=='GENERAL') | (stype=='SIMPLE'), "Data must be either 'SURFACE' | 'GENERAL' | 'SIMPLE'"
fid = open(fileName,'w')
fid.write('! ' + stype + ' FORMAT\n')
count = 0
for ii in range(DCsurvey.nSrc):
tx = np.c_[DCsurvey.srcList[ii].loc]
rx = DCsurvey.srcList[ii].rxList[0].locs
nD = DCsurvey.srcList[ii].nD
M = rx[0]
N = rx[1]
# Adapt source-receiver location for dtype and stype
if dtype=='2D':
if stype == 'SIMPLE':
#fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
A = np.repeat(tx[0,0],M.shape[0],axis=0)
B = np.repeat(tx[0,1],M.shape[0],axis=0)
M = M[:,0]
N = N[:,0]
np.savetxt(fid, np.c_[A, B, M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
else:
if stype == 'SURFACE':
fid.writelines("%e " % ii for ii in mkvc(tx[0,:]))
M = M[:,0]
N = N[:,0]
if stype == 'GENERAL':
fid.writelines("%e " % ii for ii in mkvc(tx[::2,:]))
M = M[:,0::2]
N = N[:,0::2]
fid.write('%i\n'% nD)
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
if dtype=='3D':
if stype == 'SURFACE':
fid.writelines("%e " % ii for ii in mkvc(tx[0:2,:]))
M = M[:,0:2]
N = N[:,0:2]
if stype == 'GENERAL':
fid.writelines("%e " % ii for ii in mkvc(tx))
fid.write('%i\n'% nD)
np.savetxt(fid, np.c_[ M, N , DCsurvey.dobs[count:count+nD], DCsurvey.std[count:count+nD] ], fmt='%e',delimiter=' ',newline='\n')
count += nD
fid.close()
def convertObs_DC3D_to_2D(DCsurvey,lineID):
"""
Read DC survey and data and change
coordinate system to distance along line assuming
all data is acquired along line.
First transmitter pole is assumed to be at the origin
Assumes flat topo for now...
Input:
:param Tx, Rx
Output:
:figure Tx2d, Rx2d
Edited Feb 17th, 2016
@author: dominiquef
"""
from SimPEG import np
def stn_id(v0,v1,r):
"""
Compute station ID along line
"""
dl = int(v0.dot(v1)) * r
return dl
srcLists = []
srcMat = getSrc_locs(DCsurvey)
# Find all unique line id
uniqueID = np.unique(lineID)
for jj in range(len(uniqueID)):
indx = np.where(lineID==uniqueID[jj])[0]
# Find origin of survey
r = 1e+8 # Initialize to some large number
Tx = srcMat[indx]
x0 = Tx[0][0,0:2] # Define station zero along line
vecTx, r1 = r_unit(x0,Tx[-1][1,0:2])
for ii in range(len(indx)):
# Get all receivers
Rx = DCsurvey.srcList[indx[ii]].rxList[0].locs
nrx = Rx[0].shape[0]
# Find A electrode along line
vec, r = r_unit(x0,Tx[ii][0,0:2])
A = stn_id(vecTx,vec,r)
# Find B electrode along line
vec, r = r_unit(x0,Tx[ii][1,0:2])
B = stn_id(vecTx,vec,r)
M = np.zeros(nrx)
N = np.zeros(nrx)
for kk in range(nrx):
# Find all M electrodes along line
vec, r = r_unit(x0,Rx[0][kk,0:2])
M[kk] = stn_id(vecTx,vec,r)
# Find all N electrodes along line
vec, r = r_unit(x0,Rx[1][kk,0:2])
N[kk] = stn_id(vecTx,vec,r)
Rx = DC.RxDipole(np.c_[M,np.zeros(nrx),Rx[0][:,2]],np.c_[N,np.zeros(nrx),Rx[1][:,2]])
srcLists.append( DC.SrcDipole( [Rx], np.asarray([A,0,Tx[ii][0,2]]),np.asarray([B,0,Tx[ii][1,2]]) ) )
DCsurvey2D = DC.SurveyDC(srcLists)
DCsurvey2D.dobs = np.asarray(DCsurvey.dobs)
DCsurvey2D.std = np.asarray(DCsurvey.std)
return DCsurvey2D
def readUBC_DC3Dobs(fileName):
"""
Read UBC GIF DCIP 3D observation file and generate arrays for tx-rx location
Input:
:param fileName, path to the UBC GIF 3D obs file
Output:
:param rx, tx, d, wd
:return
Created on Mon December 7th, 2015
@author: dominiquef
"""
# Load file
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
# Pre-allocate
srcLists = []
Rx = []
d = []
wd = []
zflag = True # Flag for z value provided
# Countdown for number of obs/tx
count = 0
for ii in range(obsfile.shape[0]):
if not obsfile[ii]:
continue
# First line is transmitter with number of receivers
if count==0:
temp = (np.fromstring(obsfile[ii], dtype=float,sep=' ').T)
count = int(temp[-1])
# Check if z value is provided, if False -> nan
if len(temp)==5:
tx = np.r_[temp[0:2],np.nan,temp[0:2],np.nan]
zflag = False
else:
tx = temp[:-1]
rx = []
continue
temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
if zflag:
rx.append(temp[:-2])
# Check if there is data with the location
if len(temp)==8:
d.append(temp[-2])
wd.append(temp[-1])
else:
rx.append(np.r_[temp[0:2],np.nan,temp[0:2],np.nan] )
# Check if there is data with the location
if len(temp)==6:
d.append(temp[-2])
wd.append(temp[-1])
count = count -1
# Reach the end of transmitter block
if count == 0:
rx = np.asarray(rx)
Rx = DC.RxDipole(rx[:,:3],rx[:,3:])
srcLists.append( DC.SrcDipole( [Rx], tx[:3],tx[3:]) )
# Create survey class
survey = DC.SurveyDC(srcLists)
survey.dobs = np.asarray(d)
survey.std = np.asarray(wd)
return {'DCsurvey':survey}
def readUBC_DC2Dobs(fileName):
"""
Read UBC GIF 2D observation file and generate arrays for tx-rx location
Input:
:param fileName, path to the UBC GIF 2D model file
Output:
:param rx, tx
:return
Created on Thu Nov 12 13:14:10 2015
@author: dominiquef
"""
from SimPEG import np
# Load file
obsfile = np.genfromtxt(fileName,delimiter=' \n',dtype=np.str,comments='!')
# Check first line and figure out if 2D or 3D file format
line = np.array(obsfile[0].split(),dtype=float)
tx_A = []
tx_B = []
rx_M = []
rx_N = []
d = []
wd = []
for ii in range(obsfile.shape[0]):
# If len==3, then simple format where tx-rx is listed on each line
if len(line) == 4:
temp = np.fromstring(obsfile[ii], dtype=float,sep=' ')
tx_A = np.hstack((tx_A,temp[0]))
tx_B = np.hstack((tx_B,temp[1]))
rx_M = np.hstack((rx_M,temp[2]))
rx_N = np.hstack((rx_N,temp[3]))
rx = np.transpose(np.array((rx_M,rx_N)))
tx = np.transpose(np.array((tx_A,tx_B)))
return tx, rx, d, wd
def readUBC_DC2DMesh(fileName):
"""
Read UBC GIF 2DTensor mesh and generate 2D Tensor mesh in simpeg
Input:
:param fileName, path to the UBC GIF mesh file
Output:
:param SimPEG TensorMesh 2D object
:return
Created on Thu Nov 12 13:14:10 2015
@author: dominiquef
"""
from SimPEG import np
# Open file
fopen = open(fileName,'r')
# Read down the file and unpack dx vector
def unpackdx(fid,nrows):
for ii in range(nrows):
line = fid.readline()
var = np.array(line.split(),dtype=float)
if ii==0:
x0= var[0]
xvec = np.ones(int(var[2])) * (var[1] - var[0]) / int(var[2])
xend = var[1]
else:
xvec = np.hstack((xvec,np.ones(int(var[1])) * (var[0] - xend) / int(var[1])))
xend = var[0]
return x0, xvec
#%% Start with dx block
# First line specifies the number of rows for x-cells
line = fopen.readline()
nl = np.array(line.split(),dtype=float)
[x0, dx] = unpackdx(fopen,nl)
#%% Move down the file until reaching the z-block
line = fopen.readline()
if not line:
line = fopen.readline()
#%% End with dz block
# First line specifies the number of rows for z-cells
line = fopen.readline()
nl = np.array(line.split(),dtype=float)
[z0, dz] = unpackdx(fopen,nl)
# Flip z0 to be the bottom of the mesh for SimPEG
z0 = z0 - sum(dz)
dz = dz[::-1]
#%% Make the mesh using SimPEG
from SimPEG import Mesh
tensMsh = Mesh.TensorMesh([dx,dz],(x0, z0))
return tensMsh
def xy_2_lineID(DCsurvey):
"""
Read DC survey class and append line ID.
Assumes that the locations are listed in the order
they were collected. May need to generalize for random
point locations, but will be more expensive
Input:
:param DCdict Vectors of station location
Output:
:param LineID Vector of integers
:return
Created on Thu Feb 11, 2015
@author: dominiquef
"""
# Compute unit vector between two points
nstn = DCsurvey.nSrc
# Pre-allocate space
lineID = np.zeros(nstn)
linenum = 0
indx = 0
for ii in range(nstn):
if ii == 0:
A = DCsurvey.srcList[ii].loc[0]
B = DCsurvey.srcList[ii].loc[1]
xout = np.mean([A[0:2],B[0:2]], axis = 0)
xy0 = A[:2]
xym = xout
# Deal with replicate pole location
if np.all(xy0==xym):
xym[0] = xym[0] + 1e-3
continue
A = DCsurvey.srcList[ii].loc[0]
B = DCsurvey.srcList[ii].loc[1]
xin = np.mean([A[0:2],B[0:2]], axis = 0)
# Compute vector between neighbours
vec1, r1 = r_unit(xout,xin)
# Compute vector between current stn and mid-point
vec2, r2 = r_unit(xym,xin)
# Compute vector between current stn and start line
vec3, r3 = r_unit(xy0,xin)
# Compute vector between mid-point and start line
vec4, r4 = r_unit(xym,xy0)
# Compute dot product
ang1 = np.abs(vec1.dot(vec2))
ang2 = np.abs(vec3.dot(vec4))
# If the angles are smaller then 45d, than next point is on a new line
if ((ang1 < np.cos(np.pi/4.)) | (ang2 < np.cos(np.pi/4.))) & (np.all(np.r_[r1,r2,r3,r4] > 0)):
# Re-initiate start and mid-point location
xy0 = A[:2]
xym = xin
# Deal with replicate pole location
if np.all(xy0==xym):
xym[0] = xym[0] + 1e-3
linenum += 1
indx = ii
else:
xym = np.mean([xy0,xin], axis = 0)
lineID[ii] = linenum
xout = xin
return lineID
def r_unit(p1,p2):
"""
r_unit(x,y) : Function computes the unit vector
between two points with coordinates p1(x1,y1) and p2(x2,y2)
"""
assert len(p1)==len(p2), 'locs must be the same shape.'
dx = []
for ii in range(len(p1)):
dx.append((p2[ii] - p1[ii]))
# Compute length of vector
r = np.linalg.norm(np.asarray(dx))
if r!=0:
vec = dx/r
else:
vec = np.zeros(len(p1))
return vec, r
def getSrc_locs(DCsurvey):
"""
"""
srcMat = np.zeros((DCsurvey.nSrc,2,3))
for ii in range(DCsurvey.nSrc):
print np.asarray(DCsurvey.srcList[ii].loc).shape
srcMat[ii,:,:] = np.asarray(DCsurvey.srcList[ii].loc)
return srcMat
+38
View File
@@ -0,0 +1,38 @@
import numpy as np
def WennerSrcList(nElecs, aSpacing, in2D=False, plotIt=False):
import SimPEG.DCIP as DC
elocs = np.arange(0,aSpacing*nElecs,aSpacing)
elocs -= (nElecs*aSpacing - aSpacing)/2
space = 1
WENNER = np.zeros((0,),dtype=int)
for ii in range(nElecs):
for jj in range(nElecs):
test = np.r_[jj,jj+space,jj+space*2,jj+space*3]
if np.any(test >= nElecs):
break
WENNER = np.r_[WENNER, test]
space += 1
WENNER = WENNER.reshape((-1,4))
if plotIt:
for i, s in enumerate('rbkg'):
plt.plot(elocs[WENNER[:,i]],s+'.')
plt.show()
# 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]
srcList = []
for i in range(WENNER.shape[0]):
rx = DC.RxDipole(getLoc(i,1),getLoc(i,2))
src = DC.SrcDipole([rx], getLoc(i,0),getLoc(i,3))
srcList += [src]
return srcList
+4
View File
@@ -0,0 +1,4 @@
from BaseDC import *
from BaseIP import *
from DCIPUtils import *
import Utils
+33
View File
@@ -237,6 +237,39 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
# Save the file as a npz
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
class SaveOutputDictEveryIteration(_SaveEveryIteration):
"""SaveOutputDictEveryIteration
A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
"""
def initialize(self):
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
def endIter(self):
# Save the data.
ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
phi_ms = 0.5*ms.dot(ms)
if self.reg.smoothModel == True:
mref = self.reg.mref
else:
mref = 0
mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mx = 0.5 * mx.dot(mx)
if self.prob.mesh.dim==2:
my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_my = 0.5 * my.dot(my)
else:
phi_my = 'NaN'
if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
phi_mz = 0.5 * mz.dot(mz)
else:
phi_mz = 'NaN'
# Save the file as a npz
np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
# class UpdateReferenceModel(Parameter):
+2 -2
View File
@@ -99,7 +99,7 @@ class BaseFDEMProblem(BaseEMProblem):
Df_Dm = np.array(df_dudu_dm + df_dm,dtype=complex)
P = lambda v: rx.projectFieldsDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
P = lambda v: rx.evalDeriv(src, self.mesh, u, v) # wrt u, also have wrt m
Jv[src, rx] = P(Df_Dm)
@@ -137,7 +137,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_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
+11 -12
View File
@@ -66,30 +66,29 @@ class Rx(SimPEG.Survey.BaseRx):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][2]
def projectFields(self, src, mesh, u):
def eval(self, src, mesh, f):
"""
Project fields to recievers to get data.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields u: fields object
:param Fields f: fields object
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
u_part_complex = u[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
P = self.getP(mesh) # get interpolation to recievers
u_part_complex = f[src, self.projField]
real_or_imag = self.projComp # get the real or imag component
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
def projectFieldsDeriv(self, src, mesh, u, 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.
:param Source src: FDEM source
:param Mesh mesh: mesh used
:param Fields u: fields object
:param Fields f: fields object
:param numpy.ndarray v: vector to multiply
:rtype: numpy.ndarray
:return: fields projected to recievers
@@ -171,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
@@ -181,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):
raise Exception('Use Sources to project fields deriv.')
def evalDeriv(self, u):
raise Exception('Use Receivers to project fields deriv.')
+3 -2
View File
@@ -27,6 +27,7 @@ class FieldsTDEM(Problem.TimeFields):
else:
e = np.zeros((nE,nSrc)) # if nSrc == 1 else (nE, nSrc))
u = np.concatenate((u, b, e))
return Utils.mkvc(u,nSrc)
@@ -128,7 +129,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 +156,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)
+31 -12
View File
@@ -51,12 +51,12 @@ class RxTDEM(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:
@@ -79,12 +79,32 @@ class SrcTDEM(Survey.BaseSrc):
class SrcTDEM_VMD_MVP(SrcTDEM):
def __init__(self,rxList,loc):
def __init__(self,rxList,loc,waveformType="STEPOFF"):
self.loc = loc
self.waveformType = waveformType
SrcTDEM.__init__(self,rxList)
def getInitialFields(self, mesh):
"""Vertical magnetic dipole, magnetic vector potential"""
if self.waveformType == "STEPOFF":
print ">> Step waveform: Non-zero initial condition"
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
else:
raise NotImplementedError('Non-symmetric cyl mesh not implemented yet!')
elif mesh._meshType is 'TENSOR':
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
else:
raise Exception('Unknown mesh for VMD')
return {"b": mesh.edgeCurl*MVP}
elif self.waveformType == "GENERAL":
print ">> General waveform: Zero initial condition"
return {"b": np.zeros(mesh.nF)}
else:
raise NotImplementedError("Only use STEPOFF or GENERAL")
def getMeS(self, mesh, MfMui):
if mesh._meshType is 'CYL':
if mesh.isSymmetric:
MVP = MagneticDipoleVectorPotential(self.loc, mesh, 'Ey')
@@ -93,13 +113,12 @@ class SrcTDEM_VMD_MVP(SrcTDEM):
elif mesh._meshType is 'TENSOR':
MVP = MagneticDipoleVectorPotential(self.loc, mesh, ['Ex','Ey','Ez'])
else:
raise Exception('Unknown mesh for VMD')
return {"b": mesh.edgeCurl*MVP}
raise Exception('Unknown mesh for VMD')
return mesh.edgeCurl.T*MfMui*mesh.edgeCurl*MVP
class SrcTDEM_CircularLoop_MVP(SrcTDEM):
def __init__(self,rxList,loc,radius,waveformType):
def __init__(self,rxList,loc,radius,waveformType="STEPOFF"):
self.loc = loc
self.radius = radius
self.waveformType = waveformType
@@ -149,27 +168,27 @@ class SurveyTDEM(Survey.BaseSurvey):
self.srcList = srcList
Survey.BaseSurvey.__init__(self, **kwargs)
def projectFields(self, u):
def eval(self, u):
data = Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, self.prob.timeMesh, u)
data[src, rx] = rx.eval(src, self.mesh, self.prob.timeMesh, u)
return data
def projectFieldsDeriv(self, u, v=None, adjoint=False):
def evalDeriv(self, u, v=None, adjoint=False):
assert v is not None, 'v to multiply must be provided.'
if not adjoint:
data = Survey.Data(self)
for src in self.srcList:
for rx in src.rxList:
data[src, rx] = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v)
data[src, rx] = rx.evalDeriv(src, self.mesh, self.prob.timeMesh, u, v)
return data
else:
f = FieldsTDEM(self.mesh, self)
for src in self.srcList:
for rx in src.rxList:
Ptv = rx.projectFieldsDeriv(src, self.mesh, self.prob.timeMesh, u, v, adjoint=True)
Ptv = rx.evalDeriv(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[src, rx.projField, :] = Ptv
+68
View File
@@ -0,0 +1,68 @@
from SimPEG import *
import SimPEG.DCIP as DC
def run(plotIt=False):
cs = 25.
hx = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,21),(cs,7, 1.3)]
hz = [(cs,7, -1.3),(cs,20)]
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
sighalf = 1e-2
sigma = np.ones(mesh.nC)*sighalf
xtemp = np.linspace(-150, 150, 21)
ytemp = np.linspace(-150, 150, 21)
xyz_rxP = Utils.ndgrid(xtemp-10., ytemp, np.r_[0.])
xyz_rxN = Utils.ndgrid(xtemp+10., ytemp, np.r_[0.])
xyz_rxM = Utils.ndgrid(xtemp, ytemp, np.r_[0.])
# if plotIt:
# fig, ax = plt.subplots(1,1, figsize = (5,5))
# mesh.plotSlice(sigma, grid=True, ax = ax)
# ax.plot(xyz_rxP[:,0],xyz_rxP[:,1], 'w.')
# ax.plot(xyz_rxN[:,0],xyz_rxN[:,1], 'r.', ms = 3)
rx = DC.RxDipole(xyz_rxP, xyz_rxN)
src = DC.SrcDipole([rx], [-200, 0, -12.5], [+200, 0, -12.5])
survey = DC.SurveyDC([src])
problem = DC.ProblemDC_CC(mesh)
problem.pair(survey)
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except Exception, e:
pass
data = survey.dpred(sigma)
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_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:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1,2, figsize = (12, 5))
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_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_ana)/np.linalg.norm(data_ana)
if __name__ == '__main__':
print run(plotIt=True)
+187
View File
@@ -0,0 +1,187 @@
from SimPEG import Mesh, Utils, np, sp
import SimPEG.DCIP as DC
import time
def run(loc=None, sig=None, radi=None, param=None, stype='dpdp', plotIt=True):
"""
DC Forward Simulation
=====================
Forward model conductive spheres in a half-space and plot a pseudo-section
Created by @fourndo on Mon Feb 01 19:28:06 2016
"""
assert stype in ['pdp', 'dpdp'], "Source type (stype) must be pdp or dpdp (pole dipole or dipole dipole)"
if loc is None:
loc = np.c_[[-50.,0.,-50.],[50.,0.,-50.]]
if sig is None:
sig = np.r_[1e-2,1e-1,1e-3]
if radi is None:
radi = np.r_[25.,25.]
if param is None:
param = np.r_[30.,30.,5]
# First we need to create a mesh and a model.
# This is our mesh
dx = 5.
hxind = [(dx,15,-1.3), (dx, 75), (dx,15,1.3)]
hyind = [(dx,15,-1.3), (dx, 10), (dx,15,1.3)]
hzind = [(dx,15,-1.3),(dx, 15)]
mesh = Mesh.TensorMesh([hxind, hyind, hzind], 'CCN')
# Set background conductivity
model = np.ones(mesh.nC) * sig[0]
# First anomaly
ind = Utils.ModelBuilder.getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
model[ind] = sig[1]
# Second anomaly
ind = Utils.ModelBuilder.getIndicesSphere(loc[:,1],radi[1],mesh.gridCC)
model[ind] = sig[2]
# Get index of the center
indy = int(mesh.nCy/2)
# Plot the model for reference
# Define core mesh extent
xlim = 200
zlim = 125
# Specify the survey type: "pdp" | "dpdp"
# Then specify the end points of the survey. Let's keep it simple for now and survey above the anomalies, top of the mesh
ends = [(-175,0),(175,0)]
ends = np.c_[np.asarray(ends),np.ones(2).T*mesh.vectorNz[-1]]
# Snap the endpoints to the grid. Easier to create 2D section.
indx = Utils.closestPoints(mesh, ends )
locs = np.c_[mesh.gridCC[indx,0],mesh.gridCC[indx,1],np.ones(2).T*mesh.vectorNz[-1]]
# We will handle the geometry of the survey for you and create all the combination of tx-rx along line
# [Tx, Rx] = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
survey, Tx, Rx = DC.gen_DCIPsurvey(locs, mesh, stype, param[0], param[1], param[2])
# Define some global geometry
dl_len = np.sqrt( np.sum((locs[0,:] - locs[1,:])**2) )
dl_x = ( Tx[-1][0,1] - Tx[0][0,0] ) / dl_len
dl_y = ( Tx[-1][1,1] - Tx[0][1,0] ) / dl_len
azm = np.arctan(dl_y/dl_x)
#Set boundary conditions
mesh.setCellGradBC('neumann')
# Define the differential operators needed for the DC problem
Div = mesh.faceDiv
Grad = mesh.cellGrad
Msig = Utils.sdiag(1./(mesh.aveF2CC.T*(1./model)))
A = Div*Msig*Grad
# Change one corner to deal with nullspace
A[0,0] = 1
A = sp.csc_matrix(A)
# We will solve the system iteratively, so a pre-conditioner is helpful
# This is simply a Jacobi preconditioner (inverse of the main diagonal)
dA = A.diagonal()
P = sp.spdiags(1/dA,0,A.shape[0],A.shape[0])
# Now we can solve the system for all the transmitters
# We want to store the data
data = []
# There is probably a more elegant way to do this, but we can just for-loop through the transmitters
for ii in range(len(Tx)):
start_time = time.time() # Let's time the calculations
#print("Transmitter %i / %i\r" % (ii+1,len(Tx)))
# Select dipole locations for receiver
rxloc_M = np.asarray(Rx[ii][:,0:3])
rxloc_N = np.asarray(Rx[ii][:,3:])
# For usual cases "dpdp" or "gradient"
if stype == 'pdp':
# Create an "inifinity" pole
tx = np.squeeze(Tx[ii][:,0:1])
tinf = tx + np.array([dl_x,dl_y,0])*dl_len*2
inds = Utils.closestPoints(mesh, np.c_[tx,tinf].T)
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1] / mesh.vol[inds] )
else:
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
# Iterative Solve
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
# We now have the potential everywhere
phi = Utils.mkvc(Ainvb[0])
# Solve for phi on pole locations
P1 = mesh.getInterpolationMat(rxloc_M, 'CC')
P2 = mesh.getInterpolationMat(rxloc_N, 'CC')
# Compute the potential difference
dtemp = (P1*phi - P2*phi)*np.pi
data.append( dtemp )
print '\rTransmitter {0} of {1} -> Time:{2} sec'.format(ii,len(Tx),time.time()- start_time),
print 'Transmitter {0} of {1}'.format(ii,len(Tx))
print 'Forward completed'
# Let's just convert the 3D format into 2D (distance along line) and plot
# [Tx2d, Rx2d] = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
survey2D = DC.convertObs_DC3D_to_2D(survey, np.ones(survey.nSrc))
survey2D.dobs =np.hstack(data)
# Here is an example for the first tx-rx array
if plotIt:
import matplotlib.pyplot as plt
fig = plt.figure()
ax = plt.subplot(2,1,1, aspect='equal')
mesh.plotSlice(np.log10(model), ax =ax, normal = 'Y', ind = indy,grid=True)
ax.set_title('E-W section at '+str(mesh.vectorCCy[indy])+' m')
plt.gca().set_aspect('equal', adjustable='box')
plt.scatter(Tx[0][0,:],Tx[0][2,:],s=40,c='g', marker='v')
plt.scatter(Rx[0][:,0::3],Rx[0][:,2::3],s=40,c='y')
plt.xlim([-xlim,xlim])
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
ax = plt.subplot(2,1,2, aspect='equal')
# Plot the location of the spheres for reference
circle1=plt.Circle((loc[0,0]-Tx[0][0,0],loc[2,0]),radi[0],color='w',fill=False, lw=3)
circle2=plt.Circle((loc[0,1]-Tx[0][0,0],loc[2,1]),radi[1],color='k',fill=False, lw=3)
ax.add_artist(circle1)
ax.add_artist(circle2)
# Add the speudo section
DC.plot_pseudoSection(survey2D,ax,stype)
# plt.scatter(Tx2d[0][:],Tx[0][2,:],s=40,c='g', marker='v')
# plt.scatter(Rx2d[0][:],Rx[0][:,2::3],s=40,c='y')
# plt.plot(np.r_[Tx2d[0][0],Rx2d[-1][-1,-1]],np.ones(2)*mesh.vectorNz[-1], color='k')
plt.ylim([-zlim,mesh.vectorNz[-1]+dx])
plt.show()
return fig, ax
if __name__ == '__main__':
run()
+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
+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
+2 -24
View File
@@ -10,28 +10,6 @@ def run(N=100, plotIt=True):
"""
class LinearSurvey(Survey.BaseSurvey):
def projectFields(self, u):
return u
class LinearProblem(Problem.BaseProblem):
surveyPair = LinearSurvey
def __init__(self, mesh, G, **kwargs):
Problem.BaseProblem.__init__(self, mesh, **kwargs)
self.G = G
def fields(self, m, u=None):
return self.G.dot(m)
def Jvec(self, m, v, u=None):
return self.G.dot(v)
def Jtvec(self, m, v, u=None):
return self.G.T.dot(v)
np.random.seed(1)
mesh = Mesh.TensorMesh([N])
@@ -53,8 +31,8 @@ def run(N=100, plotIt=True):
mtrue[mesh.vectorCCx > 0.45] = -0.5
mtrue[mesh.vectorCCx > 0.6] = 0
prob = LinearProblem(mesh, G)
survey = LinearSurvey()
prob = Problem.LinearProblem(mesh, G)
survey = Survey.LinearSurvey()
survey.pair(prob)
survey.makeSyntheticData(mtrue, std=0.01)
@@ -0,0 +1,129 @@
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
from scipy.constants import mu_0
import matplotlib.pyplot as plt
def run(plotIt=True):
"""
MT: 1D: Inversion
=======================
Forward model 1D MT data.
Setup and run a MT 1D inversion.
"""
## Setup the forward modeling
# Setting up 1D mesh and conductivity models to forward model data.
# Frequency
nFreq = 31
freqs = np.logspace(3,-3,nFreq)
# Set mesh parameters
ct = 20
air = simpeg.Utils.meshTensor([(ct,16,1.4)])
core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Make the model
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Setup model varibles
active = m1d.vectorCCx<0.
layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
# Set the conductivity values
sig_half = 2e-3
sig_air = 1e-8
sig_layer1 = .2
sig_layer2 = .2
# Make the true model
sigma_true = np.ones(m1d.nCx)*sig_air
sigma_true[active] = sig_half
sigma_true[layer1] = sig_layer1
sigma_true[layer2] = sig_layer2
# Extract the model
m_true = np.log(sigma_true[active])
# Make the background model
sigma_0 = np.ones(m1d.nCx)*sig_air
sigma_0[active] = sig_half
m_0 = np.log(sigma_0[active])
# Set the mapping
actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
## Setup the layout of the survey, set the sources and the connected receivers
# Receivers
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Make the survey
survey = MT.Survey(srcList)
survey.mtrue = m_true
## Set the problem
problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
problem.pair(survey)
## Forward model data
# Project the data
survey.dtrue = survey.dpred(m_true)
survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
fig.suptitle('Target - smooth true')
# Assign uncertainties
std = 0.05 # 5% std
survey.std = np.abs(survey.dobs*std)
# Assign the data weight
Wd = 1./survey.std
## Setup the inversion proceedure
# Define a counter
C = simpeg.Utils.Counter()
# Set the optimization
opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
opt.counter = C
opt.LSshorten = 0.5
opt.remember('xc')
# Data misfit
dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
dmis.Wd = Wd
# Regularization - with a regularization mesh
regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
reg = simpeg.Regularization.Tikhonov(regMesh)
reg.smoothModel = True
reg.alpha_s = 1e-7
reg.alpha_x = 1.
# Inversion problem
invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
invProb.counter = C
# Beta cooling
beta = simpeg.Directives.BetaSchedule()
beta.coolingRate = 4
betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
targmis = simpeg.Directives.TargetMisfit()
targmis.target = survey.nD
saveModel = simpeg.Directives.SaveModelEveryIteration()
saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
# Create an inversion object
inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
## Run the inversion
mopt = inv.run(m_0)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
fig.suptitle('Target - smooth true')
plt.show()
if __name__ == '__main__':
run()
+64
View File
@@ -0,0 +1,64 @@
# Test script to use SimPEG.MT platform to forward model synthetic data.
# Import
import SimPEG as simpeg
from SimPEG import MT
import numpy as np
try:
from pymatsolver import MumpsSolver as Solver
except:
from SimPEG import Solver
def run(plotIt=True, nFreq=1):
"""
MT: 3D: Forward
=======================
Forward model 3D MT data.
"""
# Make a mesh
M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
# Setup the model
conds = [1e-2,1]
sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds)
sig[M.gridCC[:,2]>0] = 1e-8
sig[M.gridCC[:,2]<-600] = 1e-1
sigBG = np.zeros(M.nC) + conds[0]
sigBG[M.gridCC[:,2]>0] = 1e-8
## Setup the the survey object
# Receiver locations
rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50))
rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1))))
# Make a receiver list
rxList = []
for loc in rx_loc:
# NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
srcList =[]
for freq in np.logspace(3,-3,nFreq):
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
## Setup the problem object
problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
problem.pair(survey)
problem.Solver = Solver
# Calculate the data
fields = problem.fields(sig)
dataVec = survey.eval(fields)
# Make the data
mtData = MT.Data(survey,dataVec)
# Add plots
if plotIt:
pass
if __name__ == '__main__':
run()
+5 -1
View File
@@ -1,6 +1,8 @@
# Run this file to add imports.
##### AUTOIMPORTS #####
import DC_Analytic_Dipole
import DC_Forward_PseudoSection
import EM_FDEM_1D_Inversion
import EM_FDEM_Analytic_MagDipoleWholespace
import EM_TDEM_1D_Inversion
@@ -14,8 +16,10 @@ import Mesh_QuadTree_Creation
import Mesh_QuadTree_FaceDiv
import Mesh_QuadTree_HangingNodes
import Mesh_Tensor_Creation
import MT_1D_ForwardAndInversion
import MT_3D_Foward
__examples__ = ["EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation"]
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
##### AUTOIMPORTS #####
+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.
+132
View File
@@ -0,0 +1,132 @@
from SimPEG import SolverLU as SimpegSolver, PropMaps, Utils, mkvc, sp, np
from SimPEG.EM.FDEM.FDEM import BaseFDEMProblem
from SurveyMT import Survey, Data
from FieldsMT import BaseMTFields
class BaseMTProblem(BaseFDEMProblem):
"""
Base class for all Natural source problems.
"""
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
Utils.setKwargs(self, **kwargs)
# Set the default pairs of the problem
surveyPair = Survey
dataPair = Data
fieldsPair = BaseMTFields
# Set the solver
Solver = SimpegSolver
solverOpts = {}
verbose = False
# Notes:
# Use the forward and devs from BaseFDEMProblem
# Might need to add more stuff here.
## NEED to clean up the Jvec and Jtvec to use Zero and Identities for None components.
def Jvec(self, m, v, u=None):
"""
Function to calculate the data sensitivities dD/dm times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nC, 1) - random vector
:param MTfields object (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
# Calculate the fields
if u is None:
u = self.fields(m)
# Set current model
self.curModel = m
# Initiate the Jv object
Jv = self.dataPair(self.survey)
# Loop all the frequenies
for freq in self.survey.freqs:
dA_du = self.getA(freq) #
dA_duI = self.Solver(dA_du, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
# We need fDeriv_m = df/du*du/dm + df/dm
# Construct du/dm, it requires a solve
# NOTE: need to account for the 2 polarizations in the derivatives.
u_src = u[src,:]
# dA_dm and dRHS_dm should be of size nE,2, so that we can multiply by dA_duI. The 2 columns are each of the polarizations.
dA_dm = self.getADeriv_m(freq, u_src, v) # Size: nE,2 (u_px,u_py) in the columns.
dRHS_dm = self.getRHSDeriv_m(freq, v) # Size: nE,2 (u_px,u_py) in the columns.
if dRHS_dm is None:
du_dm = dA_duI * ( -dA_dm )
else:
du_dm = dA_duI * ( -dA_dm + dRHS_dm )
# Calculate the projection derivatives
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
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
return mkvc(Jv)
def Jtvec(self, m, v, u=None):
"""
Function to calculate the transpose of the data sensitivities (dD/dm)^T times a vector.
:param numpy.ndarray m (nC, 1) - conductive model
:param numpy.ndarray v (nD, 1) - vector
:param MTfields object u (optional) - MT fields object, if not given it is calculated
:rtype: MTdata object
:return: Data sensitivities wrt m
"""
if u is None:
u = self.fields(m)
self.curModel = m
# Ensure v is a data object.
if not isinstance(v, self.dataPair):
v = self.dataPair(self.survey, v)
Jtv = np.zeros(m.size)
for freq in self.survey.freqs:
AT = self.getA(freq).T
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, :]
for rx in src.rxList:
# Get the adjoint evalDeriv
# PTv needs to be nE,
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)
dRHS_dmT = self.getRHSDeriv_m(freq, mkvc(dA_duIT), adjoint=True)
# Make du_dmT
if dRHS_dmT is None:
du_dmT = -dA_dmT
else:
du_dmT = -dA_dmT + dRHS_dmT
# Select the correct component
# du_dmT needs to be of size nC,
real_or_imag = rx.projComp
if real_or_imag == 'real':
Jtv += du_dmT.real
elif real_or_imag == 'imag':
Jtv += -du_dmT.real
else:
raise Exception('Must be real or imag')
# Clean the factorization, clear memory.
ATinv.clean()
return Jtv
+351
View File
@@ -0,0 +1,351 @@
from SimPEG import Survey, Utils, Problem, np, sp, mkvc
from scipy.constants import mu_0
import sys
from numpy.lib import recfunctions as recFunc
from SimPEG.EM.Utils import omega
##############
### Fields ###
##############
class BaseMTFields(Problem.Fields):
"""Field Storage for a MT survey."""
knownFields = {}
dtype = complex
class Fields1D_e(BaseMTFields):
"""
Fields storage for the 1D MT solution.
"""
knownFields = {'e_1dSolution':'F'}
aliasFields = {
'e_1d' : ['e_1dSolution','F','_e'],
'e_1dPrimary' : ['e_1dSolution','F','_ePrimary'],
'e_1dSecondary' : ['e_1dSolution','F','_eSecondary'],
'b_1d' : ['e_1dSolution','E','_b'],
'b_1dPrimary' : ['e_1dSolution','E','_bPrimary'],
'b_1dSecondary' : ['e_1dSolution','E','_bSecondary']
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
def _ePrimary(self, eSolution, srcList):
ePrimary = np.zeros_like(eSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
if ep is not None:
ePrimary[:,i] = ep[:,-1]
return ePrimary
def _eSecondary(self, eSolution, srcList):
return eSolution
def _e(self, eSolution, srcList):
return self._ePrimary(eSolution,srcList) + self._eSecondary(eSolution,srcList)
def _eDeriv_u(self, src, v, adjoint = False):
return v
def _eDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return None
def _bPrimary(self, eSolution, srcList):
bPrimary = np.zeros([self.survey.mesh.nE,eSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
if bp is not None:
bPrimary[:,i] += bp[:,-1]
return bPrimary
def _bSecondary(self, eSolution, srcList):
C = self.mesh.nodalGrad
b = (C * eSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
return b
def _b(self, eSolution, srcList):
return self._bPrimary(eSolution, srcList) + self._bSecondary(eSolution, srcList)
def _bSecondaryDeriv_u(self, src, v, adjoint = False):
C = self.mesh.nodalGrad
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _bSecondaryDeriv_m(self, src, v, adjoint = False):
# Doesn't depend on m
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
# S_eDeriv = S_eDeriv(v)
# if S_eDeriv is not None:
# return 1./(1j * omega(src.freq)) * S_eDeriv
return None
def _bDeriv_u(self, src, v, adjoint=False):
# Primary does not depend on u
return self._bSecondaryDeriv_u(src, v, adjoint)
def _bDeriv_m(self, src, v, adjoint=False):
# Assuming the primary does not depend on the model
return self._bSecondaryDeriv_m(src, v, adjoint)
def _fDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
return a vector of size (nreEle+nrbEle)
"""
de_du = v #Utils.spdiag(np.ones((self.nF,)))
db_du = self._bDeriv_u(src, v, adjoint)
# Return the stack
# This doesn't work...
return np.vstack((de_du,db_du))
def _fDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt m.
This function stacks the fields derivatives appropriately
"""
return None
class Fields3D_e(BaseMTFields):
"""
Fields storage for the 3D MT solution. Labels polarizations by px and py.
:param SimPEG object mesh: The solution mesh
:param SimPEG object survey: A survey object
"""
# Define the known the alias fields
# Assume that the solution of e on the E.
## NOTE: Need to make this more general, to allow for other solutions formats.
knownFields = {'e_pxSolution':'E','e_pySolution':'E'}
aliasFields = {
'e_px' : ['e_pxSolution','E','_e_px'],
'e_pxPrimary' : ['e_pxSolution','E','_e_pxPrimary'],
'e_pxSecondary' : ['e_pxSolution','E','_e_pxSecondary'],
'e_py' : ['e_pySolution','E','_e_py'],
'e_pyPrimary' : ['e_pySolution','E','_e_pyPrimary'],
'e_pySecondary' : ['e_pySolution','E','_e_pySecondary'],
'b_px' : ['e_pxSolution','F','_b_px'],
'b_pxPrimary' : ['e_pxSolution','F','_b_pxPrimary'],
'b_pxSecondary' : ['e_pxSolution','F','_b_pxSecondary'],
'b_py' : ['e_pySolution','F','_b_py'],
'b_pyPrimary' : ['e_pySolution','F','_b_pyPrimary'],
'b_pySecondary' : ['e_pySolution','F','_b_pySecondary']
}
def __init__(self,mesh,survey,**kwargs):
BaseMTFields.__init__(self,mesh,survey,**kwargs)
def _e_pxPrimary(self, e_pxSolution, srcList):
e_pxPrimary = np.zeros_like(e_pxSolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
if ep is not None:
e_pxPrimary[:,i] = ep[:,0]
return e_pxPrimary
def _e_pyPrimary(self, e_pySolution, srcList):
e_pyPrimary = np.zeros_like(e_pySolution)
for i, src in enumerate(srcList):
ep = src.ePrimary(self.survey.prob)
if ep is not None:
e_pyPrimary[:,i] = ep[:,1]
return e_pyPrimary
def _e_pxSecondary(self, e_pxSolution, srcList):
return e_pxSolution
def _e_pySecondary(self, e_pySolution, srcList):
return e_pySolution
def _e_px(self, e_pxSolution, srcList):
return self._e_pxPrimary(e_pxSolution,srcList) + self._e_pxSecondary(e_pxSolution,srcList)
def _e_py(self, e_pySolution, srcList):
return self._e_pyPrimary(e_pySolution,srcList) + self._e_pySecondary(e_pySolution,srcList)
#NOTE: For e_p?Deriv_u,
# v has to be u(2*nE) long for the not adjoint and nE long for adjoint.
# Returns nE long for not adjoint and 2*nE long for adjoint
def _e_pxDeriv_u(self, src, v, adjoint = False):
'''
Takes the derivative of e_px wrt u
'''
if adjoint:
# adjoint: returns a 2*nE long vector with zero's for py
return np.vstack((v,np.zeros_like(v)))
# Not adjoint: return only the px part of the vector
return v[:len(v)/2]
def _e_pyDeriv_u(self, src, v, adjoint = False):
'''
Takes the derivative of e_py wrt u
'''
if adjoint:
# adjoint: returns a 2*nE long vector with zero's for px
return np.vstack((np.zeros_like(v),v))
# Not adjoint: return only the px part of the vector
return v[len(v)/2::]
def _e_pxDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return None
def _e_pyDeriv_m(self, src, v, adjoint = False):
# assuming primary does not depend on the model
return None
def _b_pxPrimary(self, e_pxSolution, srcList):
b_pxPrimary = np.zeros([self.survey.mesh.nF,e_pxSolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
if bp is not None:
b_pxPrimary[:,i] += bp[:,0]
return b_pxPrimary
def _b_pyPrimary(self, e_pySolution, srcList):
b_pyPrimary = np.zeros([self.survey.mesh.nF,e_pySolution.shape[1]], dtype = complex)
for i, src in enumerate(srcList):
bp = src.bPrimary(self.survey.prob)
if bp is not None:
b_pyPrimary[:,i] += bp[:,1]
return b_pyPrimary
def _b_pxSecondary(self, e_pxSolution, srcList):
C = self.mesh.edgeCurl
b = (C * e_pxSolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
return b
def _b_pySecondary(self, e_pySolution, srcList):
C = self.mesh.edgeCurl
b = (C * e_pySolution)
for i, src in enumerate(srcList):
b[:,i] *= - 1./(1j*omega(src.freq))
# There is no magnetic source in the MT problem
# S_m, _ = src.eval(self.survey.prob)
# if S_m is not None:
# b[:,i] += 1./(1j*omega(src.freq)) * S_m
return b
def _b_px(self, eSolution, srcList):
return self._b_pxPrimary(eSolution, srcList) + self._b_pxSecondary(eSolution, srcList)
def _b_py(self, eSolution, srcList):
return self._b_pyPrimary(eSolution, srcList) + self._b_pySecondary(eSolution, srcList)
# NOTE: v needs to be length 2*nE to account for both polarizations
def _b_pxSecondaryDeriv_u(self, src, v, adjoint = False):
# C = sp.kron(self.mesh.edgeCurl,[[1,0],[0,0]])
C = sp.hstack((self.mesh.edgeCurl,Utils.spzeros(self.mesh.nF,self.mesh.nE))) # This works for adjoint = None
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _b_pySecondaryDeriv_u(self, src, v, adjoint = False):
# C = sp.kron(self.mesh.edgeCurl,[[0,0],[0,1]])
C = sp.hstack((Utils.spzeros(self.mesh.nF,self.mesh.nE),self.mesh.edgeCurl)) # This works for adjoint = None
if adjoint:
return - 1./(1j*omega(src.freq)) * (C.T * v)
return - 1./(1j*omega(src.freq)) * (C * v)
def _b_pxSecondaryDeriv_m(self, src, v, adjoint = False):
# Doesn't depend on m
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
# S_eDeriv = S_eDeriv(v)
# if S_eDeriv is not None:
# return 1./(1j * omega(src.freq)) * S_eDeriv
return None
def _b_pySecondaryDeriv_m(self, src, v, adjoint = False):
# Doesn't depend on m
# _, S_eDeriv = src.evalDeriv(self.survey.prob, adjoint)
# S_eDeriv = S_eDeriv(v)
# if S_eDeriv is not None:
# return 1./(1j * omega(src.freq)) * S_eDeriv
return None
def _b_pxDeriv_u(self, src, v, adjoint=False):
# Primary does not depend on u
return self._b_pxSecondaryDeriv_u(src, v, adjoint)
def _b_pyDeriv_u(self, src, v, adjoint=False):
# Primary does not depend on u
return self._b_pySecondaryDeriv_u(src, v, adjoint)
def _b_pxDeriv_m(self, src, v, adjoint=False):
# Assuming the primary does not depend on the model
return self._b_pxSecondaryDeriv_m(src, v, adjoint)
def _b_pyDeriv_m(self, src, v, adjoint=False):
# Assuming the primary does not depend on the model
return self._b_pySecondaryDeriv_m(src, v, adjoint)
def _f_pxDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
return a vector of size (nreEle+nrbEle)
"""
de_du = v #Utils.spdiag(np.ones((self.nF,)))
db_du = self._b_pxDeriv_u(src, v, adjoint)
# Return the stack
# This doesn't work...
return np.vstack((de_du,db_du))
def _f_pyDeriv_u(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt u.
:param MTsrc src: MT source
:param numpy.ndarray v: random vector of f_sol.size
This function stacks the fields derivatives appropriately
return a vector of size (nreEle+nrbEle)
"""
de_du = v #Utils.spdiag(np.ones((self.nF,)))
db_du = self._b_pyDeriv_u(src, v, adjoint)
# Return the stack
# This doesn't work...
return np.vstack((de_du,db_du))
def _f_pxDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt m.
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
def _f_pyDeriv_m(self, src, v, adjoint=False):
"""
Derivative of the fields object wrt m.
This function stacks the fields derivatives appropriately
"""
# The fields have no dependance to the model.
return None
+291
View File
@@ -0,0 +1,291 @@
from SimPEG.EM.Utils import omega
from SimPEG import mkvc
from scipy.constants import mu_0
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields1D_e
from SimPEG.MT.Utils.MT1Danalytic import getEHfields
import numpy as np
import multiprocessing, sys, time
class eForm_psField(BaseMTProblem):
"""
A MT problem soving a e formulation and primary/secondary fields decomposion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^e_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^f_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^f_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly half space ).
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e_1d'
_eqLocs = 'EF'
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
self.fieldsPair = Fields1D_e
# self._sigmaPrimary = sigmaPrimary
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A matrix.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
MeMui = self.MeMui
MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
The derivative of A wrt sigma
"""
dsig_dm = self.curModel.sigmaDeriv
MeMui = self.MeMui
#
u_src = u['e_1dSolution']
dMfSigma_dm = self.mesh.getFaceInnerProductDeriv(self.curModel.sigma)(u_src) * self.curModel.sigmaDeriv
if adjoint:
return 1j * omega(freq) * ( dMfSigma_dm.T * v )
# Note: output has to be nN/nF, not nC/nE.
# v should be nC
return 1j * omega(freq) * ( dMfSigma_dm * v )
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nF, 1), numpy.ndarray (nF, 1)
:return: RHS for 1 polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS wrt sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# NOTE: only store the e_solution(secondary), all other components calculated in the fields object
F[Src, 'e_1dSolution'] = e_s[:,-1] # Only storing the yx polarization as 1d
# Note curl e = -iwb so b = -curl e /iw
# b = -( self.mesh.nodalGrad * e )/( 1j*omega(freq) )
# F[Src, 'b_1d'] = b[:,1]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
# Note this is not fully functional.
# Missing:
# Fields class corresponding to the fields
# Update Jvec and Jtvec to include all the derivatives components
# Other things ...
class eForm_TotalField(BaseMTProblem):
"""
A MT problem solving a e formulation and a Total bondary domain decompostion.
Solves the equation:
Math:
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'EF'
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def MeMui(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MeMui', None) is None:
self._MeMui = self.mesh.getEdgeInnerProduct(1.0/mu_0)
return self._MeMui
@property
def MfSigma(self):
"""
Edge inner product matrix
"""
if getattr(self, '_MfSigma', None) is None:
self._MfSigma = self.mesh.getFaceInnerProduct(self.curModel.sigma)
return self._MfSigma
def getA(self, freq, full=False):
"""
Function to get the A matrix.
:param float freq: Frequency
:param logic full: Return full A or the inner part
:rtype: scipy.sparse.csr_matrix
:return: A
"""
MeMui = self.MeMui
MfSigma = self.MfSigma
# Note: need to use the code above since in the 1D problem I want
# e to live on Faces(nodes) and h on edges(cells). Might need to rethink this
# Possible that _fieldType and _eqLocs can fix this
# MeMui = self.MfMui
# MfSigma = self.MfSigma
C = self.mesh.nodalGrad
# Make A
A = C.T*MeMui*C + 1j*omega(freq)*MfSigma
# Either return full or only the inner part of A
if full:
return A
else:
return A[1:-1,1:-1]
def getADeriv_m(self, freq, u, v, adjoint=False):
raise NotImplementedError('getADeriv is not implemented')
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequency
# NOTE: Need to use the source information, doesn't really apply in 1D
src = self.survey.getSrcByFreq(freq)
# Get the full A
A = self.getA(freq,full=True)
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
Ed, Eu, Hd, Hu = getEHfields(self.mesh,self.curModel.sigma,freq,self.mesh.vectorNx)
Etot = (Ed + Eu)
sourceAmp = 1.0
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
eBC = np.r_[Etot[0],Etot[-1]]
# The right hand side
return -Aio*eBC, eBC
def getRHSderiv_m(self, freq, backSigma, u, v, adjoint=False):
raise NotImplementedError('getRHSDeriv not implemented yet')
return None
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
:param np.ndarray (nC,) m_back: Background conductivity model
'''
self.curModel = m
# RHS, CalcFields = self.getRHS(freq,m_back), self.calcFields
F = Fields1D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs, e_o = self.getRHS(freq)
Ainv = self.Solver(A, **self.solverOpts)
e_i = Ainv * rhs
e = mkvc(np.r_[e_o[0], e_i, e_o[1]],2)
# Store the fields
Src = self.survey.getSrcByFreq(freq)
# NOTE: only store e fields
F[Src, 'e_1dSolution'] = e[:,0]
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
return F
+1
View File
@@ -0,0 +1 @@
from Probs import eForm_TotalField, eForm_psField
View File
+1
View File
@@ -0,0 +1 @@
pass
+138
View File
@@ -0,0 +1,138 @@
from SimPEG import Survey, Problem, Utils, Models, np, sp, mkvc, SolverLU as SimpegSolver
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from SimPEG.MT.BaseMT import BaseMTProblem
from SimPEG.MT.SurveyMT import Survey, Data
from SimPEG.MT.FieldsMT import Fields3D_e
import multiprocessing, sys, time
class eForm_ps(BaseMTProblem):
"""
A MT problem solving a e formulation and a primary/secondary fields decompostion.
By eliminating the magnetic flux density using
.. math ::
\mathbf{b} = \\frac{1}{i \omega}\\left(-\mathbf{C} \mathbf{e} \\right)
we can write Maxwell's equations as a second order system in \\\(\\\mathbf{e}\\\) only:
.. math ::
\\left(\mathbf{C}^T \mathbf{M^f_{\mu^{-1}}} \mathbf{C} + i \omega \mathbf{M^e_\sigma}] \mathbf{e}_{s} =& i \omega \mathbf{M^e_{\delta \sigma}} \mathbf{e}_{p}
which we solve for \\\(\\\mathbf{e_s}\\\). The total field \\\mathbf{e}\\ = \\\mathbf{e_p}\\ + \\\mathbf{e_s}\\.
The primary field is estimated from a background model (commonly as a 1D model).
"""
# From FDEMproblem: Used to project the fields. Currently not used for MTproblem.
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields3D_e
_sigmaPrimary = None
def __init__(self, mesh, **kwargs):
BaseMTProblem.__init__(self, mesh, **kwargs)
@property
def sigmaPrimary(self):
"""
A background model, use for the calculation of the primary fields.
"""
return self._sigmaPrimary
@sigmaPrimary.setter
def sigmaPrimary(self, val):
# Note: TODO add logic for val, make sure it is the correct size.
self._sigmaPrimary = val
def getA(self, freq):
"""
Function to get the A system.
:param float freq: Frequency
:rtype: scipy.sparse.csr_matrix
:return: A
"""
Mmui = self.MfMui
Msig = self.MeSigma
C = self.mesh.edgeCurl
return C.T*Mmui*C + 1j*omega(freq)*Msig
def getADeriv_m(self, freq, u, v, adjoint=False):
"""
Calculate the derivative of A wrt m.
"""
# This considers both polarizations and returns a nE,2 matrix for each polarization
if adjoint:
dMe_dsigV = sp.hstack(( self.MeSigmaDeriv( u['e_pxSolution'] ).T, self.MeSigmaDeriv(u['e_pySolution'] ).T ))*v
else:
# Need a nE,2 matrix to be returned
dMe_dsigV = np.hstack(( mkvc(self.MeSigmaDeriv( u['e_pxSolution'] )*v,2), mkvc( self.MeSigmaDeriv(u['e_pySolution'] )*v,2) ))
return 1j * omega(freq) * dMe_dsigV
def getRHS(self, freq):
"""
Function to return the right hand side for the system.
:param float freq: Frequency
:rtype: numpy.ndarray (nE, 2), numpy.ndarray (nE, 2)
:return: RHS for both polarizations, primary fields
"""
# Get sources for the frequncy(polarizations)
Src = self.survey.getSrcByFreq(freq)[0]
S_e = Src.S_e(self)
return -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, v, adjoint=False):
"""
The derivative of the RHS with respect to sigma
"""
Src = self.survey.getSrcByFreq(freq)[0]
S_eDeriv = Src.S_eDeriv_m(self, v, adjoint)
return -1j * omega(freq) * S_eDeriv
def fields(self, m):
'''
Function to calculate all the fields for the model m.
:param np.ndarray (nC,) m: Conductivity model
'''
# Set the current model
self.curModel = m
F = Fields3D_e(self.mesh, self.survey)
for freq in self.survey.freqs:
if self.verbose:
startTime = time.time()
print 'Starting work for {:.3e}'.format(freq)
sys.stdout.flush()
A = self.getA(freq)
rhs = self.getRHS(freq)
# Solve the system
Ainv = self.Solver(A, **self.solverOpts)
e_s = Ainv * rhs
# Store the fields
Src = self.survey.getSrcByFreq(freq)[0]
# Store the fieldss
F[Src, 'e_pxSolution'] = e_s[:,0]
F[Src, 'e_pySolution'] = e_s[:,1]
# Note curl e = -iwb so b = -curl/iw
if self.verbose:
print 'Ran for {:f} seconds'.format(time.time()-startTime)
sys.stdout.flush()
Ainv.clean()
return F
+1
View File
@@ -0,0 +1 @@
from Probs import eForm_ps
+206
View File
@@ -0,0 +1,206 @@
from SimPEG import Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Utils.sourceUtils import homo1DModelSource
from Utils import rec2ndarr
import sys
#################
### Sources ###
#################
class BaseMTSrc(FDEMBaseSrc):
'''
Sources for the MT problem.
Use the SimPEG BaseSrc, since the source fields share properties with the transmitters.
:param float freq: The frequency of the source
:param list rxList: A list of receivers associated with the source
'''
freq = None #: Frequency (float)
def __init__(self, rxList, freq):
self.freq = float(freq)
FDEMBaseSrc.__init__(self, rxList)
# 1D sources
class polxy_1DhomotD(BaseMTSrc):
"""
MT source for both polarizations (x and y) for the total Domain.
It calculates fields calculated based on conditions on the boundary of the domain.
"""
def __init__(self, rxList, freq):
BaseMTSrc.__init__(self, rxList, freq)
# TODO: need to add the primary fields calc and source terms into the problem.
# Need to implement such that it works for all dims.
class polxy_1Dprimary(BaseMTSrc):
"""
MT source for both polarizations (x and y) given a 1D primary models.
It assigns fields calculated from the 1D model as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigma1d = None
BaseMTSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
def ePrimary(self,problem):
# Get primary fields for both polarizations
if self.sigma1d is None:
# Set the sigma1d as the 1st column in the background model
if len(problem._sigmaPrimary) == problem.mesh.nC:
if problem.mesh.dim == 1:
self.sigma1d = problem.mesh.r(problem._sigmaPrimary,'CC','CC','M')[:]
elif problem.mesh.dim == 3:
self.sigma1d = problem.mesh.r(problem._sigmaPrimary,'CC','CC','M')[0,0,:]
# Or as the 1D model that matches the vertical cell number
elif len(problem._sigmaPrimary) == problem.mesh.nCz:
self.sigma1d = problem._sigmaPrimary
if self._ePrimary is None:
self._ePrimary = homo1DModelSource(problem.mesh,self.freq,self.sigma1d)
return self._ePrimary
def bPrimary(self,problem):
# Project ePrimary to bPrimary
# Satisfies the primary(background) field conditions
if problem.mesh.dim == 1:
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
"""
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
# Need to deal with the edge/face discrepencies between 1d/2d/3d
if problem.mesh.dim == 1:
Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma)
Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p)
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
Mesigma = problem.MeSigma
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
return (Mesigma - Mesigma_p) * e_p
def S_eDeriv_m(self, problem, v, adjoint = False):
'''
Get the derivative of S_e wrt to sigma (m)
'''
# Need to deal with
if problem.mesh.dim == 1:
# Need to use the faceInnerProduct
MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,1]) * problem.curModel.sigmaDeriv
# MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
# Need to take the derivative of both u_px and u_py
ePri = self.ePrimary(problem)
# MsigmaDeriv = problem.MeSigmaDeriv(ePri[:,0]) + problem.MeSigmaDeriv(ePri[:,1])
# MsigmaDeriv = problem.MeSigmaDeriv(np.sum(ePri,axis=1))
if adjoint:
return sp.hstack(( problem.MeSigmaDeriv(ePri[:,0]).T, problem.MeSigmaDeriv(ePri[:,1]).T ))*v
else:
return np.hstack(( mkvc(problem.MeSigmaDeriv(ePri[:,0]) * v,2), mkvc(problem.MeSigmaDeriv(ePri[:,1])*v,2) ))
if adjoint:
#
return MsigmaDeriv.T * v
else:
# v should be nC size
return MsigmaDeriv * v
class polxy_3Dprimary(BaseMTSrc):
"""
MT source for both polarizations (x and y) given a 3D primary model. It assigns fields calculated from the 1D model
as fields in the full space of the problem.
"""
def __init__(self, rxList, freq):
# assert mkvc(self.mesh.hz.shape,1) == mkvc(sigma1d.shape,1),'The number of values in the 1D background model does not match the number of vertical cells (hz).'
self.sigmaPrimary = None
BaseMTSrc.__init__(self, rxList, freq)
# Hidden property of the ePrimary
self._ePrimary = None
def ePrimary(self,problem):
# Get primary fields for both polarizations
self.sigmaPrimary = problem._sigmaPrimary
if self._ePrimary is None:
self._ePrimary = homo3DModelSource(problem.mesh,self.sigmaPrimary,self.freq)
return self._ePrimary
def bPrimary(self,problem):
# Project ePrimary to bPrimary
# Satisfies the primary(background) field conditions
if problem.mesh.dim == 1:
C = problem.mesh.nodalGrad
elif problem.mesh.dim == 3:
C = problem.mesh.edgeCurl
bBG_bp = (- C * self.ePrimary(problem) )*(1/( 1j*omega(self.freq) ))
return bBG_bp
def S_e(self,problem):
"""
Get the electrical field source
"""
e_p = self.ePrimary(problem)
Map_sigma_p = Maps.Vertical1DMap(problem.mesh)
sigma_p = Map_sigma_p._transform(self.sigma1d)
# Make mass matrix
# Note: M(sig) - M(sig_p) = M(sig - sig_p)
# Need to deal with the edge/face discrepencies between 1d/2d/3d
if problem.mesh.dim == 1:
Mesigma = problem.mesh.getFaceInnerProduct(problem.curModel.sigma)
Mesigma_p = problem.mesh.getFaceInnerProduct(sigma_p)
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
Mesigma = problem.MeSigma
Mesigma_p = problem.mesh.getEdgeInnerProduct(sigma_p)
return (Mesigma - Mesigma_p) * e_p
def S_eDeriv_m(self, problem, v, adjoint = False):
'''
Get the derivative of S_e wrt to sigma (m)
'''
# Need to deal with
if problem.mesh.dim == 1:
# Need to use the faceInnerProduct
MsigmaDeriv = problem.mesh.getFaceInnerProductDeriv(problem.curModel.sigma)(self.ePrimary(problem)[:,1]) * problem.curModel.sigmaDeriv
# MsigmaDeriv = ( MsigmaDeriv * MsigmaDeriv.T)**2
if problem.mesh.dim == 2:
pass
if problem.mesh.dim == 3:
# Need to take the derivative of both u_px and u_py
ePri = self.ePrimary(problem)
# MsigmaDeriv = problem.MeSigmaDeriv(ePri[:,0]) + problem.MeSigmaDeriv(ePri[:,1])
# MsigmaDeriv = problem.MeSigmaDeriv(np.sum(ePri,axis=1))
if adjoint:
return sp.hstack(( problem.MeSigmaDeriv(ePri[:,0]).T, problem.MeSigmaDeriv(ePri[:,1]).T ))*v
else:
return np.hstack(( mkvc(problem.MeSigmaDeriv(ePri[:,0]) * v,2), mkvc(problem.MeSigmaDeriv(ePri[:,1])*v,2) ))
if adjoint:
#
return MsigmaDeriv.T * v
else:
# v should be nC size
return MsigmaDeriv * v
+562
View File
@@ -0,0 +1,562 @@
from SimPEG import Survey as SimPEGsurvey, Utils, Problem, Maps, np, sp, mkvc
from SimPEG.EM.FDEM.SrcFDEM import BaseSrc as FDEMBaseSrc
from SimPEG.EM.Utils import omega
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from Utils import rec2ndarr
import SrcMT
import sys
#################
### Receivers ###
#################
class Rx(SimPEGsurvey.BaseRx):
"""
Class that defines natural source receivers.
See knownRxTypes for types of allowed receivers.
:param ndArray locs: Locations of the receivers
:param str rxType: The type of receiver
"""
knownRxTypes = {
# 3D impedance
'zxxr':['Z3D', 'real'],
'zxyr':['Z3D', 'real'],
'zyxr':['Z3D', 'real'],
'zyyr':['Z3D', 'real'],
'zxxi':['Z3D', 'imag'],
'zxyi':['Z3D', 'imag'],
'zyxi':['Z3D', 'imag'],
'zyyi':['Z3D', 'imag'],
# 2D impedance
# TODO:
# 1D impedance
'z1dr':['Z1D', 'real'],
'z1di':['Z1D', 'imag'],
# Tipper
'tzxr':['T3D','real'],
'tzxi':['T3D','imag'],
'tzyr':['T3D','real'],
'tzyi':['T3D','imag']
}
# TODO: Have locs as single or double coordinates for both or numerator and denominator separately, respectively.
def __init__(self, locs, rxType):
SimPEGsurvey.BaseRx.__init__(self, locs, rxType)
@property
def projType(self):
"""
Receiver type for projection.
"""
return self.knownRxTypes[self.rxType][0]
@property
def projComp(self):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][1]
def eval(self, src, mesh, f):
'''
Project the fields to natural source data.
:param SrcMT src: The source of the fields to project
:param SimPEG.Mesh mesh:
:param FieldsMT f: Natural source fields object to project
'''
## NOTE: Assumes that e is on t
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
ex = Pex*mkvc(f[src,'e_1d'],2)
bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
# Note: Has a minus sign in front, to comply with quadrant calculations.
# Can be derived from zyx case for the 3D case.
f_part_complex = -ex/bx
# elif self.projType is 'Z2D':
elif self.projType is 'Z3D':
## NOTE: Assumes that e is on edges and b on the faces. Need to generalize that or use a prop of fields to determine that.
if self.locs.ndim == 3:
eFLocs = self.locs[:,:,0]
bFLocs = self.locs[:,:,1]
else:
eFLocs = self.locs
bFLocs = self.locs
# Get the projection
Pex = mesh.getInterpolationMat(eFLocs,'Ex')
Pey = mesh.getInterpolationMat(eFLocs,'Ey')
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
ex_px = Pex*f[src,'e_px']
ey_px = Pey*f[src,'e_px']
ex_py = Pex*f[src,'e_py']
ey_py = Pey*f[src,'e_py']
hx_px = Pbx*f[src,'b_px']/mu_0
hy_px = Pby*f[src,'b_px']/mu_0
hx_py = Pbx*f[src,'b_py']/mu_0
hy_py = Pby*f[src,'b_py']/mu_0
# Make the complex data
if 'zxx' in self.rxType:
f_part_complex = ( ex_px*hy_py - ex_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zxy' in self.rxType:
f_part_complex = (-ex_px*hx_py + ex_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyx' in self.rxType:
f_part_complex = ( ey_px*hy_py - ey_py*hy_px)/(hx_px*hy_py - hx_py*hy_px)
elif 'zyy' in self.rxType:
f_part_complex = (-ey_px*hx_py + ey_py*hx_px)/(hx_px*hy_py - hx_py*hy_px)
elif self.projType is 'T3D':
if self.locs.ndim == 3:
horLoc = self.locs[:,:,0]
vertLoc = self.locs[:,:,1]
else:
horLoc = self.locs
vertLoc = self.locs
Pbx = mesh.getInterpolationMat(horLoc,'Fx')
Pby = mesh.getInterpolationMat(horLoc,'Fy')
Pbz = mesh.getInterpolationMat(vertLoc,'Fz')
bx_px = Pbx*f[src,'b_px']
by_px = Pby*f[src,'b_px']
bz_px = Pbz*f[src,'b_px']
bx_py = Pbx*f[src,'b_py']
by_py = Pby*f[src,'b_py']
bz_py = Pbz*f[src,'b_py']
if 'tzx' in self.rxType:
f_part_complex = (- by_px*bz_py + by_py*bz_px)/(bx_px*by_py - bx_py*by_px)
if 'tzy' in self.rxType:
f_part_complex = ( bx_px*bz_py - bx_py*bz_px)/(bx_px*by_py - bx_py*by_px)
else:
NotImplementedError('Projection of {:s} receiver type is not implemented.'.format(self.rxType))
# Get the real or imag component
real_or_imag = self.projComp
f_part = getattr(f_part_complex, real_or_imag)
# print f_part
return f_part
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
The derivative of the projection wrt u
:param MTsrc src: MT source
:param TensorMesh mesh: Mesh defining the topology of the problem
:param MTfields f: MT fields object of the source
:param numpy.ndarray v: Random vector of size
"""
real_or_imag = self.projComp
if not adjoint:
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_de = -mkvc(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Pex*v),2)
dP_db = mkvc( Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)))*(Pbx*f._bDeriv_u(src,v)/mu_0),2)
PDeriv_complex = np.sum(np.hstack((dP_de,dP_db)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not been implement for 2D impedance tensor')
elif self.projType is 'Z3D':
if self.locs.ndim == 3:
eFLocs = self.locs[:,:,0]
bFLocs = self.locs[:,:,1]
else:
eFLocs = self.locs
bFLocs = self.locs
# Get the projection
Pex = mesh.getInterpolationMat(eFLocs,'Ex')
Pey = mesh.getInterpolationMat(eFLocs,'Ey')
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
ex_px = Pex*f[src,'e_px']
ey_px = Pey*f[src,'e_px']
ex_py = Pex*f[src,'e_py']
ey_py = Pey*f[src,'e_py']
hx_px = Pbx*f[src,'b_px']/mu_0
hy_px = Pby*f[src,'b_px']/mu_0
hx_py = Pbx*f[src,'b_py']/mu_0
hy_py = Pby*f[src,'b_py']/mu_0
# Derivatives as lambda functions
# The size of the diratives should be nD,nU
ex_px_u = lambda vec: Pex*f._e_pxDeriv_u(src,vec)
ey_px_u = lambda vec: Pey*f._e_pxDeriv_u(src,vec)
ex_py_u = lambda vec: Pex*f._e_pyDeriv_u(src,vec)
ey_py_u = lambda vec: Pey*f._e_pyDeriv_u(src,vec)
# NOTE: Think b_p?Deriv_u should return a 2*nF size matrix
hx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)/mu_0
hy_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)/mu_0
hx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)/mu_0
hy_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)/mu_0
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(1./(sDiag(hx_px)*hy_py - sDiag(hx_py)*hy_px))
Hd_uV = sDiag(hy_py)*hx_px_u(v) + sDiag(hx_px)*hy_py_u(v) - sDiag(hx_py)*hy_px_u(v) - sDiag(hy_px)*hx_py_u(v)
# Calculate components
if 'zxx' in self.rxType:
Zij = sDiag(Hd*( sDiag(ex_px)*hy_py - sDiag(ex_py)*hy_px ))
ZijN_uV = sDiag(hy_py)*ex_px_u(v) + sDiag(ex_px)*hy_py_u(v) - sDiag(ex_py)*hy_px_u(v) - sDiag(hy_px)*ex_py_u(v)
elif 'zxy' in self.rxType:
Zij = sDiag(Hd*(-sDiag(ex_px)*hx_py + sDiag(ex_py)*hx_px ))
ZijN_uV = -sDiag(hx_py)*ex_px_u(v) - sDiag(ex_px)*hx_py_u(v) + sDiag(ex_py)*hx_px_u(v) + sDiag(hx_px)*ex_py_u(v)
elif 'zyx' in self.rxType:
Zij = sDiag(Hd*( sDiag(ey_px)*hy_py - sDiag(ey_py)*hy_px ))
ZijN_uV = sDiag(hy_py)*ey_px_u(v) + sDiag(ey_px)*hy_py_u(v) - sDiag(ey_py)*hy_px_u(v) - sDiag(hy_px)*ey_py_u(v)
elif 'zyy' in self.rxType:
Zij = sDiag(Hd*(-sDiag(ey_px)*hx_py + sDiag(ey_py)*hx_px ))
ZijN_uV = -sDiag(hx_py)*ey_px_u(v) - sDiag(ey_px)*hx_py_u(v) + sDiag(ey_py)*hx_px_u(v) + sDiag(hx_px)*ey_py_u(v)
# Calculate the complex derivative
PDeriv_complex = Hd * (ZijN_uV - Zij * Hd_uV )
elif self.projType is 'T3D':
if self.locs.ndim == 3:
eFLocs = self.locs[:,:,0]
bFLocs = self.locs[:,:,1]
else:
eFLocs = self.locs
bFLocs = self.locs
# Get the projection
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
Pbz = mesh.getInterpolationMat(bFLocs,'Fz')
# Get the fields at location
# px: x-polaration and py: y-polaration.
bx_px = Pbx*f[src,'b_px']
by_px = Pby*f[src,'b_px']
bz_px = Pbz*f[src,'b_px']
bx_py = Pbx*f[src,'b_py']
by_py = Pby*f[src,'b_py']
bz_py = Pbz*f[src,'b_py']
# Derivatives as lambda functions
# NOTE: Think b_p?Deriv_u should return a 2*nF size matrix
bx_px_u = lambda vec: Pbx*f._b_pxDeriv_u(src,vec)
by_px_u = lambda vec: Pby*f._b_pxDeriv_u(src,vec)
bz_px_u = lambda vec: Pbz*f._b_pxDeriv_u(src,vec)
bx_py_u = lambda vec: Pbx*f._b_pyDeriv_u(src,vec)
by_py_u = lambda vec: Pby*f._b_pyDeriv_u(src,vec)
bz_py_u = lambda vec: Pbz*f._b_pyDeriv_u(src,vec)
# Update the input vector
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
# Define the components of the derivative
Hd = sDiag(1./(sDiag(bx_px)*by_py - sDiag(bx_py)*by_px))
Hd_uV = sDiag(by_py)*bx_px_u(v) + sDiag(bx_px)*by_py_u(v) - sDiag(bx_py)*by_px_u(v) - sDiag(by_px)*bx_py_u(v)
if 'tzx' in self.rxType:
Tij = sDiag(Hd*( - sDiag(by_px)*bz_py + sDiag(by_py)*bz_px ))
TijN_uV = -sDiag(by_px)*bz_py_u(v) - sDiag(bz_py)*by_px_u(v) + sDiag(by_py)*bz_px_u(v) + sDiag(bz_px)*by_py_u(v)
elif 'tzy' in self.rxType:
Tij = sDiag(Hd*( sDiag(bx_px)*bz_py - sDiag(bx_py)*bz_px ))
TijN_uV = sDiag(bz_py)*bx_px_u(v) + sDiag(bx_px)*bz_py_u(v) - sDiag(bx_py)*bz_px_u(v) - sDiag(bz_px)*bx_py_u(v)
# Calculate the complex derivative
PDeriv_complex = Hd * (TijN_uV - Tij * Hd_uV )
# Extract the real number for the real/imag components.
Pv = np.array(getattr(PDeriv_complex, real_or_imag))
elif adjoint:
# Note: The v vector is real and the return should be complex
if self.projType is 'Z1D':
Pex = mesh.getInterpolationMat(self.locs[:,-1],'Fx')
Pbx = mesh.getInterpolationMat(self.locs[:,-1],'Ex')
# ex = Pex*mkvc(f[src,'e_1d'],2)
# bx = Pbx*mkvc(f[src,'b_1d'],2)/mu_0
dP_deTv = -mkvc(Pex.T*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0)).T*v,2)
db_duv = Pbx.T/mu_0*Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))*(Utils.sdiag(1./(Pbx*mkvc(f[src,'b_1d'],2)/mu_0))).T*Utils.sdiag(Pex*mkvc(f[src,'e_1d'],2)).T*v
dP_dbTv = mkvc(f._bDeriv_u(src,db_duv,adjoint=True),2)
PDeriv_real = np.sum(np.hstack((dP_deTv,dP_dbTv)),1)
elif self.projType is 'Z2D':
raise NotImplementedError('Has not be implement for 2D impedance tensor')
elif self.projType is 'Z3D':
if self.locs.ndim == 3:
eFLocs = self.locs[:,:,0]
bFLocs = self.locs[:,:,1]
else:
eFLocs = self.locs
bFLocs = self.locs
# Get the projection
Pex = mesh.getInterpolationMat(eFLocs,'Ex')
Pey = mesh.getInterpolationMat(eFLocs,'Ey')
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
# Get the fields at location
# px: x-polaration and py: y-polaration.
aex_px = mkvc(mkvc(f[src,'e_px'],2).T*Pex.T)
aey_px = mkvc(mkvc(f[src,'e_px'],2).T*Pey.T)
aex_py = mkvc(mkvc(f[src,'e_py'],2).T*Pex.T)
aey_py = mkvc(mkvc(f[src,'e_py'],2).T*Pey.T)
ahx_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pbx.T)
ahy_px = mkvc(mkvc(f[src,'b_px'],2).T/mu_0*Pby.T)
ahx_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pbx.T)
ahy_py = mkvc(mkvc(f[src,'b_py'],2).T/mu_0*Pby.T)
# Derivatives as lambda functions
aex_px_u = lambda vec: f._e_pxDeriv_u(src,Pex.T*vec,adjoint=True)
aey_px_u = lambda vec: f._e_pxDeriv_u(src,Pey.T*vec,adjoint=True)
aex_py_u = lambda vec: f._e_pyDeriv_u(src,Pex.T*vec,adjoint=True)
aey_py_u = lambda vec: f._e_pyDeriv_u(src,Pey.T*vec,adjoint=True)
ahx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
ahx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)/mu_0
ahy_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)/mu_0
# Update the input vector
# Define shortcuts
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(1./(sDiag(ahx_px)*ahy_py - sDiag(ahx_py)*ahy_px))
aHd_uV = lambda x: ahx_px_u(sDiag(ahy_py)*x) + ahx_px_u(sDiag(ahy_py)*x) - ahy_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(ahy_px)*x)
# Need to fix this to reflect the adjoint
if 'zxx' in self.rxType:
Zij = sDiag(aHd*( sDiag(ahy_py)*aex_px - sDiag(ahy_px)*aex_py))
ZijN_uV = lambda x: aex_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aex_px)*x) - ahy_px_u(sDiag(aex_py)*x) - aex_py_u(sDiag(ahy_px)*x)
elif 'zxy' in self.rxType:
Zij = sDiag(aHd*(-sDiag(ahx_py)*aex_px + sDiag(ahx_px)*aex_py))
ZijN_uV = lambda x:-aex_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aex_px)*x) + ahx_px_u(sDiag(aex_py)*x) + aex_py_u(sDiag(ahx_px)*x)
elif 'zyx' in self.rxType:
Zij = sDiag(aHd*( sDiag(ahy_py)*aey_px - sDiag(ahy_px)*aey_py))
ZijN_uV = lambda x: aey_px_u(sDiag(ahy_py)*x) + ahy_py_u(sDiag(aey_px)*x) - ahy_px_u(sDiag(aey_py)*x) - aey_py_u(sDiag(ahy_px)*x)
elif 'zyy' in self.rxType:
Zij = sDiag(aHd*(-sDiag(ahx_py)*aey_px + sDiag(ahx_px)*aey_py))
ZijN_uV = lambda x:-aey_px_u(sDiag(ahx_py)*x) - ahx_py_u(sDiag(aey_px)*x) + ahx_px_u(sDiag(aey_py)*x) + aey_py_u(sDiag(ahx_px)*x)
# Calculate the complex derivative
PDeriv_real = ZijN_uV(aHd*v) - aHd_uV(Zij.T*aHd*v)#
# NOTE: Need to reshape the output to go from 2*nU array to a (nU,2) matrix for each polarization
# PDeriv_real = np.hstack((mkvc(PDeriv_real[:len(PDeriv_real)/2],2),mkvc(PDeriv_real[len(PDeriv_real)/2::],2)))
PDeriv_real = PDeriv_real.reshape((2,mesh.nE)).T
elif self.projType is 'T3D':
if self.locs.ndim == 3:
bFLocs = self.locs[:,:,1]
else:
bFLocs = self.locs
# Get the projection
Pbx = mesh.getInterpolationMat(bFLocs,'Fx')
Pby = mesh.getInterpolationMat(bFLocs,'Fy')
Pbz = mesh.getInterpolationMat(bFLocs,'Fz')
# Get the fields at location
# px: x-polaration and py: y-polaration.
abx_px = mkvc(mkvc(f[src,'b_px'],2).T*Pbx.T)
aby_px = mkvc(mkvc(f[src,'b_px'],2).T*Pby.T)
abz_px = mkvc(mkvc(f[src,'b_px'],2).T*Pbz.T)
abx_py = mkvc(mkvc(f[src,'b_py'],2).T*Pbx.T)
aby_py = mkvc(mkvc(f[src,'b_py'],2).T*Pby.T)
abz_py = mkvc(mkvc(f[src,'b_py'],2).T*Pbz.T)
# Derivatives as lambda functions
abx_px_u = lambda vec: f._b_pxDeriv_u(src,Pbx.T*vec,adjoint=True)
aby_px_u = lambda vec: f._b_pxDeriv_u(src,Pby.T*vec,adjoint=True)
abz_px_u = lambda vec: f._b_pxDeriv_u(src,Pbz.T*vec,adjoint=True)
abx_py_u = lambda vec: f._b_pyDeriv_u(src,Pbx.T*vec,adjoint=True)
aby_py_u = lambda vec: f._b_pyDeriv_u(src,Pby.T*vec,adjoint=True)
abz_py_u = lambda vec: f._b_pyDeriv_u(src,Pbz.T*vec,adjoint=True)
# Update the input vector
# Define shortcuts
sDiag = lambda t: Utils.sdiag(mkvc(t,2))
sVec = lambda t: Utils.sp.csr_matrix(mkvc(t,2))
# Define the components of the derivative
aHd = sDiag(1./(sDiag(abx_px)*aby_py - sDiag(abx_py)*aby_px))
aHd_uV = lambda x: abx_px_u(sDiag(aby_py)*x) + abx_px_u(sDiag(aby_py)*x) - aby_px_u(sDiag(abx_py)*x) - abx_py_u(sDiag(aby_px)*x)
# Need to fix this to reflect the adjoint
if 'tzx' in self.rxType:
Tij = sDiag(aHd*( -sDiag(abz_py)*aby_px + sDiag(abz_px)*aby_py))
TijN_uV = lambda x: -abz_py_u(sDiag(aby_px)*x) - aby_px_u(sDiag(abz_py)*x) + aby_py_u(sDiag(abz_px)*x) + abz_px_u(sDiag(aby_py)*x)
elif 'tzy' in self.rxType:
Tij = sDiag(aHd*( sDiag(abz_py)*abx_px - sDiag(abz_px)*abx_py))
TijN_uV = lambda x: abx_px_u(sDiag(abz_py)*x) + abz_py_u(sDiag(abx_px)*x) - abx_py_u(sDiag(abz_px)*x) - abz_px_u(sDiag(abx_py)*x)
# Calculate the complex derivative
PDeriv_real = TijN_uV(aHd*v) - aHd_uV(Tij.T*aHd*v)#
# NOTE: Need to reshape the output to go from 2*nU array to a (nU,2) matrix for each polarization
# PDeriv_real = np.hstack((mkvc(PDeriv_real[:len(PDeriv_real)/2],2),mkvc(PDeriv_real[len(PDeriv_real)/2::],2)))
PDeriv_real = PDeriv_real.reshape((2,mesh.nE)).T
# Extract the data
if real_or_imag == 'imag':
Pv = 1j*PDeriv_real
elif real_or_imag == 'real':
Pv = PDeriv_real.astype(complex)
return Pv
#################
### Survey ###
#################
class Survey(SimPEGsurvey.BaseSurvey):
"""
Survey class for MT. Contains all the sources associated with the survey.
:param list srcList: List of sources associated with the survey
"""
srcPair = SrcMT.BaseMTSrc
def __init__(self, srcList, **kwargs):
# Sort these by frequency
self.srcList = srcList
SimPEGsurvey.BaseSurvey.__init__(self, **kwargs)
_freqDict = {}
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])
@property
def freqs(self):
"""Frequencies"""
return self._freqs
@property
def nFreq(self):
"""Number of frequencies"""
return len(self._freqDict)
# TODO: Rename to getSources
def getSrcByFreq(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 eval(self, u):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def evalDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
#################
### Data ###
#################
class Data(SimPEGsurvey.Data):
'''
Data class for MTdata. Stores the data vector indexed by the survey.
:param SimPEG survey object survey:
:param v vector of the data in order matching of the survey
'''
def __init__(self, survey, v=None):
# Pass the variables to the "parent" method
SimPEGsurvey.Data.__init__(self, survey, v)
# # Import data
# @classmethod
# def fromEDIFiles():
# pass
def toRecArray(self,returnType='RealImag'):
'''
Function that returns a numpy.recarray for a SimpegMT impedance data object.
:param str returnType: Switches between returning a rec array where the impedance is split to real and imaginary ('RealImag') or is a complex ('Complex')
'''
# Define the record fields
dtRI = [('freq',float),('x',float),('y',float),('z',float),('zxxr',float),('zxxi',float),('zxyr',float),('zxyi',float),
('zyxr',float),('zyxi',float),('zyyr',float),('zyyi',float),('tzxr',float),('tzxi',float),('tzyr',float),('tzyi',float)]
dtCP = [('freq',float),('x',float),('y',float),('z',float),('zxx',complex),('zxy',complex),('zyx',complex),('zyy',complex),('tzx',complex),('tzy',complex)]
impList = ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']
for src in self.survey.srcList:
# Temp array for all the receivers of the source.
# Note: needs to be written more generally, using diffterent rxTypes and not all the data at the locaitons
# Assume the same locs for all RX
locs = src.rxList[0].locs
if locs.shape[1] == 1:
locs = np.hstack((np.array([[0.0,0.0]]),locs))
elif locs.shape[1] == 2:
locs = np.hstack((np.array([[0.0]]),locs))
tArrRec = np.concatenate((src.freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],12))),axis=1).view(dtRI)
# np.array([(src.freq,rx.locs[0,0],rx.locs[0,1],rx.locs[0,2],np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ,np.nan ) for rx in src.rxList],dtype=dtRI)
# Get the type and the value for the DataMT object as a list
typeList = [[rx.rxType.replace('z1d','zyx'),self[src,rx]] for rx in src.rxList]
# Insert the values to the temp array
for nr,(key,val) in enumerate(typeList):
tArrRec[key] = mkvc(val,2)
# Masked array
mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype)
# Unique freq and loc of the masked array
uniFLmarr = np.unique(mArrRec[['freq','x','y','z']]).copy()
try:
outTemp = recFunc.stack_arrays((outTemp,mArrRec))
#outTemp = np.concatenate((outTemp,dataBlock),axis=0)
except NameError as e:
outTemp = mArrRec
if 'RealImag' in returnType:
outArr = outTemp
elif 'Complex' in returnType:
# Add the real and imaginary to a complex number
outArr = np.empty(outTemp.shape,dtype=dtCP)
for comp in ['freq','x','y','z']:
outArr[comp] = outTemp[comp].copy()
for comp in ['zxx','zxy','zyx','zyy','tzx','tzy']:
outArr[comp] = outTemp[comp+'r'].copy() + 1j*outTemp[comp+'i'].copy()
else:
raise NotImplementedError('{:s} is not implemented, as to be RealImag or Complex.')
# Return
return outArr
@classmethod
def fromRecArray(cls, recArray, srcType='primary'):
"""
Class method that reads in a numpy record array to MTdata object.
Only imports the impedance data.
"""
if srcType=='primary':
src = SrcMT.polxy_1Dprimary
elif srcType=='total':
src = SrcMT.polxy_1DhomotD
else:
raise NotImplementedError('{:s} is not a valid source type for MTdata')
# Find all the frequencies in recArray
uniFreq = np.unique(recArray['freq'])
srcList = []
dataList = []
for freq in uniFreq:
# Initiate rxList
rxList = []
# Find that data for freq
dFreq = recArray[recArray['freq'] == freq].copy()
# Find the impedance rxTypes in the recArray.
rxTypes = [ comp for comp in recArray.dtype.names if (len(comp)==4 or len(comp)==3) and 'z' in comp]
for rxType in rxTypes:
# Find index of not nan values in rxType
notNaNind = ~np.isnan(dFreq[rxType])
if np.any(notNaNind): # Make sure that there is any data to add.
locs = rec2ndarr(dFreq[['x','y','z']][notNaNind].copy())
if dFreq[rxType].dtype.name in 'complex128':
rxList.append(Rx(locs,rxType+'r'))
dataList.append(dFreq[rxType][notNaNind].real.copy())
rxList.append(Rx(locs,rxType+'i'))
dataList.append(dFreq[rxType][notNaNind].imag.copy())
else:
rxList.append(Rx(locs,rxType))
dataList.append(dFreq[rxType][notNaNind].copy())
srcList.append(src(rxList,freq))
# Make a survey
survey = Survey(srcList)
dataVec = np.hstack(dataList)
return cls(survey,dataVec)
+108
View File
@@ -0,0 +1,108 @@
# Analytic solution of EM fields due to a plane wave
import numpy as np, SimPEG as simpeg
from scipy.constants import mu_0, epsilon_0 as eps_0
def getEHfields(m1d,sigma,freq,zd,scaleUD=True):
'''Analytic solution for MT 1D layered earth. Returns E and H fields.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property of conductivity corresponding with the mesh.
:param float, freq: Frequency to calculate data at.
:param numpy array, vector zd: location to calculate EH fields at
:param bollean, scaleUD: scales the output to be 1 at the top, increases numeracal stability.
Assumes a halfspace with the same conductive as the last cell below.
'''
# Note add an error check for the mesh and sigma are the same size.
# Constants: Assume constant
mu = mu_0*np.ones((m1d.nC+1))
eps = eps_0*np.ones((m1d.nC+1))
# Angular freq
w = 2*np.pi*freq
# Add the halfspace value to the property
sig = np.concatenate((np.array([sigma[0]]),sigma))
# Calculate the wave number
k = np.sqrt(eps*mu*w**2-1j*mu*sig*w)
# Initiate the propagation matrix, in the order down up.
UDp = np.zeros((2,m1d.nC+1),dtype=complex)
UDp[1,0] = 1. # Set the wave amplitude as 1 into the half-space at the bottom of the mesh
# Loop over all the layers, starting at the bottom layer
for lnr, h in enumerate(m1d.hx): # lnr-number of layer, h-thickness of the layer
# Calculate
yp1 = k[lnr]/(w*mu[lnr]) # Admittance of the layer below the current layer
zp = (w*mu[lnr+1])/k[lnr+1] # Impedance in the current layer
# Build the propagation matrix
# Convert fields to down/up going components in layer below current layer
Pj1 = np.array([[1,1],[yp1,-yp1]])
# Convert fields to down/up going components in current layer
Pjinv = 1./2*np.array([[1,zp],[1,-zp]])
# Propagate down and up components through the current layer
elamh = np.array([[np.exp(-1j*k[lnr+1]*h),0],[0,np.exp(1j*k[lnr+1]*h)]])
# The down and up component in current layer.
UDp[:,lnr+1] = elamh.dot(Pjinv.dot(Pj1)).dot(UDp[:,lnr])
if scaleUD:
UDp[:,lnr+1::-1] = UDp[:,lnr+1::-1]/UDp[1,lnr+1]
# Calculate the fields
Ed = np.empty((zd.size,),dtype=complex)
Eu = np.empty((zd.size,),dtype=complex)
Hd = np.empty((zd.size,),dtype=complex)
Hu = np.empty((zd.size,),dtype=complex)
# Loop over the layers and calculate the fields
# In the halfspace below the mesh
dup = m1d.vectorNx[0]
dind = dup >= zd
Ed[dind] = UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Eu[dind] = UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
Hd[dind] = (k[0]/(w*mu[0]))*UDp[1,0]*np.exp(-1j*k[0]*(dup-zd[dind]))
Hu[dind] = -(k[0]/(w*mu[0]))*UDp[0,0]*np.exp(1j*k[0]*(dup-zd[dind]))
for ki,mui,epsi,dlow,dup,Up,Dp in zip(k[1::],mu[1::],eps[1::],m1d.vectorNx[:-1],m1d.vectorNx[1::],UDp[0,1::],UDp[1,1::]):
dind = np.logical_and(dup >= zd, zd > dlow)
Ed[dind] = Dp*np.exp(-1j*ki*(dup-zd[dind]))
Eu[dind] = Up*np.exp(1j*ki*(dup-zd[dind]))
Hd[dind] = (ki/(w*mui))*Dp*np.exp(-1j*ki*(dup-zd[dind]))
Hu[dind] = -(ki/(w*mui))*Up*np.exp(1j*ki*(dup-zd[dind]))
# Return return the fields
return Ed, Eu, Hd, Hu
def getImpedance(m1d,sigma,freq):
"""Analytic solution for MT 1D layered earth. Returns the impedance at the surface.
:param SimPEG.mesh, object m1d: Mesh object with the 1D spatial information.
:param numpy.array, vector sigma: Physical property corresponding with the mesh.
:param numpy.array, vector freq: Frequencies to calculate data at.
"""
# Initiate the impedances
Z1d = np.empty(len(freq) , dtype='complex')
h = m1d.hx #vectorNx[:-1]
# Start the process
for nrFr, fr in enumerate(freq):
om = 2*np.pi*fr
Zall = np.empty(len(h)+1,dtype='complex')
# Calculate the impedance for the bottom layer
Zall[0] = (mu_0*om)/np.sqrt(mu_0*eps_0*(om)**2 - 1j*mu_0*sigma[0]*om)
for nr,hi in enumerate(h):
# Calculate the wave number
# print nr,sigma[nr]
k = np.sqrt(mu_0*eps_0*om**2 - 1j*mu_0*sigma[nr]*om)
Z = (mu_0*om)/k
Zall[nr+1] = Z *((Zall[nr] + Z*np.tanh(1j*k*hi))/(Z + Zall[nr]*np.tanh(1j*k*hi)))
#pdb.set_trace()
Z1d[nrFr] = Zall[-1]
return Z1d
+45
View File
@@ -0,0 +1,45 @@
import numpy as np, SimPEG as simpeg
from MT1Danalytic import getEHfields
from scipy.constants import mu_0
def get1DEfields(m1d,sigma,freq,sourceAmp=1.0):
"""Function to get 1D electrical fields"""
# Get the gradient
G = m1d.nodalGrad
# Mass matrices
# Magnetic permeability
Mmu = simpeg.Utils.sdiag(m1d.vol*(1.0/mu_0))
# Conductivity
Msig = m1d.getFaceInnerProduct(sigma)
# Set up the solution matrix
A = G.T*Mmu*G + 1j*2.*np.pi*freq*Msig
# Define the inner part of the solution matrix
Aii = A[1:-1,1:-1]
# Define the outer part of the solution matrix
Aio = A[1:-1,[0,-1]]
# Set the boundary conditions
Ed, Eu, Hd, Hu = getEHfields(m1d,sigma,freq,m1d.vectorNx)
Etot = (Ed + Eu)
if sourceAmp is not None:
Etot = ((Etot/Etot[-1])*sourceAmp) # Scale the fields to be equal to sourceAmp at the top
## Note: The analytic solution is derived with e^iwt
bc = np.r_[Etot[0],Etot[-1]]
# The right hand side
rhs = Aio*bc
# Solve the system
Aii_inv = simpeg.Solver(Aii)
eii = Aii_inv*rhs
# Assign the boundary conditions
e = np.r_[bc[0],eii,bc[1]]
# Return the electrical fields
return e
if __name__ == '__main__':
hz = [(100.,18)]
M = simpeg.Mesh.TensorMesh([hz],'C')
sig = np.zeros(M.nC) + 1e-8
sig[M.vectorCCx<=0] = sigHalf
+4
View File
@@ -0,0 +1,4 @@
from MT1Dsolutions import * # Add the names of the functions
from MT1Danalytic import *
from dataUtils import *
from ediFilesUtils import *
+245
View File
@@ -0,0 +1,245 @@
# Utils used for the data,
import numpy as np, matplotlib.pyplot as plt, sys
import SimPEG as simpeg
import numpy.lib.recfunctions as recFunc
from scipy.constants import mu_0
from scipy import interpolate as sciint
def getAppRes(MTdata):
# Make impedance
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def rotateData(MTdata,rotAngle):
'''
Function that rotates clockwist by rotAngle (- negative for a counter-clockwise rotation)
'''
recData = MTdata.toRecArray('Complex')
impData = rec2ndarr(recData[['zxx','zxy','zyx','zyy']],complex)
# Make the rotation matrix
# c,s,zxx,zxy,zyx,zyy = sympy.symbols('c,s,zxx,zxy,zyx,zyy')
# rotM = sympy.Matrix([[c,-s],[s, c]])
# zM = sympy.Matrix([[zxx,zxy],[zyx,zyy]])
# rotM*zM*rotM.T
# [c*(c*zxx - s*zyx) - s*(c*zxy - s*zyy), c*(c*zxy - s*zyy) + s*(c*zxx - s*zyx)],
# [c*(c*zyx + s*zxx) - s*(c*zyy + s*zxy), c*(c*zyy + s*zxy) + s*(c*zyx + s*zxx)]])
s = np.sin(-np.deg2rad(rotAngle))
c = np.cos(-np.deg2rad(rotAngle))
rotMat = np.array([[c,-s],[s,c]])
rotData = (rotMat.dot(impData.reshape(-1,2,2).dot(rotMat.T))).transpose(1,0,2).reshape(-1,4)
outRec = recData.copy()
for nr,comp in enumerate(['zxx','zxy','zyx','zyy']):
outRec[comp] = rotData[:,nr]
from SimPEG import MT
return MT.Data.fromRecArray(outRec)
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def skindepth(rho,freq):
''' Function to calculate the skindepth of EM waves'''
return np.sqrt( (rho*((1/(freq * mu_0 * np.pi )))))
def rec2ndarr(x,dt=float):
return x.view((dt, len(x.dtype.names)))
def makeAnalyticSolution(mesh,model,elev,freqs):
from SimPEG import MT
data1D = []
for freq in freqs:
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
anaZ = anaE/anaH
# Add to the list
data1D.append((freq,0,0,elev,anaZ[0]))
dataRec = np.array(data1D,dtype=[('freq',float),('x',float),('y',float),('z',float),('zyx',complex)])
return dataRec
def plotMT1DModelData(problem,models,symList=None):
from SimPEG import MT
# Setup the figure
fontSize = 15
fig = plt.figure(figsize=[9,7])
axM = fig.add_axes([0.075,.1,.25,.875])
axM.set_xlabel('Resistivity [Ohm*m]',fontsize=fontSize)
axM.set_xlim(1e-1,1e5)
axM.set_ylim(-10000,5000)
axM.set_ylabel('Depth [km]',fontsize=fontSize)
axR = fig.add_axes([0.42,.575,.5,.4])
axR.set_xscale('log')
axR.set_yscale('log')
axR.invert_xaxis()
# axR.set_xlabel('Frequency [Hz]')
axR.set_ylabel('Apparent resistivity [Ohm m]',fontsize=fontSize)
axP = fig.add_axes([0.42,.1,.5,.4])
axP.set_xscale('log')
axP.invert_xaxis()
axP.set_ylim(0,90)
axP.set_xlabel('Frequency [Hz]',fontsize=fontSize)
axP.set_ylabel('Apparent phase [deg]',fontsize=fontSize)
# if not symList:
# symList = ['x']*len(models)
import plotDataTypes as pDt
# Loop through the models.
modelList = [problem.survey.mtrue]
modelList.extend(models)
if False:
modelList = [problem.mapping.sigmaMap*mod for mod in modelList]
for nr, model in enumerate(modelList):
# Calculate the data
if nr==0:
data1D = problem.dataPair(problem.survey,problem.survey.dobs).toRecArray('Complex')
else:
data1D = problem.dataPair(problem.survey,problem.survey.dpred(model)).toRecArray('Complex')
# Plot the data and the model
colRat = nr/((len(modelList)-1.999)*1.)
if colRat > 1.:
col = 'k'
else:
col = plt.cm.seismic(1-colRat)
# The model - make the pts to plot
meshPts = np.concatenate((problem.mesh.gridN[0:1],np.kron(problem.mesh.gridN[1::],np.ones(2))[:-1]))
modelPts = np.kron(1./(problem.mapping.sigmaMap*model),np.ones(2,))
axM.semilogx(modelPts,meshPts,color=col)
## Data
# Appres
pDt.plotIsoStaImpedance(axR,np.array([0,0]),data1D,'zyx','res',pColor=col)
# Appphs
pDt.plotIsoStaImpedance(axP,np.array([0,0]),data1D,'zyx','phs',pColor=col)
try:
allData = np.concatenate((allData,simpeg.mkvc(data1D['zyx'],2)),1)
except:
allData = simpeg.mkvc(data1D['zyx'],2)
freq = simpeg.mkvc(data1D['freq'],2)
res, phs = appResPhs(freq,allData)
stdCol = 'gray'
axRtw = axR.twinx()
axRtw.set_ylabel('Std of log10',color=stdCol)
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
axPtw = axP.twinx()
axPtw.set_ylabel('Std ',color=stdCol)
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
axRtw.plot(freq, np.std(np.log10(res),1),'--',color=stdCol)
axPtw.plot(freq, np.std(phs,1),'--',color=stdCol)
# Fix labels and ticks
yMtick = [l/1000 for l in axM.get_yticks().tolist()]
axM.set_yticklabels(yMtick)
[ l.set_rotation(90) for l in axM.get_yticklabels()]
[ l.set_rotation(90) for l in axR.get_yticklabels()]
[(t.set_color(stdCol), t.set_rotation(-45)) for t in axRtw.get_yticklabels()]
[t.set_color(stdCol) for t in axPtw.get_yticklabels()]
for ax in [axM,axR,axP]:
ax.xaxis.set_tick_params(labelsize=fontSize)
ax.yaxis.set_tick_params(labelsize=fontSize)
return fig
def printTime():
import time
print time.strftime("%a, %d %b %Y %H:%M:%S +0000", time.localtime())
def convert3Dto1Dobject(MTdata,rxType3D='zyx'):
from SimPEG import MT
# Find the unique locations
# Need to find the locations
recDataTemp = MTdata.toRecArray()
# Check if survey.std has been assigned.
## NEED TO: write this...
# Calculte and add the DET of the tensor to the recArray
if 'det' in rxType3D:
Zon = (recDataTemp['zxxr']+1j*recDataTemp['zxxi'])*(recDataTemp['zyyr']+1j*recDataTemp['zyyi'])
Zoff = (recDataTemp['zxyr']+1j*recDataTemp['zxyi'])*(recDataTemp['zyxr']+1j*recDataTemp['zyxi'])
det = np.sqrt(Zon.data - Zoff.data)
recData = recFunc.append_fields(recDataTemp,['zdetr','zdeti'],[det.real,det.imag] )
else:
recData = recDataTemp
uniLocs = rec2ndarr(np.unique(recData[['x','y','z']])).data
mtData1DList = []
if 'zxy' in rxType3D:
corr = -1 # Shift the data to comply with the quadtrature of the 1d problem
else:
corr = 1
for loc in uniLocs:
# Make the receiver list
rx1DList = []
for rxType in ['z1dr','z1di']:
rx1DList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
locrecData = recData[np.sqrt(np.sum( (rec2ndarr(recData[['x','y','z']]).data - loc )**2,axis=1)) < 1e-5]
dat1DList = []
src1DList = []
for freq in locrecData['freq']:
src1DList.append(MT.SrcMT.src_polxy_1Dprimary(rx1DList,freq))
for comp in ['r','i']:
dat1DList.append( corr * locrecData[rxType3D+comp][locrecData['freq']== freq].data )
# Make the survey
sur1D = MT.Survey(src1DList)
# Make the data
dataVec = np.hstack(dat1DList)
dat1D = MT.Data(sur1D,dataVec)
sur1D.dobs = dataVec
# Need to take MTdata.survey.std and split it as well.
std=0.05
sur1D.std = np.abs(sur1D.dobs*std) #+ 0.01*np.linalg.norm(sur1D.dobs)
mtData1DList.append(dat1D)
# Return the the list of data.
return mtData1DList
def resampleMTdataAtFreq(MTdata,freqs):
"""
Function to resample MTdata at set of frequencies
"""
from SimPEG import MT
# Make a rec array
MTrec = MTdata.toRecArray().data
# Find unique locations
uniLoc = np.unique(MTrec[['x','y','z']])
uniFreq = MTdata.survey.freqs
# Get the comps
dNames = MTrec.dtype
# Loop over all the locations and interpolate
for loc in uniLoc:
# Find the index of the station
ind = np.sqrt(np.sum((rec2ndarr(MTrec[['x','y','z']]) - rec2ndarr(loc))**2,axis=1)) < 1. # Find dist of 1 m accuracy
# Make a temporary recArray and interpolate all the components
tArrRec = np.concatenate((simpeg.mkvc(freqs,2),np.ones((len(freqs),1))*rec2ndarr(loc),np.nan*np.ones((len(freqs),12))),axis=1).view(dNames)
for comp in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
int1d = sciint.interp1d(MTrec[ind]['freq'],MTrec[ind][comp],bounds_error=False)
tArrRec[comp] = simpeg.mkvc(int1d(freqs),2)
# Join together
try:
outRecArr = recFunc.stack_arrays((outRecArr,tArrRec))
except NameError as e:
outRecArr = tArrRec
# Make the MTdata and return
return MT.Data.fromRecArray(outRecArr)
+175
View File
@@ -0,0 +1,175 @@
# Functions to import and export MT EDI files.
from SimPEG import mkvc
from scipy.constants import mu_0
from numpy.lib import recfunctions as recFunc
from SimPEG.MT.Utils.dataUtils import rec2ndarr
# Import modules
import numpy as np
import os, sys, re
try:
import osr
except ImportError as e:
print 'Could not import osr, missing the gdal package'
pass
class EDIimporter:
"""
A class to import EDIfiles.
"""
_impUnitEDI2SI = 4*np.pi*1e-4 # Convert Z[mV/km/nT] (as in EDI)to Z[V/A] SI unit
_impUnitSI2EDI = 1./_impUnitEDI2SI # ConvertZ[V/A] SI unit to Z[mV/km/nT] (as in EDI)
# Properties
filesList = None
comps = None
# Hidden properties
_outEPSG = None
_2out = None
def __init__(self, EDIfilesList, compList=None, outEPSG=None):
# Set the fileList
self.filesList = EDIfilesList
# Set the components to import
if compList is None:
self.comps = ['ZXXR','ZXYR','ZYXR','ZYYR','ZXXI','ZXYI','ZYXI','ZYYI','ZXX.VAR','ZXY.VAR','ZYX.VAR','ZYY.VAR']
else:
self.comps = compList
if outEPSG is not None:
self._outEPSG = outEPSG
def __call__(self,comps=None):
if comps is None:
return self._data
return self._data[comps]
def importFiles(self):
"""
Function to import EDI files into a object.
"""
# Constants that are needed for convertion of units
# Temp lists
tmpStaList = []
tmpCompList = ['freq','x','y','z']
tmpCompList.extend(self.comps)
# Make the outarray
dtRI = [(compS.lower().replace('.',''),float) for compS in tmpCompList]
# Loop through all the files
for nrEDI, EDIfile in enumerate(self.filesList):
# Read the file into a list of the lines
with open(EDIfile,'r') as fid:
EDIlines = fid.readlines()
# Find the location
latD, longD, elevM = _findLatLong(EDIlines)
# Transfrom coordinates
transCoord = self._transfromPoints(longD,latD)
# Extract the name of the file (station)
EDIname = EDIfile.split(os.sep)[-1].split('.')[0]
# Arrange the data
staList = [EDIname, EDIfile, transCoord[0], transCoord[1], elevM[0]]
# Add to the station list
tmpStaList.extend(staList)
# Read the frequency data
freq = _findEDIcomp('>FREQ',EDIlines)
# Make the temporary rec array.
tArrRec = ( np.nan*np.ones( (len(freq),len(dtRI)) ) ).view(dtRI) #np.concatenate((freq*np.ones((locs.shape[0],1)),locs,np.nan*np.ones((locs.shape[0],8))),axis=1).view(dtRI)
# Add data to the array
tArrRec['freq'] = mkvc(freq,2)
tArrRec['x'] = mkvc(np.ones((len(freq),1))*transCoord[0],2)
tArrRec['y'] = mkvc(np.ones((len(freq),1))*transCoord[1],2)
tArrRec['z'] = mkvc(np.ones((len(freq),1))*elevM[0],2)
for comp in self.comps:
# Deal with converting units of the impedance tensor
if 'Z' in comp:
unitConvert = self._impUnitEDI2SI
else:
unitConvert = 1
# Rotate the data since EDI x is *north, y *east but Simpeg uses x *east, y *north (* means internal reference frame)
key = [comp.lower().replace('.','').replace(s,t) for s,t in [['xx','yy'],['xy','yx'],['yx','xy'],['yy','xx']] if s in comp.lower()][0]
tArrRec[key] = mkvc(unitConvert*_findEDIcomp('>'+comp,EDIlines),2)
# Make a masked array
mArrRec = np.ma.MaskedArray(rec2ndarr(tArrRec),mask=np.isnan(rec2ndarr(tArrRec))).view(dtype=tArrRec.dtype)
try:
outTemp = recFunc.stack_arrays((outTemp,mArrRec))
except NameError as e:
outTemp = mArrRec
# Assign the data
self._data = outTemp
# % Assign the data to the obj
# nOutData=length(obj.data);
# obj.data(nOutData+1:nOutData+length(TEMP.data),:) = TEMP.data;
def _transfromPoints(self,longD,latD):
# Coordinates convertor
if self._2out is None:
src = osr.SpatialReference()
src.ImportFromEPSG(4326)
out = osr.SpatialReference()
if self._outEPSG is None:
# Find the UTM EPSG number
Nnr = 700 if latD < 0.0 else 600
utmZ = int(1+(longD+180.0)/6.0)
self._outEPSG = 32000 + Nnr + utmZ
out.ImportFromEPSG(self._outEPSG)
self._2out = osr.CoordinateTransformation(src,out)
# Return the transfrom
return self._2out.TransformPoint(longD,latD)
# Hidden functions
def _findLatLong(fileLines):
latDMS = np.array(fileLines[_findLine('LAT=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
longDMS = np.array(fileLines[_findLine('LONG=',fileLines)[0]].split('=')[1].split()[0].split(':'),float)
elevM = np.array([fileLines[_findLine('ELEV=',fileLines)[0]].split('=')[1].split()[0]],float)
# Convert to D.ddddd values
latS = np.sign(latDMS[0])
longS = np.sign(longDMS[0])
latD = latDMS[0] + latS*latDMS[1]/60 + latS*latDMS[2]/3600
longD = longDMS[0] + longS*longDMS[1]/60 + longS*longDMS[2]/3600
return latD, longD, elevM
def _findLine(comp,fileLines):
""" Find a line number in the file"""
# Line counter
c = 0
# List of indices for found lines
found = []
# Loop through all the lines
for line in fileLines:
if comp in line:
# Append if found
found.append(c)
# Increse the counter
c += 1
# Return the found indices
return found
def _findEDIcomp(comp,fileLines,dt=float):
"""
Extract the data vector.
Returns a list of the data.
"""
# Find the data
headLine, indHead = [(st,nr) for nr,st in enumerate(fileLines) if re.search(comp,st)][0]
# Extract the data
nrVec = int(headLine.split()[-1])
c = 0
dataList = []
while c < nrVec:
indHead += 1
dataList.extend(fileLines[indHead].split())
c = len(dataList)
return np.array(dataList,dt)
+416
View File
@@ -0,0 +1,416 @@
from matplotlib import pyplot as plt, colors, numpy as np
def rec2nd(structArray):
""" Converts a structured/record array to ndarray to do operations on."""
return structArray.view((np.float,len(structArray.dtype.names)))
def plotIsoFreqNSimpedance(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
indUniFreq = np.where(freq==array['freq'])
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
if par == 'abs':
zPlot = np.abs(array[flag][indUniFreq])
cmap = plt.get_cmap('OrRd_r')#seismic')
level = np.logspace(0,-5,31)
clevel = np.logspace(0,-4,5)
plotNorm = colors.LogNorm()
elif par == 'real':
zPlot = np.real(array[flag][indUniFreq])
cmap = plt.get_cmap('RdYlBu')
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
else:
plotNorm = colors.Normalize()
elif par == 'imag':
zPlot = np.imag(array[flag][indUniFreq])
cmap = plt.get_cmap('RdYlBu')
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
if cLevel:
level = np.concatenate((-np.logspace(0,-10,31),np.logspace(-10,0,31)))
clevel = np.concatenate((-np.logspace(0,-8,5),np.logspace(-8,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-10,linscale=2)
elif colorNorm=='Lin':
plotNorm = colors.Normalize()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
return cs
def plotIsoFreqNSDiff(ax,freq,arrayList,flag,par='abs',colorbar=True,cLevel=True,mask=None,contourLine=True,useLog=False):
indUniFreq0 = np.where(freq==arrayList[0]['freq'])
indUniFreq1 = np.where(freq==arrayList[1]['freq'])
seicmap = plt.get_cmap('RdYlBu')#seismic')
x, y = arrayList[0]['x'][indUniFreq0],arrayList[0]['y'][indUniFreq0]
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arrayList[0][flag][indUniFreq0])) - np.log10(np.abs(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs(arrayList[1][flag][indUniFreq1]))
else:
zPlot = (np.abs(arrayList[0][flag][indUniFreq0]) - np.abs(arrayList[1][flag][indUniFreq1]))/np.abs(arrayList[1][flag][indUniFreq1])
if mask:
maskInd = np.logical_or(np.abs(arrayList[0][flag][indUniFreq0])< 1e-3,np.abs(arrayList[1][flag][indUniFreq1]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'real':
if useLog:
zPlot = (np.log10(np.real(arrayList[0][flag][indUniFreq0])) -np.log10(np.real(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.real(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = (np.real(arrayList[0][flag][indUniFreq0]) -np.real(arrayList[1][flag][indUniFreq1]))/np.abs((np.real(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.real(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.real(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
elif par == 'imag':
if useLog:
zPlot = (np.log10(np.imag(arrayList[0][flag][indUniFreq0])) -np.log10(np.imag(arrayList[1][flag][indUniFreq1])))/np.log10(np.abs((np.imag(arrayList[1][flag][indUniFreq1]))))
else:
zPlot = (np.imag(arrayList[0][flag][indUniFreq0]) -np.imag(arrayList[1][flag][indUniFreq1]))/np.abs((np.imag(arrayList[1][flag][indUniFreq1])))
if mask:
maskInd = np.logical_or(np.abs(np.imag(arrayList[0][flag][indUniFreq0])) < 1e-3,np.abs(np.imag(arrayList[1][flag][indUniFreq1])) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
if cLevel:
level = np.arange(-200,201,10)
clevel = np.arange(-200,201,25)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
cs = ax.tricontourf(x,y,zPlot*100,levels=level*100,cmap=seicmap,extend='both') #,norm=colors.SymLogNorm(1e-2,linscale=2))
if contourLine:
csl = ax.tricontour(x,y,zPlot*100,levels=clevel*100,colors='k')
plt.clabel(csl, fontsize=7, inline=1,fmt='%1.1e',inline_spacing=10)
if colorbar:
cb = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.1e')
for t in cb.ax.get_yticklabels():
t.set_rotation(60)
t.set_fontsize(8)
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoFreqNStipper(ax,freq,array,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=True,contour=True):
indUniFreq = np.where(freq==array['freq'])
x, y = array['x'][indUniFreq],array['y'][indUniFreq]
if par == 'abs':
cmap = plt.get_cmap('OrRd_r')#seismic')
zPlot = np.abs(array[flag][indUniFreq])
if cLevel:
level = np.logspace(-4,0,33)
clevel = np.logspace(-4,0,5)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.LogNorm()
else:
plotNorm = colors.Normalize()
elif par == 'real':
cmap = plt.get_cmap('RdYlBu')
zPlot = np.real(array[flag][indUniFreq])
if cLevel:
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
else:
plotNorm = colors.Normalize()
elif par == 'imag':
cmap = plt.get_cmap('RdYlBu')
zPlot = np.imag(array[flag][indUniFreq])
if cLevel:
level = np.concatenate((-np.logspace(0,-4,33),np.logspace(-4,0,33)))
clevel = np.concatenate((-np.logspace(0,-4,5),np.logspace(-4,0,5)))
else:
level = np.linspace(zPlot.min(),zPlot.max(),100)
clevel = np.linspace(zPlot.min(),zPlot.max(),10)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(1e-4,linscale=2)
else:
plotNorm = colors.Normalize()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),levels=level,cmap=cmap,norm=plotNorm,edgecolors='k', linewidths=0.5)
if colorbar:
plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
ax.set_title(flag+' '+par,fontsize=8)
def plotIsoStaImpedance(ax,loc,array,flag,par='abs',pSym='s',pColor=None):
appResFact = 1/(8*np.pi**2*10**(-7))
treshold = 1.0 # 1 meter
indUniSta = np.sqrt(np.sum((rec2nd(array[['x','y']])-loc)**2,axis=1)) < treshold
freq = array['freq'][indUniSta]
if par == 'abs':
zPlot = np.abs(array[flag][indUniSta])
elif par == 'real':
zPlot = np.real(array[flag][indUniSta])
elif par == 'imag':
zPlot = np.imag(array[flag][indUniSta])
elif par == 'res':
zPlot = (appResFact/freq)*np.abs(array[flag][indUniSta])**2
elif par == 'phs':
zPlot = np.arctan2(array[flag][indUniSta].imag,array[flag][indUniSta].real)*(180/np.pi)
if not pColor:
if 'xx' in flag:
lab = 'XX'
pColor = 'g'
elif 'xy' in flag:
lab = 'XY'
pColor = 'r'
elif 'yx' in flag:
lab = 'YX'
pColor = 'b'
elif 'yy' in flag:
lab = 'YY'
pColor = 'y'
ax.plot(freq,zPlot,color=pColor,marker=pSym,label=flag)
def plotPsudoSectNSimpedance(ax,sectDict,array,flag,par='abs',colorbar=True,colorNorm='None',cLevel=None,contour=True):
indSect = np.where(sectDict.values()[0]==array[sectDict.keys()[0]])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
x = array['y'][indSect]
else:
x = array['x'][indSect]
y = array['freq'][indSect]
if par == 'abs':
zPlot = np.abs(array[flag][indSect])
cmap = plt.get_cmap('OrRd_r')#seismic')
if cLevel:
level = np.logspace(0,-5,31,endpoint=True)
clevel = np.logspace(0,-4,5,endpoint=True)
else:
level = np.linspace(zPlot.min(),zPlot.max(),100,endpoint=True)
clevel = np.linspace(zPlot.min(),zPlot.max(),10,endpoint=True)
elif par == 'ares':
zPlot = np.abs(array[flag][indSect])**2/(8*np.pi**2*10**(-7)*array['freq'][indSect])
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.logspace(zMin,zMax,(zMax-zMin)*8+1,endpoint=True)
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
plotNorm = colors.LogNorm()
elif par == 'aphs':
zPlot = np.arctan2(array[flag][indSect].imag,array[flag][indSect].real)*(180/np.pi)
cmap = plt.get_cmap('RdYlBu')#seismic)
if cLevel:
zMax = cLevel[1]
zMin = cLevel[0]
else:
zMax = (np.ceil(zPlot).max())
zMin = (np.floor(zPlot).min())
level = np.arange(zMin,zMax+.1,1)
clevel = np.arange(zMin,zMax+.1,10)
plotNorm = colors.Normalize()
elif par == 'real':
zPlot = np.real(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif par == 'imag':
zPlot = np.imag(array[flag][indSect])
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
if colorNorm=='SymLog':
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
plotNorm = colors.Normalize()
elif colorNorm=='Log':
plotNorm = colors.LogNorm()
if contour:
cs = ax.tricontourf(x,y,zPlot,levels=level,cmap=cmap,norm=plotNorm)#,extend='both')
else:
uniX,uniY = np.unique(x),np.unique(y)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel,format='%1.2e')
# csB.on_mappable_changed(cs)
ax.set_title(flag+' '+par,fontsize=8)
return cs, csB
return cs,None
def plotPsudoSectNSDiff(ax,sectDict,arrayList,flag,par='abs',colorbar=True,colorNorm='SymLog',cLevel=None,contour=True,mask=None,useLog=False):
def sortInArr(arr):
return np.sort(arr,order=['freq','x','y','z'])
# Find the index for the slice
indSect0 = np.where(sectDict.values()[0]==arrayList[0][sectDict.keys()[0]])
indSect1 = np.where(sectDict.values()[0]==arrayList[1][sectDict.keys()[0]])
# Extract and sort the mats
arr0 = sortInArr(arrayList[0][indSect0])
arr1 = sortInArr(arrayList[1][indSect1])
# Define the plot axes
if 'x' in sectDict.keys()[0]:
x0 = arr0['y']
x1 = arr1['y']
else:
x0 = arr0['x']
x1 = arr1['x']
y0 = arr0['freq']
y1 = arr1['freq']
if par == 'abs':
if useLog:
zPlot = (np.log10(np.abs(arr0[flag])) - np.log10(np.abs(arr1[flag])))/np.log10(np.abs(arr1[flag]))
else:
zPlot = (np.abs(arr0[flag]) - np.abs(arr1[flag]))/np.abs(arr1[flag])
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('RdYlBu')#seismic)
elif par == 'ares':
arF = 1/(8*np.pi**2*10**(-7))
if useLog:
zPlot = (np.log10((arF/arr0['freq'])*np.abs(arr0[flag])**2) - np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2))/np.log10((arF/arr1['freq'])*np.abs(arr1[flag])**2)
else:
zPlot = ((arF/arr0['freq'])*np.abs(arr0[flag])**2 - (arF/arr1['freq'])*np.abs(arr1[flag])**2)/((arF/arr1['freq'])*np.abs(arr1[flag])**2)
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'aphs':
if useLog:
zPlot = (np.log10(np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi)) - np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi)) )/np.log10(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
else:
zPlot = ( np.arctan2(arr0[flag].imag,arr0[flag].real)*(180/np.pi) - np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi) )/(np.arctan2(arr1[flag].imag,arr1[flag].real)*(180/np.pi))
if mask:
maskInd = np.logical_or(np.abs(arr0[flag])< 1e-3,np.abs(arr1[flag]) < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral')#seismic)
elif par == 'real':
if useLog:
zPlot = (np.log10(arr0[flag].real) - np.log10(arr1[flag].real))/np.log10(arr1[flag].real)
else:
zPlot = (arr0[flag].real - arr1[flag].real)/arr1[flag].real
if mask:
maskInd = np.logical_or(arr0[flag].real< 1e-3,arr1[flag].real < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('Spectral')
elif par == 'imag':
if useLog:
zPlot = (np.log10(arr0[flag].imag) - np.log10(arr1[flag].imag))/np.log10(arr1[flag].imag)
else:
zPlot = (arr0[flag].imag - arr1[flag].imag)/arr1[flag].imag
if mask:
maskInd = np.logical_or(arr0[flag].imag< 1e-3,arr1[flag].imag < 1e-3)
zPlot = np.ma.array(zPlot)
zPlot[maskInd] = mask
cmap = plt.get_cmap('Spectral') #('RdYlBu')
if cLevel:
zMax = np.log10(cLevel[1])
zMin = np.log10(cLevel[0])
else:
zMax = (np.ceil(np.log10(np.abs(zPlot).max())))
zMin = (np.floor(np.log10(np.abs(zPlot).min())))
if colorNorm=='SymLog':
level = np.concatenate((-np.logspace(zMax,zMin-.125,(zMax-zMin)*8+1,endpoint=True),np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)))
clevel = np.concatenate((-np.logspace(zMax,zMin,(zMax-zMin)*1+1,endpoint=True),np.logspace(zMin,zMax,(zMax-zMin)*1+1,endpoint=True)))
plotNorm = colors.SymLogNorm(np.abs(level).min(),linscale=0.1)
elif colorNorm=='Lin':
if cLevel:
level = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/50.)
clevel = np.arange(cLevel[0],cLevel[1]+.1,(cLevel[1] - cLevel[0])/10.)
else:
level = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/50.)
clevel = np.arange(zPlot.min(),zPlot.max(),(zPlot.max() - zPlot.min())/10.)
plotNorm = colors.Normalize()
elif colorNorm=='Log':
level = np.logspace(zMin-.125,zMax,(zMax-zMin)*8+1,endpoint=True)
clevel = np.logspace(zMin,zMax,(zMax-zMin)*2+1,endpoint=True)
plotNorm = colors.LogNorm()
if contour:
cs = ax.tricontourf(x0,y0,zPlot*100,levels=level*100,cmap=cmap,norm=plotNorm,extend='both')#,extend='both')
else:
uniX,uniY = np.unique(x0),np.unique(y0)
X,Y = np.meshgrid(np.append(uniX-25,uniX[-1]+25),np.append(uniY-25,uniY[-1]+25))
cs = ax.pcolor(X,Y,np.reshape(zPlot,(len(uniY),len(uniX))),cmap=cmap,norm=plotNorm)
if colorbar:
csB = plt.colorbar(cs,cax=ax.cax,ticks=clevel*100,format='%1.2e')
# csB.on_mappable_changed(cs)
ax.set_title(flag+' '+par + ' diff',fontsize=8)
return cs, csB
return cs,None
+178
View File
@@ -0,0 +1,178 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,sigma_1d):
'''
Function that calculates and return background fields
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
elif mesh.dim == 2:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
elif mesh.dim == 3:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# # Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
elif mesh.dim == 2:
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
for i in np.arange(mesh.vnEx[0]):
ex_px[i,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
ey_py[i,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
elif mesh.dim == 3:
# Setup x (east) polarization (_x)
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in np.arange(mesh.vnEx[0]):
for j in np.arange(mesh.vnEx[1]):
ex_px[i,j,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
for j in np.arange(mesh.vnEy[1]):
ey_py[i,j,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
def analytic1DModelSource(mesh,freq,sigma_1d):
'''
Function that calculates and return background fields
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import getEHfields
# Get a 1d solution for a halfspace background
if mesh.dim == 1:
mesh1d = mesh
elif mesh.dim == 2:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hy],np.array([mesh.x0[1]]))
elif mesh.dim == 3:
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# # Note: Everything is using e^iwt
Eu, Ed, _, _ = getEHfields(mesh1d,sigma_1d,freq,mesh.vectorNz)
# Make the fields into a dictionary of location and the fields
e0_1d = Eu+Ed
E1dFieldDict = dict(zip(mesh.vectorNz,e0_1d))
if mesh.dim == 1:
eBG_px = simpeg.mkvc(e0_1d,2)
eBG_py = -simpeg.mkvc(e0_1d,2) # added a minus to make the results in the correct quadrents.
elif mesh.dim == 2:
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
for i in np.arange(mesh.vnEx[0]):
ex_px[i,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
ey_py[i,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
elif mesh.dim == 3:
# Setup x (east) polarization (_x)
ex_px = -np.array([E1dFieldDict[i] for i in mesh.gridEx[:,2]]).reshape(-1,1)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Construct the full fields
eBG_px = np.vstack((ex_px,ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.array([E1dFieldDict[i] for i in mesh.gridEy[:,2]]).reshape(-1,1)
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Construct the full fields
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
# def homo3DModelSource(mesh,model,freq):
# '''
# Function that estimates 1D analytic background fields from a 3D model.
# :param Simpeg mesh object mesh: Holds information on the discretization
# :param float freq: The frequency to solve at
# :param np.array sigma_1d: Background model of conductivity to base the calculations on, 1d model.
# :rtype: numpy.ndarray (mesh.nE,2)
# :return: eBG_bp, E fields for the background model at both polarizations.
# '''
# if mesh.dim < 3:
# raise IOError('Input mesh has to have 3 dimensions.')
# # Get the locations
# a = mesh.gridCC[:,0:2].copy()
# unixy = np.unique(a.view(a.dtype.descr * a.shape[1])).view(float).reshape(-1,2)
# uniz = np.unique(mesh.gridCC[:,2])
# # # Note: Everything is using e^iwt
# # Need to loop thourgh the xy locations, assess the model and calculate the fields at the phusdo cell centers.
# # Then interpolate the cc fields to the edges.
# e0_1d = get1DEfields(mesh1d,sigma_1d,freq)
# elif mesh.dim == 3:
# # Setup x (east) polarization (_x)
# ex_px = np.zeros(mesh.vnEx,dtype=complex)
# ey_px = np.zeros((mesh.nEy,1),dtype=complex)
# ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# # Assign the source to ex_x
# for i in np.arange(mesh.vnEx[0]):
# for j in np.arange(mesh.vnEx[1]):
# ex_px[i,j,:] = -e0_1d
# eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# # Setup y (north) polarization (_py)
# ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
# ey_py = np.zeros(mesh.vnEy, dtype='complex128')
# ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# # Assign the source to ey_py
# for i in np.arange(mesh.vnEy[0]):
# for j in np.arange(mesh.vnEy[1]):
# ey_py[i,j,:] = e0_1d
# # ey_py[1:-1,1:-1,1:-1] = 0
# eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# # Return the electric fields
# eBG_bp = np.hstack((eBG_px,eBG_py))
# return eBG_bp
+46
View File
@@ -0,0 +1,46 @@
import SimPEG as simpeg, numpy as np
def homo1DModelSource(mesh,freq,m_back):
'''
Function that calculates and return background fields for a 3D mesh and model.
The calculuations use 1D field solution for a vertical slice throught model (south-western most column),
which is assigned at the fields everywhere for the respective polarizations.2
:param Simpeg mesh object mesh: Holds information on the discretization
:param float freq: The frequency to solve at
:param np.array m_back: Background model of conductivity to base the calculations on.
:rtype: numpy.ndarray (mesh.nE,2)
:return: eBG_bp, E fields for the background model at both polarizations.
'''
# import
from SimPEG.MT.Utils import get1DEfields
# Get a 1d solution for a halfspace background
mesh1d = simpeg.Mesh.TensorMesh([mesh.hz],np.array([mesh.x0[2]]))
# Note: Everything is using e^iwt
e0_1d = get1DEfields(mesh1d,mesh.r(m_back,'CC','CC','M')[0,0,:],freq)
# Setup x (east) polarization (_x)
ex_px = np.zeros(mesh.vnEx,dtype=complex)
ey_px = np.zeros((mesh.nEy,1),dtype=complex)
ez_px = np.zeros((mesh.nEz,1),dtype=complex)
# Assign the source to ex_x
for i in np.arange(mesh.vnEx[0]):
for j in np.arange(mesh.vnEx[1]):
ex_px[i,j,:] = -e0_1d
eBG_px = np.vstack((simpeg.Utils.mkvc(ex_px,2),ey_px,ez_px))
# Setup y (north) polarization (_py)
ex_py = np.zeros((mesh.nEx,1), dtype='complex128')
ey_py = np.zeros(mesh.vnEy, dtype='complex128')
ez_py = np.zeros((mesh.nEz,1), dtype='complex128')
# Assign the source to ey_py
for i in np.arange(mesh.vnEy[0]):
for j in np.arange(mesh.vnEy[1]):
ey_py[i,j,:] = e0_1d
# ey_py[1:-1,1:-1,1:-1] = 0
eBG_py = np.vstack((ex_py,simpeg.Utils.mkvc(ey_py,2),ez_py))
# Return the electric fields
eBG_bp = np.hstack((eBG_px,eBG_py))
return eBG_bp
+5
View File
@@ -0,0 +1,5 @@
import Utils
from SurveyMT import Rx, Survey, Data
from FieldsMT import Fields1D_e, Fields3D_e
import Problem1D, Problem2D, Problem3D
import SrcMT
+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):
"""
-1
View File
@@ -746,4 +746,3 @@ class DiffOperators(object):
kron3(av(n[2]), speye(n[1]+1), av(n[0])),
kron3(speye(n[2]+1), av(n[1]), av(n[0]))), format="csr")
return self._aveN2F
+15
View File
@@ -213,5 +213,20 @@ class BaseTimeProblem(BaseProblem):
if hasattr(self, '_timeMesh'):
del self._timeMesh
class LinearProblem(BaseProblem):
surveyPair = Survey.LinearSurvey
def __init__(self, mesh, G, **kwargs):
BaseProblem.__init__(self, mesh, **kwargs)
self.G = G
def fields(self, m):
return self.G.dot(m)
def Jvec(self, m, v, u=None):
return self.G.dot(v)
def Jtvec(self, m, v, u=None):
return self.G.T.dot(v)
+240
View File
@@ -282,3 +282,243 @@ class Tikhonov(BaseRegularization):
out = mD.T * ( self.W.T * r )
return out
# <<<<<<< HEAD
# class Simple(BaseRegularization):
# """
# Only for tensor mesh
# """
# smoothModel = True #: SMOOTH and SMOOTH_MOD_DIF options
# alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
# alpha_x = Utils.dependentProperty('_alpha_x', 1.0, ['_W', '_Wx'], "Weight for the first derivative in the x direction")
# alpha_y = Utils.dependentProperty('_alpha_y', 1.0, ['_W', '_Wy'], "Weight for the first derivative in the y direction")
# alpha_z = Utils.dependentProperty('_alpha_z', 1.0, ['_W', '_Wz'], "Weight for the first derivative in the z direction")
# alpha_xx = Utils.dependentProperty('_alpha_xx', 0.0, ['_W', '_Wxx'], "Weight for the second derivative in the x direction")
# alpha_yy = Utils.dependentProperty('_alpha_yy', 0.0, ['_W', '_Wyy'], "Weight for the second derivative in the y direction")
# alpha_zz = Utils.dependentProperty('_alpha_zz', 0.0, ['_W', '_Wzz'], "Weight for the second derivative in the z direction")
# def __init__(self, mesh, mapping=None, **kwargs):
# BaseRegularization.__init__(self, mesh, mapping=mapping, **kwargs)
# @property
# def Ws(self):
# """Regularization matrix Ws"""
# if getattr(self,'_Ws', None) is None:
# self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s)**0.5)
# return self._Ws
# @property
# def Wx(self):
# """Regularization matrix Wx"""
# if getattr(self, '_Wx', None) is None:
# self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x)**0.5)*self.mesh.unitCellGradx
# return self._Wx
# @property
# def Wy(self):
# """Regularization matrix Wy"""
# if getattr(self, '_Wy', None) is None:
# self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y)**0.5)*self.mesh.unitCellGrady
# return self._Wy
# @property
# def Wz(self):
# """Regularization matrix Wz"""
# if getattr(self, '_Wz', None) is None:
# self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z)**0.5)*self.mesh.unitCellGradz
# return self._Wz
# @property
# def Wxx(self):
# """Regularization matrix Wxx"""
# if getattr(self, '_Wxx', None) is None:
# self._Wxx = Utils.sdiag((self.mesh.vol*self.alpha_xx)**0.5)*self.mesh.faceDivx*self.mesh.cellGradx
# return self._Wxx
# @property
# def Wyy(self):
# """Regularization matrix Wyy"""
# if getattr(self, '_Wyy', None) is None:
# self._Wyy = Utils.sdiag((self.mesh.vol*self.alpha_yy)**0.5)*self.mesh.faceDivy*self.mesh.cellGrady
# return self._Wyy
# @property
# def Wzz(self):
# """Regularization matrix Wzz"""
# if getattr(self, '_Wzz', None) is None:
# self._Wzz = Utils.sdiag((self.mesh.vol*self.alpha_zz)**0.5)*self.mesh.faceDivz*self.mesh.cellGradz
# return self._Wzz
# @property
# def Wsmooth(self):
# """Full smoothness regularization matrix W"""
# if getattr(self, '_Wsmooth', None) is None:
# wlist = (self.Wx, self.Wxx)
# if self.mesh.dim > 1:
# wlist += (self.Wy, self.Wyy)
# if self.mesh.dim > 2:
# wlist += (self.Wz, self.Wzz)
# self._Wsmooth = sp.vstack(wlist)
# return self._Wsmooth
# @property
# def W(self):
# """Full regularization matrix W"""
# if getattr(self, '_W', None) is None:
# wlist = (self.Ws, self.Wsmooth)
# self._W = sp.vstack(wlist)
# return self._W
# @Utils.timeIt
# def eval(self, m):
# if self.smoothModel == True:
# r1 = self.Wsmooth * ( self.mapping * (m) )
# r2 = self.Ws * ( self.mapping * (m - self.mref) )
# return 0.5*(r1.dot(r1)+r2.dot(r2))
# elif self.smoothModel == False:
# r = self.W * ( self.mapping * (m - self.mref) )
# return 0.5*r.dot(r)
# @Utils.timeIt
# def evalDeriv(self, m):
# """
# The regularization is:
# .. math::
# R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
# So the derivative is straight forward:
# .. math::
# R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
# """
# if self.smoothModel == True:
# mD1 = self.mapping.deriv(m)
# mD2 = self.mapping.deriv(m - self.mref)
# r1 = self.Wsmooth * ( self.mapping * (m))
# r2 = self.Ws * ( self.mapping * (m - self.mref) )
# out1 = mD1.T * ( self.Wsmooth.T * r1 )
# out2 = mD2.T * ( self.Ws.T * r2 )
# out = out1+out2
# elif self.smoothModel == False:
# mD = self.mapping.deriv(m - self.mref)
# r = self.W * ( self.mapping * (m - self.mref) )
# out = mD.T * ( self.W.T * r )
# return out
# class SparseRegularization(Simple):
# eps = 1e-1
# m = None
# gamma = 1.
# p = 0.
# qx = 2.
# qy = 2.
# qz = 2.
# def __init__(self, mesh, mapping=None, **kwargs):
# Simple.__init__(self, mesh, mapping=mapping, **kwargs)
# @property
# def Wsmooth(self):
# """Full smoothness regularization matrix W"""
# if getattr(self, '_Wsmooth', None) is None:
# wlist = (self.Wx, self.Wxx)
# if self.mesh.dim > 1:
# wlist += (self.Wy, self.Wyy)
# if self.mesh.dim > 2:
# wlist += (self.Wz, self.Wzz)
# self._Wsmooth = sp.vstack(wlist)
# return self._Wsmooth
# @property
# def W(self):
# """Full regularization matrix W"""
# if getattr(self, '_W', None) is None:
# wlist = (self.Ws, self.Wsmooth)
# self._W = sp.vstack(wlist)
# return self._W
# @property
# def Ws(self):
# """Regularization matrix Ws"""
# if getattr(self, 'm', None) is None:
# self.Rs = Utils.speye(self.mesh.nC)
# else:
# f_m = self.m
# self.rs = self.R(f_m , self.p, self.eps)
# #print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
# self.Rs = Utils.sdiag( self.rs )
# self._Ws = Utils.sdiag((self.mesh.vol*self.alpha_s*self.gamma)**0.5)*self.Rs
# return self._Ws
# @property
# def Wx(self):
# """Regularization matrix Wx"""
# if getattr(self, 'm', None) is None:
# self.Rx = Utils.speye(self.mesh.unitCellGradx.shape[0])
# else:
# f_m = self.mesh.unitCellGradx * self.m
# self.rx = self.R( f_m , self.qx, self.eps)
# self.Rx = Utils.sdiag( self.rx )
# if getattr(self, '_Wx', None) is None:
# self._Wx = Utils.sdiag((self.mesh.vol*self.alpha_x*self.gamma)**0.5)*self.Rx*self.mesh.unitCellGradx
# return self._Wx
# @property
# def Wy(self):
# """Regularization matrix Wy"""
# if getattr(self, 'm', None) is None:
# self.Ry = Utils.speye(self.mesh.unitCellGrady.shape[0])
# else:
# f_m = self.mesh.unitCellGrady * self.m
# self.ry = self.R( f_m , self.qy, self.eps)
# self.Ry = Utils.sdiag( self.ry )
# if getattr(self, '_Wy', None) is None:
# self._Wy = Utils.sdiag((self.mesh.vol*self.alpha_y*self.gamma)**0.5)*self.Ry*self.mesh.unitCellGrady
# return self._Wy
# @property
# def Wz(self):
# """Regularization matrix Wz"""
# if getattr(self, 'm', None) is None:
# self.Rz = Utils.speye(self.mesh.unitCellGradz.shape[0])
# else:
# f_m = self.mesh.unitCellGradz * self.m
# self.rz = self.R( f_m , self.qz, self.eps)
# self.Rz = Utils.sdiag( self.rz )
# if getattr(self, '_Wz', None) is None:
# self._Wz = Utils.sdiag((self.mesh.vol*self.alpha_z*self.gamma)**0.5)*self.Rz*self.mesh.unitCellGradz
# return self._Wz
# def R(self, f_m , p, dec):
# eta = (self.eps**(1-p/2.))**0.5
# r = eta / (f_m**2.+self.eps**2.)**((1-p/2.)/2.)
# return r
# =======
# >>>>>>> 834de582844e8e1eac95819fbe03eed55dbeb001
+15 -8
View File
@@ -1,6 +1,5 @@
import Utils, numpy as np, scipy.sparse as sp, uuid
class BaseRx(object):
"""SimPEG Receiver Object"""
@@ -307,12 +306,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 +319,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 +331,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):
@@ -375,3 +374,11 @@ class BaseSurvey(object):
self.dobs = self.dtrue+noise
self.std = self.dobs*0 + std
return self.dobs
class LinearSurvey(BaseSurvey):
def eval(self, u):
return u
@property
def nD(self):
return self.prob.G.shape[0]
+38
View File
@@ -118,6 +118,44 @@ def defineElipse(ccMesh, center=[0,0,0], anisotropy=[1,1,1], slope=10., theta=0.
D = np.sqrt(np.sum(G**2,axis=1))
return -np.arctan((D-1)*slope)*(2./np.pi)/2.+0.5
def getIndicesSphere(center,radius,ccMesh):
"""
Creates a vector containing the sphere indices in the cell centers mesh.
Returns a tuple
The sphere is defined by the points
p0, describe the position of the center of the cell
r, describe the radius of the sphere.
ccMesh represents the cell-centered mesh
The points p0 must live in the the same dimensional space as the mesh.
"""
# Validation: mesh and point (p0) live in the same dimensional space
dimMesh = np.size(ccMesh[0,:])
assert len(center) == dimMesh, "Dimension mismatch. len(p0) != dimMesh"
if dimMesh == 1:
# Define the reference points
ind = np.abs(center[0] - ccMesh[:,0]) < radius
elif dimMesh == 2:
# Define the reference points
ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 ) < radius
elif dimMesh == 3:
# Define the points
ind = np.sqrt( ( center[0] - ccMesh[:,0] )**2 + ( center[1] - ccMesh[:,1] )**2 + ( center[2] - ccMesh[:,2] )**2 ) < radius
# Return a tuple
return ind
def defineTwoLayers(ccMesh,depth,vals=[0,1]):
"""
Define a two layered model. Depth of the first layer must be specified.
+150
View File
@@ -0,0 +1,150 @@
.. _api_DC:
.. math::
\renewcommand{\div}{\nabla\cdot\,}
\newcommand{\grad}{\vec \nabla}
\newcommand{\curl}{{\vec \nabla}\times\,}
\newcommand{\dcurl}{{\mathbf C}}
\newcommand{\dgrad}{{\mathbf G}}
\newcommand{\Acf}{{\mathbf A_c^f}}
\newcommand{\Ace}{{\mathbf A_c^e}}
\renewcommand{\S}{{\mathbf \Sigma}}
\renewcommand{\Div}{{\mathbf {Div}}}
\renewcommand{\Grad}{{\mathbf {Grad}}}
\newcommand{\St}{{\mathbf \Sigma_\tau}}
\newcommand{\diag}{\mathbf{diag}}
\newcommand{\M}{{\mathbf M}}
\newcommand{\Me}{{\M^e}}
\newcommand{\Mes}[1]{{\M^e_{#1}}}
\newcommand{\be}{\mathbf{e}}
\newcommand{\bj}{\mathbf{j}}
\newcommand{\bphi}{\mathbf{\phi}}
\newcommand{\bq}{\mathbf{q}}
\newcommand{\bJ}{\mathbf{J}}
\newcommand{\bG}{\mathbf{G}}
\newcommand{\bP}{\mathbf{P}}
\newcommand{\bA}{\mathbf{A}}
\newcommand{\bm}{\mathbf{m}}
\newcommand{\B}{\vec{B}}
\newcommand{\D}{\vec{D}}
\renewcommand{\H}{\vec{H}}
\renewcommand {\j} { {\vec j} }
\newcommand {\h} { {\vec h} }
\renewcommand {\b} { {\vec b} }
\newcommand {\e} { {\vec e} }
\newcommand {\c} { {\vec c} }
\renewcommand {\d} { {\vec d} }
\renewcommand {\u} { {\vec u} }
\newcommand{\I}{\vec{I}}
DC resistivity survey
*********************
Electrical resistivity of subsurface materials is measured by causing an electrical current to flow in the earth between one pair of electrodes while the voltage across a second pair of electrodes is measured. The result is an "apparent" resistivity which is a value representing the weighted average resistivity over a volume of the earth. Variations in this measurement are caused by variations in the soil, rock, and pore fluid electrical resistivity. Surveys require contact with the ground, so they can be labour intensive. Results are sometimes interpreted directly, but more commonly, 1D, 2D or 3D models are estimated using inversion procedures (`GPG <http://www.eos.ubc.ca/courses/eosc350/content/>`_).
Background
==========
As direct current (DC) implies, in DC resistivity survey, we assume steady-state. We consider Maxwell's equations in steady state as
.. math::
\curl \frac{1}{\mu} \vec{b} - \j = \j_s \\
\curl \e = 0
Then by taking \\(\\curl\\) for the first equation, we have
.. math::
- \div\j = q \\
where
.. math::
\div \j_s = q = I(\delta(\vec{r}-\vec{r}_{s+})-\delta(\vec{r}-\vec{r}_{s-}))
Since \\(\\curl \\e = 0\\), we have
.. math::
\e = \grad \phi
And by Ohm's law, we have
.. math::
\j = \sigma \grad \phi
Finally, we can compute the solution of the system:
.. math::
- \div\j = q
\j = \sigma \grad \phi
\frac{\partial \phi}{\partial r}\Big|_{\partial \Omega_{BC}} = 0
Discretization
==============
By using finite volume method (FVM), we discretize our system as
.. math::
-\Div \bj = \bq
\diag(\Acf^{T}\sigma^{-1}) \bj = \Grad \bphi
Here boundary condtions are embedded in the discrete differential operators. With some linear algebra we have
.. math::
\bA\bphi = -\bq
where
.. math::
\bA = \Div (\diag(\Acf^{T}\sigma^{-1}))^{-1} \Grad
By solving this linear equation, we can compute the solution of \\(\\phi\\). Based on this discretization, we derive sensitivity in discretized space. Sensitivity matrix can be in general can be written as
.. math ::
\bJ = -\bP\bA^{-1}\bG
where
.. math ::
\bP: \text{Projection}
\bJ = \bP\frac{\partial \phi}{\partial \bm}
Here \\(\\bm\\) indicates model parameters in discretized space.
Verification
============
Comparing to the analytic function:
.. plot::
import simpegDC as DC
DC.Examples.Verification.run(plotIt=True)
API
===
.. automodule:: simpegDC.BaseDC
:show-inheritance:
:members:
:undoc-members:
:inherited-members:
+21
View File
@@ -0,0 +1,21 @@
.. _examples_DC_Analytic_Dipole:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
DC Analytic Dipole
==================
.. plot::
from SimPEG import Examples
Examples.DC_Analytic_Dipole.run()
.. literalinclude:: ../../SimPEG/Examples/DC_Analytic_Dipole.py
:language: python
:linenos:
@@ -0,0 +1,28 @@
.. _examples_DC_Forward_PseudoSection:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
DC Forward Simulation
=====================
Forward model conductive spheres in a half-space and plot a pseudo-section
Created by @fourndo on Mon Feb 01 19:28:06 2016
.. plot::
from SimPEG import Examples
Examples.DC_Forward_PseudoSection.run()
.. literalinclude:: ../../SimPEG/Examples/DC_Forward_PseudoSection.py
:language: python
:linenos:
@@ -0,0 +1,27 @@
.. _examples_MT_1D_ForwardAndInversion:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
MT: 1D: Inversion
=======================
Forward model 1D MT data.
Setup and run a MT 1D inversion.
.. plot::
from SimPEG import Examples
Examples.MT_1D_ForwardAndInversion.run()
.. literalinclude:: ../../SimPEG/Examples/MT_1D_ForwardAndInversion.py
:language: python
:linenos:
+26
View File
@@ -0,0 +1,26 @@
.. _examples_MT_3D_Foward:
.. --------------------------------- ..
.. ..
.. THIS FILE IS AUTO GENEREATED ..
.. ..
.. SimPEG/Examples/__init__.py ..
.. ..
.. --------------------------------- ..
MT: 3D: Forward
=======================
Forward model 3D MT data.
.. plot::
from SimPEG import Examples
Examples.MT_3D_Foward.run()
.. literalinclude:: ../../SimPEG/Examples/MT_3D_Foward.py
:language: python
:linenos:
+1 -3
View File
@@ -49,9 +49,7 @@ Examples
.. toctree::
:maxdepth: 2
api_Examples
Packages
********
@@ -60,9 +58,9 @@ Packages
:maxdepth: 3
em/index
mt/index
flow/index
Finite Volume
*************
+19
View File
@@ -0,0 +1,19 @@
Magnetotellurics
****************
SimPEG (Simulation and Parameter Estimation in Geophysics) is a python
package for simulation and gradient based parameter estimation in the
context of geoscience applications.
simpegMT uses SimPEG as the framework for the forward and inverse
magnetotellurics geophysical problems.
Problem
=======
.. autoclass:: SimPEG.MT.BaseMT.BaseMTProblem
:show-inheritance:
:members:
:undoc-members:
+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__':
+12
View File
@@ -0,0 +1,12 @@
import os
import glob
import unittest
if __name__ == '__main__':
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
unittest.TextTestRunner(verbosity=2).run(testSuite)
+77
View File
@@ -0,0 +1,77 @@
import unittest
from SimPEG import *
import SimPEG.DCIP as DC
class DCProblemTests(unittest.TestCase):
def setUp(self):
aSpacing=2.5
nElecs=10
surveySize = nElecs*aSpacing - aSpacing
cs = surveySize/nElecs/4
mesh = Mesh.TensorMesh([
[(cs,10, -1.3),(cs,surveySize/cs),(cs,10, 1.3)],
[(cs,3, -1.3),(cs,3,1.3)],
# [(cs,5, -1.3),(cs,10)]
],'CN')
srcList = DC.Utils.WennerSrcList(nElecs, aSpacing, in2D=True)
survey = DC.SurveyDC(srcList)
problem = DC.ProblemDC_CC(mesh)
problem.pair(survey)
mSynth = np.ones(mesh.nC)
survey.makeSyntheticData(mSynth)
# Now set up the problem to do some minimization
dmis = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Tikhonov(mesh)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
inv = Inversion.BaseInversion(invProb)
self.inv = inv
self.reg = reg
self.p = problem
self.mesh = mesh
self.m0 = mSynth
self.survey = survey
self.dmis = dmis
def test_misfit(self):
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
self.assertTrue(passed)
def test_adjoint(self):
# Adjoint Test
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))
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
passed = np.abs(wtJv - vtJtw) < 1e-10
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
self.assertTrue(passed)
def test_dataObj(self):
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
self.assertTrue(passed)
def test_massMatrices(self):
Gu = np.random.rand(self.mesh.nF)
def derChk(m):
self.p.curModel = m
return [self.p.Msig * Gu, self.p.dMdsig(Gu)]
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
self.assertTrue(passed)
if __name__ == '__main__':
unittest.main()
+65
View File
@@ -0,0 +1,65 @@
import unittest
import SimPEG.DCIP as DC
from SimPEG import *
class IPforwardTests(unittest.TestCase):
def test_IPforward(self):
cs = 12.5
nc = 200/cs+1
hx = [(cs,7, -1.3),(cs,nc),(cs,7, 1.3)]
hy = [(cs,7, -1.3),(cs,int(nc/2+1)),(cs,7, 1.3)]
hz = [(cs,7, -1.3),(cs,int(nc/2+1))]
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
sighalf = 1e-2
sigma = np.ones(mesh.nC)*sighalf
p0 = np.r_[-50., 50., -50.]
p1 = np.r_[ 50.,-50., -150.]
blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
sigma[blk_ind] = 1e-3
eta = np.zeros_like(sigma)
eta[blk_ind] = 0.1
sigmaInf = sigma.copy()
sigma0 = sigma*(1-eta)
nElecs = 11
x_temp = np.linspace(-100, 100, nElecs)
aSpacing = x_temp[1]-x_temp[0]
y_temp = 0.
xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
survey = DC.SurveyDC(srcList)
imap = Maps.IdentityMap(mesh)
problem = DC.ProblemDC_CC(mesh, mapping=imap)
try:
from pymatsolver import MumpsSolver
solver = MumpsSolver
except ImportError, e:
solver = SolverLU
problem.Solver = solver
problem.pair(survey)
phi0 = survey.dpred(sigma0)
phiInf = survey.dpred(sigmaInf)
phiIP_true = phi0-phiInf
surveyIP = DC.SurveyIP(srcList)
problemIP = DC.ProblemIP(mesh, sigma=sigma)
problemIP.pair(surveyIP)
problemIP.Solver = solver
phiIP_approx = surveyIP.dpred(eta)
err = np.linalg.norm(phiIP_true-phiIP_approx) / np.linalg.norm(phiIP_true)
self.assertTrue(err < 0.02)
if __name__ == '__main__':
unittest.main()
+82
View File
@@ -0,0 +1,82 @@
import unittest
from SimPEG import *
import SimPEG.DCIP as DC
class IPProblemTests(unittest.TestCase):
def setUp(self):
cs = 12.5
nc = 500/cs+1
hx = [(cs,0, -1.3),(cs,nc),(cs,0, 1.3)]
hy = [(cs,0, -1.3),(cs,int(nc/2+1)),(cs,0, 1.3)]
hz = [(cs,0, -1.3),(cs,int(nc/2+1))]
mesh = Mesh.TensorMesh([hx, hy, hz], 'CCN')
sighalf = 1e-2
sigma = np.ones(mesh.nC)*sighalf
p0 = np.r_[-50., 50., -50.]
p1 = np.r_[ 50.,-50., -150.]
blk_ind = Utils.ModelBuilder.getIndicesBlock(p0, p1, mesh.gridCC)
sigma[blk_ind] = 1e-3
eta = np.zeros_like(sigma)
eta[blk_ind] = 0.1
nElecs = 5
x_temp = np.linspace(-250, 250, nElecs)
aSpacing = x_temp[1]-x_temp[0]
y_temp = 0.
xyz = Utils.ndgrid(x_temp, np.r_[y_temp], np.r_[0.])
srcList = DC.Utils.WennerSrcList(nElecs,aSpacing)
survey = DC.SurveyIP(srcList)
imap = Maps.IdentityMap(mesh)
problem = DC.ProblemIP(mesh, sigma=sigma, mapping= imap)
problem.pair(survey)
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except ImportError, e:
problem.Solver = SolverLU
mSynth = eta
survey.makeSyntheticData(mSynth)
# Now set up the problem to do some minimization
dmis = DataMisfit.l2_DataMisfit(survey)
reg = Regularization.Tikhonov(mesh)
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
inv = Inversion.BaseInversion(invProb)
self.inv = inv
self.reg = reg
self.p = problem
self.mesh = mesh
self.m0 = mSynth
self.survey = survey
self.dmis = dmis
def test_misfit(self):
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
passed = Tests.checkDerivative(derChk, self.m0*0, plotIt=False)
self.assertTrue(passed)
def test_adjoint(self):
# Adjoint Test
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))
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
passed = np.abs(wtJv - vtJtw) < 1e-10
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
self.assertTrue(passed)
def test_dataObj(self):
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
self.assertTrue(passed)
if __name__ == '__main__':
unittest.main()
+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)
@@ -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)
+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.
+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.]))
+12
View File
@@ -0,0 +1,12 @@
import os
import glob
import unittest
if __name__ == '__main__':
test_file_strings = glob.glob('test_*.py')
module_strings = [str[0:len(str)-3] for str in test_file_strings]
suites = [unittest.defaultTestLoader.loadTestsFromName(str) for str
in module_strings]
testSuite = unittest.TestSuite(suites)
unittest.TextTestRunner(verbosity=2).run(testSuite)
@@ -0,0 +1,48 @@
import unittest
from SimPEG import *
from SimPEG import MT
TOL = 1e-6
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(-z.imag,z.real)*(180/np.pi)
return app_res, app_phs
def appResNorm(sigmaHalf):
nFreq = 26
m1d = Mesh.TensorMesh([[(100,5,1.5),(100.,10),(100,5,1.5)]], x0=['C'])
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC[:]>200] = 1e-8
# Calculate the analytic fields
freqs = np.logspace(4,-4,nFreq)
Z = []
for freq in freqs:
Ed, Eu, Hd, Hu = MT.Utils.getEHfields(m1d,sigma,freq,np.array([200]))
Z.append((Ed + Eu)/(Hd + Hu))
Zarr = np.concatenate(Z)
app_r, app_p = appResPhs(freqs,Zarr)
return np.linalg.norm(np.abs(app_r - np.ones(nFreq)/sigmaHalf)) / np.log10(sigmaHalf)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
def test_appRes2en1(self):self.assertLess(appResNorm(2e-1), TOL)
def test_appRes2en2(self):self.assertLess(appResNorm(2e-2), TOL)
def test_appRes2en3(self):self.assertLess(appResNorm(2e-3), TOL)
def test_appRes2en4(self):self.assertLess(appResNorm(2e-4), TOL)
def test_appRes2en5(self):self.assertLess(appResNorm(2e-5), TOL)
def test_appRes2en6(self):self.assertLess(appResNorm(2e-6), TOL)
if __name__ == '__main__':
unittest.main()
@@ -0,0 +1,162 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
TOLr = 5e-2
TOLp = 5e-2
def setupSurvey(sigmaHalf,tD=True):
# Frequency
nFreq = 33
freqs = np.logspace(3,-3,nFreq)
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],10,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC > 0 ] = 1e-8
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
if tD:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def appRes_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.eval(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
def appPhs_TotalFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf)
problem = MT.Problem1D.eForm_TotalField(mesh)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.eval(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
def appRes_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.eval(fields)
# Calculate the app res and phs
app_r = np.array(getAppResPhs(data))[:,0]
return np.linalg.norm(np.abs(app_r - np.ones(survey.nFreq)/sigmaHalf)*sigmaHalf)
def appPhs_psFieldNorm(sigmaHalf):
# Make the survey
survey, sigma, mesh = setupSurvey(sigmaHalf,False)
problem = MT.Problem1D.eForm_psField(mesh, sigmaPrimary = sigma)
problem.pair(survey)
# Get the fields
fields = problem.fields(sigma)
# Project the data
data = survey.eval(fields)
# Calculate the app phs
app_p = np.array(getAppResPhs(data))[:,1]
return np.linalg.norm(np.abs(app_p - np.ones(survey.nFreq)*45)/ 45)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
# Total Fields
# def test_appRes2en1(self):self.assertLess(appRes_TotalFieldNorm(2e-1), TOLr)
# def test_appPhs2en1(self):self.assertLess(appPhs_TotalFieldNorm(2e-1), TOLp)
# def test_appRes2en2(self):self.assertLess(appRes_TotalFieldNorm(2e-2), TOLr)
# def test_appPhs2en2(self):self.assertLess(appPhs_TotalFieldNorm(2e-2), TOLp)
# def test_appRes2en3(self):self.assertLess(appRes_TotalFieldNorm(2e-3), TOLr)
# def test_appPhs2en3(self):self.assertLess(appPhs_TotalFieldNorm(2e-3), TOLp)
# def test_appRes2en4(self):self.assertLess(appRes_TotalFieldNorm(2e-4), TOLr)
# def test_appPhs2en4(self):self.assertLess(appPhs_TotalFieldNorm(2e-4), TOLp)
# def test_appRes2en5(self):self.assertLess(appRes_TotalFieldNorm(2e-5), TOLr)
# def test_appPhs2en5(self):self.assertLess(appPhs_TotalFieldNorm(2e-5), TOLp)
# def test_appRes2en6(self):self.assertLess(appRes_TotalFieldNorm(2e-6), TOLr)
# def test_appPhs2en6(self):self.assertLess(appPhs_TotalFieldNorm(2e-6), TOLp)
# Primary/secondary
def test_appRes2en2_ps(self):self.assertLess(appRes_psFieldNorm(2e-2), TOLr)
def test_appPhs2en2_ps(self):self.assertLess(appPhs_psFieldNorm(2e-2), TOLp)
if __name__ == '__main__':
unittest.main()
@@ -0,0 +1,135 @@
import unittest
import SimPEG as simpeg
from SimPEG import MT
from SimPEG.Utils import meshTensor
import numpy as np
# Define the tolerances
TOLr = 5e-2
TOLp = 5e-2
def setupSurvey(sigmaHalf,tD=True):
# Frequency
nFreq = 33
freqs = np.logspace(3,-3,nFreq)
# Make the mesh
ct = 5
air = meshTensor([(ct,25,1.3)])
# coreT0 = meshTensor([(ct,15,1.2)])
# coreT1 = np.kron(meshTensor([(coreT0[-1],15,1.3)]),np.ones((7,)))
core = np.concatenate( ( np.kron(meshTensor([(ct,15,-1.2)]),np.ones((10,))) , meshTensor([(ct,20)]) ) )
bot = meshTensor([(core[0],15,-1.3)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
# Make the model
sigma = np.zeros(m1d.nC) + sigmaHalf
sigma[m1d.gridCC > 0 ] = 1e-8
sigmaBack = sigma.copy()
# Add structure
shallow = (m1d.gridCC < -200) * (m1d.gridCC > -600)
deep = (m1d.gridCC < -3000) * (m1d.gridCC > -5000)
sigma[shallow] = 1
sigma[deep] = 0.1
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
if tD:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1DhomotD(rxList,freq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
survey = MT.Survey(srcList)
return survey, sigma, m1d
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
zList = []
for src in MTdata.survey.srcList:
zc = [src.freq]
for rx in src.rxList:
if 'i' in rx.rxType:
m=1j
else:
m = 1
zc.append(m*MTdata[src,rx])
zList.append(zc)
return [appResPhs(zList[i][0],np.sum(zList[i][1:3])) for i in np.arange(len(zList))]
def calculateAnalyticSolution(srcList,mesh,model):
surveyAna = MT.Survey(srcList)
data1D = MT.Data(surveyAna)
for src in surveyAna.srcList:
elev = src.rxList[0].locs[0]
anaEd, anaEu, anaHd, anaHu = MT.Utils.MT1Danalytic.getEHfields(mesh,model,src.freq,elev)
anaE = anaEd+anaEu
anaH = anaHd+anaHu
# Scale the solution
# anaE = (anaEtemp/anaEtemp[-1])#.conj()
# anaH = (anaHtemp/anaEtemp[-1])#.conj()
anaZ = anaE/anaH
for rx in src.rxList:
data1D[src,rx] = getattr(anaZ, rx.projComp)
return data1D
def dataMis_AnalyticTotalDomain(sigmaHalf):
# Make the survey
# Total domain solution
surveyTD, sigma, mesh = setupSurvey(sigmaHalf)
problemTD = MT.Problem1D.eForm_TotalField(mesh)
problemTD.pair(surveyTD)
# Analytic data
dataAnaObj = calculateAnalyticSolution(surveyTD.srcList,mesh,sigma)
# dataTDObj = MT.DataMT.DataMT(surveyTD, surveyTD.dpred(sigma))
dataTD = surveyTD.dpred(sigma)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataTD - dataAna)/dataAna < 2.)
# surveyTD.dtrue = -simpeg.mkvc(dataAna,2)
# surveyTD.dobs = -simpeg.mkvc(dataAna,2)
# surveyTD.Wd = np.ones(surveyTD.dtrue.shape) #/(np.abs(surveyTD.dtrue)*0.01)
# # Setup the data misfit
# dmis = simpeg.DataMisfit.l2_DataMisfit(surveyTD)
# dmis.Wd = surveyTD.Wd
# return dmis.eval(sigma)
def dataMis_AnalyticPrimarySecondary(sigmaHalf):
# Make the survey
# Primary secondary
surveyPS, sigmaPS, mesh = setupSurvey(sigmaHalf,tD=False)
problemPS = MT.Problem1D.eForm_psField(mesh)
problemPS.sigmaPrimary = sigmaPS
problemPS.pair(surveyPS)
# Analytic data
dataAnaObj = calculateAnalyticSolution(surveyPS.srcList,mesh,sigmaPS)
dataPS = surveyPS.dpred(sigmaPS)
dataAna = simpeg.mkvc(dataAnaObj)
return np.all((dataPS - dataAna)/dataAna < 2.)
class TestNumericVsAnalytics(unittest.TestCase):
def setUp(self):
pass
# Total Fields
# def test_appRes2en2(self):self.assertTrue(dataMis_AnalyticTotalDomain(2e-2))
# Primary/secondary
def test_appRes2en2_ps(self):self.assertTrue(dataMis_AnalyticPrimarySecondary(2e-2))
if __name__ == '__main__':
unittest.main()
+271
View File
@@ -0,0 +1,271 @@
# Test functions
from glob import glob
import numpy as np, sys, os, time, scipy, subprocess
import SimPEG as simpeg
import unittest
from SimPEG import MT
from SimPEG.Utils import meshTensor
from scipy.constants import mu_0
TOLr = 5e-2
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]
addrandoms = True
def getInputs():
"""
Function that returns Mesh, freqs, rx_loc, elev.
"""
# Make a mesh
# M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
# M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,6),(1000,6,1.5)],[(1000,6,-1.5),(1000.,2),(1000,6,1.5)],[(1000,6,-1.3),(1000.,6),(1000,6,1.3)]], x0=['C','C','C'])# Setup the model
M = simpeg.Mesh.TensorMesh([[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(1000,6,-1.5),(1000.,4),(1000,6,1.5)],[(500,8,-1.3),(500.,8),(500,8,1.3)]], x0=['C','C','C'])# Setup the model
# Set the frequencies
freqs = np.logspace(1,-3,5)
elev = 0
## Setup the the survey object
# Receiver locations
rx_x, rx_y = np.meshgrid(np.arange(-1000,1001,500),np.arange(-1000,1001,500))
rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),elev+np.zeros((np.prod(rx_x.shape),1))))
return M, freqs, rx_loc, elev
def random(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Backround
sigBG = np.ones(M.nC)*conds
# Add randomness to the model (10% of the value).
sig = np.exp( np.log(sigBG) + np.random.randn(M.nC)*(conds)*1e-1 )
return (M, freqs, sig, sigBG, rx_loc)
def halfSpace(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
# conds = [1e-2]
groundInd = ccM[:,2] < elev
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds
return (M, freqs, sig, sigBG, rx_loc)
def blockInhalfSpace(conds):
''' Returns a halfspace model based on the inputs'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
# conds = [1e-2]
groundInd = ccM[:,2] < elev
sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,np.array([-1000,-1000,-1500]),np.array([1000,1000,-1000]),conds)
sig[~groundInd] = 1e-8
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds[1]
return (M, freqs, sig, sigBG, rx_loc)
def twoLayer(conds):
''' Returns a 2 layer model based on the conductivity values given'''
M, freqs, rx_loc, elev = getInputs()
# Model
ccM = M.gridCC
groundInd = ccM[:,2] < elev
botInd = ccM[:,2] < -3000
sig = np.zeros(M.nC) + 1e-8
sig[groundInd] = conds[1]
sig[botInd] = conds[0]
# Set the background, not the same as the model
sigBG = np.zeros(M.nC) + 1e-8
sigBG[groundInd] = conds[1]
return (M, freqs, sig, sigBG, rx_loc)
def setupSimpegMTfwd_eForm_ps(inputSetup,comp='Imp',singleFreq=False,expMap=True):
M,freqs,sig,sigBG,rx_loc = inputSetup
# Make a receiver list
rxList = []
if comp == 'All':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(rx_loc,rxType))
elif comp == 'Imp':
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi']:
rxList.append(MT.Rx(rx_loc,rxType))
elif comp == 'Tip':
for rxType in ['tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(rx_loc,rxType))
else:
rxList.append(MT.Rx(rx_loc,comp))
# Source list
srcList =[]
if singleFreq:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,singleFreq))
else:
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
## Setup the problem object
sigma1d = M.r(sigBG,'CC','CC','M')[0,0,:]
if expMap:
problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= np.log(sigma1d) )
problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
problem.curModel = np.log(sig)
else:
problem = MT.Problem3D.eForm_ps(M,sigmaPrimary= sigma1d)
problem.curModel = sig
problem.pair(survey)
problem.verbose = False
try:
from pymatsolver import MumpsSolver
problem.Solver = MumpsSolver
except:
pass
return (survey, problem)
def getAppResPhs(MTdata):
# Make impedance
def appResPhs(freq,z):
app_res = ((1./(8e-7*np.pi**2))/freq)*np.abs(z)**2
app_phs = np.arctan2(z.imag,z.real)*(180/np.pi)
return app_res, app_phs
recData = MTdata.toRecArray('Complex')
return appResPhs(recData['freq'],recData['zxy']), appResPhs(recData['freq'],recData['zyx'])
def JvecAdjointTest(inputSetup,comp='All',freq=False):
(M, freqs, sig, sigBG, rx_loc) = inputSetup
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp='All',singleFreq=freq)
print 'Adjoint test of eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,str(survey.freqs))
m = sig
u = problem.fields(m)
v = np.random.rand(survey.nD,)
# print problem.PropMap.PropModel.nP
w = np.random.rand(problem.mesh.nC,)
vJw = v.ravel().dot(problem.Jvec(m, w, u))
wJtv = w.ravel().dot(problem.Jtvec(m, v, u))
tol = np.max([TOL*(10**int(np.log10(np.abs(vJw)))),FLR])
print ' vJw wJtv vJw - wJtv tol abs(vJw - wJtv) < tol'
print vJw, wJtv, vJw - wJtv, tol, np.abs(vJw - wJtv) < tol
return np.abs(vJw - wJtv) < tol
# Test the Jvec derivative
def DerivJvecTest(inputSetup,comp='All',freq=False,expMap=True):
(M, freqs, sig, sigBG, rx_loc) = inputSetup
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp=comp,singleFreq=freq,expMap=expMap)
print 'Derivative test of Jvec for eForm primary/secondary for {:s} comp at {:s}\n'.format(comp,survey.freqs)
# problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
# problem.sigmaPrimary = np.log(sigBG)
x0 = np.log(sigBG)
# cond = sig[0]
# x0 = np.log(np.ones(problem.mesh.nC)*cond)
# problem.sigmaPrimary = x0
# if True:
# x0 = x0 + np.random.randn(problem.mesh.nC)*cond*1e-1
survey = problem.survey
def fun(x):
return survey.dpred(x), lambda x: problem.Jvec(x0, x)
return simpeg.Tests.checkDerivative(fun, x0, num=3, plotIt=False, eps=FLR)
def DerivProjfieldsTest(inputSetup,comp='All',freq=False):
survey, problem = setupSimpegMTfwd_eForm_ps(inputSetup,comp,freq)
print 'Derivative test of data projection for eFormulation primary/secondary\n\n'
# problem.mapping = simpeg.Maps.ExpMap(problem.mesh)
# Initate things for the derivs Test
src = survey.srcList[0]
rx = src.rxList[0]
u0x = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0y = np.random.randn(survey.mesh.nE)+np.random.randn(survey.mesh.nE)*1j
u0 = np.vstack((simpeg.mkvc(u0x,2),simpeg.mkvc(u0y,2)))
f0 = problem.fieldsPair(survey.mesh,survey)
# u0 = np.hstack((simpeg.mkvc(u0_px,2),simpeg.mkvc(u0_py,2)))
f0[src,'e_pxSolution'] = u0[:len(u0)/2]#u0x
f0[src,'e_pySolution'] = u0[len(u0)/2::]#u0y
def fun(u):
f = problem.fieldsPair(survey.mesh,survey)
f[src,'e_pxSolution'] = u[:len(u)/2]
f[src,'e_pySolution'] = u[len(u)/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)
def appResPhsHalfspace_eFrom_ps_Norm(sigmaHalf,appR=True,expMap=False):
if appR:
label = 'resistivity'
else:
label = 'phase'
# Make the survey and the problem
survey, problem = setupSimpegMTfwd_eForm_ps(halfSpace(sigmaHalf),expMap=expMap)
print 'Apperent {:s} test of eFormulation primary/secondary at {:g}\n\n'.format(label,sigmaHalf)
data = problem.dataPair(survey,survey.dpred(problem.curModel))
# Calculate the app phs
app_rpxy, app_rpyx = np.array(getAppResPhs(data))
if appR:
return np.all(np.abs(app_rpxy[0,:] - 1./sigmaHalf) * sigmaHalf < .4)
else:
return np.all(np.abs(app_rpxy[1,:] + 135) / 135 < .4)
class TestAnalytics(unittest.TestCase):
def setUp(self):
pass
# # Test apparent resistivity and phase
def test_appRes1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2))
def test_appPhs1en2(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-2,False))
def test_appRes1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3))
def test_appPhs1en3(self):self.assertTrue(appResPhsHalfspace_eFrom_ps_Norm(1e-3,False))
# Do a derivative test
def test_derivProj1(self):self.assertTrue(DerivProjfieldsTest(halfSpace(1e-2)))
# Do a derivative test of Jvec
# def test_derivJvec_zxxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxr',.1))
# def test_derivJvec_zxxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxxi',.1))
# def test_derivJvec_zxyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyr',.1))
# def test_derivJvec_zxyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zxyi',.1))
# def test_derivJvec_zyxr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxr',.1))
# def test_derivJvec_zyxi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyxi',.1))
# def test_derivJvec_zyyr(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyr',.1))
# def test_derivJvec_zyyi(self):self.assertTrue(DerivJvecTest(random(1e-2),'zyyi',.1))
def test_derivJvec_All(self):self.assertTrue(DerivJvecTest(random(1e-2),'All',.1))
# Test the adjoint of Jvec and Jtvec
# def test_JvecAdjoint_zxxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxr',.1))
# def test_JvecAdjoint_zxxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxxi',.1))
# def test_JvecAdjoint_zxyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyr',.1))
# def test_JvecAdjoint_zxyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zxyi',.1))
# def test_JvecAdjoint_zyxr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxr',.1))
# def test_JvecAdjoint_zyxi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyxi',.1))
# def test_JvecAdjoint_zyyr(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyr',.1))
# def test_JvecAdjoint_zyyi(self):self.assertTrue(JvecAdjointTest(random(1e-2),'zyyi',.1))
def test_JvecAdjoint_All(self):self.assertTrue(JvecAdjointTest(random(1e-2),'All',.1))
if __name__ == '__main__':
unittest.main()