mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 18:25:42 +08:00
Merge branch 'dcip/spectralIP' of https://github.com/simpeg/simpeg into dcip/dev
Merge spectral IP stuff, and incorporate Lindsey's comments
This commit is contained in:
seogi_macbook
+13
-13
@@ -2,7 +2,7 @@ import numpy as np
|
||||
from scipy.constants import mu_0, pi
|
||||
from scipy import special
|
||||
|
||||
def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
|
||||
def DCAnalyticHalf(txloc, rxlocs, sigma, earth_type="wholespace"):
|
||||
"""
|
||||
Analytic solution for electric potential from a postive pole
|
||||
|
||||
@@ -15,7 +15,7 @@ def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
|
||||
N: xyz locations of N (-) electrode (np.c_[xnlocs, ynlocs, znlocs])
|
||||
|
||||
sigma = conductivity (either float or complex)
|
||||
flag = "wholsespace" or "halfspace"
|
||||
earth_type = "wholsespace" or "halfspace"
|
||||
|
||||
"""
|
||||
M = rxlocs[0]
|
||||
@@ -28,7 +28,7 @@ def DCAnalyticHalf(txloc, rxlocs, sigma, flag="wholespace"):
|
||||
phiN = 1./(4*np.pi*rN*sigma)
|
||||
phi = phiM - phiN
|
||||
|
||||
if flag == "halfspace":
|
||||
if earth_type == "halfspace":
|
||||
phi *= 2
|
||||
|
||||
return phi
|
||||
@@ -37,9 +37,9 @@ deg2rad = lambda deg: deg/180.*np.pi
|
||||
rad2deg = lambda rad: rad*180./np.pi
|
||||
|
||||
def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
|
||||
flag = "sec", order=12, halfspace=False):
|
||||
field_type = "secondary", order=12, halfspace=False):
|
||||
# def DCSpherePointCurrent(txloc, rxloc, xc, radius, rho, rho1, \
|
||||
# flag = "sec", order=12):
|
||||
# field_type = "secondary", order=12):
|
||||
"""
|
||||
|
||||
Parameters:
|
||||
@@ -51,11 +51,11 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
|
||||
radius (float): radius of the sphere (m)
|
||||
rho (float) : resistivity of the background (ohm-m)
|
||||
rho1 (float) : resistivity of the sphere
|
||||
flag (string) : "sec", "total", "prim"
|
||||
(default="sec")
|
||||
"sec": secondary potential only due to sphere
|
||||
"prim": primary potential from the point source
|
||||
"total": "sec"+"prim"
|
||||
field_type (string) : "secondary", "total", "primary"
|
||||
(default="secondary")
|
||||
"secondary": secondary potential only due to sphere
|
||||
"primary": primary potential from the point source
|
||||
"total": "secondary"+"primary"
|
||||
order (float) : maximum order of Legendre polynomial
|
||||
(default=12)
|
||||
|
||||
@@ -86,7 +86,7 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
|
||||
# primary potential in a whole space
|
||||
prim = rho*1./(4*np.pi*R)
|
||||
|
||||
if flag =="prim":
|
||||
if field_type =="primary":
|
||||
return prim
|
||||
|
||||
sphind = r < radius
|
||||
@@ -105,9 +105,9 @@ def DCAnalyticSphere(txloc, rxloc, xc, radius, sigma, sigma1, \
|
||||
else:
|
||||
scale = 1
|
||||
|
||||
if flag == "sec":
|
||||
if field_type == "secondary":
|
||||
return scale*(out-prim)
|
||||
elif flag == "total":
|
||||
elif field_type == "total":
|
||||
return scale*out
|
||||
|
||||
def AnBnfun(n, radius, x0, rho, rho1, I=1.):
|
||||
|
||||
@@ -77,7 +77,7 @@ class BaseDCProblem(BaseEMProblem):
|
||||
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
Jtv += df_dmT + du_dmT
|
||||
Jtv += (df_dmT + du_dmT).astype(float)
|
||||
|
||||
return Utils.mkvc(Jtv)
|
||||
|
||||
@@ -122,13 +122,12 @@ class Problem3D_CC(BaseDCProblem):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = D MfRhoI G
|
||||
|
||||
"""
|
||||
|
||||
D = self.Div
|
||||
G = self.Grad
|
||||
# TODO: this won't work for full anisotropy
|
||||
MfRhoI = self.MfRhoI
|
||||
A = D * MfRhoI * G
|
||||
|
||||
@@ -144,13 +143,8 @@ class Problem3D_CC(BaseDCProblem):
|
||||
MfRhoIDeriv = self.MfRhoIDeriv
|
||||
|
||||
if adjoint:
|
||||
# if self._makeASymmetric is True:
|
||||
# v = V * v
|
||||
return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
|
||||
|
||||
# I think we should deprecate this for DC problem.
|
||||
# if self._makeASymmetric is True:
|
||||
# return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
|
||||
return D * (MfRhoIDeriv( G * u ) * v)
|
||||
|
||||
def getRHS(self):
|
||||
@@ -162,10 +156,6 @@ class Problem3D_CC(BaseDCProblem):
|
||||
|
||||
RHS = self.getSourceTerm()
|
||||
|
||||
# I think we should deprecate this for DC problem.
|
||||
# if self._makeASymmetric is True:
|
||||
# return self.Vol.T * RHS
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv(self, src, v, adjoint=False):
|
||||
@@ -255,11 +245,10 @@ class Problem3D_N(BaseDCProblem):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = G.T MeSigma G
|
||||
|
||||
"""
|
||||
|
||||
# TODO: this won't work for full anisotropy
|
||||
MeSigma = self.MeSigma
|
||||
Grad = self.mesh.nodalGrad
|
||||
A = Grad.T * MeSigma * Grad
|
||||
|
||||
@@ -161,14 +161,13 @@ class Problem2D_CC(BaseDCProblem_2D):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = D MfRhoI G
|
||||
|
||||
"""
|
||||
|
||||
D = self.Div
|
||||
G = self.Grad
|
||||
vol = self.mesh.vol
|
||||
# TODO: this won't work for full anisotropy
|
||||
MfRhoI = self.MfRhoI
|
||||
# Get resistivity rho
|
||||
rho = self.curModel.rho
|
||||
@@ -304,11 +303,10 @@ class Problem2D_N(BaseDCProblem_2D):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = D MfRhoI G
|
||||
|
||||
"""
|
||||
|
||||
# TODO: this won't work for full anisotropy
|
||||
MeSigma = self.MeSigma
|
||||
MnSigma = self.MnSigma
|
||||
Grad = self.mesh.nodalGrad
|
||||
|
||||
@@ -89,7 +89,7 @@ class BaseIPProblem(BaseEMProblem):
|
||||
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
Jtv += df_dmT + du_dmT
|
||||
Jtv += (df_dmT + du_dmT).astype(float)
|
||||
# Conductivity ((d u / d log sigma).T)
|
||||
if self._formulation is 'EB':
|
||||
return -Utils.mkvc(Jtv)
|
||||
@@ -180,13 +180,12 @@ class Problem3D_CC(BaseIPProblem):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = D MfRhoI G
|
||||
|
||||
"""
|
||||
|
||||
D = self.Div
|
||||
G = self.Grad
|
||||
# TODO: this won't work for full anisotropy
|
||||
MfRhoI = self.MfRhoI
|
||||
A = D * MfRhoI * G
|
||||
|
||||
@@ -313,11 +312,10 @@ class Problem3D_N(BaseIPProblem):
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI D^\\top V
|
||||
A = G.T MeSigma G
|
||||
|
||||
"""
|
||||
|
||||
# TODO: this won't work for full anisotropy
|
||||
MeSigma = self.MeSigma
|
||||
Grad = self.mesh.nodalGrad
|
||||
A = Grad.T * MeSigma * Grad
|
||||
|
||||
@@ -0,0 +1,445 @@
|
||||
from SimPEG import Problem, Utils, Maps, Mesh
|
||||
from SimPEG.EM.Base import BaseEMProblem
|
||||
from SimPEG.EM.Static.DC.FieldsDC import Fields, Fields_CC, Fields_N
|
||||
from SimPEG.Utils import sdiag
|
||||
import numpy as np
|
||||
from SimPEG.Utils import Zero
|
||||
from SimPEG.EM.Static.DC import getxBCyBC_CC
|
||||
from SurveySIP import Survey, Data
|
||||
|
||||
class ColeColePropMap(Maps.PropMap):
|
||||
"""
|
||||
Property Map for EM Problems. The electrical conductivity (\\(\\sigma\\)) is the default inversion property, and the default value of the magnetic permeability is that of free space (\\(\\mu = 4\\pi\\times 10^{-7} \\) H/m)
|
||||
"""
|
||||
|
||||
eta = Maps.Property("Electrical Conductivity", defaultInvProp=True)
|
||||
tau = Maps.Property("Electrical Conductivity", defaultVal=0.1, propertyLink=('taui', Maps.ReciprocalMap))
|
||||
taui = Maps.Property("Electrical Conductivity", defaultVal=1., propertyLink=('tau', Maps.ReciprocalMap))
|
||||
c = Maps.Property("Electrical Conductivity", defaultVal=1.)
|
||||
|
||||
|
||||
class BaseSIPProblem(BaseEMProblem):
|
||||
|
||||
surveyPair = Survey
|
||||
fieldsPair = Fields
|
||||
dataPair = Data
|
||||
PropMap = ColeColePropMap
|
||||
Ainv = None
|
||||
sigma = None
|
||||
rho = None
|
||||
f = None
|
||||
Ainv = None
|
||||
|
||||
def DebyeTime(self, t):
|
||||
peta = self.curModel.eta*np.exp(-self.curModel.taui*t)
|
||||
return peta
|
||||
|
||||
def EtaDeriv(self, t, v, adjoint=False):
|
||||
v = np.array(v, dtype=float)
|
||||
if adjoint:
|
||||
return self.curModel.etaDeriv.T * (np.exp(-self.curModel.taui*t)*v)
|
||||
else:
|
||||
return np.exp(-self.curModel.taui*t) * (self.curModel.etaDeriv*v)
|
||||
|
||||
|
||||
def TauiDeriv(self, t, v, adjoint=False):
|
||||
v = np.array(v, dtype=float)
|
||||
if adjoint:
|
||||
return -self.curModel.tauiDeriv.T * (self.curModel.eta*t*np.exp(-self.curModel.taui*t)*v)
|
||||
else:
|
||||
return -self.curModel.eta*t*np.exp(-self.curModel.taui*t) * (self.curModel.tauiDeriv*v)
|
||||
|
||||
def fields(self, m):
|
||||
self.curModel = m
|
||||
if self.f is None:
|
||||
self.f = self.fieldsPair(self.mesh, self.survey)
|
||||
if self.Ainv == None:
|
||||
A = self.getA()
|
||||
self.Ainv = self.Solver(A, **self.solverOpts)
|
||||
RHS = self.getRHS()
|
||||
u = self.Ainv * RHS
|
||||
Srcs = self.survey.srcList
|
||||
self.f[Srcs, self._solutionType] = u
|
||||
return self.f
|
||||
|
||||
def forward(self, m, f=None):
|
||||
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
Jv = self.dataPair(self.survey) #same size as the data
|
||||
# A = self.getA()
|
||||
JvAll = []
|
||||
for tind in range(len(self.survey.times)):
|
||||
#Pseudo-chareability
|
||||
t = self.survey.times[tind]
|
||||
v = self.DebyeTime(t)
|
||||
for src in self.survey.srcList:
|
||||
u_src = f[src, self._solutionType] # solution vector
|
||||
dA_dm_v = self.getADeriv(u_src, v)
|
||||
dRHS_dm_v = self.getRHSDeriv(src, v)
|
||||
du_dm_v = self.Ainv * ( - dA_dm_v + dRHS_dm_v )
|
||||
for rx in src.rxList:
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dm_v = df_dmFun(src, du_dm_v, v, adjoint=False)
|
||||
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v)
|
||||
|
||||
# Conductivity (d u / d log sigma)
|
||||
if self._formulation is 'EB':
|
||||
return -Utils.mkvc(Jv)
|
||||
# Resistivity (d u / d log rho)
|
||||
if self._formulation is 'HJ':
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
def Jvec(self, m, v, f=None):
|
||||
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
Jv = self.dataPair(self.survey) #same size as the data
|
||||
# A = self.getA()
|
||||
JvAll = []
|
||||
#Assume only eta and tau (eta first then tau)
|
||||
# v = [2*Mx1]
|
||||
v = v.reshape((int(v.size/2), 2), order='F')
|
||||
|
||||
for tind in range(len(self.survey.times)):
|
||||
t = self.survey.times[tind]
|
||||
v0 = self.EtaDeriv(t, v[:,0])
|
||||
v1 = self.TauiDeriv(t, v[:,1])
|
||||
for src in self.survey.srcList:
|
||||
u_src = f[src, self._solutionType] # solution vector
|
||||
dA_dm_v0 = self.getADeriv(u_src, v0)
|
||||
dRHS_dm_v0 = self.getRHSDeriv(src, v0)
|
||||
du_dm_v0 = self.Ainv * ( - dA_dm_v0 + dRHS_dm_v0 )
|
||||
dA_dm_v1 = self.getADeriv(u_src, v1)
|
||||
dRHS_dm_v1 = self.getRHSDeriv(src, v1)
|
||||
du_dm_v1 = self.Ainv * ( - dA_dm_v1 + dRHS_dm_v1 )
|
||||
for rx in src.rxList:
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
df_dmFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_dm_v0 = df_dmFun(src, du_dm_v0, v0, adjoint=False)
|
||||
df_dm_v1 = df_dmFun(src, du_dm_v1, v1, adjoint=False)
|
||||
Jv[src, rx, t] = rx.evalDeriv(src, self.mesh, f, df_dm_v0)
|
||||
Jv[src, rx, t] += rx.evalDeriv(src, self.mesh, f, df_dm_v1)
|
||||
# Conductivity (d u / d log sigma)
|
||||
if self._formulation is 'EB':
|
||||
return -Jv.tovec()
|
||||
# Resistivity (d u / d log rho)
|
||||
if self._formulation is 'HJ':
|
||||
return Jv.tovec()
|
||||
|
||||
def Jtvec(self, m, v, f=None):
|
||||
if f is None:
|
||||
f = self.fields(m)
|
||||
|
||||
self.curModel = m
|
||||
|
||||
# Ensure v is a data object.
|
||||
if not isinstance(v, self.dataPair):
|
||||
v = self.dataPair(self.survey, v)
|
||||
|
||||
Jtv= np.zeros(m.size)
|
||||
for tind in range(len(self.survey.times)):
|
||||
t = self.survey.times[tind]
|
||||
for src in self.survey.srcList:
|
||||
u_src = f[src, self._solutionType]
|
||||
for rx in src.rxList:
|
||||
timeindex = rx.getTimeP(self.survey.times)
|
||||
if timeindex[tind]:
|
||||
PTv = rx.evalDeriv(src, self.mesh, f, v[src, rx, t], adjoint=True) # wrt f, need possibility wrt m
|
||||
df_duTFun = getattr(f, '_%sDeriv'%rx.projField, None)
|
||||
df_duT, df_dmT = df_duTFun(src, None, PTv, adjoint=True)
|
||||
ATinvdf_duT = self.Ainv * df_duT
|
||||
dA_dmT = self.getADeriv(u_src, ATinvdf_duT, adjoint=True)
|
||||
dRHS_dmT = self.getRHSDeriv(src, ATinvdf_duT, adjoint=True)
|
||||
du_dmT = -dA_dmT + dRHS_dmT
|
||||
Jtv += np.r_[self.EtaDeriv(self.survey.times[tind], du_dmT, adjoint=True), self.TauiDeriv(self.survey.times[tind], du_dmT, adjoint=True)]
|
||||
|
||||
# Conductivity ((d u / d log sigma).T)
|
||||
if self._formulation is 'EB':
|
||||
return -Jtv
|
||||
# Conductivity ((d u / d log rho).T)
|
||||
if self._formulation is 'HJ':
|
||||
return Jtv
|
||||
|
||||
def getSourceTerm(self):
|
||||
"""
|
||||
takes concept of source and turns it into a matrix
|
||||
"""
|
||||
"""
|
||||
Evaluates the sources, and puts them in matrix form
|
||||
|
||||
:rtype: (numpy.ndarray, numpy.ndarray)
|
||||
:return: q (nC or nN, nSrc)
|
||||
"""
|
||||
|
||||
Srcs = self.survey.srcList
|
||||
|
||||
if self._formulation is 'EB':
|
||||
n = self.mesh.nN
|
||||
# return NotImplementedError
|
||||
|
||||
elif self._formulation is 'HJ':
|
||||
n = self.mesh.nC
|
||||
|
||||
q = np.zeros((n, len(Srcs)))
|
||||
|
||||
for i, src in enumerate(Srcs):
|
||||
q[:,i] = src.eval(self)
|
||||
return q
|
||||
|
||||
@property
|
||||
def deleteTheseOnModelUpdate(self):
|
||||
toDelete = []
|
||||
return toDelete
|
||||
|
||||
# assume log rho or log cond
|
||||
@property
|
||||
def MeSigma(self):
|
||||
"""
|
||||
Edge inner product matrix for \\(\\sigma\\). Used in the E-B formulation
|
||||
"""
|
||||
if getattr(self, '_MeSigma', None) is None:
|
||||
self._MeSigma = self.mesh.getEdgeInnerProduct(self.sigma)
|
||||
return self._MeSigma
|
||||
|
||||
@property
|
||||
def MfRhoI(self):
|
||||
"""
|
||||
Inverse of :code:`MfRho`
|
||||
"""
|
||||
if getattr(self, '_MfRhoI', None) is None:
|
||||
self._MfRhoI = self.mesh.getFaceInnerProduct(self.rho, invMat=True)
|
||||
return self._MfRhoI
|
||||
|
||||
def MfRhoIDeriv(self,u):
|
||||
"""
|
||||
Derivative of :code:`MfRhoI` with respect to the model.
|
||||
"""
|
||||
|
||||
dMfRhoI_dI = -self.MfRhoI**2
|
||||
dMf_drho = self.mesh.getFaceInnerProductDeriv(self.rho)(u)
|
||||
drho_dlogrho = Utils.sdiag(self.rho)
|
||||
return dMfRhoI_dI * ( dMf_drho * ( drho_dlogrho))
|
||||
|
||||
# TODO: This should take a vector
|
||||
def MeSigmaDeriv(self, u):
|
||||
"""
|
||||
Derivative of MeSigma with respect to the model
|
||||
"""
|
||||
dsigma_dlogsigma = Utils.sdiag(self.sigma)
|
||||
return self.mesh.getEdgeInnerProductDeriv(self.sigma)(u) * dsigma_dlogsigma
|
||||
|
||||
class Problem3D_CC(BaseSIPProblem):
|
||||
|
||||
_solutionType = 'phiSolution'
|
||||
_formulation = 'HJ' # CC potentials means J is on faces
|
||||
fieldsPair = Fields_CC
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseSIPProblem.__init__(self, mesh, **kwargs)
|
||||
self.setBC()
|
||||
|
||||
def getA(self):
|
||||
"""
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = D MfRhoI G
|
||||
|
||||
"""
|
||||
|
||||
D = self.Div
|
||||
G = self.Grad
|
||||
# TODO: this won't work for full anisotropy
|
||||
MfRhoI = self.MfRhoI
|
||||
A = D * MfRhoI * G
|
||||
|
||||
# I think we should deprecate this for DC problem.
|
||||
# if self._makeASymmetric is True:
|
||||
# return V.T * A
|
||||
return A
|
||||
|
||||
def getADeriv(self, u, v, adjoint= False):
|
||||
|
||||
D = self.Div
|
||||
G = self.Grad
|
||||
MfRhoIDeriv = self.MfRhoIDeriv
|
||||
|
||||
if adjoint:
|
||||
# if self._makeASymmetric is True:
|
||||
# v = V * v
|
||||
return(MfRhoIDeriv( G * u ).T) * ( D.T * v)
|
||||
|
||||
# I think we should deprecate this for DC problem.
|
||||
# if self._makeASymmetric is True:
|
||||
# return V.T * ( D * ( MfRhoIDeriv( D.T * ( V * u ) ) * v ) )
|
||||
return D * (MfRhoIDeriv( G * u ) * v)
|
||||
|
||||
def getRHS(self):
|
||||
"""
|
||||
RHS for the DC problem
|
||||
|
||||
q
|
||||
"""
|
||||
|
||||
RHS = self.getSourceTerm()
|
||||
|
||||
# I think we should deprecate this for DC problem.
|
||||
# if self._makeASymmetric is True:
|
||||
# return self.Vol.T * RHS
|
||||
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
"""
|
||||
# TODO: add qDeriv for RHS depending on m
|
||||
# qDeriv = src.evalDeriv(self, adjoint=adjoint)
|
||||
# return qDeriv
|
||||
return Zero()
|
||||
|
||||
def setBC(self):
|
||||
if self.mesh.dim==3:
|
||||
fxm,fxp,fym,fyp,fzm,fzp = self.mesh.faceBoundaryInd
|
||||
gBFxm = self.mesh.gridFx[fxm,:]
|
||||
gBFxp = self.mesh.gridFx[fxp,:]
|
||||
gBFym = self.mesh.gridFy[fym,:]
|
||||
gBFyp = self.mesh.gridFy[fyp,:]
|
||||
gBFzm = self.mesh.gridFz[fzm,:]
|
||||
gBFzp = self.mesh.gridFz[fzp,:]
|
||||
|
||||
# Setup Mixed B.C (alpha, beta, gamma)
|
||||
temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
temp_zm, temp_zp = np.ones_like(gBFzm[:,2]), np.ones_like(gBFzp[:,2])
|
||||
|
||||
alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
|
||||
alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
|
||||
alpha_zm, alpha_zp = temp_zm*0., temp_zp*0.
|
||||
|
||||
beta_xm, beta_xp = temp_xm, temp_xp
|
||||
beta_ym, beta_yp = temp_ym, temp_yp
|
||||
beta_zm, beta_zp = temp_zm, temp_zp
|
||||
|
||||
gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
|
||||
gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
|
||||
gamma_zm, gamma_zp = temp_zm*0., temp_zp*0.
|
||||
|
||||
alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp, alpha_zm, alpha_zp]
|
||||
beta = [beta_xm, beta_xp, beta_ym, beta_yp, beta_zm, beta_zp]
|
||||
gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp, gamma_zm, gamma_zp]
|
||||
|
||||
elif self.mesh.dim==2:
|
||||
|
||||
fxm,fxp,fym,fyp = self.mesh.faceBoundaryInd
|
||||
gBFxm = self.mesh.gridFx[fxm,:]
|
||||
gBFxp = self.mesh.gridFx[fxp,:]
|
||||
gBFym = self.mesh.gridFy[fym,:]
|
||||
gBFyp = self.mesh.gridFy[fyp,:]
|
||||
|
||||
# Setup Mixed B.C (alpha, beta, gamma)
|
||||
temp_xm, temp_xp = np.ones_like(gBFxm[:,0]), np.ones_like(gBFxp[:,0])
|
||||
temp_ym, temp_yp = np.ones_like(gBFym[:,1]), np.ones_like(gBFyp[:,1])
|
||||
|
||||
alpha_xm, alpha_xp = temp_xm*0., temp_xp*0.
|
||||
alpha_ym, alpha_yp = temp_ym*0., temp_yp*0.
|
||||
|
||||
beta_xm, beta_xp = temp_xm, temp_xp
|
||||
beta_ym, beta_yp = temp_ym, temp_yp
|
||||
|
||||
gamma_xm, gamma_xp = temp_xm*0., temp_xp*0.
|
||||
gamma_ym, gamma_yp = temp_ym*0., temp_yp*0.
|
||||
|
||||
alpha = [alpha_xm, alpha_xp, alpha_ym, alpha_yp]
|
||||
beta = [beta_xm, beta_xp, beta_ym, beta_yp]
|
||||
gamma = [gamma_xm, gamma_xp, gamma_ym, gamma_yp]
|
||||
|
||||
x_BC, y_BC = getxBCyBC_CC(self.mesh, alpha, beta, gamma)
|
||||
V = self.Vol
|
||||
self.Div = V * self.mesh.faceDiv
|
||||
P_BC, B = self.mesh.getBCProjWF_simple()
|
||||
M = B*self.mesh.aveCC2F
|
||||
self.Grad = self.Div.T - P_BC*Utils.sdiag(y_BC)*M
|
||||
|
||||
|
||||
class Problem3D_N(BaseSIPProblem):
|
||||
|
||||
_solutionType = 'phiSolution'
|
||||
_formulation = 'EB' # N potentials means B is on faces
|
||||
fieldsPair = Fields_N
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
BaseSIPProblem.__init__(self, mesh, **kwargs)
|
||||
|
||||
def getA(self):
|
||||
"""
|
||||
|
||||
Make the A matrix for the cell centered DC resistivity problem
|
||||
|
||||
A = G.T MeSigma G
|
||||
|
||||
"""
|
||||
|
||||
# TODO: this won't work for full anisotropy
|
||||
MeSigma = self.MeSigma
|
||||
Grad = self.mesh.nodalGrad
|
||||
A = Grad.T * MeSigma * Grad
|
||||
|
||||
# Handling Null space of A
|
||||
A[0,0] = A[0,0] + 1.
|
||||
|
||||
return A
|
||||
|
||||
def getADeriv(self, u, v, adjoint=False):
|
||||
"""
|
||||
|
||||
Product of the derivative of our system matrix with respect to the model and a vector
|
||||
|
||||
"""
|
||||
MeSigma = self.MeSigma
|
||||
Grad = self.mesh.nodalGrad
|
||||
if not adjoint:
|
||||
return Grad.T*(self.MeSigmaDeriv(Grad*u)*v)
|
||||
elif adjoint:
|
||||
return self.MeSigmaDeriv(Grad*u).T * (Grad*v)
|
||||
|
||||
|
||||
def getRHS(self):
|
||||
"""
|
||||
RHS for the DC problem
|
||||
|
||||
q
|
||||
"""
|
||||
|
||||
RHS = self.getSourceTerm()
|
||||
return RHS
|
||||
|
||||
def getRHSDeriv(self, src, v, adjoint=False):
|
||||
"""
|
||||
Derivative of the right hand side with respect to the model
|
||||
"""
|
||||
# TODO: add qDeriv for RHS depending on m
|
||||
# qDeriv = src.evalDeriv(self, adjoint=adjoint)
|
||||
# return qDeriv
|
||||
return Zero()
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
|
||||
cs = 12.5
|
||||
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],x0="CCN")
|
||||
sigma = np.ones(mesh.nC)
|
||||
prob = BaseSIPProblem(mesh, sigma=sigma)
|
||||
|
||||
|
||||
@@ -0,0 +1,204 @@
|
||||
from SimPEG import Utils, Maps, Mesh, sp, np
|
||||
from SimPEG.Regularization import BaseRegularization, Simple
|
||||
|
||||
class MultiRegularization(Simple):
|
||||
"""
|
||||
**MultiRegularization Class**
|
||||
|
||||
This is used to regularize the model space
|
||||
having multiple models [m1, m2, m3, ...] ::
|
||||
|
||||
reg = Regularization(mesh)
|
||||
|
||||
"""
|
||||
nModels = None # Number of models
|
||||
ratios = None # Ratio for different models
|
||||
crossgrad = False # Use cross gradient or not
|
||||
betacross = 1.
|
||||
wx = []
|
||||
wy = []
|
||||
wz = []
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
if self.nModels == None:
|
||||
raise Exception("Put nModels as a initial input!")
|
||||
if self.ratios == None:
|
||||
self.ratios = [1. for imodel in range(self.nModels)]
|
||||
|
||||
@property
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self,'_Wsmall', None) is None:
|
||||
vecs = []
|
||||
for imodel in range(self.nModels):
|
||||
vecs.append((self.regmesh.vol*self.alpha_s*self.wght*self.ratios[imodel])**0.5)
|
||||
self._Wsmall = Utils.sdiag(np.hstack(vecs))
|
||||
return self._Wsmall
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
if getattr(self, '_Wx', None) is None:
|
||||
mats = []
|
||||
for imodel in range(self.nModels):
|
||||
self.wx.append(Utils.sdiag((self.regmesh.aveCC2Fx * self.regmesh.vol*self.alpha_x*self.ratios[imodel]*(self.regmesh.aveCC2Fx*self.wght))**0.5))
|
||||
mats.append(self.wx[imodel]*self.regmesh.cellDiffxStencil)
|
||||
self._Wx = sp.block_diag(mats)
|
||||
return self._Wx
|
||||
|
||||
@property
|
||||
def Wy(self):
|
||||
"""Regularization matrix Wy"""
|
||||
if getattr(self, '_Wy', None) is None:
|
||||
mats = []
|
||||
for imodel in range(self.nModels):
|
||||
self.wy.append(Utils.sdiag((self.regmesh.aveCC2Fy * self.regmesh.vol*self.alpha_y*self.ratios[imodel]*(self.regmesh.aveCC2Fy*self.wght))**0.5))
|
||||
mats.append(self.wy[imodel]*self.regmesh.cellDiffyStencil)
|
||||
self._Wy = sp.block_diag(mats)
|
||||
return self._Wy
|
||||
|
||||
@property
|
||||
def Wz(self):
|
||||
"""Regularization matrix Wz"""
|
||||
if getattr(self, '_Wz', None) is None:
|
||||
mats = []
|
||||
for imodel in range(self.nModels):
|
||||
self.wz.append(Utils.sdiag((self.regmesh.aveCC2Fz * self.regmesh.vol*self.alpha_z*self.ratios[imodel]*(self.regmesh.aveCC2Fz*self.wght))**0.5))
|
||||
mats.append(self.wz[imodel]*self.regmesh.cellDiffzStencil)
|
||||
self._Wz = sp.block_diag(mats)
|
||||
return self._Wz
|
||||
|
||||
@property
|
||||
def Wsmooth(self):
|
||||
"""Full smoothness regularization matrix W"""
|
||||
if getattr(self, '_Wsmooth', None) is None:
|
||||
wlist = (self.Wx,)
|
||||
if self.regmesh.dim > 1:
|
||||
wlist += (self.Wy,)
|
||||
if self.regmesh.dim > 2:
|
||||
wlist += (self.Wz,)
|
||||
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.Wsmall, self.Wsmooth)
|
||||
self._W = sp.vstack(wlist)
|
||||
return self._W
|
||||
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
return self._evalSmall(m) + self._evalSmooth(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmooth(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wsmooth * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wsmooth * ( self.mapping * m)
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
def cross(a,b):
|
||||
ax, ay, az = a[0], a[1], a[2]
|
||||
bx, by, bz = b[0], b[1], b[2]
|
||||
cx = ay*bz - az*by
|
||||
cy = az*bx - ax*bz
|
||||
cz = ax*by - ay*bx
|
||||
return [cx, cy, cz]
|
||||
|
||||
# TODO: Implement Cross Gradients..
|
||||
@Utils.timeIt
|
||||
def _evalCross(self, m):
|
||||
if self.crossgrad == False:
|
||||
return 0.
|
||||
elif self.crossgrad == True:
|
||||
M = (self.mapping * m).reshape((self.regmesh.nC, self.nModels), order="F")
|
||||
|
||||
ax = self.regmesh.aveFx2CC*self.regmesh.wx[0]*M[:,0]
|
||||
ay = self.regmesh.aveFy2CC*self.regmesh.wy[0]*M[:,0]
|
||||
az = self.regmesh.aveFz2CC*self.regmesh.wz[0]*M[:,0]
|
||||
bx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
|
||||
by = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
|
||||
bz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
|
||||
#ab
|
||||
out_ab = cross([ax, ay, az], [bx, by, bz])
|
||||
r = np.r_[out_ab[0], out_ab[1], out_ab[2]]*np.sqrt(self.betacross)
|
||||
|
||||
if self.nModels == 3:
|
||||
cx = self.regmesh.aveFx2CC*self.regmesh.wx[1]*M[:,1]
|
||||
cy = self.regmesh.aveFy2CC*self.regmesh.wy[1]*M[:,1]
|
||||
cz = self.regmesh.aveFz2CC*self.regmesh.wz[1]*M[:,1]
|
||||
#ac
|
||||
out_ac = cross([ax, ay, az], [cx, cy, cz])
|
||||
#bc
|
||||
out_bc = cross([bx, by, bz], [cx, cy, cz])
|
||||
r = np.r_[r, np.hstack(out_ac)*np.sqrt(self.betacross), np.hstack(out_bc)*np.sqrt(self.betacross)]
|
||||
|
||||
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})}
|
||||
|
||||
"""
|
||||
deriv = self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
|
||||
if self.crossgrad==True:
|
||||
deriv += self._evalCrossDeriv(m)
|
||||
return deriv
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalCrossDeriv(self,m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
|
||||
|
||||
@Utils.timeIt
|
||||
def eval2Deriv(self, m, v=None):
|
||||
"""
|
||||
Second derivative
|
||||
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array v: vector to multiply
|
||||
:rtype: scipy.sparse.csr_matrix or numpy.ndarray
|
||||
:return: WtW or WtW*v
|
||||
|
||||
The regularization is:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \\frac{1}{2}\mathbf{(m-m_\\text{ref})^\\top W^\\top W(m-m_\\text{ref})}
|
||||
|
||||
So the second derivative is straight forward:
|
||||
|
||||
.. math::
|
||||
|
||||
R(m) = \mathbf{W^\\top W}
|
||||
|
||||
"""
|
||||
mD = self.mapping.deriv(m - self.mref)
|
||||
if v is None:
|
||||
return mD.T * self.W.T * self.W * mD
|
||||
|
||||
return mD.T * ( self.W.T * ( self.W * ( mD * v) ) )
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,88 @@
|
||||
import SimPEG
|
||||
import numpy as np
|
||||
from SimPEG.Utils import Zero, closestPoints
|
||||
|
||||
class BaseRx(SimPEG.Survey.BaseTimeRx):
|
||||
locs = None
|
||||
rxType = None
|
||||
|
||||
knownRxTypes = {
|
||||
'phi':['phi',None],
|
||||
'ex':['e','x'],
|
||||
'ey':['e','y'],
|
||||
'ez':['e','z'],
|
||||
'jx':['j','x'],
|
||||
'jy':['j','y'],
|
||||
'jz':['j','z'],
|
||||
}
|
||||
|
||||
def __init__(self, locs, times, rxType, **kwargs):
|
||||
SimPEG.Survey.BaseTimeRx.__init__(self, locs, times, rxType, **kwargs)
|
||||
|
||||
@property
|
||||
def projField(self):
|
||||
"""Field Type projection (e.g. e b ...)"""
|
||||
return self.knownRxTypes[self.rxType][0]
|
||||
|
||||
def projGLoc(self, f):
|
||||
"""Grid Location projection (e.g. Ex Fy ...)"""
|
||||
comp = self.knownRxTypes[self.rxType][1]
|
||||
if comp is not None:
|
||||
return f._GLoc(self.rxType) + comp
|
||||
return f._GLoc(self.rxType)
|
||||
|
||||
def getTimeP(self, timesall):
|
||||
"""
|
||||
Returns the time projection matrix.
|
||||
|
||||
.. note::
|
||||
|
||||
This is not stored in memory, but is created on demand.
|
||||
"""
|
||||
time_inds = np.in1d(timesall, self.times)
|
||||
return time_inds
|
||||
|
||||
def evalDeriv(self, src, mesh, f, v, adjoint=False):
|
||||
P = self.getP(mesh, self.projGLoc(f))
|
||||
if not adjoint:
|
||||
return P*v
|
||||
elif adjoint:
|
||||
return P.T*v
|
||||
|
||||
|
||||
# DC.Rx.Dipole(locs)
|
||||
class Dipole(BaseRx):
|
||||
|
||||
def __init__(self, locsM, locsN, times, rxType = 'phi', **kwargs):
|
||||
assert locsM.shape == locsN.shape, 'locsM and locsN need to be the same size'
|
||||
locs = [locsM, locsN]
|
||||
# We may not need this ...
|
||||
BaseRx.__init__(self, locs, times, rxType)
|
||||
|
||||
@property
|
||||
def nD(self):
|
||||
"""Number of data in the receiver."""
|
||||
# return self.locs[0].shape[0] * len(self.times)
|
||||
return self.locs[0].shape[0]
|
||||
|
||||
@property
|
||||
def nRx(self):
|
||||
"""Number of data in the receiver."""
|
||||
return self.locs[0].shape[0]
|
||||
|
||||
# Not sure why ...
|
||||
# return int(self.locs[0].size / 2)
|
||||
|
||||
|
||||
def getP(self, mesh, Gloc):
|
||||
if mesh in self._Ps:
|
||||
return self._Ps[mesh]
|
||||
|
||||
P0 = mesh.getInterpolationMat(self.locs[0], Gloc)
|
||||
P1 = mesh.getInterpolationMat(self.locs[1], Gloc)
|
||||
P = P0 - P1
|
||||
|
||||
if self.storeProjections:
|
||||
self._Ps[mesh] = P
|
||||
|
||||
return P
|
||||
@@ -0,0 +1,64 @@
|
||||
import SimPEG
|
||||
# from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG.Utils import Zero, closestPoints, mkvc
|
||||
import numpy as np
|
||||
|
||||
class BaseSrc(SimPEG.Survey.BaseSrc):
|
||||
|
||||
current = 1.0
|
||||
loc = None
|
||||
|
||||
def __init__(self, rxList, **kwargs):
|
||||
SimPEG.Survey.BaseSrc.__init__(self, rxList, **kwargs)
|
||||
|
||||
def eval(self, prob):
|
||||
raise NotImplementedError
|
||||
|
||||
def evalDeriv(self, prob):
|
||||
return Zero()
|
||||
|
||||
@property
|
||||
def nD(self):
|
||||
"""Number of data"""
|
||||
return self.vnD.sum()
|
||||
|
||||
@property
|
||||
def vnD(self):
|
||||
"""Vector number of data"""
|
||||
return np.array([rx.nD*len(rx.times) for rx in self.rxList])
|
||||
|
||||
|
||||
|
||||
class Dipole(BaseSrc):
|
||||
|
||||
def __init__(self, rxList, locA, locB, **kwargs):
|
||||
assert locA.shape == locB.shape, 'Shape of locA and locB should be the same'
|
||||
self.loc = [locA, locB]
|
||||
BaseSrc.__init__(self, rxList, **kwargs)
|
||||
|
||||
def eval(self, prob):
|
||||
if prob._formulation == 'HJ':
|
||||
inds = closestPoints(prob.mesh, self.loc, gridLoc='CC')
|
||||
q = np.zeros(prob.mesh.nC)
|
||||
q[inds] = self.current * np.r_[1., -1.]
|
||||
elif prob._formulation == 'EB':
|
||||
qa = prob.mesh.getInterpolationMat(self.loc[0], locType='N').todense()
|
||||
qb = -prob.mesh.getInterpolationMat(self.loc[1], locType='N').todense()
|
||||
q = self.current * mkvc(qa+qb)
|
||||
return q
|
||||
|
||||
class Pole(BaseSrc):
|
||||
|
||||
def __init__(self, rxList, loc, **kwargs):
|
||||
BaseSrc.__init__(self, rxList, loc=loc, **kwargs)
|
||||
|
||||
def eval(self, prob):
|
||||
if prob._formulation == 'HJ':
|
||||
inds = closestPoints(prob.mesh, self.loc)
|
||||
q = np.zeros(prob.mesh.nC)
|
||||
q[inds] = self.current * np.r_[1.]
|
||||
elif prob._formulation == 'EB':
|
||||
q = prob.mesh.getInterpolationMat(self.loc, locType='N').todense()
|
||||
q = self.current * mkvc(q)
|
||||
return q
|
||||
|
||||
@@ -0,0 +1,102 @@
|
||||
import SimPEG
|
||||
from SimPEG.EM.Base import BaseEMSurvey
|
||||
from SimPEG import np, sp, Survey, Utils
|
||||
from SimPEG.Utils import Zero, Identity
|
||||
from SimPEG.EM.Static.SIP.SrcSIP import BaseSrc
|
||||
from SimPEG.EM.Static.SIP.RxSIP import BaseRx
|
||||
import uuid
|
||||
|
||||
|
||||
class Survey(BaseEMSurvey):
|
||||
rxPair = BaseRx
|
||||
srcPair = BaseSrc
|
||||
times = None
|
||||
|
||||
def __init__(self, srcList, **kwargs):
|
||||
self.srcList = srcList
|
||||
BaseEMSurvey.__init__(self, srcList, **kwargs)
|
||||
self.getUniqueTimes()
|
||||
|
||||
def getUniqueTimes(self):
|
||||
time_rx = []
|
||||
for src in self.srcList:
|
||||
for rx in src.rxList:
|
||||
time_rx.append(rx.times)
|
||||
self.times = np.unique(np.hstack(time_rx))
|
||||
|
||||
def dpred(self, m, f=None):
|
||||
"""
|
||||
Predicted data.
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = Pf(m)
|
||||
"""
|
||||
return self.prob.forward(m, f=f)
|
||||
|
||||
|
||||
class Data(SimPEG.Survey.Data):
|
||||
"""Fancy data storage by Src and Rx"""
|
||||
|
||||
def __init__(self, survey, v=None):
|
||||
self.uid = str(uuid.uuid4())
|
||||
self.survey = survey
|
||||
self._dataDict = {}
|
||||
for src in self.survey.srcList:
|
||||
self._dataDict[src] = {}
|
||||
for rx in src.rxList:
|
||||
self._dataDict[src][rx] = {}
|
||||
|
||||
if v is not None:
|
||||
self.fromvec(v)
|
||||
|
||||
def _ensureCorrectKey(self, key):
|
||||
if type(key) is tuple:
|
||||
if len(key) is not 3:
|
||||
raise KeyError('Key must be [Src, Rx, tInd]')
|
||||
if key[0] not in self.survey.srcList:
|
||||
raise KeyError('Src Key must be a source in the survey.')
|
||||
if key[1] not in key[0].rxList:
|
||||
raise KeyError('Rx Key must be a receiver for the source.')
|
||||
return key
|
||||
elif isinstance(key, self.survey.srcPair):
|
||||
if key not in self.survey.srcList:
|
||||
raise KeyError('Key must be a source in the survey.')
|
||||
return key, None, None
|
||||
else:
|
||||
raise KeyError('Key must be [Src] or [Src,Rx] or [Src, Rx, tInd]')
|
||||
|
||||
def __setitem__(self, key, value):
|
||||
src, rx, t = self._ensureCorrectKey(key)
|
||||
assert rx is not None, 'set data using [Src, Rx]'
|
||||
assert isinstance(value, np.ndarray), 'value must by ndarray'
|
||||
assert value.size == rx.nD, "value must have the same number of data as the source."
|
||||
self._dataDict[src][rx][t] = Utils.mkvc(value)
|
||||
|
||||
def __getitem__(self, key):
|
||||
src, rx, t = self._ensureCorrectKey(key)
|
||||
if rx is not None:
|
||||
if rx not in self._dataDict[src]:
|
||||
raise Exception('Data for receiver has not yet been set.')
|
||||
return self._dataDict[src][rx][t]
|
||||
|
||||
return np.concatenate([self[src,rx, t] for rx in src.rxList])
|
||||
|
||||
def tovec(self):
|
||||
val = []
|
||||
for src in self.survey.srcList:
|
||||
for rx in src.rxList:
|
||||
for t in rx.times:
|
||||
val.append(self[src, rx, t])
|
||||
return np.concatenate(val)
|
||||
|
||||
|
||||
def fromvec(self, v):
|
||||
v = Utils.mkvc(v)
|
||||
assert v.size == self.survey.nD, 'v must have the correct number of data.'
|
||||
indBot, indTop = 0, 0
|
||||
for src in self.survey.srcList:
|
||||
for rx in src.rxList:
|
||||
for t in rx.times:
|
||||
indTop += rx.nRx
|
||||
self[src, rx, t] = v[indBot:indTop]
|
||||
indBot += rx.nRx
|
||||
@@ -0,0 +1,5 @@
|
||||
from ProblemSIP import Problem3D_CC, Problem3D_N
|
||||
from SurveySIP import Survey, Data
|
||||
import SrcSIP as Src #Pole
|
||||
import RxSIP as Rx
|
||||
from Regularization import MultiRegularization
|
||||
@@ -1,2 +1,3 @@
|
||||
import DC
|
||||
import IP
|
||||
import SIP
|
||||
|
||||
@@ -0,0 +1,232 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import SimPEG
|
||||
from SimPEG import Mesh, Utils, EM, Maps, np, Survey
|
||||
from SimPEG.EM.Static import SIP, DC, IP
|
||||
from pymatsolver import MumpsSolver
|
||||
|
||||
|
||||
class IPProblemTestsCC(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
|
||||
problem = SIP.Problem3D_CC(mesh, rho=1./sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta, 1./tau]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC*2)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-10
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class IPProblemTestsN(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCN")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)), ("taui", Maps.IdentityMap(mesh))]
|
||||
problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta, 1./tau]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC*2)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
class IPProblemTestsN_air(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
|
||||
cs = 25.
|
||||
hx = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hy = [(cs,0, -1.3),(cs,21),(cs,0, 1.3)]
|
||||
hz = [(cs,0, -1.3),(cs,20),(cs,0, 1.3)]
|
||||
mesh = Mesh.TensorMesh([hx, hy, hz],x0="CCC")
|
||||
blkind0 = Utils.ModelBuilder.getIndicesSphere(np.r_[-100., -100., -200.], 75., mesh.gridCC)
|
||||
blkind1 = Utils.ModelBuilder.getIndicesSphere(np.r_[100., 100., -200.], 75., mesh.gridCC)
|
||||
sigma = np.ones(mesh.nC)*1e-2
|
||||
airind = mesh.gridCC[:,2]>0.
|
||||
sigma[airind] = 1e-8
|
||||
eta = np.zeros(mesh.nC)
|
||||
tau = np.ones_like(sigma)*1.
|
||||
eta[blkind0] = 0.1
|
||||
eta[blkind1] = 0.1
|
||||
tau[blkind0] = 0.1
|
||||
tau[blkind1] = 0.01
|
||||
|
||||
actmapeta = Maps.InjectActiveCells(mesh, ~airind, 0.)
|
||||
actmaptau = Maps.InjectActiveCells(mesh, ~airind, 1.)
|
||||
|
||||
x = mesh.vectorCCx[(mesh.vectorCCx>-155.)&(mesh.vectorCCx<155.)]
|
||||
y = mesh.vectorCCx[(mesh.vectorCCy>-155.)&(mesh.vectorCCy<155.)]
|
||||
Aloc = np.r_[-200., 0., 0.]
|
||||
Bloc = np.r_[200., 0., 0.]
|
||||
M = Utils.ndgrid(x-25.,y, np.r_[0.])
|
||||
N = Utils.ndgrid(x+25.,y, np.r_[0.])
|
||||
|
||||
times = np.arange(10)*1e-3 + 1e-3
|
||||
rx = SIP.Rx.Dipole(M, N, times)
|
||||
src = SIP.Src.Dipole([rx], Aloc, Bloc)
|
||||
survey = SIP.Survey([src])
|
||||
colemap = [("eta", Maps.IdentityMap(mesh)*actmapeta), ("taui", Maps.IdentityMap(mesh)*actmaptau)]
|
||||
problem = SIP.Problem3D_N(mesh, sigma=sigma, mapping=colemap)
|
||||
problem.Solver = MumpsSolver
|
||||
problem.pair(survey)
|
||||
mSynth = np.r_[eta[~airind], 1./tau[~airind]]
|
||||
survey.makeSyntheticData(mSynth)
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
regmap = Maps.IdentityMap(nP=int(mSynth[~airind].size*2))
|
||||
reg = SIP.MultiRegularization(mesh, mapping=regmap, nModels=2, indActive=~airind)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.survey.nSrc)
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.survey.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.Jvec(self.m0, v))
|
||||
vtJtw = v.dot(self.p.Jtvec(self.m0, w))
|
||||
passed = np.abs(wtJv - vtJtw) < 1e-8
|
||||
print 'Adjoint Test', np.abs(wtJv - vtJtw), passed
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_dataObj(self):
|
||||
derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False, num=3)
|
||||
self.assertTrue(passed)
|
||||
|
||||
if __name__ == '__main__':
|
||||
unittest.main()
|
||||
Reference in New Issue
Block a user