EM: TDEM: 1D: Inversion

Here we will create and run a TDEM 1D inversion.

SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.
SimPEG.InvProblem will set Regularization.reference_model to m0.

                        SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
                        ***Done using same Solver, and solver_opts as the Simulation3DElectricField problem***

/home/vsts/conda/envs/simpeg-test/lib/python3.8/site-packages/pymatsolver/direct.py:23: PardisoTypeConversionWarning:

Converting csc_matrix matrix to CSR format.

/home/vsts/conda/envs/simpeg-test/lib/python3.8/site-packages/pymatsolver/direct.py:73: PardisoTypeConversionWarning:

Converting csc_matrix matrix to CSR format.

model has any nan: 0
============================ Inexact Gauss Newton ============================
  #     beta     phi_d     phi_m       f      |proj(x-g)-x|  LS    Comment
-----------------------------------------------------------------------------
x0 has any nan: 0
   0  9.59e+02  5.26e+03  0.00e+00  5.26e+03    6.56e+03      0
   1  9.59e+02  3.94e+02  1.65e-01  5.52e+02    1.01e+03      0
   2  1.92e+02  4.51e+01  3.00e-01  1.03e+02    2.71e+02      0   Skip BFGS
   3  1.92e+02  3.92e+00  3.72e-01  7.52e+01    3.29e+01      0   Skip BFGS
   4  3.83e+01  3.65e+00  3.66e-01  1.77e+01    6.33e+01      0   Skip BFGS
   5  3.83e+01  4.31e-01  3.91e-01  1.54e+01    3.03e+00      0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 2.2481e+00 <= tolF*(1+|f0|) = 5.2595e+02
1 : |xc-x_last| = 8.9019e-02 <= tolX*(1+|x0|) = 3.0149e+00
0 : |proj(x-g)-x|    = 3.0311e+00 <= tolG          = 1.0000e-01
0 : |proj(x-g)-x|    = 3.0311e+00 <= 1e3*eps       = 1.0000e-02
1 : maxIter   =       5    <= iter          =      5
------------------------- DONE! -------------------------

import numpy as np
from SimPEG.electromagnetics import time_domain
from SimPEG import (
    optimization,
    discretize,
    maps,
    data_misfit,
    regularization,
    inverse_problem,
    inversion,
    directives,
    utils,
)
import matplotlib.pyplot as plt


def run(plotIt=True):
    cs, ncx, ncz, npad = 5.0, 25, 15, 15
    hx = [(cs, ncx), (cs, npad, 1.3)]
    hz = [(cs, npad, -1.3), (cs, ncz), (cs, npad, 1.3)]
    mesh = discretize.CylindricalMesh([hx, 1, hz], "00C")

    active = mesh.cell_centers_z < 0.0
    layer = (mesh.cell_centers_z < 0.0) & (mesh.cell_centers_z >= -100.0)
    actMap = maps.InjectActiveCells(mesh, active, np.log(1e-8), nC=mesh.shape_cells[2])
    mapping = maps.ExpMap(mesh) * maps.SurjectVertical1D(mesh) * actMap
    sig_half = 2e-3
    sig_air = 1e-8
    sig_layer = 1e-3
    sigma = np.ones(mesh.shape_cells[2]) * sig_air
    sigma[active] = sig_half
    sigma[layer] = sig_layer
    mtrue = np.log(sigma[active])

    rxOffset = 1e-3
    rx = time_domain.Rx.PointMagneticFluxTimeDerivative(
        np.array([[rxOffset, 0.0, 30]]), np.logspace(-5, -3, 31), "z"
    )
    src = time_domain.Src.MagDipole([rx], location=np.array([0.0, 0.0, 80]))
    survey = time_domain.Survey([src])
    time_steps = [(1e-06, 20), (1e-05, 20), (0.0001, 20)]
    simulation = time_domain.Simulation3DElectricField(
        mesh, sigmaMap=mapping, survey=survey, time_steps=time_steps
    )
    # d_true = simulation.dpred(mtrue)

    # create observed data
    rel_err = 0.05
    data = simulation.make_synthetic_data(mtrue, relative_error=rel_err)

    dmisfit = data_misfit.L2DataMisfit(simulation=simulation, data=data)
    regMesh = discretize.TensorMesh([mesh.h[2][mapping.maps[-1].indActive]])
    reg = regularization.WeightedLeastSquares(regMesh, alpha_s=1e-2, alpha_x=1.0)
    opt = optimization.InexactGaussNewton(maxIter=5, LSshorten=0.5)
    invProb = inverse_problem.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)
    simulation.counter = opt.counter = utils.Counter()
    opt.remember("xc")

    mopt = inv.run(m0)

    if plotIt:
        fig, ax = plt.subplots(1, 2, figsize=(10, 6))
        ax[0].loglog(rx.times, -invProb.dpred, "b.-")
        ax[0].loglog(rx.times, -data.dobs, "r.-")
        ax[0].legend(("Noisefree", "$d^{obs}$"), fontsize=16)
        ax[0].set_xlabel("Time (s)", fontsize=14)
        ax[0].set_ylabel("$B_z$ (T)", fontsize=16)
        ax[0].set_xlabel("Time (s)", fontsize=14)
        ax[0].grid(color="k", alpha=0.5, linestyle="dashed", linewidth=0.5)

        plt.semilogx(sigma[active], mesh.cell_centers_z[active])
        plt.semilogx(np.exp(mopt), mesh.cell_centers_z[active])
        ax[1].set_ylim(-600, 0)
        ax[1].set_xlim(1e-4, 1e-2)
        ax[1].set_xlabel("Conductivity (S/m)", fontsize=14)
        ax[1].set_ylabel("Depth (m)", fontsize=14)
        ax[1].grid(color="k", alpha=0.5, linestyle="dashed", linewidth=0.5)
        plt.legend([r"$\sigma_{true}$", r"$\sigma_{pred}$"])


if __name__ == "__main__":
    run()
    plt.show()

Estimated memory usage: 9 MB