mirror of
https://github.com/wassname/simpeg.git
synced 2026-08-04 13:13:40 +08:00
Merge branch 'master' of https://github.com/simpeg/simpeg into Examples
Conflicts: SimPEG/Examples/__init__.py
This commit is contained in:
seogi_macbook
@@ -0,0 +1,116 @@
|
||||
from SimPEG import *
|
||||
import SimPEG.EM as EM
|
||||
from SimPEG.EM import mu_0
|
||||
|
||||
|
||||
def run(plotIt=True):
|
||||
"""
|
||||
EM: FDEM: 1D: Inversion
|
||||
=======================
|
||||
|
||||
Here we will create and run a FDEM 1D inversion.
|
||||
|
||||
"""
|
||||
|
||||
cs, ncx, ncz, npad = 5., 25, 15, 15
|
||||
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')
|
||||
|
||||
layerz = -100.
|
||||
|
||||
active = mesh.vectorCCz<0.
|
||||
layer = (mesh.vectorCCz<0.) & (mesh.vectorCCz>=layerz)
|
||||
actMap = Maps.ActiveCells(mesh, active, np.log(1e-8), nC=mesh.nCz)
|
||||
mapping = Maps.ExpMap(mesh) * Maps.Vertical1DMap(mesh) * actMap
|
||||
sig_half = 2e-2
|
||||
sig_air = 1e-8
|
||||
sig_layer = 1e-2
|
||||
sigma = np.ones(mesh.nCz)*sig_air
|
||||
sigma[active] = sig_half
|
||||
sigma[layer] = sig_layer
|
||||
mtrue = np.log(sigma[active])
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
fig, ax = plt.subplots(1,1, figsize = (3, 6))
|
||||
plt.semilogx(sigma[active], mesh.vectorCCz[active])
|
||||
ax.set_ylim(-500, 0)
|
||||
ax.set_xlim(1e-3, 1e-1)
|
||||
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
|
||||
ax.set_ylabel('Depth (m)', fontsize = 14)
|
||||
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
|
||||
|
||||
|
||||
rxOffset=10.
|
||||
bzi = EM.FDEM.Rx(np.array([[rxOffset, 0., 1e-3]]), 'bzi')
|
||||
|
||||
freqs = np.logspace(1,3,10)
|
||||
srcLoc = np.array([0., 0., 10.])
|
||||
|
||||
srcList = []
|
||||
[srcList.append(EM.FDEM.Src.MagDipole([bzi],freq, srcLoc,orientation='Z')) for freq in freqs]
|
||||
|
||||
survey = EM.FDEM.Survey(srcList)
|
||||
prb = EM.FDEM.Problem_b(mesh, mapping=mapping)
|
||||
|
||||
try:
|
||||
from pymatsolver import MumpsSolver
|
||||
prb.Solver = MumpsSolver
|
||||
except ImportError, e:
|
||||
prb.Solver = SolverLU
|
||||
|
||||
prb.pair(survey)
|
||||
|
||||
std = 0.05
|
||||
survey.makeSyntheticData(mtrue, std)
|
||||
|
||||
survey.std = std
|
||||
survey.eps = np.linalg.norm(survey.dtrue)*1e-5
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
fig, ax = plt.subplots(1,1, figsize = (6, 6))
|
||||
ax.semilogx(freqs,survey.dtrue[:freqs.size], 'b.-')
|
||||
ax.semilogx(freqs,survey.dobs[:freqs.size], 'r.-')
|
||||
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
|
||||
ax.set_xlabel('Time (s)', fontsize = 14)
|
||||
ax.set_ylabel('$B_z$ (T)', fontsize = 16)
|
||||
ax.set_xlabel('Time (s)', fontsize = 14)
|
||||
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
|
||||
|
||||
dmisfit = DataMisfit.l2_DataMisfit(survey)
|
||||
regMesh = Mesh.TensorMesh([mesh.hz[mapping.maps[-1].indActive]])
|
||||
reg = Regularization.Tikhonov(regMesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIter = 6)
|
||||
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
|
||||
|
||||
# Create an inversion object
|
||||
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
|
||||
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
|
||||
inv = Inversion.BaseInversion(invProb, directiveList=[beta,betaest])
|
||||
m0 = np.log(np.ones(mtrue.size)*sig_half)
|
||||
reg.alpha_s = 1e-3
|
||||
reg.alpha_x = 1.
|
||||
prb.counter = opt.counter = Utils.Counter()
|
||||
opt.LSshorten = 0.5
|
||||
opt.remember('xc')
|
||||
|
||||
mopt = inv.run(m0)
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
fig, ax = plt.subplots(1,1, figsize = (3, 6))
|
||||
plt.semilogx(sigma[active], mesh.vectorCCz[active])
|
||||
plt.semilogx(np.exp(mopt), mesh.vectorCCz[active])
|
||||
ax.set_ylim(-500, 0)
|
||||
ax.set_xlim(1e-3, 1e-1)
|
||||
ax.set_xlabel('Conductivity (S/m)', fontsize = 14)
|
||||
ax.set_ylabel('Depth (m)', fontsize = 14)
|
||||
ax.grid(color='k', alpha=0.5, linestyle='dashed', linewidth=0.5)
|
||||
plt.legend(['$\sigma_{true}$', '$\sigma_{pred}$'],loc='best')
|
||||
plt.show()
|
||||
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
@@ -0,0 +1,145 @@
|
||||
from SimPEG import *
|
||||
from SimPEG import EM
|
||||
from pymatsolver import MumpsSolver
|
||||
from scipy.constants import mu_0
|
||||
|
||||
def run(plotIt=True):
|
||||
"""
|
||||
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 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 only real part, since that is our interest.
|
||||
|
||||
"""
|
||||
# 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)
|
||||
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()
|
||||
vizfields(1, primsec="primary", realimag="real")
|
||||
vizfields(1, primsec="secondary", realimag="real")
|
||||
|
||||
if __name__ == '__main__':
|
||||
run()
|
||||
@@ -1,6 +1,6 @@
|
||||
from SimPEG import *
|
||||
import SimPEG.EM as EM
|
||||
from scipy.constants import mu_0
|
||||
from SimPEG.EM import mu_0
|
||||
|
||||
|
||||
def run(plotIt=True):
|
||||
@@ -50,20 +50,18 @@ def run(plotIt=True):
|
||||
prb.Solver = SolverLU
|
||||
prb.timeSteps = [(1e-06, 20),(1e-05, 20), (0.0001, 20)]
|
||||
prb.pair(survey)
|
||||
dtrue = survey.dpred(mtrue)
|
||||
|
||||
|
||||
survey.dtrue = dtrue
|
||||
# create observed data
|
||||
std = 0.05
|
||||
noise = std*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
|
||||
survey.dobs = survey.dtrue+noise
|
||||
survey.std = survey.dobs*0 + std
|
||||
survey.Wd = 1/(abs(survey.dobs)*std)
|
||||
|
||||
survey.dobs = survey.makeSyntheticData(mtrue,std)
|
||||
survey.std = std
|
||||
survey.eps = 1e-5*np.linalg.norm(survey.dobs)
|
||||
|
||||
if plotIt:
|
||||
import matplotlib.pyplot as plt
|
||||
fig, ax = plt.subplots(1,1, figsize = (10, 6))
|
||||
ax.loglog(rx.times, dtrue, 'b.-')
|
||||
ax.loglog(rx.times, survey.dtrue, 'b.-')
|
||||
ax.loglog(rx.times, survey.dobs, 'r.-')
|
||||
ax.legend(('Noisefree', '$d^{obs}$'), fontsize = 16)
|
||||
ax.set_xlabel('Time (s)', fontsize = 14)
|
||||
@@ -76,6 +74,7 @@ def run(plotIt=True):
|
||||
reg = Regularization.Tikhonov(regMesh)
|
||||
opt = Optimization.InexactGaussNewton(maxIter = 5)
|
||||
invProb = InvProblem.BaseInvProblem(dmisfit, reg, opt)
|
||||
|
||||
# Create an inversion object
|
||||
beta = Directives.BetaSchedule(coolingFactor=5, coolingRate=2)
|
||||
betaest = Directives.BetaEstimate_ByEig(beta0_ratio=1e0)
|
||||
|
||||
Reference in New Issue
Block a user