Merge branch 'dev' of https://github.com/simpeg/simpeg into dcip/dev

# Conflicts: # SimPEG/Examples/__init__.py
2026-07-16 11:21:38 +08:00 · 2016-02-16 22:14:22 -08:00
parent 9bf4a6228d fb1973cfe4
commit 16239aa414
48 changed files with 5327 additions and 565 deletions
@@ -0,0 +1,129 @@
+import SimPEG as simpeg
+import numpy as np
+import SimPEG.MT as MT
+from scipy.constants import mu_0
+import matplotlib.pyplot as plt
+
+def run(plotIt=True):
+    """
+        MT: 1D: Inversion
+        =======================
+
+        Forward model 1D MT data.
+        Setup and run a MT 1D inversion.
+
+    """
+
+    ## Setup the forward modeling
+    # Setting up 1D mesh and conductivity models to forward model data.
+    # Frequency
+    nFreq = 31
+    freqs = np.logspace(3,-3,nFreq)
+    # Set mesh parameters
+    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)])
+    x0 = -np.array([np.sum(np.concatenate((core,bot)))])
+    # Make the model
+    m1d = simpeg.Mesh.TensorMesh([np.concatenate((bot,core,air))], x0=x0)
+
+    # Setup model varibles
+    active = m1d.vectorCCx<0.
+    layer1 = (m1d.vectorCCx<-500.) & (m1d.vectorCCx>=-800.)
+    layer2 = (m1d.vectorCCx<-3500.) & (m1d.vectorCCx>=-5000.)
+    # Set the conductivity values
+    sig_half = 2e-3
+    sig_air = 1e-8
+    sig_layer1 = .2
+    sig_layer2 = .2
+    # Make the true model
+    sigma_true = np.ones(m1d.nCx)*sig_air
+    sigma_true[active] = sig_half
+    sigma_true[layer1] = sig_layer1
+    sigma_true[layer2] = sig_layer2
+    # Extract the model
+    m_true = np.log(sigma_true[active])
+    # Make the background model
+    sigma_0 = np.ones(m1d.nCx)*sig_air
+    sigma_0[active] = sig_half
+    m_0 = np.log(sigma_0[active])
+
+    # Set the mapping
+    actMap = simpeg.Maps.ActiveCells(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))
+    # Source list
+    srcList =[]
+    for freq in freqs:
+            srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+    # Make the survey
+    survey = MT.Survey(srcList)
+    survey.mtrue = m_true
+
+    ## Set the problem
+    problem = MT.Problem1D.eForm_psField(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)
+
+    if plotIt:
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem)
+        fig.suptitle('Target - smooth true')
+
+
+    # Assign uncertainties
+    std = 0.05 # 5% std
+    survey.std = np.abs(survey.dobs*std)
+    # Assign the data weight
+    Wd = 1./survey.std
+
+    ## Setup the inversion proceedure
+    # Define a counter
+    C =  simpeg.Utils.Counter()
+    # Set the optimization
+    opt = simpeg.Optimization.InexactGaussNewton(maxIter = 30)
+    opt.counter = C
+    opt.LSshorten = 0.5
+    opt.remember('xc')
+    # Data misfit
+    dmis = simpeg.DataMisfit.l2_DataMisfit(survey)
+    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.smoothModel = True
+    reg.alpha_s = 1e-7
+    reg.alpha_x = 1.
+    # Inversion problem
+    invProb = simpeg.InvProblem.BaseInvProblem(dmis, reg, opt)
+    invProb.counter = C
+    # Beta cooling
+    beta = simpeg.Directives.BetaSchedule()
+    beta.coolingRate = 4
+    betaest = simpeg.Directives.BetaEstimate_ByEig(beta0_ratio=0.75)
+    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])
+
+    ## Run the inversion
+    mopt = inv.run(m_0)
+
+    if plotIt:
+        fig = MT.Utils.dataUtils.plotMT1DModelData(problem,[mopt])
+        fig.suptitle('Target - smooth true')
+        plt.show()
+
+if __name__ == '__main__':
+    run()
@@ -0,0 +1,64 @@
+# Test script to use SimPEG.MT platform to forward model synthetic data.
+
+# Import
+import SimPEG as simpeg
+from SimPEG import MT
+import numpy as np
+try:
+    from pymatsolver import MumpsSolver as Solver
+except:
+    from SimPEG import Solver
+
+def run(plotIt=True, nFreq=1):
+    """
+        MT: 3D: Forward
+        =======================
+
+        Forward model 3D MT data.
+
+    """
+
+    # Make a mesh
+    M = simpeg.Mesh.TensorMesh([[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,-1.5),(100.,10),(100,5,1.5)],[(100,5,1.6),(100.,10),(100,3,2)]], x0=['C','C',-3529.5360])
+    # Setup the model
+    conds = [1e-2,1]
+    sig = simpeg.Utils.ModelBuilder.defineBlock(M.gridCC,[-1000,-1000,-400],[1000,1000,-200],conds)
+    sig[M.gridCC[:,2]>0] = 1e-8
+    sig[M.gridCC[:,2]<-600] = 1e-1
+    sigBG = np.zeros(M.nC) + conds[0]
+    sigBG[M.gridCC[:,2]>0] = 1e-8
+
+    ## Setup the the survey object
+    # Receiver locations
+    rx_x, rx_y = np.meshgrid(np.arange(-500,501,50),np.arange(-500,501,50))
+    rx_loc = np.hstack((simpeg.Utils.mkvc(rx_x,2),simpeg.Utils.mkvc(rx_y,2),np.zeros((np.prod(rx_x.shape),1))))
+    # Make a receiver list
+    rxList = []
+    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))
+    # Source list
+    srcList =[]
+    for freq in np.logspace(3,-3,nFreq):
+        srcList.append(MT.SrcMT.polxy_1Dprimary(rxList,freq))
+    # Survey MT
+    survey = MT.Survey(srcList)
+
+    ## Setup the problem object
+    problem = MT.Problem3D.eForm_ps(M, sigmaPrimary=sigBG)
+    problem.pair(survey)
+    problem.Solver = Solver
+
+    # Calculate the data
+    fields = problem.fields(sig)
+    dataVec = survey.projectFields(fields)
+
+    # Make the data
+    mtData = MT.Data(survey,dataVec)
+    # Add plots
+    if plotIt:
+        pass
+
+if __name__ == '__main__':
+    run()
@@ -16,8 +16,10 @@ import Mesh_QuadTree_Creation
 import Mesh_QuadTree_FaceDiv
 import Mesh_QuadTree_HangingNodes
 import Mesh_Tensor_Creation
+import MT_1D_ForwardAndInversion
+import MT_3D_Foward

-__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "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"]
+__examples__ = ["DC_Analytic_Dipole", "DC_Forward_PseudoSection", "EM_FDEM_1D_Inversion", "EM_FDEM_Analytic_MagDipoleWholespace", "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"]

 ##### AUTOIMPORTS #####