mirror of
https://github.com/wassname/simpeg.git
synced 2026-07-12 12:19:11 +08:00
Fixes to ModelBuilder. Start of the DCProblem.
This commit is contained in:
Rowan Cockett
@@ -1,6 +1,6 @@
|
||||
import numpy as np
|
||||
from scipy import sparse as sp
|
||||
from utils import mkvc, sdiag, speye, kron3, spzeros
|
||||
from SimPEG.utils import mkvc, sdiag, speye, kron3, spzeros
|
||||
|
||||
|
||||
def ddx(n):
|
||||
@@ -287,15 +287,19 @@ class DiffOperators(object):
|
||||
nodalVectorAve = property(**nodalVectorAve())
|
||||
|
||||
def getEdgeMass(self, materialProp=None):
|
||||
"""mass matix for products of edge functions w'*M(materialProp)*e"""
|
||||
"""mass matrix for products of edge functions w'*M(materialProp)*e"""
|
||||
if(materialProp is None):
|
||||
materialProp = np.ones(self.nC)
|
||||
Av = self.edgeAve
|
||||
return sdiag(Av.T * (self.vol * mkvc(materialProp)))
|
||||
|
||||
def getFaceMass(self, materialProp=None):
|
||||
"""mass matix for products of edge functions w'*M(materialProp)*e"""
|
||||
"""mass matrix for products of face functions w'*M(materialProp)*f"""
|
||||
if(materialProp is None):
|
||||
materialProp = np.ones(self.nC)
|
||||
Av = self.faceAve
|
||||
return sdiag(Av.T*(self.vol*mkvc(materialProp)))
|
||||
return sdiag(Av.T * (self.vol * mkvc(materialProp)))
|
||||
|
||||
def getFaceMassDeriv(self):
|
||||
Av = self.faceAve
|
||||
return Av.T * sdiag(self.vol)
|
||||
|
||||
@@ -1,5 +1,5 @@
|
||||
from scipy import sparse as sp
|
||||
from utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
|
||||
from SimPEG.utils import sub2ind, ndgrid, mkvc, getSubArray, sdiag, inv3X3BlockDiagonal, inv2X2BlockDiagonal
|
||||
import numpy as np
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,93 @@
|
||||
from SimPEG import TensorMesh
|
||||
from SimPEG.forward import Problem, SyntheticProblem
|
||||
from SimPEG.utils import ModelBuilder
|
||||
import numpy as np
|
||||
import scipy.sparse.linalg as linalg
|
||||
import DCutils
|
||||
|
||||
class DCProblem(Problem):
|
||||
"""docstring for DCProblem"""
|
||||
def __init__(self, mesh):
|
||||
super(DCProblem, self).__init__(mesh)
|
||||
self.mesh.setCellGradBC('neumann')
|
||||
|
||||
def createMatrix(self, m):
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
sigma = self.modelTransform(m)
|
||||
Msig = self.mesh.getFaceMass(sigma)
|
||||
A = D*Msig*G
|
||||
return A.tocsc()
|
||||
|
||||
def field(self, m):
|
||||
A = self.createMatrix(m)
|
||||
solve = linalg.factorized(A)
|
||||
|
||||
nRHSs = self.RHS.shape[1] # Number of RHSs
|
||||
phi = np.zeros((self.mesh.nC, nRHSs)) + np.nan
|
||||
for ii in range(nRHSs):
|
||||
phi[:,ii] = solve(self.RHS[:,ii])
|
||||
|
||||
return phi
|
||||
|
||||
def J(self, m, v, u=None, RHSii=0, solve=None):
|
||||
P = self.P
|
||||
D = self.mesh.faceDiv
|
||||
G = self.mesh.cellGrad
|
||||
A = self.createMatrix(m)
|
||||
Av_dm = self.mesh.getFaceMassDeriv()
|
||||
mT_dm = self.modelTransform(m)
|
||||
|
||||
dCdu = A
|
||||
dCdm = - D * ( sdiag( G * u[:, RHSii] ) * ( Av_dm * ( mT_dm * v ) ) )
|
||||
|
||||
if solve is None:
|
||||
solve = linalg.factorized(dCdu)
|
||||
|
||||
return - P * solve(dCdm)
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
# Create the mesh
|
||||
h1 = np.ones(100)
|
||||
h2 = np.ones(100)
|
||||
mesh = TensorMesh([h1,h2])
|
||||
|
||||
# Create some parameters for the model
|
||||
sig1 = 1
|
||||
sig2 = 0.01
|
||||
|
||||
# Create a synthetic model from a block in a half-space
|
||||
p0 = [20, 20]
|
||||
p1 = [50, 50]
|
||||
condVals = [sig1, sig2]
|
||||
mSynth = ModelBuilder.defineBlockConductivity(p0,p1,mesh.gridCC,condVals)
|
||||
mesh.plotImage(mSynth, showIt=False)
|
||||
|
||||
|
||||
# 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 = DCutils.genTxRxmat(nelec, spacelec, surfloc, elecini, mesh)
|
||||
|
||||
|
||||
# Create some data
|
||||
class syntheticDCProblem(DCProblem, SyntheticProblem):
|
||||
pass
|
||||
|
||||
synthetic = syntheticDCProblem(mesh);
|
||||
synthetic.P = Q.T
|
||||
synthetic.RHS = q
|
||||
dobs, Wd = synthetic.createData(mSynth)
|
||||
|
||||
# Now set up the problem to do some minimization
|
||||
problem = DCProblem(mesh)
|
||||
|
||||
|
||||
|
||||
@@ -0,0 +1,29 @@
|
||||
import numpy as np
|
||||
import scipy.sparse as sp
|
||||
|
||||
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
|
||||
+65
-11
@@ -84,8 +84,15 @@ class Problem(object):
|
||||
self._P = value
|
||||
|
||||
|
||||
def J(self, u):
|
||||
def J(self, m, v, u=None, RHSii=0):
|
||||
"""
|
||||
:param numpy.array m: model
|
||||
:param numpy.array v: vector to multiply
|
||||
:param numpy.array u: fields
|
||||
:param int RHSii: which RHS to calculate sensitivity too
|
||||
:rtype: numpy.array
|
||||
:return: Jv
|
||||
|
||||
Working with the general PDE, c(m, u) = 0, where m is the model and u is the field,
|
||||
the sensitivity is defined as:
|
||||
|
||||
@@ -107,15 +114,26 @@ class Problem(object):
|
||||
"""
|
||||
pass
|
||||
|
||||
def Jt(self, v):
|
||||
def Jt(self, m, v, u=None, RHSii=0):
|
||||
"""
|
||||
:param numpy.array m: model
|
||||
:param numpy.array v: vector to multiply
|
||||
:param numpy.array u: fields
|
||||
:param int RHSii: which RHS to calculate sensitivity too
|
||||
:rtype: numpy.array
|
||||
:return: JTv
|
||||
|
||||
Transpose of J
|
||||
"""
|
||||
pass
|
||||
|
||||
def field(self, m):
|
||||
"""
|
||||
The fields.
|
||||
The field given the model.
|
||||
|
||||
.. math::
|
||||
u(m)
|
||||
|
||||
"""
|
||||
pass
|
||||
|
||||
@@ -179,10 +197,10 @@ class Problem(object):
|
||||
m = np.random.rand(5)
|
||||
return checkDerivative(lambda m : [self.modelTransform(m), self.modelTransformDeriv(m)], m)
|
||||
|
||||
def misfit(self, m, R=None):
|
||||
def misfit(self, m, u=None):
|
||||
"""
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array R: residual, R = W o (dpred - dobs)
|
||||
:param numpy.array u: fields
|
||||
:rtype: float
|
||||
:return: data misfit
|
||||
|
||||
@@ -195,15 +213,15 @@ class Problem(object):
|
||||
Where P is a projection matrix that brings the field on the full domain to the data measurement locations;
|
||||
u is the field of interest; d_obs is the observed data; and W is the weighting matrix.
|
||||
"""
|
||||
if R is None:
|
||||
R = self.W*(self.dpred(m) - self.dobs)
|
||||
|
||||
R = self.W*(self.dpred(m, u=u) - self.dobs)
|
||||
R = mkvc(R)
|
||||
return 0.5*R.inner(R)
|
||||
|
||||
def misfitDeriv(self, m, R=None, u=None):
|
||||
def misfitDeriv(self, m, u=None):
|
||||
"""
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array u: fields
|
||||
:rtype: numpy.array
|
||||
:return: data misfit derivative
|
||||
|
||||
@@ -213,6 +231,12 @@ class Problem(object):
|
||||
|
||||
\mu_\\text{data} = {1\over 2}\left| \mathbf{W} \circ (\mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}) \\right|_2^2
|
||||
|
||||
If the field, u, is provided, the calculation of the data is fast:
|
||||
|
||||
.. math::
|
||||
|
||||
\mathbf{d}_\\text{pred} = \mathbf{Pu(m)}
|
||||
|
||||
\mathbf{R} = \mathbf{d}_\\text{pred} - \mathbf{d}_\\text{obs}
|
||||
|
||||
\mu_\\text{data} = {1\over 2}\left| \mathbf{W \circ R} \\right|_2^2
|
||||
@@ -230,8 +254,7 @@ class Problem(object):
|
||||
if u is None:
|
||||
u = self.field(m)
|
||||
|
||||
if R is None:
|
||||
R = self.W*(self.dpred(m, u=u) - self.dobs)
|
||||
R = self.W*(self.dpred(m, u=u) - self.dobs)
|
||||
|
||||
dmisfit = 0
|
||||
for i in range(self.RHS.shape[1]): # Loop over each right hand side
|
||||
@@ -240,9 +263,40 @@ class Problem(object):
|
||||
return dmisfit
|
||||
|
||||
|
||||
class SyntheticProblem(object):
|
||||
"""
|
||||
Has helpful functions when dealing with synthetic problems
|
||||
|
||||
To use this class, inherit to your problem::
|
||||
|
||||
class mySyntheticExample(Problem, SyntheticProblem):
|
||||
pass
|
||||
"""
|
||||
def createData(self, m, std=0.05):
|
||||
"""
|
||||
:param numpy.array m: geophysical model
|
||||
:param numpy.array std: standard deviation
|
||||
:rtype: numpy.array, numpy.array
|
||||
:return: dobs, Wd
|
||||
|
||||
Create synthetic data given a model, and a standard deviation.
|
||||
|
||||
Returns the observed data with random Gaussian noise
|
||||
and Wd which is the same size as data, and can be used to weight the inversion.
|
||||
"""
|
||||
dobs = self.dpred(m)
|
||||
dobs = dobs
|
||||
noise = std*abs(dobs)*np.random.randn(*dobs.shape)
|
||||
dobs = dobs+noise
|
||||
eps = np.linalg.norm(mkvc(dobs),2)*1e-5
|
||||
Wd = 1/(abs(dobs)*std+eps)
|
||||
return dobs, Wd
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
from SimPEG.inverse import checkDerivative
|
||||
|
||||
p = Problem(None)
|
||||
m = np.random.rand(5)
|
||||
checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m)
|
||||
checkDerivative(lambda m : [p.modelTransform(m), p.modelTransformDeriv(m)], m, plotIt=False)
|
||||
|
||||
@@ -15,10 +15,6 @@ def getIndecesBlock(p0,p1,ccMesh):
|
||||
The points p0 and p1 must live in the the same dimensional space as the mesh.
|
||||
"""
|
||||
|
||||
# Validation of the input
|
||||
assert type(p0) == np.ndarray, "Vector must be a numpy array"
|
||||
assert type(p1) == np.ndarray, "Vector must be a numpy array"
|
||||
|
||||
# Validation: p0 and p1 live in the same dimensional space
|
||||
assert len(p0) == len(p1), "Dimension mismatch. len(p0) != len(p1)"
|
||||
|
||||
@@ -47,7 +43,7 @@ def getIndecesBlock(p0,p1,ccMesh):
|
||||
|
||||
ind = np.where(indX & indY)
|
||||
|
||||
else:
|
||||
elif dimMesh == 3:
|
||||
# Define the points
|
||||
x1 = p0[0]
|
||||
y1 = p0[1]
|
||||
@@ -98,13 +94,16 @@ def defineTwoLayeredConductivity(depth,ccMesh,condVals):
|
||||
|
||||
# Identify 1st cell centered reference point
|
||||
p0[0] = ccMesh[0,0]
|
||||
p0[1] = ccMesh[0,1]
|
||||
p0[2] = ccMesh[0,2]
|
||||
if dim>1: p0[1] = ccMesh[0,1]
|
||||
if dim>2: p0[2] = ccMesh[0,2]
|
||||
|
||||
# Identify the last cell-centered reference point
|
||||
p1[0] = ccMesh[-1,0]
|
||||
p1[1] = ccMesh[-1,1]
|
||||
p1[2] = ccMesh[-1,2] - depth;
|
||||
if dim>1: p1[1] = ccMesh[-1,1]
|
||||
if dim>2: p1[2] = ccMesh[-1,2]
|
||||
|
||||
# The depth is always defined on the last one.
|
||||
p1[len(p1)-1] -= depth
|
||||
|
||||
ind = getIndecesBlock(p0,p1,ccMesh)
|
||||
|
||||
@@ -117,23 +116,24 @@ def scalarConductivity(ccMesh,pFunction):
|
||||
Define the distribution conductivity in the mesh according to the
|
||||
analytical expression given in pFunction
|
||||
"""
|
||||
xCC = ccMesh[:,0]
|
||||
yCC = ccMesh[:,1]
|
||||
zCC = ccMesh[:,2]
|
||||
dim = np.size(ccMesh[0,:])
|
||||
CC = [ccMesh[:,0]]
|
||||
if dim>1: CC.append(ccMesh[:,1])
|
||||
if dim>2: CC.append(ccMesh[:,2])
|
||||
|
||||
sigma = pFunction(xCC,yCC,zCC)
|
||||
|
||||
sigma = pFunction(*CC)
|
||||
|
||||
return sigma
|
||||
|
||||
if __name__ == '__main__':
|
||||
|
||||
import sys
|
||||
sys.path.append('../')
|
||||
from TensorMesh import TensorMesh
|
||||
from SimPEG import TensorMesh
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
# Define the mesh
|
||||
|
||||
testDim = 3
|
||||
testDim = 2
|
||||
h1 = 0.3*np.ones(7)
|
||||
h1[0] = 0.5
|
||||
h1[-1] = 0.6
|
||||
@@ -157,8 +157,8 @@ if __name__ == '__main__':
|
||||
# ------------------- Test conductivities! --------------------------
|
||||
print('Testing 1 block conductivity')
|
||||
|
||||
p0 = np.array([0.5,0.5,0.5])
|
||||
p1 = np.array([1.0,1.0,1.0])
|
||||
p0 = np.array([0.5,0.5,0.5])[:testDim]
|
||||
p1 = np.array([1.0,1.0,1.0])[:testDim]
|
||||
condVals = np.array([100,1e-6])
|
||||
|
||||
sigma = defineBlockConductivity(p0,p1,ccMesh,condVals)
|
||||
@@ -167,6 +167,7 @@ if __name__ == '__main__':
|
||||
print sigma.shape
|
||||
M.plotImage(sigma)
|
||||
print 'Done with block! :)'
|
||||
plt.show()
|
||||
|
||||
# -----------------------------------------
|
||||
print('Testing the two layered model')
|
||||
@@ -178,11 +179,17 @@ if __name__ == '__main__':
|
||||
M.plotImage(sigma)
|
||||
print sigma
|
||||
print 'layer model!'
|
||||
plt.show()
|
||||
|
||||
# -----------------------------------------
|
||||
print('Testing scalar conductivity')
|
||||
|
||||
pFunction = lambda x,y,z: np.exp(x+y+z)
|
||||
if testDim == 1:
|
||||
pFunction = lambda x: np.exp(x)
|
||||
elif testDim == 2:
|
||||
pFunction = lambda x,y: np.exp(x+y)
|
||||
elif testDim == 3:
|
||||
pFunction = lambda x,y,z: np.exp(x+y+z)
|
||||
|
||||
sigma = scalarConductivity(ccMesh,pFunction)
|
||||
|
||||
@@ -190,5 +197,6 @@ if __name__ == '__main__':
|
||||
M.plotImage(sigma)
|
||||
print sigma
|
||||
print 'Scalar conductivity defined!'
|
||||
plt.show()
|
||||
|
||||
# -----------------------------------------
|
||||
|
||||
Reference in New Issue
Block a user