Working on examples

This commit is contained in:
GudniRos
2016-06-01 01:05:21 -07:00
parent f34314bba2
commit fb7f5a53d4
2 changed files with 24 additions and 21 deletions
+18 -15
View File
@@ -1,9 +1,11 @@
import SimPEG as simpeg
import numpy as np
import SimPEG.MT as MT
from SimPEG import NSEM
from scipy.constants import mu_0
import matplotlib.pyplot as plt
np.random.seed(1983)
def run(plotIt=True):
"""
MT: 1D: Inversion
@@ -23,7 +25,7 @@ def run(plotIt=True):
ct = 20
air = simpeg.Utils.meshTensor([(ct,16,1.4)])
core = np.concatenate( ( np.kron(simpeg.Utils.meshTensor([(ct,10,-1.3)]),np.ones((5,))) , simpeg.Utils.meshTensor([(ct,5)]) ) )
bot = simpeg.Utils.meshTensor([(core[0],10,-1.4)])
bot = simpeg.Utils.meshTensor([(core[0],12,-1.4)])
x0 = -np.array([np.sum(np.concatenate((core,bot)))])
# Make the model
m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
@@ -47,41 +49,43 @@ def run(plotIt=True):
# Make the background model
sigma_0 = np.ones(m1d.nCx)*sig_air
sigma_0[active] = sig_half
sigma_0[layer1] = sig_layer1
sigma_0[layer2] = .002
m_0 = np.log(sigma_0[active])
# Set the mapping
actMap = simpeg.Maps.ActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
actMap = simpeg.Maps.InjectActiveCells(m1d, active, np.log(1e-8), nC=m1d.nCx)
mappingExpAct = simpeg.Maps.ExpMap(m1d) * actMap
## Setup the layout of the survey, set the sources and the connected receivers
# Receivers
rxList = []
for rxType in ['z1dr','z1di']:
rxList.append(MT.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
rxList.append(NSEM.Rx(simpeg.mkvc(np.array([0.0]),2).T,rxType))
# Source list
srcList =[]
for freq in freqs:
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Make the survey
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
survey.mtrue = m_true
## Set the problem
problem = MT.Problem1D.eForm_psField(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
problem = NSEM.Problem1D_ePrimSec(m1d,sigmaPrimary=sigma_0,mapping=mappingExpAct)
problem.pair(survey)
## Forward model data
# Project the data
survey.dtrue = survey.dpred(m_true)
survey.dobs = survey.dtrue + 0.025*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
survey.dobs = survey.dtrue + 0.01*abs(survey.dtrue)*np.random.randn(*survey.dtrue.shape)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[])
fig.suptitle('Target - smooth true')
# Assign uncertainties
std = 0.05 # 5% std
std = 0.025 # 5% std
survey.std = np.abs(survey.dobs*std)
# Assign the data weight
Wd = 1./survey.std
@@ -92,6 +96,7 @@ def run(plotIt=True):
# Set the optimization
opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
opt.counter = C
opt.maxStep = m1d.nC * survey.nD
opt.LSshorten = 0.5
opt.remember('xc')
# Data misfit
@@ -99,7 +104,7 @@ def run(plotIt=True):
dmis.Wd = Wd
# Regularization - with a regularization mesh
regMesh = simpeg.Mesh.TensorMesh([m1d.hx[problem.mapping.sigmaMap.maps[-1].indActive]],m1d.x0)
reg = simpeg.Regularization.Tikhonov(regMesh)
reg = simpeg.Regularization.Tikhonov(m1d,indActive=active)
reg.mrefInSmooth = True
reg.alpha_s = 1e-7
reg.alpha_x = 1.
@@ -109,11 +114,9 @@ def run(plotIt=True):
# Beta cooling
beta = simpeg.Directives.BetaSchedule()
beta.coolingRate = 4
betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=1.)
targmis = simpeg.Directives.TargetMisfit()
targmis.target = survey.nD
saveModel = simpeg.Directives.SaveModelEveryIteration()
saveModel.fileName = 'Inversion_TargMisEqnD_smoothTrue'
# Create an inversion object
inv = simpeg.Inversion.BaseInversion(invProb, directiveList=[beta,betaest,targmis])
@@ -121,7 +124,7 @@ def run(plotIt=True):
mopt = inv.run(m_0)
if plotIt:
fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
fig = NSEM.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
fig.suptitle('Target - smooth true')
plt.show()
+6 -6
View File
@@ -2,7 +2,7 @@
# Import
import SimPEG as simpeg
from SimPEG import MT
from SimPEG import NSEM
import numpy as np
try:
from pymatsolver import MumpsSolver as Solver
@@ -37,16 +37,16 @@ def run(plotIt=True, nFreq=1):
for loc in rx_loc:
# NOTE: loc has to be a (1,3) np.ndarray otherwise errors accure
for rxType in ['zxxr','zxxi','zxyr','zxyi','zyxr','zyxi','zyyr','zyyi','tzxr','tzxi','tzyr','tzyi']:
rxList.append(MT.Rx(simpeg.mkvc(loc,2).T,rxType))
rxList.append(NSEM.Rx(simpeg.mkvc(loc,2).T,rxType))
# Source list
srcList =[]
for freq in np.logspace(3,-3,nFreq):
srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
srcList.append(NSEM.SrcNSEM.polxy_1Dprimary(rxList,freq))
# Survey MT
survey = MT.Survey(srcList)
survey = NSEM.Survey(srcList)
## Setup the problem object
problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
problem = NSEM.Problem3D_ePrimSec(M, sigmaPrimary=sigBG)
problem.pair(survey)
problem.Solver = Solver
@@ -55,7 +55,7 @@ def run(plotIt=True, nFreq=1):
dataVec = survey.eval(fields)
# Make the data
mtData = MT.Data(survey,dataVec)
mtData = NSEM.Data(survey,dataVec)
# Add plots
if plotIt:
pass