.. 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 /home/vsts/work/1/s/examples/05-fdem/plot_inv_fdem_loop_loop_2Dinversion.py:104: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) /home/vsts/work/1/s/examples/05-fdem/plot_inv_fdem_loop_loop_2Dinversion.py:105: DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.) 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 = 15.94 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 Running inversion with SimPEG v0.22.2.dev13+g048ef809f 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 3.20e+00 1.42e+04 0.00e+00 1.42e+04 3.21e+03 0 1 3.20e+00 1.48e+03 5.42e+00 1.50e+03 4.61e+02 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 1.4191e+03 1 : |xc-x_last| = 6.2345e+00 <= tolX*(1+|x0|) = 1.3056e+01 0 : |proj(x-g)-x| = 4.6099e+02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 4.6099e+02 <= 1e3*eps = 1.0000e-02 0 : maxIter = 20 <= iter = 2 ------------------------- DONE! ------------------------- Inversion Complete. Elapsed Time = 203.76 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
Thu Sep 19 08:50:52 2024 UTC
OS Linux (Ubuntu 22.04) CPU(s) 2 Machine x86_64
Architecture 64bit RAM 6.8 GiB Environment Python
File system ext4
Python 3.10.14 | packaged by conda-forge | (main, Mar 20 2024, 12:45:18) [GCC 12.3.0]
simpeg 0.22.2.dev13+g048ef809f discretize 0.10.0 pymatsolver 0.2.0
numpy 1.26.4 scipy 1.14.1 matplotlib 3.9.2
empymod 2.3.1 geoana 0.6.0 pydiso 0.0.5
numba 0.60.0 dask 2024.9.0 sklearn 1.5.2
pandas 2.2.2 sympy 1.13.2 plotly 5.24.1
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 47.134 seconds) **Estimated memory usage:** 3156 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 ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_fdem_loop_loop_2Dinversion.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_