mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-13 02:21:03 +08:00
Major updates.
This commit is contained in:
@@ -0,0 +1,9 @@
|
||||
[run]
|
||||
source = simpegDC
|
||||
omit =
|
||||
*/python?.?/*
|
||||
*/lib-python/?.?/*.py
|
||||
*/lib_pypy/_*.py
|
||||
*/site-packages/ordereddict.py
|
||||
*/site-packages/nose/*
|
||||
*/unittest2/*
|
||||
+28
-17
@@ -1,24 +1,35 @@
|
||||
language: python
|
||||
python:
|
||||
- "2.7"
|
||||
virtualenv:
|
||||
system_site_packages: true
|
||||
- 2.7
|
||||
|
||||
# Setup anaconda
|
||||
before_install:
|
||||
- sudo apt-get install -qq gcc gfortran libblas-dev liblapack-dev python-numpy python-scipy python-matplotlib python-pip
|
||||
- sudo pip install scipy --upgrade
|
||||
- sudo pip install numpy --upgrade
|
||||
- cd ../
|
||||
- if [ ${TRAVIS_PYTHON_VERSION:0:1} == "2" ]; then wget http://repo.continuum.io/miniconda/Miniconda-3.3.0-Linux-x86_64.sh -O miniconda.sh; else wget http://repo.continuum.io/miniconda/Miniconda3-3.3.0-Linux-x86_64.sh -O miniconda.sh; fi
|
||||
- chmod +x miniconda.sh
|
||||
- ./miniconda.sh -b
|
||||
- export PATH=/home/travis/anaconda/bin:/home/travis/miniconda/bin:$PATH
|
||||
- conda update --yes conda
|
||||
# The next couple lines fix a crash with multiprocessing on Travis and are not specific to using Miniconda
|
||||
- sudo rm -rf /dev/shm
|
||||
- sudo ln -s /run/shm /dev/shm
|
||||
|
||||
# Install packages
|
||||
install:
|
||||
- conda install --yes pip python=$TRAVIS_PYTHON_VERSION numpy scipy matplotlib cython
|
||||
- pip install nose-cov python-coveralls
|
||||
# Remove this when SimPEG is on pip
|
||||
- git clone https://github.com/simpeg/simpeg.git
|
||||
- cd simpeg/SimPEG/
|
||||
- python setup.py
|
||||
- cd ../../
|
||||
- echo export PYTHONPATH=$PYTHONPATH:/home/travis/build/simpeg/simpeg >> .bashrc
|
||||
- source .bashrc
|
||||
- cd simpegdc
|
||||
# command to install dependencies
|
||||
install: "pip install -r requirements.txt --use-mirrors"
|
||||
# command to run tests
|
||||
script: nosetests -v
|
||||
- cd simpeg/
|
||||
- python setup.py install
|
||||
- cd ../
|
||||
|
||||
# Run test
|
||||
script:
|
||||
- nosetests --with-cov --cov simpegDC --cov-config .coveragerc
|
||||
|
||||
# Calculate coverage
|
||||
after_success:
|
||||
- coveralls --config_file .coveragerc
|
||||
|
||||
notifications:
|
||||
email:
|
||||
|
||||
@@ -0,0 +1,211 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class DipoleTx(Survey.BaseTx):
|
||||
"""A dipole transmitter, locA and locB are moved to the closest cell-centers"""
|
||||
def __init__(self, locA, locB, rxList, **kwargs):
|
||||
super(DipoleTx, self).__init__((locA, locB), 'dipole', rxList, **kwargs)
|
||||
self._rhsDict = {}
|
||||
|
||||
def getRhs(self, mesh):
|
||||
if mesh not in self._rhsDict:
|
||||
pts = [self.loc[0], self.loc[1]]
|
||||
inds = Utils.closestPoints(mesh, pts)
|
||||
q = np.zeros(mesh.nC)
|
||||
q[inds] = [1., -1.]
|
||||
self._rhsDict[mesh] = q
|
||||
return self._rhsDict[mesh]
|
||||
|
||||
|
||||
class DipoleRx(Survey.BaseRx):
|
||||
"""A dipole transmitter, 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(DipoleRx, 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.
|
||||
|
||||
"""
|
||||
|
||||
def __init__(self, txList, **kwargs):
|
||||
self.txList = txList
|
||||
Survey.BaseSurvey.__init__(self, **kwargs)
|
||||
self._rhsDict = {}
|
||||
self._Ps = {}
|
||||
|
||||
def projectFields(self, u):
|
||||
"""
|
||||
Predicted data.
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = Pu(m)
|
||||
"""
|
||||
P = self.getP(self.prob.mesh)
|
||||
return P*mkvc(u)
|
||||
|
||||
def getRhs(self, mesh):
|
||||
if mesh not in self._rhsDict:
|
||||
RHSlist = [tx.getRhs(mesh) for tx in self.txList]
|
||||
RHS = np.array(RHSlist).T
|
||||
self._rhsDict[mesh] = RHS
|
||||
return self._rhsDict[mesh]
|
||||
|
||||
def getP(self, mesh):
|
||||
if mesh in self._Ps:
|
||||
return self._Ps[mesh]
|
||||
|
||||
P_tx = [sp.vstack([rx.getP(mesh) for rx in tx.rxList]) for tx in self.txList]
|
||||
|
||||
self._Ps[mesh] = sp.block_diag(P_tx)
|
||||
return self._Ps[mesh]
|
||||
|
||||
|
||||
|
||||
|
||||
class ProblemDC(Problem.BaseProblem):
|
||||
"""
|
||||
**ProblemDC**
|
||||
|
||||
Geophysical DC resistivity problem.
|
||||
|
||||
"""
|
||||
|
||||
surveyPair = SurveyDC
|
||||
Solver = Solver
|
||||
|
||||
def __init__(self, mesh, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh)
|
||||
self.mesh.setCellGradBC('neumann')
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
|
||||
deleteTheseOnModelUpdate = ['_A']
|
||||
|
||||
@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
|
||||
sigma = self.curModel.transform
|
||||
Msig = self.mesh.getFaceInnerProduct(sigma, invProp=True, invMat=True)
|
||||
self._A = D*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 fields(self, m):
|
||||
self.curModel = m
|
||||
A = self.A
|
||||
Ainv = self.Solver(A)
|
||||
Q = self.survey.getRhs(self.mesh)
|
||||
Phi = Ainv * Q
|
||||
for ii in range(Phi.shape[1]):
|
||||
# Remove the static shift for each phi column.
|
||||
Phi[:,ii] -= Phi[-1, ii]
|
||||
return Phi
|
||||
|
||||
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)
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
# Derivative of inner product, $\left(\mathbf{M}_{1/\sigma}^f\right)^{-1}$
|
||||
dMdsig = self.mesh.getFaceInnerProductDeriv(sigma, invProp=True, invMat=True)
|
||||
# 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 transmitter
|
||||
for i in range(self.survey.nTx):
|
||||
dAdsig = D * dMdsig( G * u[:,i] )
|
||||
dCdm_x_v[:, i] = dAdsig * dsigdm_x_v
|
||||
|
||||
# Take derivative of $C(m,u)$ w.r.t. $u$
|
||||
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( dCdu_inv * dCdm_x_v )
|
||||
return J_x_v # Make $\mathbf{Jv}$ a vector.
|
||||
|
||||
def Jtvec(self, m, v, u=None):
|
||||
"""Takes data, turns it into a model..ish"""
|
||||
|
||||
if u is None:
|
||||
u = self.fields(m)
|
||||
|
||||
u = self.survey.reshapeFields(u)
|
||||
v = self.survey.reshapeFields(v)
|
||||
|
||||
P = self.survey.getP(self.mesh)
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.getA(m)
|
||||
Av_dm = self.mesh.getFaceInnerProductDeriv(m)
|
||||
mT_dm = self.mapping.deriv(m)
|
||||
|
||||
dCdu = A.T
|
||||
Ainv = self.Solver(dCdu)
|
||||
w = Ainv * (P.T*v)
|
||||
|
||||
Jtv = 0
|
||||
for i, ui in enumerate(u.T): # loop over each column
|
||||
Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] )
|
||||
|
||||
Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
|
||||
return Jtv
|
||||
-248
@@ -1,248 +0,0 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
class DCData(Data.BaseData):
|
||||
"""
|
||||
**DCData**
|
||||
|
||||
Geophysical DC resistivity data.
|
||||
|
||||
"""
|
||||
|
||||
P = None #: projection
|
||||
|
||||
def __init__(self, **kwargs):
|
||||
Data.BaseData.__init__(self, **kwargs)
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
def reshapeFields(self, u):
|
||||
if len(u.shape) == 1:
|
||||
u = u.reshape([-1, self.RHS.shape[1]], order='F')
|
||||
return u
|
||||
|
||||
def projectField(self, u):
|
||||
"""
|
||||
Predicted data.
|
||||
|
||||
.. math::
|
||||
d_\\text{pred} = Pu(m)
|
||||
"""
|
||||
u = self.reshapeFields(u)
|
||||
return Utils.mkvc(self.P*u)
|
||||
|
||||
|
||||
|
||||
class DCProblem(Problem.BaseProblem):
|
||||
"""
|
||||
**DCProblem**
|
||||
|
||||
Geophysical DC resistivity problem.
|
||||
|
||||
"""
|
||||
|
||||
dataPair = DCData
|
||||
|
||||
def __init__(self, mesh, model, **kwargs):
|
||||
Problem.BaseProblem.__init__(self, mesh, model)
|
||||
self.mesh.setCellGradBC('neumann')
|
||||
Utils.setKwargs(self, **kwargs)
|
||||
|
||||
|
||||
def createMatrix(self, m):
|
||||
"""
|
||||
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.
|
||||
"""
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
sigma = self.model.transform(m)
|
||||
Msig = self.mesh.getFaceMass(sigma)
|
||||
A = D*Msig*G
|
||||
return A.tocsc()
|
||||
|
||||
def field(self, m):
|
||||
A = self.createMatrix(m)
|
||||
solve = Solver(A)
|
||||
phi = solve.solve(self.data.RHS)
|
||||
return Utils.mkvc(phi)
|
||||
|
||||
def J(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 ) )
|
||||
"""
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
|
||||
u = self.data.reshapeFields(u)
|
||||
|
||||
P = self.data.P
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.createMatrix(m)
|
||||
Av_dm = self.mesh.getFaceMassDeriv()
|
||||
mT_dm = self.model.transformDeriv(m)
|
||||
|
||||
dCdu = A
|
||||
|
||||
dCdm = np.empty_like(u)
|
||||
for i, ui in enumerate(u.T): # loop over each column
|
||||
dCdm[:, i] = D * ( Utils.sdiag( G * ui ) * ( Av_dm * ( mT_dm * v ) ) )
|
||||
|
||||
solve = Solver(dCdu)
|
||||
Jv = - P * solve.solve(dCdm)
|
||||
return Utils.mkvc(Jv)
|
||||
|
||||
def Jt(self, m, v, u=None):
|
||||
"""Takes data, turns it into a model..ish"""
|
||||
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
|
||||
u = self.data.reshapeFields(u)
|
||||
v = self.data.reshapeFields(v)
|
||||
|
||||
P = self.data.P
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.createMatrix(m)
|
||||
Av_dm = self.mesh.getFaceMassDeriv()
|
||||
mT_dm = self.model.transformDeriv(m)
|
||||
|
||||
dCdu = A.T
|
||||
solve = Solver(dCdu)
|
||||
|
||||
w = solve.solve(P.T*v)
|
||||
|
||||
Jtv = 0
|
||||
for i, ui in enumerate(u.T): # loop over each column
|
||||
Jtv += Utils.sdiag( G * ui ) * ( D.T * w[:,i] )
|
||||
|
||||
Jtv = - mT_dm.T * ( Av_dm.T * Jtv )
|
||||
return Jtv
|
||||
|
||||
|
||||
|
||||
def genTxRxmat(nelec, spacelec, surfloc, elecini, mesh):
|
||||
""" Generate projection matrix (Q) and """
|
||||
elecend = 0.5+spacelec*(nelec-1)
|
||||
elecLocR = np.linspace(elecini, elecend, nelec)
|
||||
elecLocT = elecLocR+1
|
||||
nrx = nelec-1
|
||||
ntx = nelec-1
|
||||
q = np.zeros((mesh.nC, ntx))
|
||||
Q = np.zeros((mesh.nC, nrx))
|
||||
|
||||
for i in range(nrx):
|
||||
|
||||
rxind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i]))
|
||||
rxind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocR[i+1]))
|
||||
|
||||
txind1 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i]))
|
||||
txind2 = np.argwhere((mesh.gridCC[:,0]==surfloc) & (mesh.gridCC[:,1]==elecLocT[i+1]))
|
||||
|
||||
q[txind1,i] = 1
|
||||
q[txind2,i] = -1
|
||||
Q[rxind1,i] = 1
|
||||
Q[rxind2,i] = -1
|
||||
|
||||
Q = sp.csr_matrix(Q)
|
||||
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
|
||||
return q, Q, rxmidLoc
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
# Create the mesh
|
||||
h1 = np.ones(20)
|
||||
h2 = np.ones(100)
|
||||
M = Mesh.TensorMesh([h1,h2])
|
||||
|
||||
# Create some parameters for the model
|
||||
sig1 = np.log(1)
|
||||
sig2 = np.log(0.01)
|
||||
|
||||
# Create a synthetic model from a block in a half-space
|
||||
p0 = [5, 10]
|
||||
p1 = [15, 50]
|
||||
condVals = [sig1, sig2]
|
||||
mSynth = Utils.ModelBuilder.defineBlockConductivity(M.gridCC,p0,p1,condVals)
|
||||
plt.colorbar(M.plotImage(mSynth))
|
||||
# plt.show()
|
||||
|
||||
# Set up the projection
|
||||
nelec = 50
|
||||
spacelec = 2
|
||||
surfloc = 0.5
|
||||
elecini = 0.5
|
||||
elecend = 0.5+spacelec*(nelec-1)
|
||||
elecLocR = np.linspace(elecini, elecend, nelec)
|
||||
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
|
||||
q, Q, rxmidloc = genTxRxmat(nelec, spacelec, surfloc, elecini, M)
|
||||
P = Q.T
|
||||
|
||||
model = Model.LogModel(M)
|
||||
prob = DCProblem(M, model)
|
||||
|
||||
# Create some data
|
||||
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=q)
|
||||
|
||||
u = prob.field(mSynth)
|
||||
u = data.reshapeFields(u)
|
||||
M.plotImage(u[:,10])
|
||||
plt.show()
|
||||
|
||||
# Now set up the prob to do some minimization
|
||||
# prob.dobs = dobs
|
||||
# prob.std = dobs*0 + 0.05
|
||||
m0 = M.gridCC[:,0]*0+sig2
|
||||
|
||||
reg = Regularization.Tikhonov(model)
|
||||
objFunc = ObjFunction.BaseObjFunction(data, reg)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=3, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
inv = Inversion.BaseInversion(objFunc, opt)
|
||||
|
||||
# Check Derivative
|
||||
derChk = lambda m: [objFunc.dataObj(m), objFunc.dataObjDeriv(m)]
|
||||
# Tests.checkDerivative(derChk, mSynth)
|
||||
|
||||
print objFunc.dataObj(m0)
|
||||
print objFunc.dataObj(mSynth)
|
||||
|
||||
m = inv.run(m0)
|
||||
|
||||
plt.colorbar(M.plotImage(m))
|
||||
print m
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,48 @@
|
||||
from SimPEG import *
|
||||
import simpegDC as DC
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
|
||||
def example(cs=2.5, nElecs=10, plotIt=False):
|
||||
|
||||
mesh = Mesh.TensorMesh([
|
||||
[(cs,10, -1.3),(cs,200/cs),(cs,10, 1.3)],
|
||||
[(cs,3, -1.3),(cs,3,1.3)],
|
||||
# [(cs,5, -1.3),(cs,10)]
|
||||
],'CN')
|
||||
if plotIt:
|
||||
mesh.plotGrid(showIt=True)
|
||||
|
||||
space = 1
|
||||
elocs = np.linspace(-100, 100, nElecs)
|
||||
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 transmitters and receivers
|
||||
i = 0
|
||||
getLoc = lambda ii, abmn: np.r_[elocs[WENNER[ii,abmn]],0]
|
||||
txList = []
|
||||
for i in range(WENNER.shape[0]):
|
||||
rx = DC.DipoleRx(getLoc(i,1),getLoc(i,2))
|
||||
tx = DC.DipoleTx(getLoc(i,0),getLoc(i,3), [rx])
|
||||
txList += [tx]
|
||||
|
||||
|
||||
survey = DC.SurveyDC(txList)
|
||||
problem = DC.ProblemDC(mesh)
|
||||
problem.pair(survey)
|
||||
|
||||
return mesh, survey, problem
|
||||
@@ -0,0 +1 @@
|
||||
import WennerArray
|
||||
@@ -1,76 +1,53 @@
|
||||
import unittest
|
||||
from SimPEG import *
|
||||
import simpegDC
|
||||
import simpegDC as DC
|
||||
|
||||
|
||||
class DCProblemTests(unittest.TestCase):
|
||||
|
||||
def setUp(self):
|
||||
# Create the mesh
|
||||
h1 = np.ones(20)
|
||||
h2 = np.ones(20)
|
||||
mesh = Mesh.TensorMesh([h1,h2])
|
||||
model = Model.BaseModel(mesh)
|
||||
|
||||
# Create some parameters for the model
|
||||
sig1 = 1
|
||||
sig2 = 0.01
|
||||
mesh, survey, problem = DC.Examples.WennerArray.example()
|
||||
|
||||
# Create a synthetic model from a block in a half-space
|
||||
p0 = [2, 2]
|
||||
p1 = [5, 5]
|
||||
condVals = [sig1, sig2]
|
||||
mSynth = Utils.ModelBuilder.defineBlockConductivity(mesh.gridCC,p0,p1,condVals)
|
||||
|
||||
# Set up the projection
|
||||
nelec = 10
|
||||
spacelec = 2
|
||||
surfloc = 0.5
|
||||
elecini = 0.5
|
||||
elecend = 0.5+spacelec*(nelec-1)
|
||||
elecLocR = np.linspace(elecini, elecend, nelec)
|
||||
rxmidLoc = (elecLocR[0:nelec-1]+elecLocR[1:nelec])*0.5
|
||||
q, Q, rxmidloc = simpegDC.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
|
||||
P = Q.T
|
||||
Q = Q.toarray()
|
||||
|
||||
# Create some data
|
||||
|
||||
prob = simpegDC.DCProblem(mesh, model)
|
||||
data = prob.createSyntheticData(mSynth, std=0.05, P=P, RHS=Q)
|
||||
mSynth = np.ones(mesh.nC)
|
||||
survey.makeSyntheticData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
reg = Regularization.Tikhonov(mesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIterLS=20, maxIter=10, tolF=1e-6, tolX=1e-6, tolG=1e-6, maxIterCG=6)
|
||||
reg = Regularization.Tikhonov(model)
|
||||
objFunc = ObjFunction.BaseObjFunction(data, reg, beta=1e4)
|
||||
inv = Inversion.BaseInversion(objFunc, opt)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta=1e4)
|
||||
inv = Inversion.BaseInversion(invProb)
|
||||
|
||||
self.inv = inv
|
||||
self.reg = reg
|
||||
self.p = prob
|
||||
self.p = problem
|
||||
self.mesh = mesh
|
||||
self.m0 = mSynth
|
||||
self.data = data
|
||||
self.objFunc = objFunc
|
||||
self.survey = survey
|
||||
self.dmis = dmis
|
||||
|
||||
def test_misfit(self):
|
||||
derChk = lambda m: [self.data.dpred(m), lambda mx: self.p.J(self.m0, mx)]
|
||||
derChk = lambda m: [self.survey.dpred(m), lambda mx: self.p.Jvec(self.m0, mx)]
|
||||
passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
|
||||
self.assertTrue(passed)
|
||||
|
||||
def test_adjoint(self):
|
||||
# Adjoint Test
|
||||
u = np.random.rand(self.mesh.nC*self.data.RHS.shape[1])
|
||||
v = np.random.rand(self.mesh.nC)
|
||||
w = np.random.rand(self.data.dobs.shape[0])
|
||||
wtJv = w.dot(self.p.J(self.m0, v, u=u))
|
||||
vtJtw = v.dot(self.p.Jt(self.m0, w, u=u))
|
||||
passed = (wtJv - vtJtw) < 1e-10
|
||||
self.assertTrue(passed)
|
||||
# def test_adjoint(self):
|
||||
# # Adjoint Test
|
||||
# u = np.random.rand(self.mesh.nC*self.survey.RHS.shape[1])
|
||||
# v = np.random.rand(self.mesh.nC)
|
||||
# w = np.random.rand(self.survey.dobs.shape[0])
|
||||
# wtJv = w.dot(self.p.Jvec(self.m0, v, u=u))
|
||||
# vtJtw = v.dot(self.p.Jtvec(self.m0, w, u=u))
|
||||
# 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.objFunc.dataObj(m), self.objFunc.dataObjDeriv(m)]
|
||||
Tests.checkDerivative(derChk, self.m0, plotIt=False)
|
||||
# def test_dataObj(self):
|
||||
# derChk = lambda m: [self.dmis.eval(m), self.dmis.evalDeriv(m)]
|
||||
# passed = Tests.checkDerivative(derChk, self.m0, plotIt=False)
|
||||
# self.assertTrue(passed)
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
@@ -1 +1,2 @@
|
||||
from DC import *
|
||||
from BaseDC import *
|
||||
import Examples
|
||||
|
||||
Reference in New Issue
Block a user