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

Merge remote-tracking branch 'origin/dev' into mt/dev

Browse Source
This commit is contained in:
GudniRos
2016-03-29 16:52:01 -07:00
parent 127c51974f 7fd6ddf62a
commit 009806f2ba
52 changed files with 4181 additions and 812 deletions
Download Patch File Download Diff File Expand all files Collapse all files
SimPEG/DCIP/BaseDC.py
+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
SimPEG/DCIP/BaseIP.py
+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
SimPEG/DCIP/DCIPUtils.py
+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
SimPEG/DCIP/Utils.py
+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
SimPEG/DCIP/__init__.py
+4
View File
@@ -0,0 +1,4 @@
from BaseDC import *
from BaseIP import *
from DCIPUtils import *
import Utils
SimPEG/EM/FDEM/FDEM.py
+49 -65
View File
@@ -54,8 +54,7 @@ class BaseFDEMProblem(BaseEMProblem):
Ainv = self.Solver(A, **self.solverOpts)
sol = Ainv * rhs
Srcs = self.survey.getSrcByFreq(freq)
ftype = self._fieldType + 'Solution'
F[Srcs, ftype] = sol
F[Srcs, self._solutionType] = sol
Ainv.clean()
return F
@@ -78,30 +77,19 @@ class BaseFDEMProblem(BaseEMProblem):
Jv = self.dataPair(self.survey)
for freq in self.survey.freqs:
A = self.getA(freq) #
A = self.getA(freq)
Ainv = self.Solver(A, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, ftype]
dA_dm = self.getADeriv_m(freq, u_src, v)
dRHS_dm = self.getRHSDeriv_m(freq, src, v)
du_dm = Ainv * ( - dA_dm + dRHS_dm )
u_src = u[src, self._solutionType]
dA_dm_v = self.getADeriv(freq, u_src, v)
dRHS_dm_v = self.getRHSDeriv(freq, src, v)
du_dm_v = Ainv * ( - dA_dm_v + dRHS_dm_v )
for rx in src.rxList:
df_duFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_dudu_dm = df_duFun(src, du_dm, adjoint=False)
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
df_dm = df_dmFun(src, v, adjoint=False)
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
Jv[src, rx] = P(Df_Dm)
df_dmFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
Jv[src, rx] = rx.evalDeriv(src, self.mesh, u, df_dm_v)
Ainv.clean()
return Utils.mkvc(Jv)
@@ -132,32 +120,28 @@ class BaseFDEMProblem(BaseEMProblem):
ATinv = self.Solver(AT, **self.solverOpts)
for src in self.survey.getSrcByFreq(freq):
ftype = self._fieldType + 'Solution'
u_src = u[src, ftype]
u_src = u[src, self._solutionType]
for rx in src.rxList:
PTv = rx.projectFieldsDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, u, v[src, rx], adjoint=True) # wrt u, need possibility wrt m
df_duTFun = getattr(u, '_%sDeriv'%rx.projField, None)
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
df_duTFun = getattr(u, '_%sDeriv_u'%rx.projField, None)
df_duT = df_duTFun(src, PTv, adjoint=True)
ATinvdf_duT = ATinv * df_duT
dA_dmT = self.getADeriv_m(freq, u_src, ATinvdf_duT, adjoint=True)
dRHS_dmT = self.getRHSDeriv_m(freq,src, ATinvdf_duT, adjoint=True)
dA_dmT = self.getADeriv(freq, u_src, ATinvdf_duT, adjoint=True)
dRHS_dmT = self.getRHSDeriv(freq, src, ATinvdf_duT, adjoint=True)
du_dmT = -dA_dmT + dRHS_dmT
df_dmFun = getattr(u, '_%sDeriv_m'%rx.projField, None)
dfT_dm = df_dmFun(src, PTv, adjoint=True)
df_dmT = df_dmT + du_dmT
du_dmT += dfT_dm
# TODO: this should be taken care of by the reciever
# TODO: this should be taken care of by the reciever?
real_or_imag = rx.projComp
if real_or_imag is 'real':
Jtv += np.array(du_dmT,dtype=complex).real
Jtv += np.array(df_dmT, dtype=complex).real
elif real_or_imag is 'imag':
Jtv += - np.array(du_dmT,dtype=complex).real
Jtv += - np.array(df_dmT, dtype=complex).real
else:
raise Exception('Must be real or imag')
@@ -174,10 +158,10 @@ class BaseFDEMProblem(BaseEMProblem):
:return: S_m, S_e (nE or nF, nSrc)
"""
Srcs = self.survey.getSrcByFreq(freq)
if self._eqLocs is 'FE':
if self._formulation is 'EB':
S_m = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
elif self._eqLocs is 'EF':
elif self._formulation is 'HJ':
S_m = np.zeros((self.mesh.nE,len(Srcs)), dtype=complex)
S_e = np.zeros((self.mesh.nF,len(Srcs)), dtype=complex)
@@ -213,9 +197,9 @@ class Problem_e(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'e'
_eqLocs = 'FE'
fieldsPair = Fields_e
_solutionType = 'eSolution'
_formulation = 'EB'
fieldsPair = Fields_e
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -239,7 +223,7 @@ class Problem_e(BaseFDEMProblem):
return C.T*MfMui*C + 1j*omega(freq)*MeSigma
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -280,7 +264,7 @@ class Problem_e(BaseFDEMProblem):
return C.T * (MfMui * S_m) -1j * omega(freq) * S_e
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -324,9 +308,9 @@ class Problem_b(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'b'
_eqLocs = 'FE'
fieldsPair = Fields_b
_solutionType = 'bSolution'
_formulation = 'EB'
fieldsPair = Fields_b
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -354,7 +338,7 @@ class Problem_b(BaseFDEMProblem):
return MfMui.T*A
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -411,7 +395,7 @@ class Problem_b(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -472,9 +456,9 @@ class Problem_j(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'j'
_eqLocs = 'EF'
fieldsPair = Fields_j
_solutionType = 'jSolution'
_formulation = 'HJ'
fieldsPair = Fields_j
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -503,7 +487,7 @@ class Problem_j(BaseFDEMProblem):
return A
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -524,16 +508,16 @@ class Problem_j(BaseFDEMProblem):
MeMuI = self.MeMuI
MfRho = self.MfRho
C = self.mesh.edgeCurl
MfRhoDeriv_m = self.MfRhoDeriv(u)
MfRhoDeriv = self.MfRhoDeriv(u)
if adjoint:
if self._makeASymmetric is True:
v = MfRho * v
return MfRhoDeriv_m.T * (C * (MeMuI.T * (C.T * v)))
return MfRhoDeriv.T * (C * (MeMuI.T * (C.T * v)))
if self._makeASymmetric is True:
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv_m * v) )))
return C * (MeMuI * (C.T * (MfRhoDeriv_m * v)))
return MfRho.T * (C * ( MeMuI * (C.T * (MfRhoDeriv * v) )))
return C * (MeMuI * (C.T * (MfRhoDeriv * v)))
def getRHS(self, freq):
@@ -560,7 +544,7 @@ class Problem_j(BaseFDEMProblem):
return RHS
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
@@ -610,9 +594,9 @@ class Problem_h(BaseFDEMProblem):
:param SimPEG.Mesh mesh: mesh
"""
_fieldType = 'h'
_eqLocs = 'EF'
fieldsPair = Fields_h
_solutionType = 'hSolution'
_formulation = 'HJ'
fieldsPair = Fields_h
def __init__(self, mesh, **kwargs):
BaseFDEMProblem.__init__(self, mesh, **kwargs)
@@ -635,7 +619,7 @@ class Problem_h(BaseFDEMProblem):
return C.T * (MfRho * C) + 1j*omega(freq)*MeMu
def getADeriv_m(self, freq, u, v, adjoint=False):
def getADeriv(self, freq, u, v, adjoint=False):
"""
Product of the derivative of our system matrix with respect to the model and a vector
@@ -652,11 +636,11 @@ class Problem_h(BaseFDEMProblem):
MeMu = self.MeMu
C = self.mesh.edgeCurl
MfRhoDeriv_m = self.MfRhoDeriv(C*u)
MfRhoDeriv = self.MfRhoDeriv(C*u)
if adjoint:
return MfRhoDeriv_m.T * (C * v)
return C.T * (MfRhoDeriv_m * v)
return MfRhoDeriv.T * (C * v)
return C.T * (MfRhoDeriv * v)
def getRHS(self, freq):
"""
@@ -677,7 +661,7 @@ class Problem_h(BaseFDEMProblem):
return S_m + C.T * ( MfRho * S_e )
def getRHSDeriv_m(self, freq, src, v, adjoint=False):
def getRHSDeriv(self, freq, src, v, adjoint=False):
"""
Derivative of the right hand side with respect to the model
SimPEG/EM/FDEM/FieldsFDEM.py
+698 -313
View File
File diff suppressed because it is too large Load Diff
SimPEG/EM/FDEM/SrcFDEM.py
+117 -79
View File
@@ -1,7 +1,7 @@
from SimPEG import Survey, Problem, Utils, np, sp
from scipy.constants import mu_0
from SimPEG.EM.Utils import *
from SimPEG.Utils import Zero
from SimPEG.Utils import Zero
class BaseSrc(Survey.BaseSrc):
"""
@@ -14,7 +14,7 @@ class BaseSrc(Survey.BaseSrc):
def eval(self, prob):
"""
Evaluate the source terms.
Evaluate the source terms.
- :math:`S_m` : magnetic source term
- :math:`S_e` : electric source term
@@ -36,12 +36,12 @@ class BaseSrc(Survey.BaseSrc):
:param numpy.ndarray v: vector to take product with
:param bool adjoint: adjoint?
:rtype: (numpy.ndarray, numpy.ndarray)
:return: tuple with magnetic source term and electric source term derivatives times a vector
:return: tuple with magnetic source term and electric source term derivatives times a vector
"""
if v is not None:
return self.S_mDeriv(prob,v,adjoint), self.S_eDeriv(prob,v,adjoint)
if v is not None:
return self.S_mDeriv(prob, v, adjoint), self.S_eDeriv(prob, v, adjoint)
else:
return lambda v: self.S_mDeriv(prob,v,adjoint), lambda v: self.S_eDeriv(prob,v,adjoint)
return lambda v: self.S_mDeriv(prob, v, adjoint), lambda v: self.S_eDeriv(prob, v, adjoint)
def bPrimary(self, prob):
"""
@@ -49,7 +49,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic flux density
:return: primary magnetic flux density
"""
return Zero()
@@ -59,7 +59,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
return Zero()
@@ -69,7 +69,7 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary electric field
:return: primary electric field
"""
return Zero()
@@ -79,13 +79,13 @@ class BaseSrc(Survey.BaseSrc):
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: primary current density
:return: primary current density
"""
return Zero()
def S_m(self, prob):
"""
Magnetic source term
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -95,7 +95,7 @@ class BaseSrc(Survey.BaseSrc):
def S_e(self, prob):
"""
Electric source term
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
@@ -136,15 +136,26 @@ class RawVec_e(BaseSrc):
:param list rxList: receiver list
:param float freq: frequency
:param numpy.array S_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_e): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e,dtype=complex)
def __init__(self, rxList, freq, S_e, integrate=True): #, ePrimary=None, bPrimary=None, hPrimary=None, jPrimary=None):
self._S_e = np.array(S_e, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_e(self, prob):
"""
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
@@ -155,10 +166,11 @@ class RawVec_m(BaseSrc):
:param float freq: frequency
:param rxList: receiver list
:param numpy.array S_m: magnetic source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_m, integrate = True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m,dtype=complex)
def __init__(self, rxList, freq, S_m, integrate=True): #ePrimary=Zero(), bPrimary=Zero(), hPrimary=Zero(), jPrimary=Zero()):
self._S_m = np.array(S_m, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
@@ -166,12 +178,14 @@ class RawVec_m(BaseSrc):
def S_m(self, prob):
"""
Magnetic source term
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._S_m
return self._S_m
@@ -183,36 +197,51 @@ class RawVec(BaseSrc):
:param float freq: frequency
:param numpy.array S_m: magnetic source term
:param numpy.array S_e: electric source term
:param bool integrate: Integrate the source term (multiply by Me) [True]
"""
def __init__(self, rxList, freq, S_m, S_e, integrate = True):
self._S_m = np.array(S_m,dtype=complex)
self._S_e = np.array(S_e,dtype=complex)
def __init__(self, rxList, freq, S_m, S_e, integrate=True):
self._S_m = np.array(S_m, dtype=complex)
self._S_e = np.array(S_e, dtype=complex)
self.freq = float(freq)
self.integrate = integrate
BaseSrc.__init__(self, rxList)
def S_m(self, prob):
if prob._eqLocs is 'EF' and self.integrate is True:
"""
Magnetic source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: magnetic source term on mesh
"""
if prob._formulation is 'HJ' and self.integrate is True:
return prob.Me * self._S_m
return self._S_m
def S_e(self, prob):
if prob._eqLocs is 'FE' and self.integrate is True:
"""
Electric source term
:param Problem prob: FDEM Problem
:rtype: numpy.ndarray
:return: electric source term on mesh
"""
if prob._formulation is 'EB' and self.integrate is True:
return prob.Me * self._S_e
return self._S_e
class MagDipole(BaseSrc):
"""
"""
Point magnetic dipole source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency. Here we show the
derivation for E-B formulation noting that similar steps are followed for
the H-J formulation.
.. math::
.. math::
\mathbf{C} \mathbf{e} + i \omega \mathbf{b} = \mathbf{s_m} \\\\
{\mathbf{C}^T \mathbf{M_{\mu^{-1}}^f} \mathbf{b} - \mathbf{M_{\sigma}^e} \mathbf{e} = \mathbf{s_e}}
@@ -225,17 +254,17 @@ class MagDipole(BaseSrc):
and define a zero-frequency primary problem, noting that the source is
generated by a divergence free electric current
.. math::
.. math::
\mathbf{C} \mathbf{e^P} = \mathbf{s_m^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} - \mathbf{M_{\sigma}^e} \mathbf{e^P} = \mathbf{M^e} \mathbf{s_e^P}}
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
Since :math:`\mathbf{e^P}` is curl-free, divergence-free, we assume that there is no constant field background, the :math:`\mathbf{e^P} = 0`, so our primary problem is
.. math::
.. math::
\mathbf{e^P} = 0 \\\\
{\mathbf{C}^T \mathbf{{M_{\mu^{-1}}^f}^P} \mathbf{b^P} = \mathbf{s_e^P}}
Our secondary problem is then
Our secondary problem is then
.. math::
\mathbf{C} \mathbf{e^S} + i \omega \mathbf{b^S} = - i \omega \mathbf{b^P} \\\\
@@ -245,15 +274,15 @@ class MagDipole(BaseSrc):
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu=mu_0):
self.freq = float(freq)
self.loc = loc
self.orientation = orientation
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.moment = moment
self.mu = mu
self.integrate = False
@@ -265,17 +294,17 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
@@ -303,10 +332,10 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob,b)
return 1./self.mu * b
def S_m(self, prob):
"""
@@ -314,10 +343,12 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b_p = self.bPrimary(prob)
if prob._formulation is 'HJ':
b_p = prob.Me * b_p
return -1j*omega(self.freq)*b_p
def S_e(self, prob):
@@ -326,21 +357,21 @@ class MagDipole(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
@@ -353,21 +384,20 @@ class MagDipole_Bfield(BaseSrc):
fields from a magnetic dipole. No discrete curl is taken, so the magnetic
flux density may not be strictly divergence free.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
#TODO: neither does moment
def __init__(self, rxList, freq, loc, orientation='Z', moment=1., mu = mu_0):
self.freq = float(freq)
self.loc = loc
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.orientation = orientation
self.moment = moment
self.mu = mu
@@ -379,18 +409,18 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
@@ -418,10 +448,10 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return h_from_b(prob, b)
return 1/self.mu * b
def S_m(self, prob):
"""
@@ -429,9 +459,11 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
if prob._formulation is 'HJ':
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
@@ -440,20 +472,20 @@ class MagDipole_Bfield(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
@@ -463,22 +495,22 @@ class CircularLoop(BaseSrc):
"""
Circular loop magnetic source calculated by taking the curl of a magnetic
vector potential. By taking the discrete curl, we ensure that the magnetic
flux density is divergence free (no magnetic monopoles!).
flux density is divergence free (no magnetic monopoles!).
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
This approach uses a primary-secondary in frequency in the same fashion as the MagDipole.
:param list rxList: receiver list
:param float freq: frequency
:param numpy.ndarray loc: source location (ie: :code:`np.r_[xloc,yloc,zloc]`)
:param string orientation: 'X', 'Y', 'Z'
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
:param float moment: magnetic dipole moment
:param float mu: background magnetic permeability
"""
#TODO: right now, orientation doesn't actually do anything! The methods in SrcUtils should take care of that
def __init__(self, rxList, freq, loc, orientation='Z', radius = 1., mu=mu_0):
def __init__(self, rxList, freq, loc, orientation='Z', radius=1., mu=mu_0):
self.freq = float(freq)
self.orientation = orientation
assert orientation in ['X','Y','Z'], "Orientation (right now) doesn't actually do anything! The methods in SrcUtils should take care of this..."
self.radius = radius
self.mu = mu
self.loc = loc
@@ -491,17 +523,17 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
gridX = prob.mesh.gridEx
gridY = prob.mesh.gridEy
gridZ = prob.mesh.gridEz
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
gridX = prob.mesh.gridFx
gridY = prob.mesh.gridFy
gridZ = prob.mesh.gridFz
@@ -528,7 +560,7 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
return 1./self.mu*b
@@ -539,9 +571,11 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
b = self.bPrimary(prob)
if prob._formulation is 'HJ':
b = prob.Me * b
return -1j*omega(self.freq)*b
def S_e(self, prob):
@@ -550,22 +584,26 @@ class CircularLoop(BaseSrc):
:param Problem prob: FDEM problem
:rtype: numpy.ndarray
:return: primary magnetic field
:return: primary magnetic field
"""
if all(np.r_[self.mu] == np.r_[prob.curModel.mu]):
return Zero()
else:
eqLocs = prob._eqLocs
formulation = prob._formulation
if eqLocs is 'FE':
if formulation is 'EB':
mui_s = prob.curModel.mui - 1./self.mu
MMui_s = prob.mesh.getFaceInnerProduct(mui_s)
C = prob.mesh.edgeCurl
elif eqLocs is 'EF':
elif formulation is 'HJ':
mu_s = prob.curModel.mu - self.mu
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s,invMat=True)
MMui_s = prob.mesh.getEdgeInnerProduct(mu_s, invMat=True)
C = prob.mesh.edgeCurl.T
return -C.T * (MMui_s * self.bPrimary(prob))
SimPEG/EM/FDEM/SurveyFDEM.py
+44 -38
View File
@@ -3,6 +3,7 @@ from SimPEG.EM.Utils import *
from scipy.constants import mu_0
from SimPEG.Utils import Zero, Identity
import SrcFDEM as Src
from SimPEG import sp
####################################################
@@ -18,33 +19,33 @@ class Rx(SimPEG.Survey.BaseRx):
"""
knownRxTypes = {
'exr':['e', 'Ex', 'real'],
'eyr':['e', 'Ey', 'real'],
'ezr':['e', 'Ez', 'real'],
'exi':['e', 'Ex', 'imag'],
'eyi':['e', 'Ey', 'imag'],
'ezi':['e', 'Ez', 'imag'],
'exr':['e', 'x', 'real'],
'eyr':['e', 'y', 'real'],
'ezr':['e', 'z', 'real'],
'exi':['e', 'x', 'imag'],
'eyi':['e', 'y', 'imag'],
'ezi':['e', 'z', 'imag'],
'bxr':['b', 'Fx', 'real'],
'byr':['b', 'Fy', 'real'],
'bzr':['b', 'Fz', 'real'],
'bxi':['b', 'Fx', 'imag'],
'byi':['b', 'Fy', 'imag'],
'bzi':['b', 'Fz', 'imag'],
'bxr':['b', 'x', 'real'],
'byr':['b', 'y', 'real'],
'bzr':['b', 'z', 'real'],
'bxi':['b', 'x', 'imag'],
'byi':['b', 'y', 'imag'],
'bzi':['b', 'z', 'imag'],
'jxr':['j', 'Fx', 'real'],
'jyr':['j', 'Fy', 'real'],
'jzr':['j', 'Fz', 'real'],
'jxi':['j', 'Fx', 'imag'],
'jyi':['j', 'Fy', 'imag'],
'jzi':['j', 'Fz', 'imag'],
'jxr':['j', 'x', 'real'],
'jyr':['j', 'y', 'real'],
'jzr':['j', 'z', 'real'],
'jxi':['j', 'x', 'imag'],
'jyi':['j', 'y', 'imag'],
'jzi':['j', 'z', 'imag'],
'hxr':['h', 'Ex', 'real'],
'hyr':['h', 'Ey', 'real'],
'hzr':['h', 'Ez', 'real'],
'hxi':['h', 'Ex', 'imag'],
'hyi':['h', 'Ey', 'imag'],
'hzi':['h', 'Ez', 'imag'],
'hxr':['h', 'x', 'real'],
'hyr':['h', 'y', 'real'],
'hzr':['h', 'z', 'real'],
'hxi':['h', 'x', 'imag'],
'hyi':['h', 'y', 'imag'],
'hzi':['h', 'z', 'imag'],
}
radius = None
@@ -56,34 +57,37 @@ class Rx(SimPEG.Survey.BaseRx):
"""Field Type projection (e.g. e b ...)"""
return self.knownRxTypes[self.rxType][0]
@property
def projGLoc(self):
"""Grid Location projection (e.g. Ex Fy ...)"""
return self.knownRxTypes[self.rxType][1]
@property
def projComp(self):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][2]
def projectFields(self, src, mesh, u):
def projGLoc(self, u):
"""Grid Location projection (e.g. Ex Fy ...)"""
return u._GLoc(self.rxType[0]) + self.knownRxTypes[self.rxType][1]
def eval(self, src, mesh, u):
"""
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)
# projGLoc = u._GLoc(self.knownRxTypes[self.rxType][0])
# projGLoc += self.knownRxTypes[self.rxType][1]
P = self.getP(mesh, self.projGLoc(u))
u_part_complex = u[src, self.projField]
# get the real or imag component
real_or_imag = self.projComp
u_part = getattr(u_part_complex, real_or_imag)
return P*u_part
def projectFieldsDeriv(self, src, mesh, u, v, adjoint=False):
def evalDeriv(self, src, mesh, u, v, adjoint=False):
"""
Derivative of projected fields with respect to the inversion model times a vector.
@@ -94,7 +98,8 @@ class Rx(SimPEG.Survey.BaseRx):
:rtype: numpy.ndarray
:return: fields projected to recievers
"""
P = self.getP(mesh)
P = self.getP(mesh, self.projGLoc(u))
if not adjoint:
Pv_complex = P * v
@@ -171,7 +176,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 +186,9 @@ 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.')
SimPEG/EM/TDEM/BaseTDEM.py
+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)
SimPEG/EM/TDEM/SurveyTDEM.py
+7 -7
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:
@@ -168,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
SimPEG/EM/Utils/EMUtils.py
-33
View File
@@ -13,37 +13,4 @@ def k(freq, sigma, mu=mu_0, eps=epsilon_0):
beta = w * np.sqrt( mu*eps/2 * ( np.sqrt(1. + (sigma / (eps*w))**2 ) - 1) )
return alp - 1j*beta
# Constitutive relations
def e_from_j(prob,j):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MSigmaI = prob.MeSigmaI
elif eqLocs is 'EF':
MSigmaI = prob.MfRho
return MSigmaI*j
def j_from_e(prob,e):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MSigma = prob.MeSigma
elif eqLocs is 'EF':
MSigma = prob.MfRhoI
return MSigma*e
def b_from_h(prob,h):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MMu = prob.MfMuiI
elif eqLocs is 'EF':
MMu = prob.MeMu
return MMu*h
def h_from_b(prob,b):
eqLocs = prob._eqLocs
if eqLocs is 'FE':
MMuI = prob.MfMui
elif eqLocs is 'EF':
MMuI = prob.MeMuI
return MMuI*b
SimPEG/EM/Utils/__init__.py
+1 -4
View File
@@ -1,5 +1,2 @@
# import Sources
# import Ana
# import Solver
from EMUtils import omega, e_from_j, j_from_e, b_from_h, h_from_b
from EMUtils import omega, k
from AnalyticUtils import MagneticDipoleFields, MagneticDipoleVectorPotential, MagneticLoopVectorPotential
SimPEG/EM/Utils/testingUtils.py
+64 -13
View File
@@ -4,19 +4,28 @@ from SimPEG import EM
import sys
from scipy.constants import mu_0
def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
cs = 5.
ncx, ncy, ncz = 6, 6, 6
npad = 3
FLR = 1e-20 # "zero", so if residual below this --> pass regardless of order
CONDUCTIVITY = 1e1
MU = mu_0
freq = 5e-1
def getFDEMProblem(fdemType, comp, SrcList, freq, useMu=False, verbose=False):
cs = 10.
ncx, ncy, ncz = 0, 0, 0
npad = 8
hx = [(cs,npad,-1.3), (cs,ncx), (cs,npad,1.3)]
hy = [(cs,npad,-1.3), (cs,ncy), (cs,npad,1.3)]
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
mesh = Mesh.TensorMesh([hx,hy,hz],['C','C','C'])
mapping = Maps.ExpMap(mesh)
if useMu is True:
mapping = [('sigma', Maps.ExpMap(mesh)), ('mu', Maps.IdentityMap(mesh))]
else:
mapping = Maps.ExpMap(mesh)
x = np.array([np.linspace(-30,-15,3),np.linspace(15,30,3)]) #don't sample right by the source
XYZ = Utils.ndgrid(x,x,np.r_[0.])
x = np.array([np.linspace(-5.*cs,-2.*cs,3),np.linspace(5.*cs,2.*cs,3)]) + cs/4. #don't sample right by the source, slightly off alignment from either staggered grid
XYZ = Utils.ndgrid(x,x,np.linspace(-2.*cs,2.*cs,5))
Rx0 = EM.FDEM.Rx(XYZ, comp)
Src = []
@@ -32,15 +41,15 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
if fdemType is 'e' or fdemType is 'b':
S_m = np.zeros(mesh.nF)
S_e = np.zeros(mesh.nE)
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
elif fdemType is 'h' or fdemType is 'j':
S_m = np.zeros(mesh.nE)
S_e = np.zeros(mesh.nF)
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1.
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1.
S_m[Utils.closestPoints(mesh,[0.,0.,0.],'Ez') + np.sum(mesh.vnE[:1])] = 1e-3
S_e[Utils.closestPoints(mesh,[0.,0.,0.],'Fz') + np.sum(mesh.vnF[:1])] = 1e-3
Src.append(EM.FDEM.Src.RawVec([Rx0], freq, S_m, S_e))
if verbose:
@@ -70,6 +79,48 @@ def getFDEMProblem(fdemType, comp, SrcList, freq, verbose=False):
from pymatsolver import MumpsSolver
prb.Solver = MumpsSolver
except ImportError, e:
pass
prb.Solver = SolverLU
return prb
return prb
def crossCheckTest(SrcList, fdemType1, fdemType2, comp, addrandoms = False, useMu=False, TOL=1e-5, verbose=False):
l2norm = lambda r: np.sqrt(r.dot(r))
prb1 = getFDEMProblem(fdemType1, comp, SrcList, freq, useMu, verbose)
mesh = prb1.mesh
print 'Cross Checking Forward: %s, %s formulations - %s' % (fdemType1, fdemType2, comp)
logsig = np.log(np.ones(mesh.nC)*CONDUCTIVITY)
mu = np.ones(mesh.nC)*MU
if addrandoms is True:
logsig += np.random.randn(mesh.nC)*np.log(CONDUCTIVITY)*1e-1
mu += np.random.randn(mesh.nC)*MU*1e-1
if useMu is True:
m = np.r_[logsig, mu]
else:
m = logsig
survey1 = prb1.survey
d1 = survey1.dpred(m)
if verbose:
print ' Problem 1 solved'
prb2 = getFDEMProblem(fdemType2, comp, SrcList, freq, useMu, verbose)
survey2 = prb2.survey
d2 = survey2.dpred(m)
if verbose:
print ' Problem 2 solved'
r = d2-d1
l2r = l2norm(r)
tol = np.max([TOL*(10**int(np.log10(0.5* (l2norm(d1) + l2norm(d2)) ))),FLR])
print l2norm(d1), l2norm(d2), l2r , tol, l2r < tol
return l2r < tol
SimPEG/Examples/DC_Analytic_Dipole.py
+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)
SimPEG/Examples/DC_Forward_PseudoSection.py
+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()
SimPEG/Examples/EM_FDEM_1D_Inversion.py
+2 -2
View File
@@ -21,8 +21,8 @@ def run(plotIt=True):
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
sig_half = 2e-2
sig_air = 1e-8
sig_layer = 1e-2
SimPEG/Examples/EM_TDEM_1D_Inversion.py
+2 -2
View File
@@ -19,8 +19,8 @@ def run(plotIt=True):
active = mesh.vectorCCz<0.
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=-100.)
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
actMap = Maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
mapping = Maps.ExpMap(mesh) * Maps.SurjectVertical1D(mesh) * actMap
sig_half = 2e-3
sig_air = 1e-8
sig_layer = 1e-3
SimPEG/Examples/Inversion_Linear.py
+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)
SimPEG/Examples/MT_3D_Foward.py
+1 -1
View File
@@ -52,7 +52,7 @@ def run(plotIt=True, nFreq=1):
# Calculate the data
fields = problem.fields(sig)
dataVec = survey.projectFields(fields)
dataVec = survey.eval(fields)
# Make the data
mtData = MT.Data(survey,dataVec)
SimPEG/Examples/__init__.py
+3 -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
@@ -17,7 +19,7 @@ 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", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
__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 #####
SimPEG/FLOW/Richards/RichardsProblem.py
+10 -10
View File
@@ -8,7 +8,7 @@ class RichardsRx(Survey.BaseTimeRx):
knownRxTypes = ['saturation','pressureHead']
def projectFields(self, U, m, mapping, mesh, timeMesh):
def eval(self, U, m, mapping, mesh, timeMesh):
if self.rxType == 'pressureHead':
u = np.concatenate(U)
@@ -17,7 +17,7 @@ class RichardsRx(Survey.BaseTimeRx):
return self.getP(mesh, timeMesh) * u
def projectFieldsDeriv(self, U, m, mapping, mesh, timeMesh):
def evalDeriv(self, U, m, mapping, mesh, timeMesh):
P = self.getP(mesh, timeMesh)
if self.rxType == 'pressureHead':
@@ -57,13 +57,13 @@ class RichardsSurvey(Survey.BaseSurvey):
Where P is a projection of the fields onto the data space.
"""
if u is None: u = self.prob.fields(m)
return Utils.mkvc(self.projectFields(u, m))
return Utils.mkvc(self.eval(u, m))
@Utils.requires('prob')
def projectFields(self, U, m):
def eval(self, U, m):
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFields(U, m,
Ds[ii] = rx.eval(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -71,11 +71,11 @@ class RichardsSurvey(Survey.BaseSurvey):
return np.concatenate(Ds)
@Utils.requires('prob')
def projectFieldsDeriv(self, U, m):
def evalDeriv(self, U, m):
"""The Derivative with respect to the fields."""
Ds = range(len(self.rxList))
for ii, rx in enumerate(self.rxList):
Ds[ii] = rx.projectFieldsDeriv(U, m,
Ds[ii] = rx.evalDeriv(U, m,
self.prob.mapping,
self.prob.mesh,
self.prob.timeMesh)
@@ -251,7 +251,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
B = np.array(sp.vstack(Bs).todense())
Ainv = self.Solver(A, **self.solverOpts)
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
AinvB = Ainv * B
z = np.zeros((self.mesh.nC, B.shape[1]))
zAinvB = np.vstack((z, AinvB))
@@ -277,7 +277,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
Adiaginv = self.Solver(Adiag, **self.solverOpts)
JvC[ii] = Adiaginv * (B*v - Asub*JvC[ii-1])
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
return P * np.concatenate([np.zeros(self.mesh.nC)] + JvC)
@Utils.timeIt
@@ -285,7 +285,7 @@ class RichardsProblem(Problem.BaseTimeProblem):
if u is None:
u = self.field(m)
P = self.survey.projectFieldsDeriv(u, m)
P = self.survey.evalDeriv(u, m)
PTv = P.T*v
# This is done via backward substitution.
SimPEG/MT/BaseMT.py
+3 -3
View File
@@ -68,7 +68,7 @@ class BaseMTProblem(BaseFDEMProblem):
for rx in src.rxList:
# Get the projection derivative
# v should be of size 2*nE (for 2 polarizations)
PDeriv_u = lambda t: rx.projectFieldsDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
PDeriv_u = lambda t: rx.evalDeriv(src, self.mesh, u, t) # wrt u, we don't have have PDeriv wrt m
Jv[src, rx] = PDeriv_u(mkvc(du_dm))
dA_duI.clean()
# Return the vectorized sensitivities
@@ -106,9 +106,9 @@ class BaseMTProblem(BaseFDEMProblem):
u_src = u[src, :]
for rx in src.rxList:
# Get the adjoint projectFieldsDeriv
# Get the adjoint evalDeriv
# PTv needs to be nE,
PTv = rx.projectFieldsDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
PTv = rx.evalDeriv(src, self.mesh, u, mkvc(v[src, rx],2), adjoint=True) # wrt u, need possibility wrt m
# Get the
dA_duIT = ATinv * PTv
dA_dmT = self.getADeriv_m(freq, u_src, mkvc(dA_duIT), adjoint=True)
SimPEG/MT/SurveyMT.py
+5 -5
View File
@@ -59,7 +59,7 @@ class Rx(SimPEGsurvey.BaseRx):
"""Component projection (real/imag)"""
return self.knownRxTypes[self.rxType][1]
def projectFields(self, src, mesh, f):
def eval(self, src, mesh, f):
'''
Project the fields to natural source data.
@@ -139,7 +139,7 @@ class Rx(SimPEGsurvey.BaseRx):
# print f_part
return f_part
def projectFieldsDeriv(self, src, mesh, f, v, adjoint=False):
def evalDeriv(self, src, mesh, f, v, adjoint=False):
"""
The derivative of the projection wrt u
@@ -427,15 +427,15 @@ class Survey(SimPEGsurvey.BaseSurvey):
assert freq in self._freqDict, "The requested frequency is not in this survey."
return self._freqDict[freq]
def projectFields(self, u):
def eval(self, u):
data = Data(self)
for src in self.srcList:
sys.stdout.flush()
for rx in src.rxList:
data[src, rx] = rx.projectFields(src, self.mesh, u)
data[src, rx] = rx.eval(src, self.mesh, u)
return data
def projectFieldsDeriv(self, u):
def evalDeriv(self, u):
raise Exception('Use Transmitters to project fields deriv.')
#################
SimPEG/Maps.py
+41 -8
View File
@@ -4,6 +4,7 @@ from Tests import checkDerivative
from PropMaps import PropMap, Property
from numpy.polynomial import polynomial
from scipy.interpolate import UnivariateSpline
import warnings
class IdentityMap(object):
"""
@@ -296,11 +297,11 @@ class LogMap(IdentityMap):
def inverse(self, m):
return np.exp(Utils.mkvc(m))
class FullMap(IdentityMap):
class SurjectFull(IdentityMap):
"""
FullMap
SurjectFull
Given a scalar, the FullMap maps the value to the
Given a scalar, the SurjectFull maps the value to the
full model space.
"""
@@ -327,9 +328,15 @@ class FullMap(IdentityMap):
"""
return np.ones([self.mesh.nC,1])
class FullMap(SurjectFull):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`FullMap` is deprecated and will be removed in future versions. Use `SurjectFull` instead",
FutureWarning)
SurjectFull.__init__(self,mesh,**kwargs)
class Vertical1DMap(IdentityMap):
"""Vertical1DMap
class SurjectVertical1D(IdentityMap):
"""SurjectVertical1DMap
Given a 1D vector through the last dimension
of the mesh, this will extend to the full
@@ -369,8 +376,14 @@ class Vertical1DMap(IdentityMap):
), shape=(repNum, 1))
return sp.kron(sp.identity(self.nP), repVec)
class Vertical1DMap(SurjectVertical1D):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Vertical1DMap` is deprecated and will be removed in future versions. Use `SurjectVertical1D` instead",
FutureWarning)
SurjectVertical1D.__init__(self,mesh,**kwargs)
class Map2Dto3D(IdentityMap):
class Surject2Dto3D(IdentityMap):
"""Map2Dto3D
Given a 2D vector, this will extend to the full
@@ -425,6 +438,13 @@ class Map2Dto3D(IdentityMap):
), shape=(nC, nP))
return P
class Map2Dto3D(Surject2Dto3D):
def __init__(self,mesh,**kwargs):
warnings.warn(
"`Map2Dto3D` is deprecated and will be removed in future versions. Use `Surject2Dto3D` instead",
FutureWarning)
Surject2Dto3D.__init__(self,mesh,**kwargs)
class Mesh2Mesh(IdentityMap):
"""
Takes a model on one mesh are translates it to another mesh.
@@ -458,7 +478,7 @@ class Mesh2Mesh(IdentityMap):
return self.P
class ActiveCells(IdentityMap):
class InjectActiveCells(IdentityMap):
"""
Active model parameters.
@@ -506,7 +526,14 @@ class ActiveCells(IdentityMap):
def deriv(self, m):
return self.P
class ActiveCellsTopo(IdentityMap):
class ActiveCells(InjectActiveCells):
def __init__(self, mesh, indActive, valInactive, nC=None):
warnings.warn(
"`ActiveCells` is deprecated and will be removed in future versions. Use `InjectActiveCells` instead",
FutureWarning)
InjectActiveCells.__init__(self, mesh, indActive, valInactive, nC)
class InjectActiveCellsTopo(IdentityMap):
"""
Active model parameters. Extend for cells on topography to air cell (only works for tensor mesh)
@@ -577,6 +604,12 @@ class ActiveCellsTopo(IdentityMap):
def deriv(self, m):
return self.P
class ActiveCellsTopo(InjectActiveCellsTopo):
def __init__(self, mesh, indActive, valInactive, nC=None):
warnings.warn(
"`ActiveCellsTopo` is deprecated and will be removed in future versions. Use `InjectActiveCellsTopo` instead",
FutureWarning)
InjectActiveCellsTopo.__init__(self, mesh, indActive, valInactive, nC)
class Weighting(IdentityMap):
"""
SimPEG/Mesh/DiffOperators.py
-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
SimPEG/Mesh/TensorMesh.py
+13
View File
@@ -234,6 +234,9 @@ class BaseTensorMesh(BaseMesh):
'Fz' -> z-component of field defined on faces
'N' -> scalar field defined on nodes
'CC' -> scalar field defined on cell centers
'CCVx' -> x-component of vector field defined on cell centers
'CCVy' -> y-component of vector field defined on cell centers
'CCVz' -> z-component of vector field defined on cell centers
"""
if self._meshType == 'CYL' and self.isSymmetric and locType in ['Ex','Ez','Fy']:
raise Exception('Symmetric CylMesh does not support %s interpolation, as this variable does not exist.' % locType)
@@ -257,6 +260,16 @@ class BaseTensorMesh(BaseMesh):
Q = sp.hstack(components)
elif locType in ['CC', 'N']:
Q = Utils.interpmat(loc, *self.getTensor(locType))
elif locType in ['CCVx', 'CCVy', 'CCVz']:
Q = Utils.interpmat(loc, *self.getTensor('CC'))
Z = Utils.spzeros(loc.shape[0],self.nC)
if locType == 'CCVx':
Q = sp.hstack([Q,Z,Z])
elif locType == 'CCVy':
Q = sp.hstack([Z,Q,Z])
elif locType == 'CCVz':
Q = sp.hstack([Z,Z,Q])
else:
raise NotImplementedError('getInterpolationMat: locType=='+locType+' and mesh.dim=='+str(self.dim))
SimPEG/Problem.py
+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)
SimPEG/Regularization.py
+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
SimPEG/Survey.py
+20 -10
View File
@@ -1,6 +1,5 @@
import Utils, numpy as np, scipy.sparse as sp, uuid
class BaseRx(object):
"""SimPEG Receiver Object"""
@@ -35,7 +34,7 @@ class BaseRx(object):
"""Number of data in the receiver."""
return self.locs.shape[0]
def getP(self, mesh):
def getP(self, mesh, projGLoc=None):
"""
Returns the projection matrices as a
list for all components collected by
@@ -48,7 +47,10 @@ class BaseRx(object):
if mesh in self._Ps:
return self._Ps[mesh]
P = mesh.getInterpolationMat(self.locs, self.projGLoc)
if projGLoc is None:
projGLoc = self.projGLoc
P = mesh.getInterpolationMat(self.locs, projGLoc)
if self.storeProjections:
self._Ps[mesh] = P
return P
@@ -307,12 +309,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 +322,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 +334,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 +377,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]
SimPEG/Utils/ModelBuilder.py
+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.