.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/05-fdem/plot_inv_fdem_loop_loop_2Dinversion.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_examples_05-fdem_plot_inv_fdem_loop_loop_2Dinversion.py: 2D inversion of Loop-Loop EM Data ================================= In this example, we consider a single line of loop-loop EM data at 30kHz with 3 different coil separations [0.32m, 0.71m, 1.18m]. We will use only Horizontal co-planar orientations (vertical magnetic dipole), and look at the real and imaginary parts of the secondary magnetic field. We use the :class:`SimPEG.maps.Surject2Dto3D` mapping to invert for a 2D model and perform the forward modelling in 3D. .. GENERATED FROM PYTHON SOURCE LINES 14-37 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import time try: from pymatsolver import Pardiso as Solver except ImportError: from SimPEG import SolverLU as Solver import discretize from SimPEG import ( maps, optimization, data_misfit, regularization, inverse_problem, inversion, directives, Report, ) from SimPEG.electromagnetics import frequency_domain as FDEM .. GENERATED FROM PYTHON SOURCE LINES 38-43 Setup ----- Define the survey and model parameters .. GENERATED FROM PYTHON SOURCE LINES 43-54 .. code-block:: Python sigma_surface = 10e-3 sigma_deep = 40e-3 sigma_air = 1e-8 coil_separations = [0.32, 0.71, 1.18] freq = 30e3 print("skin_depth: {:1.2f}m".format(500 / np.sqrt(sigma_deep * freq))) .. rst-class:: sphx-glr-script-out .. code-block:: none skin_depth: 14.43m .. GENERATED FROM PYTHON SOURCE LINES 55-57 Define a dipping interface between the surface layer and the deeper layer .. GENERATED FROM PYTHON SOURCE LINES 58-80 .. code-block:: Python z_interface_shallow = -0.25 z_interface_deep = -1.5 x_dip = np.r_[0.0, 8.0] def interface(x): interface = np.zeros_like(x) interface[x < x_dip[0]] = z_interface_shallow dipping_unit = (x >= x_dip[0]) & (x <= x_dip[1]) x_dipping = (-(z_interface_shallow - z_interface_deep) / x_dip[1]) * ( x[dipping_unit] ) + z_interface_shallow interface[dipping_unit] = x_dipping interface[x > x_dip[1]] = z_interface_deep return interface .. GENERATED FROM PYTHON SOURCE LINES 81-91 Forward Modelling Mesh ---------------------- Here, we set up a 3D tensor mesh which we will perform the forward simulations on. .. note:: In practice, a smaller horizontal discretization should be used to improve accuracy, particularly for the shortest offset (eg. you can try 0.25m). .. GENERATED FROM PYTHON SOURCE LINES 91-123 .. code-block:: Python csx = 0.5 # cell size for the horizontal direction csz = 0.125 # cell size for the vertical direction pf = 1.3 # expansion factor for the padding cells npadx = 7 # number of padding cells in the x-direction npady = 7 # number of padding cells in the y-direction npadz = 11 # number of padding cells in the z-direction core_domain_x = np.r_[-11.5, 11.5] # extent of uniform cells in the x-direction core_domain_z = np.r_[-2.0, 0.0] # extent of uniform cells in the z-direction # number of cells in the core region ncx = int(np.diff(core_domain_x) / csx) ncz = int(np.diff(core_domain_z) / csz) # create a 3D tensor mesh mesh = discretize.TensorMesh( [ [(csx, npadx, -pf), (csx, ncx), (csx, npadx, pf)], [(csx, npady, -pf), (csx, 1), (csx, npady, pf)], [(csz, npadz, -pf), (csz, ncz), (csz, npadz, pf)], ] ) # set the origin mesh.x0 = np.r_[ -mesh.h[0].sum() / 2.0, -mesh.h[1].sum() / 2.0, -mesh.h[2][: npadz + ncz].sum() ] print("the mesh has {} cells".format(mesh.nC)) mesh.plot_grid() .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_001.png :alt: plot inv fdem loop loop 2Dinversion :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none the mesh has 34200 cells .. GENERATED FROM PYTHON SOURCE LINES 124-129 Inversion Mesh -------------- Here, we set up a 2D tensor mesh which we will represent the inversion model on .. GENERATED FROM PYTHON SOURCE LINES 129-134 .. code-block:: Python inversion_mesh = discretize.TensorMesh([mesh.h[0], mesh.h[2][mesh.cell_centers_z <= 0]]) inversion_mesh.x0 = [-inversion_mesh.h[0].sum() / 2.0, -inversion_mesh.h[1].sum()] inversion_mesh.plot_grid() .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_002.png :alt: plot inv fdem loop loop 2Dinversion :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 135-141 Mappings --------- Mappings are used to take the inversion model and represent it as electrical conductivity on the inversion mesh. We will invert for log-conductivity below the surface, fixing the conductivity of the air cells to 1e-8 S/m .. GENERATED FROM PYTHON SOURCE LINES 141-157 .. code-block:: Python # create a 2D mesh that includes air cells mesh2D = discretize.TensorMesh([mesh.h[0], mesh.h[2]], x0=mesh.x0[[0, 2]]) active_inds = mesh2D.gridCC[:, 1] < 0 # active indices are below the surface mapping = ( maps.Surject2Dto3D(mesh) * maps.InjectActiveCells( # populates 3D space from a 2D model mesh2D, active_inds, sigma_air ) * maps.ExpMap( # adds air cells nP=inversion_mesh.nC ) # takes the exponential (log(sigma) --> sigma) ) .. GENERATED FROM PYTHON SOURCE LINES 158-162 True Model ---------- Create our true model which we will use to generate synthetic data for .. GENERATED FROM PYTHON SOURCE LINES 162-174 .. code-block:: Python m_true = np.log(sigma_deep) * np.ones(inversion_mesh.nC) interface_depth = interface(inversion_mesh.gridCC[:, 0]) m_true[inversion_mesh.gridCC[:, 1] > interface_depth] = np.log(sigma_surface) fig, ax = plt.subplots(1, 1) cb = plt.colorbar(inversion_mesh.plot_image(m_true, ax=ax, grid=True)[0], ax=ax) cb.set_label(r"$\log(\sigma)$") ax.set_title("true model") ax.set_xlim([-10, 10]) ax.set_ylim([-2, 0]) .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_003.png :alt: true model :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none (-2.0, 0.0) .. GENERATED FROM PYTHON SOURCE LINES 175-179 Survey ------ Create our true model which we will use to generate synthetic data for .. GENERATED FROM PYTHON SOURCE LINES 179-215 .. code-block:: Python src_locations = np.arange(-11, 11, 0.5) src_z = 0.25 # src is 0.25m above the surface orientation = "z" # z-oriented dipole for horizontal co-planar loops # reciever offset in 3D space rx_offsets = np.vstack([np.r_[sep, 0.0, 0.0] for sep in coil_separations]) # create our source list - one source per location source_list = [] for x in src_locations: src_loc = np.r_[x, 0.0, src_z] rx_locs = src_loc - rx_offsets rx_real = FDEM.Rx.PointMagneticFluxDensitySecondary( locations=rx_locs, orientation=orientation, component="real" ) rx_imag = FDEM.Rx.PointMagneticFluxDensitySecondary( locations=rx_locs, orientation=orientation, component="imag" ) src = FDEM.Src.MagDipole( receiver_list=[rx_real, rx_imag], location=src_loc, orientation=orientation, frequency=freq, ) source_list.append(src) # create the survey and problem objects for running the forward simulation survey = FDEM.Survey(source_list) prob = FDEM.Simulation3DMagneticFluxDensity( mesh, survey=survey, sigmaMap=mapping, solver=Solver ) .. GENERATED FROM PYTHON SOURCE LINES 216-221 Set up data for inversion ------------------------- Generate clean, synthetic data. Later we will invert the clean data, and assign a standard deviation of 0.05, and a floor of 1e-11. .. GENERATED FROM PYTHON SOURCE LINES 221-269 .. code-block:: Python t = time.time() data = prob.make_synthetic_data( m_true, relative_error=0.05, noise_floor=1e-11, add_noise=False ) dclean = data.dclean print("Done forward simulation. Elapsed time = {:1.2f} s".format(time.time() - t)) def plot_data(data, ax=None, color="C0", label=""): if ax is None: fig, ax = plt.subplots(1, 3, figsize=(15, 5)) # data is [re, im, re, im, ...] data_real = data[0::2] data_imag = data[1::2] for i, offset in enumerate(coil_separations): ax[i].plot( src_locations, data_real[i :: len(coil_separations)], color=color, label="{} real".format(label), ) ax[i].plot( src_locations, data_imag[i :: len(coil_separations)], "--", color=color, label="{} imag".format(label), ) ax[i].set_title("offset = {:1.2f}m".format(offset)) ax[i].legend() ax[i].grid(which="both") ax[i].set_ylim(np.r_[data.min(), data.max()] + 1e-11 * np.r_[-1, 1]) ax[i].set_xlabel("source location x (m)") ax[i].set_ylabel("Secondary B-Field (T)") plt.tight_layout() return ax ax = plot_data(dclean) .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_004.png :alt: offset = 0.32m, offset = 0.71m, offset = 1.18m :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_004.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Done forward simulation. Elapsed time = 8.33 s .. GENERATED FROM PYTHON SOURCE LINES 270-285 Set up the inversion -------------------- We create the data misfit, simple regularization (a least-squares-style regularization, :class:`SimPEG.regularization.LeastSquareRegularization`) The smoothness and smallness contributions can be set by including `alpha_s, alpha_x, alpha_y` as input arguments when the regularization is created. The default reference model in the regularization is the starting model. To set something different, you can input an `mref` into the regularization. We estimate the trade-off parameter, beta, between the data misfit and regularization by the largest eigenvalue of the data misfit and the regularization. Here, we use a fixed beta, but could alternatively employ a beta-cooling schedule using :class:`SimPEG.directives.BetaSchedule` .. GENERATED FROM PYTHON SOURCE LINES 285-299 .. code-block:: Python dmisfit = data_misfit.L2DataMisfit(simulation=prob, data=data) reg = regularization.WeightedLeastSquares(inversion_mesh) opt = optimization.InexactGaussNewton(maxIterCG=10, remember="xc") invProb = inverse_problem.BaseInvProblem(dmisfit, reg, opt) betaest = directives.BetaEstimate_ByEig(beta0_ratio=0.05, n_pw_iter=1, seed=1) target = directives.TargetMisfit() directiveList = [betaest, target] inv = inversion.BaseInversion(invProb, directiveList=directiveList) print("The target misfit is {:1.2f}".format(target.target)) .. rst-class:: sphx-glr-script-out .. code-block:: none The target misfit is 264.00 .. GENERATED FROM PYTHON SOURCE LINES 300-304 Run the inversion ------------------ We start from a half-space equal to the deep conductivity. .. GENERATED FROM PYTHON SOURCE LINES 304-311 .. code-block:: Python m0 = np.log(sigma_deep) * np.ones(inversion_mesh.nC) t = time.time() mrec = inv.run(m0) print("\n Inversion Complete. Elapsed Time = {:1.2f} s".format(time.time() - t)) .. rst-class:: sphx-glr-script-out .. code-block:: none 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 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 Simulation3DMagneticFluxDensity problem*** 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 2.53e+00 1.42e+04 0.00e+00 1.42e+04 3.21e+03 0 1 2.53e+00 1.48e+03 5.42e+00 1.49e+03 4.61e+02 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.4191e+03 1 : |xc-x_last| = 6.2505e+00 <= tolX*(1+|x0|) = 1.3056e+01 0 : |proj(x-g)-x| = 4.6148e+02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 4.6148e+02 <= 1e3*eps = 1.0000e-02 0 : maxIter = 20 <= iter = 2 ------------------------- DONE! ------------------------- Inversion Complete. Elapsed Time = 183.86 s .. GENERATED FROM PYTHON SOURCE LINES 312-315 Plot the predicted and observed data ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 315-320 .. code-block:: Python fig, ax = plt.subplots(1, 3, figsize=(15, 5)) plot_data(dclean, ax=ax, color="C0", label="true") plot_data(invProb.dpred, ax=ax, color="C1", label="predicted") .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_005.png :alt: offset = 0.32m, offset = 0.71m, offset = 1.18m :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_005.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none array([, , ], dtype=object) .. GENERATED FROM PYTHON SOURCE LINES 321-324 Plot the recovered model ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 324-356 .. code-block:: Python fig, ax = plt.subplots(1, 2, figsize=(12, 5)) # put both plots on the same colorbar clim = np.r_[np.log(sigma_surface), np.log(sigma_deep)] # recovered model cb = plt.colorbar( inversion_mesh.plot_image(mrec, ax=ax[0], clim=clim)[0], ax=ax[0], ) ax[0].set_title("recovered model") cb.set_label(r"$\log(\sigma)$") # true model cb = plt.colorbar( inversion_mesh.plot_image(m_true, ax=ax[1], clim=clim)[0], ax=ax[1], ) ax[1].set_title("true model") cb.set_label(r"$\log(\sigma)$") # # uncomment to plot the true interface # x = np.linspace(-10, 10, 50) # [a.plot(x, interface(x), 'k') for a in ax] [a.set_xlim([-10, 10]) for a in ax] [a.set_ylim([-2, 0]) for a in ax] plt.tight_layout() plt.show() .. image-sg:: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_006.png :alt: recovered model, true model :srcset: /content/examples/05-fdem/images/sphx_glr_plot_inv_fdem_loop_loop_2Dinversion_006.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 357-360 Print the version of SimPEG and dependencies -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 360-363 .. code-block:: Python Report() .. raw:: html
Wed Jun 19 15:02:22 2024 PDT
OS Linux CPU(s) 48 Machine x86_64
Architecture 64bit RAM 125.8 GiB Environment Python
File system ext4
Python 3.8.19 | packaged by conda-forge | (default, Mar 20 2024, 12:47:35) [GCC 12.3.0]
SimPEG 0.21.1.dev3+ge1c28e94c discretize 0.10.0 pymatsolver 0.2.0
numpy 1.24.4 scipy 1.10.1 sklearn 1.3.2
matplotlib 3.7.3 empymod 2.2.2 geoana 0.6.0
pandas 2.0.3 pydiso 0.0.5 numba 0.58.1
dask 2023.5.0 sympy 1.12 IPython 8.12.2
ipywidgets 8.1.2 plotly 5.22.0 vtk 9.2.6
memory_profiler 0.61.0 choclo 0.2.0


.. GENERATED FROM PYTHON SOURCE LINES 364-375 Moving Forward -------------- If you have suggestions for improving this example, please create a `pull request on the example in SimPEG `_ You might try: - improving the discretization - changing beta - changing the noise model - playing with the regulariztion parameters - ... .. rst-class:: sphx-glr-timing **Total running time of the script:** (3 minutes 22.270 seconds) **Estimated memory usage:** 3227 MB .. _sphx_glr_download_content_examples_05-fdem_plot_inv_fdem_loop_loop_2Dinversion.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_inv_fdem_loop_loop_2Dinversion.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_fdem_loop_loop_2Dinversion.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_