.. 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-33 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import time 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 34-39 Setup ----- Define the survey and model parameters .. GENERATED FROM PYTHON SOURCE LINES 39-50 .. 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 51-53 Define a dipping interface between the surface layer and the deeper layer .. GENERATED FROM PYTHON SOURCE LINES 54-76 .. 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 77-87 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 87-119 .. 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:100: 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:101: 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 120-125 Inversion Mesh -------------- Here, we set up a 2D tensor mesh which we will represent the inversion model on .. GENERATED FROM PYTHON SOURCE LINES 125-130 .. 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 131-137 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 137-153 .. 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 154-158 True Model ---------- Create our true model which we will use to generate synthetic data for .. GENERATED FROM PYTHON SOURCE LINES 158-170 .. 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 171-175 Survey ------ Create our true model which we will use to generate synthetic data for .. GENERATED FROM PYTHON SOURCE LINES 175-209 .. 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) .. GENERATED FROM PYTHON SOURCE LINES 210-215 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 215-263 .. 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 /home/vsts/work/1/s/simpeg/simulation.py:197: DefaultSolverWarning: Using the default solver: Pardiso. If you would like to suppress this notification, add warnings.filterwarnings('ignore', simpeg.utils.solver_utils.DefaultSolverWarning) to your script. Done forward simulation. Elapsed time = 15.82 s .. GENERATED FROM PYTHON SOURCE LINES 264-279 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 279-293 .. 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, random_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 294-298 Run the inversion ------------------ We start from a half-space equal to the deep conductivity. .. GENERATED FROM PYTHON SOURCE LINES 298-305 .. 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.23.0 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 = 200.52 s .. GENERATED FROM PYTHON SOURCE LINES 306-309 Plot the predicted and observed data ------------------------------------ .. GENERATED FROM PYTHON SOURCE LINES 309-314 .. 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 315-318 Plot the recovered model ------------------------ .. GENERATED FROM PYTHON SOURCE LINES 318-350 .. 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 351-354 Print the version of SimPEG and dependencies -------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 354-357 .. code-block:: Python Report() .. raw:: html
Sun Nov 03 19:22:56 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.15 | packaged by conda-forge | (main, Oct 16 2024, 01:24:24) [GCC 13.3.0]
simpeg 0.23.0 discretize 0.11.0 pymatsolver 0.3.1
numpy 2.0.2 scipy 1.14.1 matplotlib 3.9.2
geoana 0.7.2 libdlf 0.3.0 pydiso 0.1.2
numba 0.60.0 dask 2024.10.0 sklearn 1.5.2
pandas 2.2.3 sympy 1.13.3 plotly 5.24.1
memory_profiler 0.61.0 choclo 0.3.0


.. GENERATED FROM PYTHON SOURCE LINES 358-369 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 44.645 seconds) **Estimated memory usage:** 3135 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 `_