mirror of
https://github.com/wassname/simpeg.git
synced 2026-06-27 19:48:52 +08:00
Merge branch 'feat/sparse-regularization' into patch/sparse-dcip
# Conflicts: # SimPEG/Examples/__init__.py # SimPEG/Optimization.py
This commit is contained in:
Lindsey Heagy
+84
-46
@@ -237,39 +237,41 @@ class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
# Save the file as a npz
|
||||
np.savez('{:03d}-{:s}'.format(self.opt.iter,self.fileName), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
|
||||
|
||||
class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
"""SaveOutputDictEveryIteration
|
||||
A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
|
||||
"""
|
||||
|
||||
def initialize(self):
|
||||
print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
|
||||
|
||||
def endIter(self):
|
||||
# Save the data.
|
||||
ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
|
||||
phi_ms = 0.5*ms.dot(ms)
|
||||
if self.reg.mrefInSmooth == True:
|
||||
mref = self.reg.mref
|
||||
else:
|
||||
mref = 0
|
||||
mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
phi_mx = 0.5 * mx.dot(mx)
|
||||
if self.prob.mesh.dim==2:
|
||||
my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
phi_my = 0.5 * my.dot(my)
|
||||
else:
|
||||
phi_my = 'NaN'
|
||||
if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
|
||||
mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
phi_mz = 0.5 * mz.dot(mz)
|
||||
else:
|
||||
phi_mz = 'NaN'
|
||||
|
||||
|
||||
# Save the file as a npz
|
||||
np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
|
||||
|
||||
#==============================================================================
|
||||
# class SaveOutputDictEveryIteration(_SaveEveryIteration):
|
||||
# """SaveOutputDictEveryIteration
|
||||
# A directive that saves some relevant information from the inversion run to a numpy .npz dictionary file (see numpy.savez function for further info).
|
||||
# """
|
||||
#
|
||||
# def initialize(self):
|
||||
# print "SimPEG.SaveOutputDictEveryIteration will save your inversion progress as dictionary: '%s-###.npz'"%self.fileName
|
||||
#
|
||||
# def endIter(self):
|
||||
# # Save the data.
|
||||
# ms = self.reg.Ws * ( self.reg.mapping * (self.invProb.curModel - self.reg.mref) )
|
||||
# phi_ms = 0.5*ms.dot(ms)
|
||||
# if self.reg.mrefInSmooth == True:
|
||||
# mref = self.reg.mref
|
||||
# else:
|
||||
# mref = 0
|
||||
# mx = self.reg.Wx * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
# phi_mx = 0.5 * mx.dot(mx)
|
||||
# if self.prob.mesh.dim==2:
|
||||
# my = self.reg.Wy * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
# phi_my = 0.5 * my.dot(my)
|
||||
# else:
|
||||
# phi_my = 'NaN'
|
||||
# if self.prob.mesh.dim==3 and 'CYL' not in self.prob.mesh._meshType:
|
||||
# mz = self.reg.Wz * ( self.reg.mapping * (self.invProb.curModel - mref) )
|
||||
# phi_mz = 0.5 * mz.dot(mz)
|
||||
# else:
|
||||
# phi_mz = 'NaN'
|
||||
#
|
||||
#
|
||||
# # Save the file as a npz
|
||||
# np.savez('{:s}-{:03d}'.format(self.fileName,self.opt.iter), iter=self.opt.iter, beta=self.invProb.beta, phi_d=self.invProb.phi_d, phi_m=self.invProb.phi_m, phi_ms=phi_ms, phi_mx=phi_mx, phi_my=phi_my, phi_mz=phi_mz,f=self.opt.f, m=self.invProb.curModel,dpred=self.invProb.dpred)
|
||||
#
|
||||
#==============================================================================
|
||||
|
||||
# class UpdateReferenceModel(Parameter):
|
||||
|
||||
@@ -295,6 +297,8 @@ class update_IRLS(InversionDirective):
|
||||
|
||||
# Scale the regularization for changes in norm
|
||||
if getattr(self, 'phi_m_last', None) is not None:
|
||||
|
||||
self.reg.curModel = self.invProb.curModel
|
||||
self.reg.gamma = 1.
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
self.gamma = self.phi_m_last / phim_new
|
||||
@@ -320,25 +324,59 @@ class update_IRLS(InversionDirective):
|
||||
|
||||
# Update the model used for the IRLS weights
|
||||
self.reg.curModel = 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.reg.curModel.size))**2.
|
||||
PC = Utils.sdiag(diagA**-1.)
|
||||
self.opt.approxHinv = PC
|
||||
|
||||
# Temporarely set gamma to 1.
|
||||
self.reg.gamma = 1.
|
||||
|
||||
# Compute change in model objective function and update scaling
|
||||
phim_new = self.reg.eval(self.invProb.curModel)
|
||||
|
||||
self.reg.gamma = self.phi_m_last / phim_new
|
||||
|
||||
# TO DO: Re-scale beta if too much change in misfit
|
||||
self.invProb.beta = self.invProb.beta * self.phi_d_last / self.invProb.phi_d
|
||||
|
||||
#==============================================================================
|
||||
# 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))
|
||||
#==============================================================================
|
||||
# TO DO: Check optimization class, data misfit not matching reality
|
||||
#dpred = self.prob.fields(self.invProb.curModel)
|
||||
#phid = self.invProb.dmisfit.eval(self.invProb.curModel)
|
||||
#print self.survey.std[0]
|
||||
#print phid
|
||||
#print self.invProb.phi_d
|
||||
#print self.invProb.phi_d_last
|
||||
|
||||
self.invProb.beta = self.invProb.beta * self.survey.nD*0.5 / self.invProb.phi_d
|
||||
|
||||
class update_lin_PreCond(InversionDirective):
|
||||
|
||||
|
||||
def endIter(self):
|
||||
# Cool the threshold parameter
|
||||
|
||||
if getattr(self.opt, 'approxHinv', None) is not None:
|
||||
# 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.reg.curModel.size))**2.
|
||||
PC = Utils.sdiag(diagA**-1.)
|
||||
self.opt.approxHinv = PC
|
||||
print 'Updated pre-cond'
|
||||
|
||||
class update_Wj(InversionDirective):
|
||||
"""
|
||||
Create approx-sensitivity base weighting
|
||||
"""
|
||||
k = None
|
||||
|
||||
def endIter(self):
|
||||
|
||||
if self.opt.iter == 2:
|
||||
m = self.invProb.curModel
|
||||
if self.k is None:
|
||||
self.k = int(self.survey.nD/10)
|
||||
|
||||
def JtJv(v):
|
||||
|
||||
Jv = self.prob.Jvec(m, v)
|
||||
|
||||
return self.prob.Jtvec(m,Jv)
|
||||
|
||||
JtJdiag = Utils.diagEst(JtJv,len(m),k=self.k)
|
||||
JtJdiag = JtJdiag / max(JtJdiag)
|
||||
|
||||
self.reg.wght = JtJdiag
|
||||
|
||||
@@ -0,0 +1,220 @@
|
||||
from SimPEG import *
|
||||
# import simpegDCIP as DC
|
||||
from SimPEG import DCIP 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
|
||||
|
||||
"""
|
||||
|
||||
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
|
||||
# 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 = getIndicesSphere(loc[:,0],radi[0],mesh.gridCC)
|
||||
model[ind] = sig[1]
|
||||
|
||||
# Second anomaly
|
||||
ind = 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()
|
||||
@@ -0,0 +1,148 @@
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
from pymatsolver import MumpsSolver
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def run(plotIt=True):
|
||||
"""
|
||||
EM: FDEM: Effects of susceptibility
|
||||
===================================
|
||||
|
||||
When airborne freqeuncy domain EM (AFEM) survey is flown over
|
||||
the earth including significantly susceptible bodies (magnetite-rich rocks),
|
||||
negative data is often observed in the real part of the lowest frequency
|
||||
(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
|
||||
occurs due to a susceptible body when the magnetic field is applied.
|
||||
|
||||
To clarify what is happening in the earth when we are exciting the earth with
|
||||
a loop source in the frequency domain we run three forward modelling:
|
||||
|
||||
- F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
|
||||
- F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
|
||||
- F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
|
||||
|
||||
We plot vector magnetic fields in the earth. For secondary fields we provide
|
||||
F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
|
||||
figure show both real and parts.
|
||||
|
||||
"""
|
||||
# Generate Cylindrical mesh
|
||||
cs, ncx, ncz, npad = 5, 25, 24, 20.
|
||||
hx = [(cs,ncx), (cs,npad,1.3)]
|
||||
hz = [(cs,npad,-1.3), (cs,ncz), (cs,npad,1.3)]
|
||||
mesh = Mesh.CylMesh([hx,1,hz], '00C')
|
||||
sighalf = 1e-3
|
||||
sigma = np.ones(mesh.nC)*1e-8
|
||||
sigmahomo = sigma.copy()
|
||||
mu = np.ones(mesh.nC)*mu_0
|
||||
sigma[mesh.gridCC[:,-1]<0.] = sighalf
|
||||
blkind = np.logical_and(mesh.gridCC[:,0]<30., (mesh.gridCC[:,2]<0)&(mesh.gridCC[:,2]>-150)&(mesh.gridCC[:,2]<-50))
|
||||
sigma[blkind] = 1e-1
|
||||
mu[blkind] = mu_0*1.1
|
||||
offset = 0.
|
||||
frequency = np.r_[10., 100., 1000.]
|
||||
rx0 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzr')
|
||||
rx1 = EM.FDEM.Rx(np.array([[8., 0., 30.]]), 'bzi')
|
||||
srcLists = []
|
||||
nfreq = frequency.size
|
||||
for ifreq in range(nfreq):
|
||||
src = EM.FDEM.Src.CircularLoop([rx0, rx1], frequency[ifreq], np.array([[0., 0., 30.]]), radius=5.)
|
||||
srcLists.append(src)
|
||||
survey = EM.FDEM.Survey(srcLists)
|
||||
iMap = Maps.IdentityMap(nP=int(mesh.nC))
|
||||
# Use PhysPropMap
|
||||
maps = [('sigma', iMap), ('mu', iMap)]
|
||||
prob = EM.FDEM.Problem_b(mesh, mapping=maps)
|
||||
prob.Solver = MumpsSolver
|
||||
survey.pair(prob)
|
||||
m = np.r_[sigma, mu]
|
||||
survey0 = EM.FDEM.Survey(srcLists)
|
||||
prob0 = EM.FDEM.Problem_b(mesh, mapping=maps)
|
||||
prob0.Solver = MumpsSolver
|
||||
survey0.pair(prob0)
|
||||
m = np.r_[sigma, mu]
|
||||
m0 = np.r_[sigma, np.ones(mesh.nC)*mu_0]
|
||||
m00 = np.r_[np.ones(mesh.nC)*1e-8, np.ones(mesh.nC)*mu_0]
|
||||
# Anomalous conductivity and susceptibility
|
||||
F = prob.fields(m)
|
||||
# Only anomalous conductivity
|
||||
F0 = prob.fields(m0)
|
||||
# Primary field
|
||||
F00 = prob.fields(m00)
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
def vizfields(ifreq=0, primsec="secondary",realimag="real"):
|
||||
|
||||
titles = ["F[$\sigma$, $\mu$]", "F[$\sigma$, $\mu_0$]", "F[$\sigma$, $\mu$]-F[$\sigma$, $\mu_0$]"]
|
||||
actind = np.logical_and(mesh.gridCC[:,0]<200., (mesh.gridCC[:,2]>-400)&(mesh.gridCC[:,2]<200))
|
||||
|
||||
if primsec=="secondary":
|
||||
bCCprim = (mesh.aveF2CCV*F00[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
|
||||
bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
|
||||
bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')-bCCprim
|
||||
elif primsec=="primary":
|
||||
bCC = (mesh.aveF2CCV*F[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
|
||||
bCC0 = (mesh.aveF2CCV*F0[:,'b'][:,ifreq]).reshape(mesh.nC, 2, order='F')
|
||||
|
||||
XYZ = mesh.gridCC[actind,:]
|
||||
X = XYZ[:,0].reshape((31,43), order='F')
|
||||
Z = XYZ[:,2].reshape((31,43), order='F')
|
||||
bx = bCC[actind,0].reshape((31,43), order='F')
|
||||
bz = bCC[actind,1].reshape((31,43), order='F')
|
||||
bx0 = bCC0[actind,0].reshape((31,43), order='F')
|
||||
bz0 = bCC0[actind,1].reshape((31,43), order='F')
|
||||
|
||||
bxsec = (bCC[actind,0]-bCC0[actind,0]).reshape((31,43), order='F')
|
||||
bzsec = (bCC[actind,1]-bCC0[actind,1]).reshape((31,43), order='F')
|
||||
|
||||
absbreal = np.sqrt(bx.real**2+bz.real**2)
|
||||
absbimag = np.sqrt(bx.imag**2+bz.imag**2)
|
||||
absb0real = np.sqrt(bx0.real**2+bz0.real**2)
|
||||
absb0imag = np.sqrt(bx0.imag**2+bz0.imag**2)
|
||||
|
||||
absbrealsec = np.sqrt(bxsec.real**2+bzsec.real**2)
|
||||
absbimagsec = np.sqrt(bxsec.imag**2+bzsec.imag**2)
|
||||
|
||||
fig = plt.figure(figsize=(15,5))
|
||||
ax1 = plt.subplot(131)
|
||||
ax2 = plt.subplot(132)
|
||||
ax3 = plt.subplot(133)
|
||||
typefield="real"
|
||||
scale=20
|
||||
if realimag=="real":
|
||||
ax1.contourf(X, Z,np.log10(absbreal), 100)
|
||||
ax1.quiver(X, Z,bx.real/absbreal,bz.real/absbreal,scale=scale,width=0.005, alpha = 0.5)
|
||||
ax2.contourf(X, Z,np.log10(absb0real), 100)
|
||||
ax2.quiver(X, Z,bx0.real/absb0real,bz0.real/absb0real,scale=scale,width=0.005, alpha = 0.5)
|
||||
ax3.contourf(X, Z,np.log10(absbrealsec), 100)
|
||||
ax3.quiver(X, Z,bxsec.real/absbrealsec,bzsec.real/absbrealsec,scale=scale,width=0.005, alpha = 0.5)
|
||||
elif realimag=="imag":
|
||||
ax1.contourf(X, Z,np.log10(absbimag), 100)
|
||||
ax1.quiver(X, Z,bx.imag/absbimag,bz.imag/absbimag,scale=scale,width=0.005, alpha = 0.5)
|
||||
ax2.contourf(X, Z,np.log10(absb0imag), 100)
|
||||
ax2.quiver(X, Z,bx0.imag/absb0imag,bz0.imag/absb0imag,scale=scale,width=0.005, alpha = 0.5)
|
||||
ax3.contourf(X, Z,np.log10(absbimagsec), 100)
|
||||
ax3.quiver(X, Z,bxsec.imag/absbimagsec,bzsec.imag/absbimagsec,scale=scale,width=0.005, alpha = 0.5)
|
||||
|
||||
ax = [ax1, ax2, ax3]
|
||||
ax3.text(30, 50, ("Frequency=%5.2f Hz")%(frequency[ifreq]), color="k", fontsize=18)
|
||||
ax2.text(30, 50, primsec, color="k", fontsize=18)
|
||||
ax1.text(30, 50, realimag, color="k", fontsize=18)
|
||||
for i, axtemp in enumerate(ax):
|
||||
axtemp.plot(np.r_[0, 29.75], np.r_[-50, -50], 'w', lw=3)
|
||||
|
||||
axtemp.plot(np.r_[29.5, 29.5], np.r_[-50, -142.5], 'w', lw=3)
|
||||
axtemp.plot(np.r_[0, 29.5], np.r_[-142.5, -142.5], 'w', lw=3)
|
||||
axtemp.plot(np.r_[0, 100.], np.r_[0, 0], 'w', lw=3)
|
||||
axtemp.set_ylim(-200, 100.)
|
||||
axtemp.set_xlim(10, 100.)
|
||||
axtemp.set_title(titles[i])
|
||||
plt.show()
|
||||
return fig, ax
|
||||
fig1, ax1 = vizfields(1, primsec="primary", realimag="real")
|
||||
fig2, ax2 = vizfields(1, primsec="secondary", realimag="real")
|
||||
fig4, ax4 = vizfields(1, primsec="secondary", realimag="imag")
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
@@ -0,0 +1,132 @@
|
||||
from SimPEG import *
|
||||
|
||||
|
||||
def run(N=200, plotIt=True):
|
||||
"""
|
||||
Inversion: Linear Problem
|
||||
=========================
|
||||
|
||||
Here we go over the basics of creating a linear problem and inversion.
|
||||
|
||||
"""
|
||||
|
||||
|
||||
np.random.seed(1)
|
||||
|
||||
std_noise = 1e-2
|
||||
|
||||
mesh = Mesh.TensorMesh([N])
|
||||
|
||||
m0 = np.ones(mesh.nC) * 1e-4
|
||||
nk = 10
|
||||
jk = np.linspace(1.,nk,nk)
|
||||
p = -2.
|
||||
q = 1.
|
||||
|
||||
g = lambda k: np.exp(p*jk[k]*mesh.vectorCCx)*np.cos(np.pi*q*jk[k]*mesh.vectorCCx)
|
||||
|
||||
G = np.empty((nk, mesh.nC))
|
||||
|
||||
for i in range(nk):
|
||||
G[i,:] = g(i)
|
||||
|
||||
mtrue = np.zeros(mesh.nC)
|
||||
mtrue[mesh.vectorCCx > 0.3] = 1.
|
||||
mtrue[mesh.vectorCCx > 0.45] = -0.5
|
||||
mtrue[mesh.vectorCCx > 0.6] = 0
|
||||
|
||||
|
||||
prob = Problem.LinearProblem(mesh, G)
|
||||
survey = Survey.LinearSurvey()
|
||||
survey.pair(prob)
|
||||
survey.dobs = prob.fields(mtrue) + std_noise * np.random.randn(nk)
|
||||
#survey.makeSyntheticData(mtrue, std=std_noise)
|
||||
|
||||
wd = np.ones(nk) * std_noise
|
||||
|
||||
#print survey.std[0]
|
||||
#M = prob.mesh
|
||||
# Distance weighting
|
||||
wr = np.sum(prob.G**2.,axis=0)**0.5
|
||||
wr = ( wr/np.max(wr) )
|
||||
|
||||
reg = Regularization.Simple(mesh)
|
||||
reg.wght = wr
|
||||
|
||||
dmis = DataMisfit.l2_DataMisfit(survey)
|
||||
dmis.Wd = 1./wd
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=30,lower=-2.,upper=2., maxIterCG= 20, tolCG = 1e-4)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt)
|
||||
invProb.curModel = m0
|
||||
|
||||
beta = Directives.BetaSchedule(coolingFactor=2, coolingRate=1)
|
||||
target = Directives.TargetMisfit()
|
||||
|
||||
betaest = Directives.BetaEstimate_ByEig()
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta, betaest, target])
|
||||
|
||||
|
||||
mrec = inv.run(m0)
|
||||
ml2 = mrec
|
||||
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
|
||||
|
||||
# Switch regularization to sparse
|
||||
phim = invProb.phi_m_last
|
||||
phid = invProb.phi_d
|
||||
|
||||
reg = Regularization.Sparse(mesh)
|
||||
|
||||
#==============================================================================
|
||||
# fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
# dmdx = reg.mesh.cellDiffxStencil * mrec
|
||||
# plt.plot(np.sort(dmdx))
|
||||
#==============================================================================
|
||||
|
||||
#reg.recModel = mrec
|
||||
reg.wght = np.ones(mesh.nC)
|
||||
reg.mref = np.zeros(mesh.nC)
|
||||
reg.eps_p = 2e-3
|
||||
reg.eps_q = 2e-3
|
||||
reg.norms = [0., 0., 2., 2.]
|
||||
reg.wght = wr
|
||||
|
||||
opt = Optimization.ProjectedGNCG(maxIter=5 ,lower=-2.,upper=2., maxIterCG= 100, tolCG = 1e-3)
|
||||
invProb = InvProblem.BaseInvProblem(dmis, reg, opt, beta = invProb.beta*2.)
|
||||
beta = Directives.BetaSchedule(coolingFactor=1, coolingRate=1)
|
||||
#betaest = Directives.BetaEstimate_ByEig()
|
||||
target = Directives.TargetMisfit()
|
||||
IRLS =Directives.update_IRLS( phi_m_last = phim, phi_d_last = phid )
|
||||
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,IRLS])
|
||||
|
||||
m0 = mrec
|
||||
|
||||
# Run inversion
|
||||
mrec = inv.run(m0)
|
||||
|
||||
print "Final misfit:" + str(invProb.dmisfit.eval(mrec))
|
||||
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
|
||||
fig, axes = plt.subplots(1,2,figsize=(12*1.2,4*1.2))
|
||||
for i in range(prob.G.shape[0]):
|
||||
axes[0].plot(prob.G[i,:])
|
||||
axes[0].set_title('Columns of matrix G')
|
||||
|
||||
axes[1].plot(mesh.vectorCCx, mtrue, 'b-')
|
||||
axes[1].plot(mesh.vectorCCx, ml2, 'r-')
|
||||
#axes[1].legend(('True Model', 'Recovered Model'))
|
||||
axes[1].set_ylim(-1.0,1.25)
|
||||
|
||||
axes[1].plot(mesh.vectorCCx, mrec, 'k-',lw = 2)
|
||||
axes[1].legend(('True Model', 'Smooth l2-l2',
|
||||
'Sparse lp:' + str(reg.norms[0]) + ', lqx:' + str(reg.norms[1]) ), fontsize = 12)
|
||||
plt.show()
|
||||
|
||||
return prob, survey, mesh, mrec
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
@@ -0,0 +1,443 @@
|
||||
from scipy.constants import epsilon_0, mu_0
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from ipywidgets import *
|
||||
#from SimPEG.EM.Utils import k, omega
|
||||
|
||||
"""
|
||||
MT1D: n layered earth problem
|
||||
*****************************
|
||||
|
||||
Author: Thibaut Astic
|
||||
Contact: thast@eos.ubc.ca
|
||||
Date: January 2016
|
||||
|
||||
This code compute the analytic response of a n-layered Earth to a plane wave (Magneto-Tellurics).
|
||||
|
||||
We start by looking at Maxwell's equations in the electric
|
||||
field \\\(\\\mathbf{E}\\) and the magnetic flux
|
||||
\\\(\\\mathbf{H}\\) to write the wave equations
|
||||
\\(\\ \nabla ^2 \mathbf{E_x} + k^2 \mathbf{E_x} = 0 \\) &
|
||||
\\(\\ \nabla ^2 \mathbf{H_y} + k^2 \mathbf{H_y} = 0 \\)
|
||||
|
||||
Then solving the equations in each layer "j" between z_{j-1} and z_j in the form of
|
||||
\\(\\ E_{x,j} (z) = U_j e^{i k (z-z_{j-1})} + D_j e^{-i k (z-z_{j-1})} \\)
|
||||
\\(\\ H_{y,j} (z) = \frac{1}{Z_j} (D_j e^{-i k (z-z_{j-1})} - U_j e^{i k (z-z_{j-1})}) \\)
|
||||
|
||||
With U and D the Up and Down components of the E-field.
|
||||
|
||||
The iteration from one layer to another is ensure by:
|
||||
|
||||
\\(\\ \left(\begin{matrix} E_{x,j} \\ H_{y,j} \end{matrix} \right) =
|
||||
P_j T_j P^{-1}_J \left(\begin{matrix} E_{x,j+1} \\ H_{y,j+1} \end{matrix} \right) \\)
|
||||
|
||||
And the Boundary Condition is set for the E-field in the last layer, with no Up component (=0)
|
||||
and only a down component (=1 then normalized by the highest amplitude to ensure numeric stability)
|
||||
|
||||
The layer 0 is assumed to be the air layer.
|
||||
|
||||
"""
|
||||
|
||||
#Frequency conversion
|
||||
omega = lambda f: 2.*np.pi*f
|
||||
|
||||
#Evaluate k wavenumber
|
||||
k = lambda mu,sig,eps,f: np.sqrt(mu*mu_0*eps*epsilon_0*(2.*np.pi*f)**2.-1.j*mu*mu_0*sig*omega(f))
|
||||
|
||||
#Define a frquency range for a survey
|
||||
frange = lambda minfreq, maxfreq, step: np.logspace(minfreq,maxfreq,num = step, base = 10.)
|
||||
|
||||
#Functions to create random physical Perties for a n-layered earth
|
||||
thick = lambda minthick, maxthick, nlayer: np.append(np.array([1.2*10.**5]),
|
||||
np.ndarray.round(minthick + (maxthick-minthick)* np.random.rand(nlayer-1,1)
|
||||
,decimals =1))
|
||||
|
||||
sig = lambda minsig, maxsig, nlayer: np.append(np.array([0.]),
|
||||
np.ndarray.round(10.**minsig + (10.**maxsig-10.**minsig)* np.random.rand(nlayer,1)
|
||||
,decimals=3))
|
||||
|
||||
mu = lambda minmu, maxmu, nlayer: np.append(np.array([1.]),
|
||||
np.ndarray.round(minmu + (maxmu-minmu)* np.random.rand(nlayer,1)
|
||||
,decimals=1))
|
||||
|
||||
eps = lambda mineps, maxeps, nlayer: np.append(np.array([1.]),
|
||||
np.ndarray.round(mineps + (maxeps-mineps)* np.random.rand(nlayer,1)
|
||||
,decimals=1))
|
||||
|
||||
#Evaluate Impedance Z of a layer
|
||||
ImpZ = lambda f, mu, k: omega(f)*mu*mu_0/k
|
||||
|
||||
#Complex Cole-Cole Conductivity - EM utils
|
||||
PCC= lambda siginf,m,t,c,f: siginf*(1.-(m/(1.+(1j*omega(f)*t)**c)))
|
||||
|
||||
|
||||
#Converted thickness array into top of layer array
|
||||
def top(thick):
|
||||
topv= np.zeros(len(thick)+1)
|
||||
|
||||
topv[0]=-thick[0]
|
||||
|
||||
for i in range(1,len(topv),1):
|
||||
topv[i] = topv[i-1] + thick[i-1]
|
||||
|
||||
return topv
|
||||
|
||||
#Propagation Matrix and theirs inverses
|
||||
|
||||
#matrix T for transition of Up and Down components accross a layer
|
||||
T = lambda h,k: np.matrix([[np.exp(1j*k*h),0.],[0.,np.exp(-1j*k*h)]],dtype='complex_')
|
||||
|
||||
Tinv = lambda h,k: np.matrix([[np.exp(-1j*k*h),0.],[0.,np.exp(1j*k*h)]],dtype='complex_')
|
||||
|
||||
#transition of Up and Down components accross a layer
|
||||
UD_Z = lambda UD,z,zj,k : T((z-zj),k)*UD
|
||||
|
||||
|
||||
#matrix P relating Up and Down components with E and H fields
|
||||
P = lambda z: np.matrix([[1.,1,],[-1./z,1./z]],dtype='complex_')
|
||||
|
||||
Pinv = lambda z: np.matrix([[1.,-z],[1.,z]],dtype='complex_')/2.
|
||||
|
||||
|
||||
#Time Variation of E and H
|
||||
E_ZT = lambda U,D,f,t : np.exp(1j*omega(f)*t)*(U+D)
|
||||
H_ZT = lambda U,D,Z,f,t : (1./Z)*np.exp(1j*omega(f)*t)*(D-U)
|
||||
|
||||
#Plot the configuration of the problem
|
||||
def PlotConfiguration(thick,sig,eps,mu,ax,widthg,z):
|
||||
|
||||
topn = top(thick)
|
||||
widthn = np.arange(-widthg,widthg+widthg/10.,widthg/10.)
|
||||
|
||||
ax.set_ylim([z.min(),z.max()])
|
||||
ax.set_xlim([-widthg,widthg])
|
||||
|
||||
ax.set_ylabel("Depth (m)", fontsize=16.)
|
||||
ax.yaxis.tick_right()
|
||||
ax.yaxis.set_label_position("right")
|
||||
|
||||
#define filling for the different layers
|
||||
hatches=['/' , '+', 'x', '|' , '\\', '-' , 'o' , 'O' , '.' , '*' ]
|
||||
|
||||
#Write the physical properties of air
|
||||
ax.annotate(("Air, $\sigma$ =%1.0f mS/m")%(sig[0]*10**(3)),
|
||||
xy=(-widthg/2., -np.abs(z.max())/2.), xycoords='data',
|
||||
xytext=(-widthg/2., -np.abs(z.max())/2.), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\epsilon_r$= %1i")%(eps[0]),
|
||||
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
|
||||
xytext=(-widthg/2., -np.abs(z.max())/3.), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\mu_r$= %1i")%(mu[0]),
|
||||
xy=(-widthg/2., -np.abs(z.max())/3.), xycoords='data',
|
||||
xytext=(0, -np.abs(z.max())/3.), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
#Write the physical properties of the differents layers up to the (n-1)-th and fill it with pattern
|
||||
for i in range(1,len(topn)-1,1):
|
||||
if topn[i] == topn[i+1]:
|
||||
pass
|
||||
else:
|
||||
ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[i]*10**(3)),
|
||||
xy=(0., (2.*topn[i]+topn[i+1])/3), xycoords='data',
|
||||
xytext=(0., (2.*topn[i]+topn[i+1])/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\epsilon_r$= %1i")%(eps[i]),
|
||||
xy=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), xycoords='data',
|
||||
xytext=(-widthg/1.1, (2.*topn[i]+topn[i+1])/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\mu_r$= %1.2f")%(mu[i]),
|
||||
xy=(-widthg/2., (2.*topn[i]+topn[i+1])/3), xycoords='data',
|
||||
xytext=(-widthg/2., (2.*topn[i]+topn[i+1])/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.plot(widthn,topn[i]*np.ones_like(widthn),color='black')
|
||||
ax.fill_between(widthn,topn[i],topn[i+1],alpha=0.3,color="none",edgecolor='black', hatch=hatches[(i-1)%10])
|
||||
|
||||
#Write the physical properties of the n-th layer and fill it with pattern
|
||||
ax.plot(widthn,topn[-1]*np.ones_like(widthn),color='black')
|
||||
ax.fill_between(widthn,topn[-1],z.max(),alpha=0.3,color="none",edgecolor='black', hatch=hatches[(len(topn)-2)%10])
|
||||
|
||||
ax.annotate(("$\sigma$ =%3.3f mS/m")%(sig[-1]*10**(3)),
|
||||
xy=(0., (2.*topn[-1]+z.max())/3), xycoords='data',
|
||||
xytext=(0., (2.*topn[-1]+z.max())/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\epsilon_r$= %1i")%(eps[-1]),
|
||||
xy=(-widthg/1.1, (2.*topn[-1]+z.max())/3), xycoords='data',
|
||||
xytext=(-widthg/1.1, (2.*topn[-1]+z.max())/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
ax.annotate(("$\mu_r$= %1.2f")%(mu[-1]),
|
||||
xy=(-widthg/2., (2.*topn[-1]+z.max())/3), xycoords='data',
|
||||
xytext=(-widthg/2., (2.*topn[-1]+z.max())/3), textcoords='data',
|
||||
fontsize=14.)
|
||||
|
||||
#plot Trees!
|
||||
ax.annotate("",
|
||||
xy=(widthg/2., -1.*z.max()/5.), xycoords='data',
|
||||
xytext=(widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(widthg/2., -3./4.*z.max()/5.), xycoords='data',
|
||||
xytext=(widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(widthg/2., -1./2.*z.max()/5.), xycoords='data',
|
||||
xytext=(widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.2*widthg/2., -1.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.2*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.2*widthg/2., -3./4.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.2*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.2*widthg/2., -1./2.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.2*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.5*widthg/2., -1.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.5*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.2,head_length=1.2',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.5*widthg/2., -3./4.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.5*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.4,head_length=1.4',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
ax.annotate("",
|
||||
xy=(1.5*widthg/2., -1./2.*z.max()/5.), xycoords='data',
|
||||
xytext=(1.5*widthg/2., 0.), textcoords='data',
|
||||
arrowprops=dict(arrowstyle='->, head_width=1.6,head_length=1.6',color='green',linewidth=2.)
|
||||
)
|
||||
|
||||
|
||||
ax.invert_yaxis()
|
||||
|
||||
return ax
|
||||
|
||||
#Propagate Up and Down component for a certain frequency & evaluate E and H field
|
||||
|
||||
def Propagate(f,H,sig,chg,taux,c,mu,eps,n):
|
||||
|
||||
sigcm = np.zeros_like(sig,dtype='complex_')
|
||||
|
||||
for j in range(1,len(sig)):
|
||||
sigcm[j]=PCC(sig[j],chg[j],taux[j],c[j],f)
|
||||
|
||||
K = k(mu,sigcm,eps,f)
|
||||
Z = ImpZ(f,mu,K)
|
||||
|
||||
EH = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
|
||||
UD = np.matrix(np.zeros((2,n+1),dtype = 'complex_'),dtype = 'complex_')
|
||||
|
||||
UD[1,-1] = 1.
|
||||
|
||||
for i in range(-2,-(n+2),-1):
|
||||
|
||||
UD[:,i] = Tinv(H[i+1],K[i])*Pinv(Z[i])*P(Z[i+1])*UD[:,i+1]
|
||||
UD = UD/((np.abs(UD[0,:]+UD[1,:])).max())
|
||||
|
||||
for j in range(0,n+1):
|
||||
EH[:,j] = np.matrix([[1.,1,],[-1./Z[j],1./Z[j]]])*UD[:,j]
|
||||
|
||||
return UD, EH, Z ,K
|
||||
|
||||
|
||||
#Evaluate the apparent resistivity and phase for a frequency range
|
||||
def appres(F,H,sig,chg,taux,c,mu,eps,n):
|
||||
|
||||
Res = np.zeros_like(F)
|
||||
Phase = np.zeros_like(F)
|
||||
App_ImpZ= np.zeros_like(F,dtype='complex_')
|
||||
|
||||
for i in range(0,len(F)):
|
||||
|
||||
UD,EH,Z ,K = Propagate(F[i],H,sig,chg,taux,c,mu,eps,n)
|
||||
|
||||
App_ImpZ[i] = EH[0,1]/EH[1,1]
|
||||
|
||||
Res[i] = np.abs(App_ImpZ[i])**2./(mu_0*omega(F[i]))
|
||||
Phase[i] = np.angle(App_ImpZ[i], deg = True)
|
||||
|
||||
return Res,Phase
|
||||
|
||||
#Evaluate Up, Down components, E and H field, for a frequency range,
|
||||
#a discretized depth range and a time range (use to calculate envelope)
|
||||
def calculateEHzt(F,H,sig,chg,taux,c,mu,eps,n,zsample,tsample):
|
||||
|
||||
topc = top(H)
|
||||
|
||||
layer = np.zeros(len(zsample),dtype=np.int)-1
|
||||
|
||||
Exzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
|
||||
Hyzt = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
|
||||
Uz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
|
||||
Dz = np.matrix(np.zeros((len(zsample),len(tsample)),dtype = 'complex_'),dtype = 'complex_')
|
||||
UDaux = np.matrix(np.zeros((2,len(zsample)),dtype = 'complex_'),dtype = 'complex_')
|
||||
|
||||
for i in range(0,n+1,1):
|
||||
layer = layer+(zsample>=topc[i])*1
|
||||
|
||||
for j in range(0,len(F)):
|
||||
|
||||
UD,EH,Z ,K = Propagate(F[j],H,sig,chg,taux,c,mu,eps,n)
|
||||
|
||||
for p in range(0,len(zsample)):
|
||||
|
||||
UDaux[:,p] = UD_Z(UD[:,layer[p]],zsample[p],topc[layer[p]],K[layer[p]])
|
||||
|
||||
for q in range(0,len(tsample)):
|
||||
|
||||
Exzt[p,q] = Exzt[p,q] + E_ZT(UDaux[0,p],UDaux[1,p],F[j],tsample[q])/len(F)
|
||||
Hyzt[p,q] = Hyzt[p,q] + H_ZT(UDaux[0,p],UDaux[1,p],Z[layer[p]],F[j],tsample[q])/len(F)
|
||||
Uz[p,q] = Uz[p,q] + UDaux[0,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
|
||||
Dz[p,q] = Dz[p,q] + UDaux[1,p]*np.exp(1j*omega(F[j])*tsample[q])/len(F)
|
||||
|
||||
return Exzt,Hyzt,Uz,Dz,UDaux,layer
|
||||
|
||||
|
||||
#Function to Plot Apparent Resistivity and Phase
|
||||
def PlotAppRes(F,H,sig,chg,taux,c,mu,eps,n,fenvelope,PlotEnvelope):
|
||||
|
||||
Res, Phase = appres(F,H,sig,chg,taux,c,mu,eps,n)
|
||||
|
||||
fig,ax = plt.subplots(1,2,figsize=(16,10))
|
||||
|
||||
ax[0].scatter(Res,F,color='black')
|
||||
ax[0].set_xscale('Log')
|
||||
ax[0].set_yscale('Log')
|
||||
ax[0].set_xlim([10.**(np.log10(Res.min())-1.),10.**(np.log10(Res.max())+1.)])
|
||||
ax[0].set_ylim([F.min(),F.max()])
|
||||
ax[0].set_xlabel('Apparent Resistivity (Ohm*m)',fontsize=16.,color="black")
|
||||
ax[0].set_ylabel('Frequency (Hz)',fontsize=16.)
|
||||
ax[0].grid(which='major')
|
||||
|
||||
ax0 = ax[0].twiny()
|
||||
|
||||
ax0.set_xlim([0.,90.])
|
||||
ax0.set_ylim([F.min(),F.max()])
|
||||
ax0.scatter(Phase,F,color='purple')
|
||||
ax0.set_xlabel('Phase (Degrees)',fontsize=16.,color="purple")
|
||||
|
||||
zc=np.arange(-(H[1:].max()+10)*n,(H[1:].max()+10)*n,10.)
|
||||
|
||||
ax[0].tick_params(labelsize=16)
|
||||
ax[1].tick_params(labelsize=16)
|
||||
ax0.tick_params(labelsize=16)
|
||||
|
||||
if PlotEnvelope:
|
||||
|
||||
widthn=np.logspace(np.log10(Res.min())-1., np.log10(Res.max())+1., num=100, endpoint=True, base=10.0)
|
||||
fenvelope1n=np.ones(100)*fenvelope
|
||||
ax[0].plot(widthn,fenvelope1n,linestyle='dashed',color='black')
|
||||
|
||||
tc=np.arange(0.,1./fenvelope,0.01/(fenvelope))
|
||||
Exzt,Hyzt,Uz,Dz,UDaux,layer = calculateEHzt(np.array([fenvelope]),H,sig,chg,taux,c,mu,eps,n,zc,tc)
|
||||
|
||||
ax1=ax[1].twiny()
|
||||
|
||||
ax[1].tick_params(labelsize=16)
|
||||
ax1.tick_params(labelsize=16)
|
||||
|
||||
ax[1].set_xlabel('Amplitude Electric Field E (V/m)',color='blue',fontsize=16)
|
||||
|
||||
ax1.set_xlabel('Amplitude Magnetic Field H (A/m)',color='red',fontsize=16)
|
||||
|
||||
ax[1].fill_betweenx(zc,np.squeeze(np.asarray(np.real(Exzt.min(axis=1)))),
|
||||
np.squeeze(np.asarray(np.real(Exzt.max(axis=1)))),
|
||||
color='blue', alpha=0.1)
|
||||
|
||||
ax1.fill_betweenx(zc,np.squeeze(np.asarray(np.real(Hyzt.min(axis=1)))),
|
||||
np.squeeze(np.asarray(np.real(Hyzt.max(axis=1)))),
|
||||
color='red', alpha=0.1)
|
||||
|
||||
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],(1.5*np.abs(Exzt).max()),zc)
|
||||
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
|
||||
ax1.set_xlim([-1.5*np.abs(Hyzt).max(),1.5*np.abs(Hyzt).max()])
|
||||
else:
|
||||
print 'No envelop (if True, might be slow)'
|
||||
ax[1] = PlotConfiguration(H,sig,eps,mu,ax[1],1.,zc)
|
||||
ax[1].get_xaxis().set_ticks([])
|
||||
|
||||
plt.show()
|
||||
|
||||
#Interactive MT for Notebook
|
||||
def PlotAppRes3LayersInteract(h1,h2,sigl1,sigl2,sigl3,mul1,mul2,mul3,epsl1,epsl2,epsl3,PlotEnvelope,F_Envelope):
|
||||
|
||||
frangn=frange(-5,5,100.)
|
||||
sig3= np.array([0.,0.001,0.1, 0.001])
|
||||
thick3 = np.array([120000.,50.,50.])
|
||||
eps3=np.array([1.,1.,1.,1])
|
||||
mu3=np.array([1.,1.,1.,1])
|
||||
chg3=np.array([0.,0.1,0.,0.2])
|
||||
chg3_0=np.array([0.,0.1,0.,0.])
|
||||
taux3=np.array([0.,0.1,0.,0.1])
|
||||
c3=np.array([1.,1.,1.,1.])
|
||||
|
||||
sig3[1]=sigl1
|
||||
sig3[1]=10.**sig3[1]
|
||||
sig3[2]=sigl2
|
||||
sig3[2]=10.**sig3[2]
|
||||
sig3[3]=sigl3
|
||||
sig3[3]=10.**sig3[3]
|
||||
mu3[1]=mul1
|
||||
mu3[2]=mul2
|
||||
mu3[3]=mul3
|
||||
eps3[1]=epsl1
|
||||
eps3[2]=epsl2
|
||||
eps3[3]=epsl3
|
||||
thick3[1]=h1
|
||||
thick3[2]=h2
|
||||
|
||||
PlotAppRes(frangn,thick3,sig3,chg3_0,taux3,c3,mu3,eps3,3,F_Envelope,PlotEnvelope)
|
||||
|
||||
|
||||
def run(n,plotIt=True):
|
||||
# something to make a plot
|
||||
|
||||
F = frange(-5.,5.,20)
|
||||
H = thick(50.,100.,n)
|
||||
sign = sig(-5.,0.,n)
|
||||
mun = mu(1.,2.,n)
|
||||
epsn = eps(1.,9.,n)
|
||||
chg = np.zeros_like(sign)
|
||||
taux = np.zeros_like(sign)
|
||||
c = np.zeros_like(sign)
|
||||
|
||||
Res, Phase = appres(F,H,sign,chg,taux,c,mun,epsn,n)
|
||||
|
||||
if plotIt:
|
||||
|
||||
PlotAppRes(F, H, sign, chg, taux, c, mun, epsn, n, fenvelope=1000., PlotEnvelope=True)
|
||||
|
||||
return Res, Phase
|
||||
|
||||
if __name__ == '__main__':
|
||||
run(3)
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
@@ -3,12 +3,15 @@
|
||||
##### AUTOIMPORTS #####
|
||||
import DC_Analytic_Dipole
|
||||
import DC_Forward_PseudoSection
|
||||
import DC_PseudoSection_Simulation
|
||||
import EM_FDEM_1D_Inversion
|
||||
import EM_FDEM_Analytic_MagDipoleWholespace
|
||||
import EM_FDEM_SusEffects
|
||||
import EM_Schenkel_Morrison_Casing
|
||||
import EM_TDEM_1D_Inversion
|
||||
import FLOW_Richards_1D_Celia1990
|
||||
import Forward_BasicDirectCurrent
|
||||
import Inversion_IRLS
|
||||
import Inversion_Linear
|
||||
import Mesh_Basic_PlotImage
|
||||
import Mesh_Basic_Types
|
||||
@@ -17,10 +20,12 @@ import Mesh_QuadTree_Creation
|
||||
import Mesh_QuadTree_FaceDiv
|
||||
import Mesh_QuadTree_HangingNodes
|
||||
import Mesh_Tensor_Creation
|
||||
import MT_1D_analytic_nlayer_Earth
|
||||
import MT_1D_ForwardAndInversion
|
||||
import MT_3D_Foward
|
||||
import sphereElectrostatic_example
|
||||
|
||||
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_ForwardAndInversion", "MT_3D_Foward"]
|
||||
__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "DC_PseudoSection_Simulation", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "EM_FDEM_SusEffects", "EM_Schenkel_Morrison_Casing", "EM_TDEM_1D_Inversion", "FLOW_Richards_1D_Celia1990", "Forward_BasicDirectCurrent", "Inversion_IRLS", "Inversion_Linear", "Mesh_Basic_PlotImage", "Mesh_Basic_Types", "Mesh_Operators_CahnHilliard", "Mesh_QuadTree_Creation", "Mesh_QuadTree_FaceDiv", "Mesh_QuadTree_HangingNodes", "Mesh_Tensor_Creation", "MT_1D_analytic_nlayer_Earth", "MT_1D_ForwardAndInversion", "MT_3D_Foward", "sphereElectrostatic_example"]
|
||||
|
||||
##### AUTOIMPORTS #####
|
||||
|
||||
|
||||
@@ -0,0 +1,775 @@
|
||||
from scipy.constants import epsilon_0
|
||||
import matplotlib.pyplot as plt
|
||||
import matplotlib.colors as colors
|
||||
import numpy as np
|
||||
from SimPEG.Utils import ndgrid, mkvc
|
||||
|
||||
'''
|
||||
Authors: Thibaut Astic, Lindsey Heagy, Sanna Tyrvainen, Ronghua Peng
|
||||
Date: December 2015
|
||||
|
||||
This code defines function to resolve analytically the electrostatic sphere problem.
|
||||
We first define a problem configuration, with a conductive or resistive sphere in a
|
||||
wholespace background.
|
||||
We then calculate the potential, then the electric field, then the current density and
|
||||
finally the charges accumulation.
|
||||
|
||||
Several plotting functions are defined for data visualisation.
|
||||
|
||||
|
||||
'''
|
||||
|
||||
# Plot options
|
||||
ftsize_title = 18 #font size for titles
|
||||
ftsize_axis = 14 #font size for axis ticks
|
||||
ftsize_label = 14 #font size for axis labels
|
||||
|
||||
# Radius function, useful sigma ratio, and log scale converter
|
||||
r = lambda x,y,z: np.sqrt(x**2.+y**2.+z**2.)
|
||||
sigf = lambda sig0,sig1: (sig1-sig0)/(sig1+2.*sig0)
|
||||
|
||||
#tools to convert log conductivity in conductivity
|
||||
def conductivity_log_wrapper(log_sig0,log_sig1):
|
||||
sig0 = 10.**log_sig0
|
||||
sig1 = 10.**log_sig1
|
||||
|
||||
return sig0,sig1
|
||||
|
||||
# Examples
|
||||
#Plot the configuration. Label=False is used to generate a general case figure
|
||||
def get_Setup(XYZ,sig0,sig1,R,E0,ax,label,colorsphere):
|
||||
'''
|
||||
XYZ: ndgrid
|
||||
sig0: conductivity of the background
|
||||
sig1: conductivity of the sphere
|
||||
R: radius of the sphere
|
||||
E0: Amplitude of the uniform electrostatic field
|
||||
ax: ax where to plot the configuration
|
||||
label: True: plot real values, False: plot general case
|
||||
colorsphere: color of the sphere, format [x,x,x]
|
||||
'''
|
||||
|
||||
xplt = np.linspace(-R, R, num=100)
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
dx = xr[1]-xr[0]
|
||||
top = np.sqrt(R**2-xplt**2)
|
||||
bot = -np.sqrt(R**2-xplt**2)
|
||||
|
||||
if R != 0:
|
||||
ax.plot(xplt, top, xplt, bot, color=colorsphere,linewidth=1.5)
|
||||
ax.fill_between(xplt,bot,top,color=colorsphere,alpha=0.5 )
|
||||
ax.arrow(0.,0.,np.sqrt(2.)*R/2.,np.sqrt(2.)*R/2.,head_width=0.,head_length=0.)
|
||||
|
||||
if label:
|
||||
ax.annotate(("$\sigma_1$=%3.3f mS/m")%(sig1*10.**(3.)),
|
||||
xy=(0.,-R/2.), xycoords='data',
|
||||
xytext=(0.,-R/2.), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
|
||||
xy=(0.,-1.5*R), xycoords='data',
|
||||
xytext=(0.,-1.5*R), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0),
|
||||
xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
|
||||
xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.annotate(('$R$ = %1i m')%(R),
|
||||
xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
|
||||
xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
|
||||
else:
|
||||
ax.set_xticklabels([])
|
||||
ax.set_yticklabels([])
|
||||
ax.text(-1.,-np.sqrt(R)/2.-10.,'$\sigma_1$',fontsize=14)
|
||||
ax.text(-0.05,-R-10,'$\sigma_0$',fontsize=14)
|
||||
ax.annotate(('$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m'),
|
||||
xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
|
||||
xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.annotate(('$R$'),
|
||||
xy=(R/4.+(xr[1]-xr[0]),R/4.), xycoords='data',
|
||||
xytext=(R/4.+(xr[1]-xr[0]),R/4.), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.set_xlabel('x',fontsize=12)
|
||||
ax.set_ylabel('y',fontsize=12)
|
||||
|
||||
else:
|
||||
if label:
|
||||
ax.annotate(("$\sigma_0$= %3.3f mS/m")%(sig0*10.**(3.)),
|
||||
xy=(0.,-1.5*R), xycoords='data',
|
||||
xytext=(0.,-1.5*R), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.annotate(('$\mathbf{E_0} = %1i \mathbf{\hat{x}}$ V/m')%(E0),
|
||||
xy=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), xycoords='data',
|
||||
xytext=(xr.min()+np.abs(xr.max()-xr.min())/20.,0), textcoords='data',
|
||||
fontsize=14.)
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
|
||||
else:
|
||||
ax.set_xticklabels([])
|
||||
ax.set_yticklabels([])
|
||||
ax.text(-0.05,-10,'$\sigma_0$',fontsize=14)
|
||||
ax.text(xr.min()+np.abs(xr.max()-xr.min())/20., 0, '$\mathbf{E_0} = E_0 \mathbf{\hat{x}}$ V/m', fontsize=14)
|
||||
ax.set_xlabel('x',fontsize=12)
|
||||
ax.set_ylabel('y',fontsize=12)
|
||||
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
[ax.arrow(xr.min(),_,np.abs(xr.max()-xr.min())/20.,0.,head_width=5.,head_length=2.,color='k') for _ in np.linspace(yr.min(),yr.max(),num=10)]
|
||||
ax.patch.set_facecolor([0.4,0.7,0.4])
|
||||
ax.patch.set_alpha(0.2)
|
||||
|
||||
ax.set_aspect('equal')
|
||||
|
||||
|
||||
|
||||
return ax
|
||||
|
||||
def get_Conductivity(XYZ,sig0,sig1,R):
|
||||
'''
|
||||
Define the conductivity for each point of the space
|
||||
'''
|
||||
x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
|
||||
r_view=r(x,y,z)
|
||||
|
||||
ind0= (r_view>R)
|
||||
ind1= (r_view<=R)
|
||||
|
||||
assert (ind0 + ind1).all(), 'Some indicies not included'
|
||||
|
||||
Sigma = np.zeros_like(x)
|
||||
|
||||
Sigma[ind0] = sig0
|
||||
Sigma[ind1] = sig1
|
||||
|
||||
return Sigma
|
||||
|
||||
|
||||
def get_Potential(XYZ,sig0,sig1,R,E0):
|
||||
|
||||
'''
|
||||
Function that returns the total, the primary and the secondary potentials, assumes an x-oriented inducing field and that the sphere is at the origin
|
||||
:input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
|
||||
'''
|
||||
|
||||
x,y,z = XYZ[:,0],XYZ[:,1],XYZ[:,2]
|
||||
|
||||
sig_cur = sigf(sig0,sig1)
|
||||
|
||||
r_cur = r(x,y,z) # current radius
|
||||
|
||||
ind0 = (r_cur > R)
|
||||
ind1 = (r_cur <= R)
|
||||
|
||||
assert (ind0 + ind1).all(), 'Some indicies not included'
|
||||
|
||||
Vt = np.zeros_like(x)
|
||||
Vp = np.zeros_like(x)
|
||||
Vs = np.zeros_like(x)
|
||||
|
||||
Vt[ind0] = -E0*x[ind0]*(1.-sig_cur*R**3./r_cur[ind0]**3.) # total potential outside the sphere
|
||||
Vt[ind1] = -E0*x[ind1]*3.*sig0/(sig1+2.*sig0) # inside the sphere
|
||||
|
||||
|
||||
Vp = - E0*x # primary potential
|
||||
|
||||
Vs = Vt - Vp # secondary potential
|
||||
|
||||
return Vt,Vp,Vs
|
||||
|
||||
#plot the primary potential on ax
|
||||
def Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
Pplot = ax.pcolor(xr,yr,Vp.reshape(xr.size,yr.size))
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_title('Primary Potential',fontsize=ftsize_title)
|
||||
cb = plt.colorbar(Pplot,ax=ax)
|
||||
cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_aspect('equal')
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
|
||||
return ax
|
||||
|
||||
#plot the total potential on ax
|
||||
def Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
|
||||
Pplot = ax.pcolor(xr,yr,Vt.reshape(xr.size,yr.size))
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_title('Total Potential',fontsize=ftsize_title)
|
||||
cb = plt.colorbar(Pplot,ax=ax)
|
||||
cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_aspect('equal')
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
|
||||
return ax
|
||||
|
||||
#plot the secondary potential on ax
|
||||
def Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Vt,Vp,Vs = get_Potential(XYZ,sig0,sig1,R,E0)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
Pplot = ax.pcolor(xr,yr,Vs.reshape(xr.size,yr.size))
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_title('Secondary Potential',fontsize=ftsize_title)
|
||||
cb = plt.colorbar(Pplot,ax=ax)
|
||||
cb.set_label(label= 'Potential ($V$)',size=ftsize_label)
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_aspect('equal')
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
|
||||
return ax
|
||||
|
||||
|
||||
def get_ElectricField(XYZ,sig0,sig1,R,E0):
|
||||
'''
|
||||
Function that returns the total, the primary and the secondary electric fields,
|
||||
input: grid, outer sigma, inner sigma, radius of the sphere, strength of the electric field
|
||||
'''
|
||||
|
||||
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
|
||||
|
||||
r_cur=r(x,y,z) # current radius
|
||||
|
||||
ind0= (r_cur>R)
|
||||
ind1= (r_cur<=R)
|
||||
|
||||
assert (ind0 + ind1).all(), 'Some indicies not included'
|
||||
|
||||
Ep = np.zeros(shape=(len(x),3))
|
||||
Ep[:,0] = E0
|
||||
|
||||
Et = np.zeros(shape=(len(x),3))
|
||||
|
||||
Et[ind0,0] = E0 + E0*R**3./(r_cur[ind0]**5.)*sigf(sig0,sig1)*(2.*x[ind0]**2.-y[ind0]**2.-z[ind0]**2.);
|
||||
Et[ind0,1] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*y[ind0]*sigf(sig0,sig1);
|
||||
Et[ind0,2] = E0*R**3./(r_cur[ind0]**5.)*3.*x[ind0]*z[ind0]*sigf(sig0,sig1);
|
||||
|
||||
Et[ind1,0] = 3.*sig0/(sig1+2.*sig0)*E0;
|
||||
Et[ind1,1] = 0.;
|
||||
Et[ind1,2] = 0.;
|
||||
|
||||
Es = Et - Ep
|
||||
|
||||
return Et, Ep, Es
|
||||
|
||||
#plot the total electric field on ax
|
||||
def Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
EtXr = Et[:,0].reshape(xr.size, yr.size)
|
||||
EtYr = Et[:,1].reshape(xr.size, yr.size)
|
||||
EtAmp = np.sqrt(Et[:,0]**2+Et[:,1]**2 + Et[:,2]**2).reshape(xr.size, yr.size)
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_aspect('equal')
|
||||
|
||||
Eplot = ax.pcolor(xr,yr,EtAmp)
|
||||
cb = plt.colorbar(Eplot,ax=ax)
|
||||
cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.streamplot(xr,yr,EtXr,EtYr,color='gray',linewidth=2.,density=0.75)#angles='xy',scale_units='xy',scale=0.05)
|
||||
ax.set_title('Total Field',fontsize=ftsize_title)
|
||||
|
||||
|
||||
return ax
|
||||
|
||||
#plot the secondary electric field on ax
|
||||
def Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
EsXr = Es[:,0].reshape(xr.size, yr.size)
|
||||
EsYr = Es[:,1].reshape(xr.size, yr.size)
|
||||
EsAmp = np.sqrt(Es[:,0]**2+Es[:,1]**2+Es[:,2]**2).reshape(xr.size, yr.size)
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize = ftsize_label)
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_aspect('equal')
|
||||
|
||||
Eplot = ax.pcolor(xr,yr,EsAmp)
|
||||
cb = plt.colorbar(Eplot,ax=ax)
|
||||
cb.set_label(label= 'Amplitude ($V/m$)',size=ftsize_label) #weight='bold')
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.streamplot(xr,yr,EsXr,EsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=0.05)
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_title('Secondary Field',fontsize=ftsize_title)
|
||||
|
||||
return ax
|
||||
|
||||
|
||||
def get_Current(XYZ,sig0,sig1,R,Et,Ep,Es):
|
||||
'''
|
||||
Function that returns the total, the primary and the secondary current densities,
|
||||
:input: grid, outer sigma, inner sigma, radius of the sphere, total, the primary and the seconadry electric fields,
|
||||
'''
|
||||
|
||||
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
|
||||
|
||||
r_cur=r(x,y,z)
|
||||
|
||||
ind0= (r_cur>R)
|
||||
ind1= (r_cur<=R)
|
||||
|
||||
assert (ind0 + ind1).all(), 'Some indicies not included'
|
||||
|
||||
Jt = np.zeros(shape=(len(x),3))
|
||||
J0 = np.zeros(shape=(len(x),3))
|
||||
Js = np.zeros(shape=(len(x),3))
|
||||
|
||||
|
||||
Jp = sig0*Ep
|
||||
|
||||
Jt[ind0,:] = sig0*Et[ind0,:]
|
||||
Jt[ind1,:] = sig1*Et[ind1,:]
|
||||
|
||||
Js[ind0,:] = sig0*(Et[ind0,:]-Ep[ind0,:])
|
||||
Js[ind1,:] = sig1*Et[ind1,:]-sig0*Ep[ind1,:]
|
||||
|
||||
return Jt,Jp,Js
|
||||
|
||||
#plot the total currents density on ax
|
||||
def Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
|
||||
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
JtXr = Jt[:,0].reshape(xr.size, yr.size)
|
||||
JtYr = Jt[:,1].reshape(xr.size, yr.size)
|
||||
JtAmp = np.sqrt(Jt[:,0]**2+Jt[:,1]**2+Jt[:,2]**2).reshape(xr.size, yr.size)
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_aspect('equal')
|
||||
|
||||
Jplot = ax.pcolor(xr,yr,JtAmp.reshape(xr.size,yr.size))
|
||||
cb = plt.colorbar(Jplot,ax=ax)
|
||||
cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.streamplot(xr,yr,JtXr,JtYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
|
||||
ax.set_title('Total Current Density',fontsize=ftsize_title)
|
||||
|
||||
return ax
|
||||
|
||||
|
||||
#plot the secondary currents density on ax
|
||||
def Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
Et,Ep,Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
|
||||
Jt,Jp,Js = get_Current(XYZ,sig0,sig1,R,Et,Ep,Es)
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
JsXr = Js[:,0].reshape(xr.size, yr.size)
|
||||
JsYr = Js[:,1].reshape(xr.size, yr.size)
|
||||
JsAmp = np.sqrt(Js[:,1]**2+Js[:,0]**2+Jt[:,2]**2).reshape(xr.size,yr.size)
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_aspect('equal')
|
||||
|
||||
Jplot = ax.pcolor(xr,yr,JsAmp.reshape(xr.size,yr.size))
|
||||
cb = plt.colorbar(Jplot,ax=ax)
|
||||
cb.set_label(label= 'Current Density ($A/m^2$)',size=ftsize_label) #weight='bold')
|
||||
cb.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.streamplot(xr,yr,JsXr,JsYr,color='gray',linewidth=2.,density=0.75)#,angles='xy',scale_units='xy',scale=1)
|
||||
ax.set_title('Secondary Current Density',fontsize=ftsize_title)
|
||||
|
||||
return ax
|
||||
|
||||
|
||||
def get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep):
|
||||
'''
|
||||
Function that returns the charges accumulation at the background/sphere interface,
|
||||
:input: grid, outer sigma, inner sigma, radius of the sphere, total and the primary electric fields,
|
||||
'''
|
||||
|
||||
x,y,z= XYZ[:,0], XYZ[:,1], XYZ[:,2]
|
||||
|
||||
dx = x[1]-x[0]
|
||||
|
||||
r_cur=r(x,y,z)
|
||||
|
||||
ind0 = (r_cur > R)
|
||||
ind1 = (r_cur < R)
|
||||
ind2 = ((r_cur < (R+dx/2)) & (r_cur > (R-dx/2)) )
|
||||
|
||||
assert (ind0 + ind1 + ind2).all(), 'Some indicies not included'
|
||||
|
||||
rho = np.zeros_like(x)
|
||||
|
||||
rho[ind0] = 0
|
||||
rho[ind1] = 0
|
||||
rho[ind2] = epsilon_0*3.*Ep[ind2,0]*sigf(sig0,sig1)*x[ind2]/(np.sqrt(x[ind2]**2.+y[ind2]**2.))
|
||||
|
||||
return rho
|
||||
|
||||
#Plot charges density on ax
|
||||
def Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax):
|
||||
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
xcirc = xr[np.abs(xr) <= R]
|
||||
|
||||
Et, Ep, Es = get_ElectricField(XYZ,sig0,sig1,R,E0)
|
||||
rho = get_ChargesDensity(XYZ,sig0,sig1,R,Et,Ep)
|
||||
|
||||
ax.set_xlim([xr.min(),xr.max()])
|
||||
ax.set_ylim([yr.min(),yr.max()])
|
||||
ax.set_aspect('equal')
|
||||
Cplot = ax.pcolor(xr,yr,rho.reshape(xr.size, yr.size))
|
||||
cb1 = plt.colorbar(Cplot,ax=ax)
|
||||
cb1.set_label(label= 'Charge Density ($C/m^2$)',size=ftsize_label) #weight='bold')
|
||||
cb1.ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.plot(xcirc,np.sqrt(R**2-xcirc**2),'--k',xcirc,-np.sqrt(R**2-xcirc**2),'--k')
|
||||
ax.set_ylabel('Y coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.set_xlabel('X coordinate ($m$)',fontsize=ftsize_label)
|
||||
ax.tick_params(labelsize=ftsize_axis)
|
||||
ax.set_title('Charges Density', fontsize=ftsize_title)
|
||||
|
||||
return ax
|
||||
|
||||
def MN_Potential_total(sig0,sig1,R,E0,start,end,nbmp,mn):
|
||||
|
||||
'''
|
||||
Function that return array of midpoints electrodes, electrodes positions,
|
||||
potentials differences for total and secondary potentials fields, unormalized and
|
||||
normalized to electrodes distances.
|
||||
sig0: background conductivity
|
||||
sig1: sphere conductivity
|
||||
R: Sphere's radius
|
||||
E0: uniform E field value
|
||||
start: start point for the profile start.shape = (2,)
|
||||
end: end point for the profile end.shape = (2,)
|
||||
nbmp: number of dipoles
|
||||
mn: Space between the M and N electrodes
|
||||
'''
|
||||
|
||||
#D: total distance from start to end
|
||||
D = np.sqrt((start[0]-end[0])**2.+(start[1]-end[1])**2.)
|
||||
|
||||
#MP: dipoles'midpoint positions (x,y)
|
||||
MP = np.zeros(shape=(nbmp,2))
|
||||
MP[:,0] = np.linspace(start[0],end[0],nbmp)
|
||||
MP[:,1] = np.linspace(start[1],end[1],nbmp)
|
||||
|
||||
#Dipoles'Electrodes positions around each midpoints
|
||||
EL = np.zeros(shape=(2*nbmp,2))
|
||||
for n in range(0,len(EL),2):
|
||||
EL[n,0] = MP[n/2,0] - ((end[0]-start[0])/D)*mn/2.
|
||||
EL[n+1,0] = MP[n/2,0] + ((end[0]-start[0])/D)*mn/2.
|
||||
EL[n,1] = MP[n/2,1] - ((end[1]-start[1])/D)*mn/2.
|
||||
EL[n+1,1] = MP[n/2,1] + ((end[1]-start[1])/D)*mn/2.
|
||||
|
||||
VtEL = np.zeros(2*nbmp) #Total Potential (Vt-) at each electrode (-EL)
|
||||
VsEL = np.zeros(2*nbmp) #Secondary Potential (Vt-) at each electrode (-EL)
|
||||
dVtMP = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP)
|
||||
dVtMPn = np.zeros(nbmp) #Diffence (d-) of Total Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
|
||||
dVsMP = np.zeros(nbmp) #Diffence (d-) of Secondaty Potential (Vt-) at each dipole (-MP)
|
||||
dVsMPn = np.zeros(nbmp) #Diffence (d-) of Secondary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
|
||||
dVpMP = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP)
|
||||
dVpMPn = np.zeros(nbmp) #Diffence (d-) of Primary Potential (Vt-) at each dipole (-MP) normalized for the mn spacing (n)
|
||||
|
||||
#Computing VtEL
|
||||
for m in range(0,2*nbmp):
|
||||
if (r(EL[m,0],EL[m,1],0) > R):
|
||||
VtEL[m] = -E0*EL[m,0]*(1.-sigf(sig0,sig1)*R**3./r(EL[m,0],EL[m,1],0)**3.)
|
||||
else:
|
||||
VtEL[m] = -E0*EL[m,0]*3.*sig0/(sig1+2.*sig0)
|
||||
|
||||
#Computing VsEL
|
||||
VsEL = VtEL + E0*EL[:,0]
|
||||
|
||||
#Computing dVtMP, dVsMP
|
||||
for p in range(0,nbmp):
|
||||
dVtMP[p] = VtEL[2*p]-VtEL[2*p+1]
|
||||
dVtMPn[p] = dVtMP[p]/mn
|
||||
dVsMP[p] = VsEL[2*p]-VsEL[2*p+1]
|
||||
dVsMPn[p] = dVsMP[p]/mn
|
||||
|
||||
return MP,EL,dVtMP,dVtMPn,dVsMP,dVsMPn
|
||||
|
||||
#Compare the DC response of two configurations
|
||||
def two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,nb_dipole,electrode_spacing,PlotOpt):#,linearcolor):
|
||||
|
||||
#Define the mesh
|
||||
xr,yr,zr = np.unique(XYZ[:,0]),np.unique(XYZ[:,1]),np.unique(XYZ[:,2])
|
||||
|
||||
#Defining the Profile
|
||||
start = np.array([xstart,ystart])
|
||||
end = np.array([xend,yend])
|
||||
|
||||
#Calculating the data from the defined survey line for Configuration 0 and 1
|
||||
MP0,EL0,VtdMP0,VtdMPn0,VsdMP0,VsdMPn0 = MN_Potential_total(sig0,sig1,R0,E0,start,end,nb_dipole,electrode_spacing)
|
||||
MP1,EL1,VtdMP1,VtdMPn1,VsdMP1,VsdMPn1 = MN_Potential_total(sig0,sig2,R1,E0,start,end,nb_dipole,electrode_spacing)
|
||||
|
||||
|
||||
# Initializing the figure
|
||||
fig = plt.figure(figsize=(20,20))
|
||||
ax0 = plt.subplot2grid((20,12), (0, 0),colspan=6,rowspan=6)
|
||||
ax1 = plt.subplot2grid((20,12), (0, 6),colspan=6,rowspan=6)
|
||||
ax2 = plt.subplot2grid((20,12), (16, 2), colspan=9,rowspan=4)
|
||||
ax3 = plt.subplot2grid((20,12), (8, 0),colspan=6,rowspan=6)
|
||||
ax4 = plt.subplot2grid((20,12), (8, 6),colspan=6,rowspan=6)
|
||||
|
||||
#Plotting the Configuration 0
|
||||
ax0 = get_Setup(XYZ,sig0,sig1,R0,E0,ax0,True,[0.6,0.1,0.1])
|
||||
|
||||
#Plotting the Configuration 1
|
||||
ax1 = get_Setup(XYZ,sig0,sig2,R1,E0,ax1,True,[0.1,0.1,0.6])
|
||||
|
||||
#Plotting the Data (Legends)
|
||||
ax2.set_title('Potential Differences',fontsize=ftsize_title)
|
||||
ax2.set_ylabel('Potential difference ($V$)',fontsize=ftsize_label)
|
||||
ax2.set_xlabel('Distance from start point ($m$)',fontsize=ftsize_label)
|
||||
ax2.tick_params(labelsize=ftsize_axis)
|
||||
ax2.grid()
|
||||
|
||||
if PlotOpt == 'Total':
|
||||
ax3= Plot_Total_Potential(XYZ,sig0,sig1,R0,E0,ax3)
|
||||
ax4= Plot_Total_Potential(XYZ,sig0,sig2,R1,E0,ax4)
|
||||
|
||||
#Plot the Data (from Configuration 0)
|
||||
gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VtdMP0
|
||||
,marker='o',color='blue',linewidth=3.,label ='Left Model Response' )
|
||||
|
||||
#Plot the Data (from Configuration 1)
|
||||
gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VtdMP1
|
||||
,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
|
||||
ax2.legend(('Left Model Response','Right Model Response'),loc=4)
|
||||
|
||||
elif PlotOpt == 'Secondary':
|
||||
#plot the secondary potentials
|
||||
ax3= Plot_Secondary_Potential(XYZ,sig0,sig1,R0,E0,ax3)
|
||||
ax4= Plot_Secondary_Potential(XYZ,sig0,sig2,R1,E0,ax4)
|
||||
|
||||
#Plot the data(from configuration 0)
|
||||
gphy0 = ax2.plot(np.sqrt((MP0[0,0]-MP0[:,0])**2+(MP0[:,1]-MP0[0,1])**2),VsdMP0,color='blue'
|
||||
,marker='o',linewidth=3.,label ='Left Model Response' )
|
||||
|
||||
|
||||
#Plot the Data (from Configuration 1)
|
||||
gphy1 = ax2.plot(np.sqrt((MP1[0,0]-MP1[:,0])**2+(MP1[:,1]-MP1[0,1])**2),VsdMP1
|
||||
,marker='o',color='red',linewidth=2.,label ='Right Model Response' )
|
||||
ax2.legend(('Left Model Response','Right Model Response'),loc=4 )
|
||||
|
||||
else:
|
||||
print('What dont you get? Total or Secondary?')
|
||||
|
||||
#Legends
|
||||
ax3.plot(MP0[:,0],MP0[:,1],color='gray')
|
||||
Dip_Midpoint0 = ax3.scatter(MP0[:,0],MP0[:,1],color='black')
|
||||
Electrodes0 = ax3.scatter(EL0[:,0],EL0[:,1],color='red')
|
||||
ax3.legend([Dip_Midpoint0,Electrodes0], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
|
||||
|
||||
ax4.plot(MP1[:,0],MP1[:,1],color='gray')
|
||||
Dip_Midpoint1 = ax4.scatter(MP1[:,0],MP1[:,1],color='black')
|
||||
Electrodes1 = ax4.scatter(EL1[:,0],EL1[:,1],color='red')
|
||||
ax4.legend([Dip_Midpoint1,Electrodes1], ["Dipole Midpoint", "Electrodes"],scatterpoints=1)
|
||||
|
||||
return fig
|
||||
|
||||
#Function to visualise and compare any two meaningful plots for the sphere in a uniform backgound with an unifom Electric Field
|
||||
def interact_conductiveSphere(R,log_sig0,log_sig1,Figure1a,Figure1b,Figure2a,Figure2b):
|
||||
|
||||
sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
|
||||
E0 = 1. # inducing field strength in V/m
|
||||
n = 100 #level of discretisation
|
||||
xr = np.linspace(-200., 200., n) # X-axis discretization
|
||||
yr = xr.copy() # Y-axis discretization
|
||||
zr = np.r_[0] # identical to saying `zr = np.array([0])`
|
||||
XYZ = ndgrid(xr,yr,zr) # Space Definition
|
||||
|
||||
fig, ax = plt.subplots(1,2,figsize=(18,6))
|
||||
|
||||
#Setup figure 1 with options Configuration, Total or Secondary,
|
||||
#then Potential, ElectricField, Current Density or Charges Density
|
||||
if Figure1a == 'Configuration':
|
||||
ax[0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0],True,[0.1,0.1,0.6])
|
||||
|
||||
elif Figure1a == 'Total':
|
||||
|
||||
if Figure1b == 'Potential':
|
||||
ax[0] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'ElectricField':
|
||||
ax[0] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'CurrentDensity':
|
||||
ax[0] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'ChargesDensity':
|
||||
ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1a == 'Secondary':
|
||||
|
||||
if Figure1b == 'Potential':
|
||||
ax[0] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'ElectricField':
|
||||
ax[0] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'CurrentDensity':
|
||||
ax[0] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
elif Figure1b == 'ChargesDensity':
|
||||
ax[0] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[0])
|
||||
|
||||
|
||||
if Figure1a== 'Configuration':
|
||||
ax[1] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
print 'While figure1 is plotting Configuration, figure2 plots the primary field'
|
||||
|
||||
elif Figure2a == 'Total':
|
||||
if Figure2b == 'Potential':
|
||||
ax[1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'ElectricField':
|
||||
ax[1] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'CurrentDensity':
|
||||
ax[1]=Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'ChargesDensity':
|
||||
ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
|
||||
elif Figure2a == 'Secondary':
|
||||
if Figure2b == 'Potential':
|
||||
ax[1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'ElectricField':
|
||||
ax[1] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'CurrentDensity':
|
||||
ax[1] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
elif Figure2b == 'ChargesDensity':
|
||||
ax[1] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1])
|
||||
|
||||
plt.tight_layout(True)
|
||||
plt.show()
|
||||
|
||||
#Interactive Visualisation of the responses of two configurations to a (pseudo) DC resistivity survey
|
||||
def interactive_two_configurations_comparison(log_sig0,log_sig1,log_sig2,R0,R1,xstart,ystart,xend,yend,dipole_number,electrode_spacing,matching_spheres_example):
|
||||
|
||||
sig0,sig1 = conductivity_log_wrapper(log_sig0,log_sig1)
|
||||
sig2 = 10.**log_sig2
|
||||
E0 = 1. # inducing field strength in V/m
|
||||
n = 100 #level of discretisation
|
||||
xr = np.linspace(-200., 200., n) # X-axis discretization
|
||||
yr = xr.copy() # Y-axis discretization
|
||||
zr = np.r_[0] # identical to saying `zr = np.array([0])`
|
||||
XYZ = ndgrid(xr,yr,zr) # Space Definition
|
||||
PlotOpt = 'Total'
|
||||
|
||||
if matching_spheres_example:
|
||||
sig0 = 10.**(-3)
|
||||
sig1 = 10.**(-2)
|
||||
sig2 = 1.310344828 * 10**(-3)
|
||||
R0 = 20.
|
||||
R1 = 40.
|
||||
|
||||
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
|
||||
|
||||
else:
|
||||
two_configurations_comparison(XYZ,sig0,sig1,sig2,R0,R1,E0,xstart,ystart,xend,yend,dipole_number,electrode_spacing,PlotOpt)
|
||||
|
||||
plt.tight_layout(True)
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
sig0 = -3. # conductivity of the wholespace
|
||||
sig1 = -1. # conductivity of the sphere
|
||||
sig0, sig1 = conductivity_log_wrapper(sig0,sig1)
|
||||
R = 50. # radius of the sphere
|
||||
E0 = 1. # inducing field strength
|
||||
n = 100 #level of discretisation
|
||||
xr = np.linspace(-2.*R, 2.*R, n) # X-axis discretization
|
||||
yr = xr.copy() # Y-axis discretization
|
||||
zr = np.r_[0] # identical to saying `zr = np.array([0])`
|
||||
XYZ = ndgrid(xr,yr,zr) # Space Definition
|
||||
|
||||
fig, ax = plt.subplots(2,5,figsize=(50,10))
|
||||
ax[0,0] = get_Setup(XYZ,sig0,sig1,R,E0,ax[0,0],True,[0.6,0.1,0.1])
|
||||
ax[1,0] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[1,0])
|
||||
ax[0,1] = Plot_Total_Potential(XYZ,sig0,sig1,R,E0,ax[0,1])
|
||||
ax[1,1] = Plot_Secondary_Potential(XYZ,sig0,sig1,R,E0,ax[1,1])
|
||||
ax[0,2] = Plot_Total_ElectricField(XYZ,sig0,sig1,R,E0,ax[0,2])
|
||||
ax[1,2] = Plot_Secondary_ElectricField(XYZ,sig0,sig1,R,E0,ax[1,2])
|
||||
ax[0,3] = Plot_Total_Currents(XYZ,sig0,sig1,R,E0,ax[0,3])
|
||||
ax[1,3] = Plot_Secondary_Currents(XYZ,sig0,sig1,R,E0,ax[1,3])
|
||||
ax[0,4] = Plot_Primary_Potential(XYZ,sig0,sig1,R,E0,ax[0,4])
|
||||
ax[1,4] = Plot_ChargesDensity(XYZ,sig0,sig1,R,E0,ax[1,4])
|
||||
|
||||
|
||||
plt.show()
|
||||
|
||||
+29
-16
@@ -759,15 +759,27 @@ class PolyMap(IdentityMap):
|
||||
|
||||
m = [\sigma_1, \sigma_2, c]
|
||||
|
||||
"""
|
||||
def __init__(self, mesh, order, logSigma=True, normal='X'):
|
||||
"""
|
||||
def __init__(self, mesh, order, logSigma=True, normal='X', actInd = None):
|
||||
IdentityMap.__init__(self, mesh)
|
||||
self.logSigma = logSigma
|
||||
self.order = order
|
||||
self.normal = normal
|
||||
|
||||
slope = 1e4
|
||||
|
||||
self.normal = normal
|
||||
self.actInd = actInd
|
||||
|
||||
if getattr(self, 'actInd', None) is None:
|
||||
self.actInd = range(self.mesh.nC)
|
||||
self.nC = self.mesh.nC
|
||||
|
||||
else:
|
||||
self.nC = len(self.actInd)
|
||||
|
||||
slope = 1e4
|
||||
|
||||
@property
|
||||
def shape(self):
|
||||
return (self.nC, self.nP)
|
||||
|
||||
@property
|
||||
def nP(self):
|
||||
if np.isscalar(self.order):
|
||||
@@ -785,8 +797,8 @@ class PolyMap(IdentityMap):
|
||||
sig1, sig2 = np.exp(sig1), np.exp(sig2)
|
||||
#2D
|
||||
if self.mesh.dim == 2:
|
||||
X = self.mesh.gridCC[:,0]
|
||||
Y = self.mesh.gridCC[:,1]
|
||||
X = self.mesh.gridCC[self.actInd,0]
|
||||
Y = self.mesh.gridCC[self.actInd,1]
|
||||
if self.normal =='X':
|
||||
f = polynomial.polyval(Y, c) - X
|
||||
elif self.normal =='Y':
|
||||
@@ -795,9 +807,9 @@ class PolyMap(IdentityMap):
|
||||
raise(Exception("Input for normal = X or Y or Z"))
|
||||
#3D
|
||||
elif self.mesh.dim == 3:
|
||||
X = self.mesh.gridCC[:,0]
|
||||
Y = self.mesh.gridCC[:,1]
|
||||
Z = self.mesh.gridCC[:,2]
|
||||
X = self.mesh.gridCC[self.actInd,0]
|
||||
Y = self.mesh.gridCC[self.actInd,1]
|
||||
Z = self.mesh.gridCC[self.actInd,2]
|
||||
if self.normal =='X':
|
||||
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
|
||||
elif self.normal =='Y':
|
||||
@@ -806,6 +818,7 @@ class PolyMap(IdentityMap):
|
||||
f = polynomial.polyval2d(X, Y, c.reshape((self.order[0]+1,self.order[1]+1))) - Z
|
||||
else:
|
||||
raise(Exception("Input for normal = X or Y or Z"))
|
||||
|
||||
else:
|
||||
raise(Exception("Only supports 2D"))
|
||||
|
||||
@@ -819,8 +832,8 @@ class PolyMap(IdentityMap):
|
||||
sig1, sig2 = np.exp(sig1), np.exp(sig2)
|
||||
#2D
|
||||
if self.mesh.dim == 2:
|
||||
X = self.mesh.gridCC[:,0]
|
||||
Y = self.mesh.gridCC[:,1]
|
||||
X = self.mesh.gridCC[self.actInd,0]
|
||||
Y = self.mesh.gridCC[self.actInd,1]
|
||||
|
||||
if self.normal =='X':
|
||||
f = polynomial.polyval(Y, c) - X
|
||||
@@ -832,9 +845,9 @@ class PolyMap(IdentityMap):
|
||||
raise(Exception("Input for normal = X or Y or Z"))
|
||||
#3D
|
||||
elif self.mesh.dim == 3:
|
||||
X = self.mesh.gridCC[:,0]
|
||||
Y = self.mesh.gridCC[:,1]
|
||||
Z = self.mesh.gridCC[:,2]
|
||||
X = self.mesh.gridCC[self.actInd,0]
|
||||
Y = self.mesh.gridCC[self.actInd,1]
|
||||
Z = self.mesh.gridCC[self.actInd,2]
|
||||
|
||||
if self.normal =='X':
|
||||
f = polynomial.polyval2d(Y, Z, c.reshape((self.order[0]+1,self.order[1]+1))) - X
|
||||
|
||||
@@ -1003,6 +1003,7 @@ class ProjectedGNCG(BFGS, Minimize, Remember):
|
||||
# perturb inactive set off of bounds so that they are included in the step
|
||||
delx = delx + self.stepOffBoundsFact * (rhs_a * dm_i / dm_a)
|
||||
|
||||
|
||||
# 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.
|
||||
|
||||
+109
-113
@@ -344,7 +344,6 @@ class BaseRegularization(object):
|
||||
def W(self):
|
||||
"""Full regularization weighting matrix W."""
|
||||
return sp.identity(self.regmesh.nC)
|
||||
# self.regmesh._Pac.T * sp.identity(self.regmesh.nC) * self.regmesh._Pac # or do we want sp.identity(self.mesh.nC) or even just Utils.Identity() ?
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
@@ -375,11 +374,12 @@ class BaseRegularization(object):
|
||||
@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
|
||||
: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:
|
||||
|
||||
@@ -402,25 +402,48 @@ class BaseRegularization(object):
|
||||
|
||||
class Tikhonov(BaseRegularization):
|
||||
"""
|
||||
L2 Tikhonov regularization with both smallness and smoothness (first order
|
||||
derivative) contributions.
|
||||
|
||||
.. math::
|
||||
\phi_m(\mathbf{m}) = \\alpha_s \| W_s (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
|
||||
+ \\alpha_x \| W_x \\frac{\partial}{\partial x} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
|
||||
+ \\alpha_y \| W_y \\frac{\partial}{\partial y} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
|
||||
+ \\alpha_z \| W_z \\frac{\partial}{\partial z} (\mathbf{m} - \mathbf{m_{ref}} ) \|^2
|
||||
|
||||
Note if the key word argument `mrefInSmooth` is False, then mref is not
|
||||
included in the smoothness contribution.
|
||||
|
||||
:param Mesh mesh: SimPEG mesh
|
||||
:param Maps mapping: regularization mapping, takes the model from model space to the thing you want to regularize
|
||||
:param numpy.ndarray indActive: active cell indices for reducing the size of differential operators in the definition of a regularization mesh
|
||||
:param bool mrefInSmooth: (default = False) put mref in the smoothness component?
|
||||
:param float alpha_s: (default 1e-6) smallness weight
|
||||
:param float alpha_x: (default 1) smoothness weight for first derivative in the x-direction
|
||||
:param float alpha_y: (default 1) smoothness weight for first derivative in the y-direction
|
||||
:param float alpha_z: (default 1) smoothness weight for first derivative in the z-direction
|
||||
:param float alpha_xx: (default 1) smoothness weight for second derivative in the x-direction
|
||||
:param float alpha_yy: (default 1) smoothness weight for second derivative in the y-direction
|
||||
:param float alpha_zz: (default 1) smoothness weight for second derivative in the z-direction
|
||||
"""
|
||||
mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_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")
|
||||
mrefInSmooth = False # put mref in the smoothness contribution
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1e-6, ['_W', '_Wsmall'], "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, indActive = None, **kwargs):
|
||||
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
@property
|
||||
def Ws(self):
|
||||
"""Regularization matrix Ws"""
|
||||
if getattr(self,'_Ws', None) is None:
|
||||
self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
|
||||
return self._Ws
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self,'_Wsmall', None) is None:
|
||||
self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s)**0.5)
|
||||
return self._Wsmall
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
@@ -483,25 +506,44 @@ class Tikhonov(BaseRegularization):
|
||||
def W(self):
|
||||
"""Full regularization matrix W"""
|
||||
if getattr(self, '_W', None) is None:
|
||||
wlist = (self.Ws, self.Wsmooth)
|
||||
wlist = (self.Wsmall, self.Wsmooth)
|
||||
self._W = sp.vstack(wlist)
|
||||
return self._W
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
if self.mrefInSmooth == 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.mrefInSmooth == False:
|
||||
r = self.W * ( self.mapping * (m - self.mref) )
|
||||
return 0.5*r.dot(r)
|
||||
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)
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
return self._evalSmall(m) + self._evalSmooth(m)
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmallDeriv(self,m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return r.T * ( self.Wsmall * self.mapping.deriv(m - self.mref) )
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmoothDeriv(self,m):
|
||||
if self.mrefInSmooth == True:
|
||||
r = self.Wsmooth * ( self.mapping * ( m - self.mref ) )
|
||||
return r.T * ( self.Wsmooth * self.mapping.deriv(m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.Wsmooth * ( self.mapping * m )
|
||||
return r.T * ( self.Wsmooth * self.mapping.deriv(m) )
|
||||
|
||||
@Utils.timeIt
|
||||
def evalDeriv(self, m):
|
||||
"""
|
||||
|
||||
The regularization is:
|
||||
|
||||
.. math::
|
||||
@@ -515,45 +557,33 @@ class Tikhonov(BaseRegularization):
|
||||
R(m) = \mathbf{W^\\top W (m-m_\\text{ref})}
|
||||
|
||||
"""
|
||||
if self.mrefInSmooth == 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.mrefInSmooth == 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
|
||||
return self._evalSmallDeriv(m) + self._evalSmoothDeriv(m)
|
||||
|
||||
|
||||
class Simple(BaseRegularization):
|
||||
class Simple(Tikhonov):
|
||||
"""
|
||||
Only for tensor mesh
|
||||
Simple regularization that does not include length scales in the derivatives.
|
||||
"""
|
||||
|
||||
mrefInSmooth = True #: SMOOTH and SMOOTH_MOD_DIF options
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Ws'], "Smallness weight")
|
||||
mrefInSmooth = False #: SMOOTH and SMOOTH_MOD_DIF options
|
||||
alpha_s = Utils.dependentProperty('_alpha_s', 1.0, ['_W', '_Wsmall'], "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")
|
||||
wght = 1.
|
||||
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
BaseRegularization.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
|
||||
if isinstance(self.wght,float):
|
||||
self.wght = np.ones(self.regmesh.nC) * self.wght
|
||||
|
||||
@property
|
||||
def Ws(self):
|
||||
"""Regularization matrix Ws"""
|
||||
if getattr(self,'_Ws', None) is None:
|
||||
self._Ws = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
|
||||
return self._Ws
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self,'_Wsmall', None) is None:
|
||||
self._Wsmall = Utils.sdiag((self.regmesh.vol*self.alpha_s*self.wght)**0.5)
|
||||
return self._Wsmall
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
@@ -592,99 +622,65 @@ class Simple(BaseRegularization):
|
||||
def W(self):
|
||||
"""Full regularization matrix W"""
|
||||
if getattr(self, '_W', None) is None:
|
||||
wlist = (self.Ws, self.Wsmooth)
|
||||
wlist = (self.Wsmall, self.Wsmooth)
|
||||
self._W = sp.vstack(wlist)
|
||||
return self._W
|
||||
|
||||
@Utils.timeIt
|
||||
def _evalSmall(self, m):
|
||||
r = self.Wsmall * ( self.mapping * (m - self.mref) )
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
@Utils.timeIt
|
||||
def eval(self, m):
|
||||
def _evalSmooth(self, m):
|
||||
if self.mrefInSmooth == True:
|
||||
r1 = self.Wsmooth * ( self.mapping * (m) )
|
||||
r2 = self.Ws * ( self.mapping * (m - self.mref) )
|
||||
return 0.5*(r1.dot(r1)+r2.dot(r2))
|
||||
r = self.Wsmooth * ( self.mapping * (m - self.mref) )
|
||||
elif self.mrefInSmooth == False:
|
||||
r = self.W * ( self.mapping * (m - self.mref) )
|
||||
return 0.5*r.dot(r)
|
||||
return phim
|
||||
|
||||
|
||||
|
||||
@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.mrefInSmooth == 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.mrefInSmooth == 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
|
||||
r = self.Wsmooth * ( self.mapping * m)
|
||||
return 0.5 * r.dot(r)
|
||||
|
||||
|
||||
class Sparse(Simple):
|
||||
|
||||
# set default values
|
||||
eps = 1e-1
|
||||
eps_p = 1e-1
|
||||
eps_q = 1e-1
|
||||
curModel = None # use a model to compute the weights
|
||||
gamma = 1.
|
||||
p = 0.
|
||||
qx = 2.
|
||||
qy = 2.
|
||||
qz = 2.
|
||||
norms = [0., 2., 2., 2.]
|
||||
wght = 1.
|
||||
|
||||
def __init__(self, mesh, mapping=None, indActive=None, **kwargs):
|
||||
Simple.__init__(self, mesh, mapping=mapping, indActive=indActive, **kwargs)
|
||||
|
||||
|
||||
if isinstance(self.wght,float):
|
||||
self.wght = np.ones(self.regmesh.nC) * self.wght
|
||||
|
||||
@property
|
||||
def Ws(self):
|
||||
"""Regularization matrix Ws"""
|
||||
def Wsmall(self):
|
||||
"""Regularization matrix Wsmall"""
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rs = Utils.speye(self.regmesh.nC)
|
||||
|
||||
else:
|
||||
f_m = self.curModel - self.reg.mref
|
||||
self.rs = self.R(f_m , self.p)
|
||||
self.rs = self.R(f_m , self.eps_p, self.norms[0])
|
||||
#print "Min rs: " + str(np.max(self.rs)) + "Max rs: " + str(np.min(self.rs))
|
||||
self.Rs = Utils.sdiag( self.rs )
|
||||
|
||||
|
||||
return Utils.sdiag((self.regmesh.vol*self.alpha_s*self.gamma*self.wght)**0.5)*self.Rs
|
||||
|
||||
|
||||
@property
|
||||
def Wx(self):
|
||||
"""Regularization matrix Wx"""
|
||||
|
||||
|
||||
if getattr(self, 'curModel', None) is None:
|
||||
self.Rx = Utils.speye(self.regmesh.cellDiffxStencil.shape[0])
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffxStencil * self.curModel
|
||||
self.rx = self.R( f_m , self.qx)
|
||||
self.rx = self.R( f_m , self.eps_q, self.norms[1])
|
||||
self.Rx = Utils.sdiag( self.rx )
|
||||
|
||||
return Utils.sdiag(( (self.regmesh.aveCC2Fx * self.regmesh.vol) *self.alpha_x*self.gamma*(self.regmesh.aveCC2Fx*self.wght))**0.5)*self.Rx*self.regmesh.cellDiffxStencil
|
||||
@@ -698,9 +694,9 @@ class Sparse(Simple):
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffyStencil * self.curModel
|
||||
self.ry = self.R( f_m , self.qy)
|
||||
self.ry = self.R( f_m , self.eps_q, self.norms[2])
|
||||
self.Ry = Utils.sdiag( self.ry )
|
||||
|
||||
|
||||
return Utils.sdiag(((self.regmesh.aveCC2Fy * self.regmesh.vol)*self.alpha_y*self.gamma*(self.regmesh.aveCC2Fy*self.wght))**0.5)*self.Ry*self.regmesh.cellDiffyStencil
|
||||
|
||||
@property
|
||||
@@ -712,7 +708,7 @@ class Sparse(Simple):
|
||||
|
||||
else:
|
||||
f_m = self.regmesh.cellDiffzStencil * self.curModel
|
||||
self.rz = self.R( f_m , self.qz)
|
||||
self.rz = self.R( f_m , self.eps_q, self.norms[3])
|
||||
self.Rz = Utils.sdiag( self.rz )
|
||||
|
||||
return Utils.sdiag(((self.regmesh.aveCC2Fz * self.regmesh.vol)*self.alpha_z*self.gamma*(self.regmesh.aveCC2Fz*self.wght))**0.5)*self.Rz*self.regmesh.cellDiffzStencil
|
||||
@@ -733,13 +729,13 @@ class Sparse(Simple):
|
||||
def W(self):
|
||||
"""Full regularization matrix W"""
|
||||
#if getattr(self, '_W', None) is None:
|
||||
wlist = (self.Ws, self.Wsmooth)
|
||||
wlist = (self.Wsmall, self.Wsmooth)
|
||||
#self._W = sp.vstack(wlist)
|
||||
return sp.vstack(wlist)
|
||||
|
||||
def R(self, f_m , exponent):
|
||||
|
||||
eta = (self.eps**(1-exponent/2.))**0.5
|
||||
r = eta / (f_m**2.+self.eps**2.)**((1-exponent/2.)/2.)
|
||||
def R(self, f_m , eps, exponent):
|
||||
|
||||
eta = (eps**(1-exponent/2.))**0.5
|
||||
r = eta / (f_m**2.+ eps**2.)**((1-exponent/2.)/2.)
|
||||
|
||||
return r
|
||||
|
||||
+1
-1
@@ -375,7 +375,7 @@ class BaseSurvey(object):
|
||||
self.dtrue = self.dpred(m, f=f)
|
||||
noise = std*abs(self.dtrue)*np.random.randn(*self.dtrue.shape)
|
||||
self.dobs = self.dtrue+noise
|
||||
self.std = self.dobs*0 + std
|
||||
self.std = self.dobs*0. + std
|
||||
return self.dobs
|
||||
|
||||
class LinearSurvey(BaseSurvey):
|
||||
|
||||
@@ -0,0 +1,31 @@
|
||||
.. _examples_DC_PseudoSection_Simulation:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
|
||||
|
||||
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
|
||||
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.DC_PseudoSection_Simulation.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/DC_PseudoSection_Simulation.py
|
||||
:language: python
|
||||
:linenos:
|
||||
@@ -0,0 +1,41 @@
|
||||
.. _examples_EM_FDEM_SusEffects:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
|
||||
EM: FDEM: Effects of susceptibility
|
||||
===================================
|
||||
|
||||
When airborne freqeuncy domain EM (AFEM) survey is flown over
|
||||
the earth including significantly susceptible bodies (magnetite-rich rocks),
|
||||
negative data is often observed in the real part of the lowest frequency
|
||||
(e.g. Dighem system 900 Hz). This phenomenon mostly based upon magnetization
|
||||
occurs due to a susceptible body when the magnetic field is applied.
|
||||
|
||||
To clarify what is happening in the earth when we are exciting the earth with
|
||||
a loop source in the frequency domain we run three forward modelling:
|
||||
|
||||
- F[:math:`\sigma`, :math:`\mu`]: Anomalous conductivity and susceptibility
|
||||
- F[:math:`\sigma`, :math:`\mu_0`]: Anomalous conductivity
|
||||
- F[:math:`\sigma_{air}`, :math:`\mu_0`]: primary field
|
||||
|
||||
We plot vector magnetic fields in the earth. For secondary fields we provide
|
||||
F[:math:`\sigma`, :math:`\mu`]-F[:math:`\sigma`, :math:`\mu_0`]. Following
|
||||
figure show both real and parts.
|
||||
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.EM_FDEM_SusEffects.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/EM_FDEM_SusEffects.py
|
||||
:language: python
|
||||
:linenos:
|
||||
@@ -0,0 +1,26 @@
|
||||
.. _examples_Inversion_IRLS:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
|
||||
Inversion: Linear Problem
|
||||
=========================
|
||||
|
||||
Here we go over the basics of creating a linear problem and inversion.
|
||||
|
||||
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.Inversion_IRLS.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/Inversion_IRLS.py
|
||||
:language: python
|
||||
:linenos:
|
||||
@@ -0,0 +1,21 @@
|
||||
.. _examples_MT_1D_analytic_nlayer_Earth:
|
||||
|
||||
.. --------------------------------- ..
|
||||
.. ..
|
||||
.. THIS FILE IS AUTO GENEREATED ..
|
||||
.. ..
|
||||
.. SimPEG/Examples/__init__.py ..
|
||||
.. ..
|
||||
.. --------------------------------- ..
|
||||
|
||||
MT 1D analytic nlayer Earth
|
||||
===========================
|
||||
|
||||
.. plot::
|
||||
|
||||
from SimPEG import Examples
|
||||
Examples.MT_1D_analytic_nlayer_Earth.run()
|
||||
|
||||
.. literalinclude:: ../../SimPEG/Examples/MT_1D_analytic_nlayer_Earth.py
|
||||
:language: python
|
||||
:linenos:
|
||||
Reference in New Issue
Block a user