.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/user-guide/examples/03-magnetics/plot_inv_mag_MVI_VectorAmplitude.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_user-guide_examples_03-magnetics_plot_inv_mag_MVI_VectorAmplitude.py: Magnetic inversion on a TreeMesh ================================ In this example, we demonstrate the use of a Magnetic Vector Inversion on 3D TreeMesh for the inversion of magnetic data. The inverse problem uses the :class:`simpeg.regularization.VectorAmplitude` regularization. .. GENERATED FROM PYTHON SOURCE LINES 12-35 .. code-block:: Python from simpeg import ( data, data_misfit, directives, maps, inverse_problem, optimization, inversion, regularization, ) from simpeg import utils from simpeg.utils import mkvc, sdiag from discretize.utils import mesh_builder_xyz, refine_tree_xyz, active_from_xyz from simpeg.potential_fields import magnetics import numpy as np import matplotlib.pyplot as plt # sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 36-45 Setup ----- Define the survey and model parameters First we need to define the direction of the inducing field As a simple case, we pick a vertical inducing field of magnitude 50,000 nT. .. GENERATED FROM PYTHON SOURCE LINES 45-74 .. code-block:: Python np.random.seed(1) # We will assume a vertical inducing field h0_amplitude, h0_inclination, h0_declination = (50000.0, 90.0, 0.0) # Create grid of points for topography # Lets create a simple Gaussian topo and set the active cells [xx, yy] = np.meshgrid(np.linspace(-200, 200, 50), np.linspace(-200, 200, 50)) b = 100 A = 50 zz = A * np.exp(-0.5 * ((xx / b) ** 2.0 + (yy / b) ** 2.0)) topo = np.c_[utils.mkvc(xx), utils.mkvc(yy), utils.mkvc(zz)] # Create an array of observation points xr = np.linspace(-100.0, 100.0, 20) yr = np.linspace(-100.0, 100.0, 20) X, Y = np.meshgrid(xr, yr) Z = A * np.exp(-0.5 * ((X / b) ** 2.0 + (Y / b) ** 2.0)) + 5 # Create a MAGsurvey xyzLoc = np.c_[mkvc(X.T), mkvc(Y.T), mkvc(Z.T)] rxLoc = magnetics.receivers.Point(xyzLoc) srcField = magnetics.sources.UniformBackgroundField( receiver_list=[rxLoc], amplitude=h0_amplitude, inclination=h0_inclination, declination=h0_declination, ) survey = magnetics.survey.Survey(srcField) .. GENERATED FROM PYTHON SOURCE LINES 75-80 Inversion Mesh -------------- Here, we create a TreeMesh with base cell size of 5 m. .. GENERATED FROM PYTHON SOURCE LINES 80-97 .. code-block:: Python # Create a mesh h = [5, 5, 5] padDist = np.ones((3, 2)) * 100 mesh = mesh_builder_xyz( xyzLoc, h, padding_distance=padDist, depth_core=100, mesh_type="tree" ) mesh = refine_tree_xyz( mesh, topo, method="surface", octree_levels=[2, 6], finalize=True ) # Define an active cells from topo actv = active_from_xyz(mesh, topo) nC = int(actv.sum()) .. rst-class:: sphx-glr-script-out .. code-block:: none /usr/share/miniconda/envs/simpeg-test/lib/python3.11/site-packages/discretize/utils/mesh_utils.py:528: FutureWarning: In discretize v1.0 the TreeMesh will change the default value of diagonal_balance to True, which will likely slightly change meshes you have previously created. If you need to keep the current behavior, explicitly set diagonal_balance=False. /home/vsts/work/1/s/examples/03-magnetics/plot_inv_mag_MVI_VectorAmplitude.py:88: DeprecationWarning: The surface option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_surface` functionality. It will be removed in a future version of discretize. .. GENERATED FROM PYTHON SOURCE LINES 98-103 Forward modeling data --------------------- We can now create a magnetization model and generate data. .. GENERATED FROM PYTHON SOURCE LINES 103-142 .. code-block:: Python model_azm_dip = np.zeros((mesh.nC, 2)) model_amp = np.ones(mesh.nC) * 1e-8 ind = utils.model_builder.get_indices_block( np.r_[-30, -20, -10], np.r_[30, 20, 25], mesh.gridCC, ) model_amp[ind] = 0.05 model_azm_dip[ind, 0] = 45.0 model_azm_dip[ind, 1] = 90.0 # Remove air cells model_azm_dip = model_azm_dip[actv, :] model_amp = model_amp[actv] model = sdiag(model_amp) * utils.mat_utils.dip_azimuth2cartesian( model_azm_dip[:, 0], model_azm_dip[:, 1] ) # Create reduced identity map idenMap = maps.IdentityMap(nP=nC * 3) # Create the simulation simulation = magnetics.simulation.Simulation3DIntegral( survey=survey, mesh=mesh, chiMap=idenMap, active_cells=actv, model_type="vector" ) # Compute some data and add some random noise d = simulation.dpred(mkvc(model)) std = 10 # nT synthetic_data = d + np.random.randn(len(d)) * std wd = np.ones(len(d)) * std # Assign data and uncertainties to the survey data_object = data.Data(survey, dobs=synthetic_data, standard_deviation=wd) # Create a projection matrix for plotting later actv_plot = maps.InjectActiveCells(mesh, actv, np.nan) .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/examples/03-magnetics/plot_inv_mag_MVI_VectorAmplitude.py:105: BreakingChangeWarning: Since SimPEG v0.25.0, the 'get_indices_block' function returns a single array with the cell indices, instead of a tuple with a single element. This means that we don't need to unpack the tuple anymore to access to the cell indices. If you were using this function as in: ind = get_indices_block(p0, p1, mesh.cell_centers)[0] Make sure you update it to: ind = get_indices_block(p0, p1, mesh.cell_centers) To hide this warning, add this to your script or notebook: import warnings from simpeg.utils import BreakingChangeWarning warnings.filterwarnings(action='ignore', category=BreakingChangeWarning) .. GENERATED FROM PYTHON SOURCE LINES 143-148 Inversion --------- We can now attempt the inverse calculations. .. GENERATED FROM PYTHON SOURCE LINES 148-202 .. code-block:: Python # Create sensitivity weights from our linear forward operator rxLoc = survey.source_field.receiver_list[0].locations # This Mapping connects the regularizations for the three-component # vector model wires = maps.Wires(("p", nC), ("s", nC), ("t", nC)) m0 = np.ones(3 * nC) * 1e-4 # Starting model # Create the regularization on the amplitude of magnetization reg = regularization.VectorAmplitude( mesh, mapping=idenMap, active_cells=actv, reference_model_in_smooth=True, norms=[0.0, 2.0, 2.0, 2.0], gradient_type="total", ) # Data misfit function dmis = data_misfit.L2DataMisfit(simulation=simulation, data=data_object) dmis.W = 1.0 / data_object.standard_deviation # The optimization scheme opt = optimization.ProjectedGNCG( maxIter=20, lower=-10, upper=10.0, maxIterLS=20, cg_maxiter=20, cg_rtol=1e-3 ) # The inverse problem invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) # Estimate the initial beta factor betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e1) # Add sensitivity weights sensitivity_weights = directives.UpdateSensitivityWeights() # Here is where the norms are applied IRLS = directives.UpdateIRLS( f_min_change=1e-3, max_irls_iterations=10, misfit_tolerance=5e-1 ) # Pre-conditioner update_Jacobi = directives.UpdatePreconditioner() inv = inversion.BaseInversion( invProb, directiveList=[sensitivity_weights, IRLS, update_Jacobi, betaest] ) # Run the inversion mrec = inv.run(m0) .. rst-class:: sphx-glr-script-out .. code-block:: none Running inversion with SimPEG v0.25.0 ================================================= Projected GNCG ================================================= # beta phi_d phi_m f |proj(x-g)-x| LS iter_CG CG |Ax-b|/|b| CG |Ax-b| Comment ----------------------------------------------------------------------------------------------------------------- 0 1.34e+05 1.19e+04 0.00e+00 1.19e+04 0 inf inf 1 1.34e+05 7.00e+03 1.50e-02 9.01e+03 2.43e+03 0 20 2.13e-01 4.89e+04 2 6.69e+04 4.77e+03 3.83e-02 7.33e+03 2.02e+03 0 13 5.87e-04 5.19e+01 3 3.35e+04 2.79e+03 8.04e-02 5.48e+03 1.86e+03 0 14 6.19e-04 3.79e+01 4 1.67e+04 1.39e+03 1.39e-01 3.72e+03 1.74e+03 0 15 9.53e-04 4.31e+01 5 8.36e+03 6.18e+02 2.03e-01 2.32e+03 1.62e+03 0 18 8.18e-04 2.49e+01 6 4.18e+03 2.67e+02 2.61e-01 1.36e+03 1.49e+03 0 20 1.58e-03 3.01e+01 Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 0.008897945940906758 7 4.18e+03 4.93e+02 3.71e-01 2.04e+03 1.56e+03 0 18 5.24e-04 1.02e+01 8 4.18e+03 7.49e+02 4.52e-01 2.64e+03 1.56e+03 0 13 9.41e-04 1.89e+01 9 2.69e+03 6.57e+02 6.57e-01 2.42e+03 1.20e+03 0 14 5.10e-04 3.57e+00 10 1.83e+03 6.17e+02 8.92e-01 2.24e+03 1.53e+03 0 12 9.01e-04 1.47e+01 11 1.28e+03 5.71e+02 1.08e+00 1.95e+03 1.57e+03 0 11 9.32e-04 1.90e+01 12 1.28e+03 6.49e+02 1.06e+00 2.00e+03 1.70e+03 0 9 9.21e-04 2.79e+01 13 8.72e+02 5.29e+02 1.15e+00 1.53e+03 9.50e+02 0 11 9.02e-04 5.47e+00 14 8.72e+02 5.37e+02 1.07e+00 1.47e+03 1.47e+03 0 10 4.60e-04 5.17e+00 15 8.72e+02 5.35e+02 1.00e+00 1.41e+03 1.37e+03 0 10 3.35e-04 2.90e+00 16 8.72e+02 5.32e+02 9.51e-01 1.36e+03 1.33e+03 0 9 9.80e-04 7.17e+00 Reach maximum number of IRLS cycles: 10 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 1.6611e+01 <= tolF*(1+|f0|) = 1.1929e+03 1 : |xc-x_last| = 2.9300e-02 <= tolX*(1+|x0|) = 1.0290e-01 0 : |proj(x-g)-x| = 1.3253e+03 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 1.3253e+03 <= 1e3*eps = 1.0000e-02 0 : maxIter = 20 <= iter = 16 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 203-210 Final Plot ---------- Let's compare the smooth and compact model .. GENERATED FROM PYTHON SOURCE LINES 210-256 .. code-block:: Python plt.figure(figsize=(12, 6)) ax = plt.subplot(2, 2, 1) im = utils.plot_utils.plot2Ddata(xyzLoc, synthetic_data, ax=ax) plt.colorbar(im[0]) ax.set_title("Predicted data.") plt.gca().set_aspect("equal", adjustable="box") for ii, (title, mvec) in enumerate( [("True model", model), ("Smooth model", invProb.l2model), ("Sparse model", mrec)] ): ax = plt.subplot(2, 2, ii + 2) mesh.plot_slice( actv_plot * mvec.reshape((-1, 3), order="F"), v_type="CCv", view="vec", ax=ax, normal="Y", grid=True, quiver_opts={ "pivot": "mid", "scale": 8 * np.abs(mvec).max(), "scale_units": "inches", }, ) ax.set_xlim([-200, 200]) ax.set_ylim([-100, 75]) ax.set_title(title) ax.set_xlabel("x") ax.set_ylabel("z") plt.gca().set_aspect("equal", adjustable="box") plt.show() print("END") # Plot the final predicted data and the residual # plt.figure() # ax = plt.subplot(1, 2, 1) # utils.plot_utils.plot2Ddata(xyzLoc, invProb.dpred, ax=ax) # ax.set_title("Predicted data.") # plt.gca().set_aspect("equal", adjustable="box") # # ax = plt.subplot(1, 2, 2) # utils.plot_utils.plot2Ddata(xyzLoc, synthetic_data - invProb.dpred, ax=ax) # ax.set_title("Data residual.") # plt.gca().set_aspect("equal", adjustable="box") .. image-sg:: /content/user-guide/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_VectorAmplitude_001.png :alt: Predicted data., True model, Smooth model, Sparse model :srcset: /content/user-guide/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_VectorAmplitude_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none END .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 24.643 seconds) **Estimated memory usage:** 448 MB .. _sphx_glr_download_content_user-guide_examples_03-magnetics_plot_inv_mag_MVI_VectorAmplitude.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_mag_MVI_VectorAmplitude.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_mag_MVI_VectorAmplitude.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_inv_mag_MVI_VectorAmplitude.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_