mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-30 18:16:21 +08:00
Merge branch 'Dom_Dev' of https://github.com/simpeg/simpeg into dcip/dev
This commit is contained in:
Rowan Cockett
+46
-9
@@ -239,14 +239,51 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
|
||||
|
||||
|
||||
# class UpdateReferenceModel(Parameter):
|
||||
class update_IRLS(InversionDirective):
|
||||
|
||||
# mref0 = None
|
||||
m = None
|
||||
eps_min = None
|
||||
factor = None
|
||||
gamma = None
|
||||
phi_m_last = None
|
||||
|
||||
def initialize(self):
|
||||
|
||||
# Scale the regularization for changes in norm
|
||||
if getattr(self, 'phi_m_last', None) is not None:
|
||||
self.reg.gamma = 1.
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.gamma = self.phi_m_last / phim_new
|
||||
|
||||
self.reg.gamma = self.gamma
|
||||
|
||||
def endIter(self):
|
||||
# Cool the threshold parameter
|
||||
if getattr(self, 'factor', None) is not None:
|
||||
eps = self.reg.eps / self.factor
|
||||
|
||||
if getattr(self, 'eps_min', None) is not None:
|
||||
self.reg.eps = np.max([self.eps_min,eps])
|
||||
else:
|
||||
self.reg.eps = eps
|
||||
|
||||
|
||||
# Update the model used for the IRLS weights
|
||||
if getattr(self, 'm', None) is None:
|
||||
self.reg.m = self.invProb.curModel
|
||||
|
||||
# Update the pre-conditioner
|
||||
diagA = np.sum(self.prob.G**2.,axis=0) + self.invProb.beta*(self.reg.W.T*self.reg.W).diagonal() * (self.reg.mapping * np.ones(self.prob.mesh.nC))**2.
|
||||
PC = Utils.sdiag(diagA**-1.)
|
||||
|
||||
# def nextIter(self):
|
||||
# mref = getattr(self, 'm_prev', None)
|
||||
# if mref is None:
|
||||
# if self.debug: print 'UpdateReferenceModel is using mref0'
|
||||
# mref = self.mref0
|
||||
# self.m_prev = self.invProb.m_current
|
||||
# return mref
|
||||
self.opt.approxHinv = PC
|
||||
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.reg.gamma = self.reg.gamma * self.invProb.phi_m_last / phim_new
|
||||
|
||||
#==============================================================================
|
||||
# import pylab as plt
|
||||
# plt.figure()
|
||||
# ax = plt.subplot(221)
|
||||
# self.prob.mesh.plotSlice(self.invProb.curModel, ax = ax, normal = 'Z', ind=-5, clim = (0, 0.005))
|
||||
#==============================================================================
|
||||
|
||||
@@ -0,0 +1,179 @@
|
||||
from SimPEG import *
|
||||
import simpegDCIP as DC
|
||||
import scipy.interpolate as interpolation
|
||||
import matplotlib.pyplot as plt
|
||||
import time
|
||||
import re
|
||||
|
||||
def run(loc=np.c_[[-50.,0.,-50.],[50.,0.,-50.]], sig=np.r_[1e-2,1e-1,1e-3], radi=np.r_[25.,25.], param = np.r_[30.,30.,5], stype = 'dpdp', plotIt=True):
|
||||
"""
|
||||
DC Forward Simulation
|
||||
|
||||
Forward model conductive spheres in a half-space and plot a pseudo-section
|
||||
|
||||
Created on Mon Feb 01 19:28:06 2016
|
||||
|
||||
@fourndo
|
||||
"""
|
||||
|
||||
# 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])
|
||||
|
||||
# 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 not re.match(stype,'pdp'):
|
||||
inds = Utils.closestPoints(mesh, np.asarray(Tx[ii]).T )
|
||||
RHS = mesh.getInterpolationMat(np.asarray(Tx[ii]).T, 'CC').T*( [-1,1] / mesh.vol[inds] )
|
||||
|
||||
else:
|
||||
|
||||
# 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] )
|
||||
|
||||
|
||||
# Iterative Solve
|
||||
Ainvb = sp.linalg.bicgstab(P*A,P*RHS, tol=1e-5)
|
||||
|
||||
# We now have the potential everywhere
|
||||
phi = 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(Tx,Rx)
|
||||
|
||||
|
||||
# Here is an example for the first tx-rx array
|
||||
if plotIt:
|
||||
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(Tx2d,Rx2d,data,mesh.vectorNz[-1],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()
|
||||
@@ -565,7 +565,58 @@ class DiffOperators(object):
|
||||
|
||||
return Pbc, Pin, Pout
|
||||
|
||||
|
||||
|
||||
def unitCellGradx():
|
||||
doc = """Cell centered Gradient in the x dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
def fget(self):
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_unitCellGradx', None) is None:
|
||||
|
||||
n = self.vnC
|
||||
gx = ddx(n[0]-1)
|
||||
gx_square = sp.vstack((gx,gx[-1,:]*-1), format="csr")
|
||||
|
||||
self._unitCellGradx = kron3(speye(n[2]), speye(n[1]), gx_square)
|
||||
|
||||
return self._unitCellGradx
|
||||
return locals()
|
||||
unitCellGradx = property(**unitCellGradx())
|
||||
|
||||
def unitCellGrady():
|
||||
doc = """Cell centered Gradient in they dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
def fget(self):
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_unitCellGrady', None) is None:
|
||||
|
||||
n = self.vnC
|
||||
gy = ddx(n[1]-1)
|
||||
gy_square = sp.vstack((gy,gy[-1,:]*-1), format="csr")
|
||||
|
||||
self._unitCellGrady = kron3(speye(n[2]), gy_square, speye(n[0]))
|
||||
|
||||
return self._unitCellGrady
|
||||
return locals()
|
||||
unitCellGrady = property(**unitCellGrady())
|
||||
|
||||
def unitCellGradz():
|
||||
doc = """Cell centered Gradient in they dimension used for
|
||||
regularization. The gradient operator is square (nC-by-nC)"""
|
||||
def fget(self):
|
||||
if self.dim < 3: return None
|
||||
if getattr(self, '_unitCellGradz', None) is None:
|
||||
|
||||
n = self.vnC
|
||||
gz = ddx(n[2]-1)
|
||||
gz_square = sp.vstack((gz,gz[-1,:]*-1), format="csr")
|
||||
|
||||
self._unitCellGradz = kron3( gz_square , speye(n[1]), speye(n[0]))
|
||||
|
||||
return self._unitCellGradz
|
||||
return locals()
|
||||
unitCellGradz = property(**unitCellGradz())
|
||||
|
||||
# --------------- Averaging ---------------------
|
||||
|
||||
@property
|
||||
|
||||
+15
-1
@@ -989,5 +989,19 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
|
||||
if np.logical_or(norm(resid)/normResid0 <= self.tolCG, cgiter == self.maxIterCG):
|
||||
cgFlag = 1
|
||||
# End CG Iterations
|
||||
|
||||
# Take a gradient step on the active cells if exist
|
||||
if temp != self.xc.size:
|
||||
|
||||
rhs_a = (Active) * -self.g
|
||||
|
||||
dm_i = max( abs( delx ) )
|
||||
dm_a = max( abs(rhs_a) )
|
||||
|
||||
delx = delx + rhs_a * dm_i / dm_a /10.
|
||||
|
||||
# Only keep gradients going in the right direction on the active set
|
||||
indx = ((self.xc<=self.lower) & (delx < 0)) | ((self.xc>=self.upper) & (delx > 0))
|
||||
delx[indx] = 0.
|
||||
|
||||
return delx
|
||||
return delx
|
||||
@@ -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):
|
||||
Problem.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)
|
||||
|
||||
|
||||
@@ -282,3 +282,239 @@ class Tikhonov(BaseRegularization):
|
||||
out = mD.T * ( self.W.T * r )
|
||||
return out
|
||||
|
||||
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
|
||||
|
||||
+8
-1
@@ -1,6 +1,5 @@
|
||||
import Utils, numpy as np, scipy.sparse as sp, uuid
|
||||
|
||||
|
||||
class BaseRx(object):
|
||||
"""SimPEG Receiver Object"""
|
||||
|
||||
@@ -375,3 +374,11 @@ class BaseSurvey(object):
|
||||
self.dobs = self.dtrue+noise
|
||||
self.std = self.dobs*0 + std
|
||||
return self.dobs
|
||||
|
||||
class LinearSurvey(BaseSurvey):
|
||||
def projectFields(self, u):
|
||||
return u
|
||||
|
||||
@property
|
||||
def nD(self):
|
||||
return self.prob.G.shape[0]
|
||||
|
||||
@@ -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.
|
||||
|
||||
Reference in New Issue
Block a user