.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.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_03-magnetics_plot_inv_mag_MVI_Sparse_TreeMesh.py: Magnetic inversion on a TreeMesh ================================ In this example, we demonstrate the use of a Magnetic Vector Inverison on 3D TreeMesh for the inversion of magnetics affected by remanence. The mesh is auto-generated based on the position of the observation locations and topography. We invert the data twice, first for a smooth starting model using the Cartesian coordinate system, and second for a compact model using the Spherical formulation. The inverse problem uses the :class:'SimPEG.regularization.Sparse' that .. GENERATED FROM PYTHON SOURCE LINES 18-42 .. 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 from discretize.utils import active_from_xyz, mesh_builder_xyz, refine_tree_xyz from SimPEG.potential_fields import magnetics import scipy as sp import numpy as np import matplotlib.pyplot as plt # sphinx_gallery_thumbnail_number = 3 .. GENERATED FROM PYTHON SOURCE LINES 43-52 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 52-95 .. 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) # The magnetization is set along a different direction (induced + remanence) M = np.array([45.0, 90.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) # Here how the topography looks with a quick interpolation, just a Gaussian... tri = sp.spatial.Delaunay(topo) fig = plt.figure() ax = fig.add_subplot(1, 1, 1, projection="3d") ax.plot_trisurf( topo[:, 0], topo[:, 1], topo[:, 2], triangles=tri.simplices, cmap=plt.cm.Spectral ) ax.scatter3D(xyzLoc[:, 0], xyzLoc[:, 1], xyzLoc[:, 2], c="k") plt.show() .. image-sg:: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_001.png :alt: plot inv mag MVI Sparse TreeMesh :srcset: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 96-107 Inversion Mesh -------------- Here, we create a TreeMesh with base cell size of 5 m. We created a small utility function to center the mesh around points and to figure out the outermost dimension for adequate padding distance. The second stage allows us to refine the mesh around points or surfaces (point assumed to follow some horizontal trend) The refinement process is repeated twice to allow for a finer level around the survey locations. .. GENERATED FROM PYTHON SOURCE LINES 107-124 .. 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=[4, 4], 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 /home/vsts/work/1/s/examples/03-magnetics/plot_inv_mag_MVI_Sparse_TreeMesh.py:115: 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 125-131 Forward modeling data --------------------- We can now create a magnetization model and generate data Lets start with a block below topography .. GENERATED FROM PYTHON SOURCE LINES 131-207 .. code-block:: Python model = np.zeros((mesh.nC, 3)) # Convert the inclination declination to vector in Cartesian M_xyz = utils.mat_utils.dip_azimuth2cartesian(M[0], M[1]) # Get the indicies of the magnetized block ind = utils.model_builder.get_indices_block( np.r_[-20, -20, -10], np.r_[20, 20, 25], mesh.gridCC, )[0] # Assign magnetization values model[ind, :] = np.kron(np.ones((ind.shape[0], 1)), M_xyz * 0.05) # Remove air cells model = model[actv, :] # Create active map to go from reduce set to full actvMap = maps.InjectActiveCells(mesh, actv, np.nan) # Creat reduced identity map idenMap = maps.IdentityMap(nP=nC * 3) # Create the simulation simulation = magnetics.simulation.Simulation3DIntegral( survey=survey, mesh=mesh, chiMap=idenMap, ind_active=actv, model_type="vector" ) # Compute some data and add some random noise d = simulation.dpred(mkvc(model)) std = 5 # 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 an projection matrix for plotting later actv_plot = maps.InjectActiveCells(mesh, actv, np.nan) # Plot the model and data plt.figure() ax = plt.subplot(2, 1, 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") # Plot the vector model ax = plt.subplot(2, 1, 2) mesh.plot_slice( actv_plot * model.reshape((-1, 3), order="F"), v_type="CCv", view="vec", ax=ax, normal="Y", ind=66, grid=True, quiver_opts={ "pivot": "mid", "scale": 5 * np.abs(model).max(), "scale_units": "inches", }, ) ax.set_xlim([-200, 200]) ax.set_ylim([-100, 75]) ax.set_xlabel("x") ax.set_ylabel("y") plt.gca().set_aspect("equal", adjustable="box") plt.show() .. image-sg:: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_002.png :alt: Predicted data., Slice 66, Y = 12.50 :srcset: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 208-218 Inversion --------- We can now attempt the inverse calculations. We put great care into designing an inversion methology that would yield a geologically reasonable solution for the non-induced problem. The inversion is done in two stages. First we compute a smooth solution using a Cartesian coordinate system, then a sparse inversion in the Spherical domain. .. GENERATED FROM PYTHON SOURCE LINES 218-275 .. 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 three regularizations for the different components # of magnetization reg_p = regularization.Sparse(mesh, active_cells=actv, mapping=wires.p) reg_p.reference_model = np.zeros(3 * nC) reg_s = regularization.Sparse(mesh, active_cells=actv, mapping=wires.s) reg_s.reference_model = np.zeros(3 * nC) reg_t = regularization.Sparse(mesh, active_cells=actv, mapping=wires.t) reg_t.reference_model = np.zeros(3 * nC) reg = reg_p + reg_s + reg_t reg.reference_model = np.zeros(3 * nC) # Data misfit function dmis = data_misfit.L2DataMisfit(simulation=simulation, data=data_object) dmis.W = 1.0 / data_object.standard_deviation # Add directives to the inversion opt = optimization.ProjectedGNCG( maxIter=10, lower=-10, upper=10.0, maxIterLS=20, maxIterCG=20, tolCG=1e-4 ) invProb = inverse_problem.BaseInvProblem(dmis, reg, opt) # A list of directive to control the inverson betaest = directives.BetaEstimate_ByEig(beta0_ratio=1e1) # Add sensitivity weights sensitivity_weights = directives.UpdateSensitivityWeights() # Here is where the norms are applied # Use a threshold parameter empirically based on the distribution of # model parameters IRLS = directives.Update_IRLS(f_min_change=1e-3, max_irls_iterations=2, beta_tol=5e-1) # Pre-conditioner update_Jacobi = directives.UpdatePreconditioner() inv = inversion.BaseInversion( invProb, directiveList=[sensitivity_weights, IRLS, update_Jacobi, betaest] ) # Run the inversion mrec_MVIC = inv.run(m0) .. rst-class:: sphx-glr-script-out .. code-block:: none SimPEG.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using the default solver Pardiso and no solver_opts.*** model has any nan: 0 =============================== Projected GNCG =============================== # beta phi_d phi_m f |proj(x-g)-x| LS Comment ----------------------------------------------------------------------------- x0 has any nan: 0 0 1.14e+06 2.48e+04 2.40e-03 2.76e+04 3.02e+03 0 1 5.72e+05 1.66e+04 3.99e-03 1.88e+04 2.82e+03 0 2 2.86e+05 1.16e+04 8.61e-03 1.41e+04 2.37e+03 0 Skip BFGS 3 1.43e+05 7.22e+03 1.97e-02 1.00e+04 2.24e+03 0 Skip BFGS 4 7.15e+04 3.72e+03 3.70e-02 6.36e+03 2.10e+03 0 Skip BFGS 5 3.57e+04 1.62e+03 5.74e-02 3.67e+03 1.94e+03 0 Skip BFGS 6 1.79e+04 6.27e+02 7.64e-02 1.99e+03 1.77e+03 0 Skip BFGS Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 0.005087485697837089 irls_threshold 0.003431354664705422 irls_threshold 0.006119051004170579 7 8.94e+03 2.33e+02 1.21e-01 1.31e+03 1.40e+03 0 Skip BFGS 8 2.35e+04 1.23e+02 1.57e-01 3.81e+03 1.95e+03 0 Skip BFGS Reach maximum number of IRLS cycles: 2 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 2.7591e+03 1 : |xc-x_last| = 1.6948e-02 <= tolX*(1+|x0|) = 1.0336e-01 0 : |proj(x-g)-x| = 1.9458e+03 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 1.9458e+03 <= 1e3*eps = 1.0000e-02 0 : maxIter = 10 <= iter = 9 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 276-283 Sparse Vector Inversion ----------------------- Re-run the MVI in the spherical domain so we can impose sparsity in the vectors. .. GENERATED FROM PYTHON SOURCE LINES 283-370 .. code-block:: Python spherical_map = maps.SphericalSystem() m_start = utils.mat_utils.cartesian2spherical(mrec_MVIC.reshape((nC, 3), order="F")) beta = invProb.beta dmis.simulation.chiMap = spherical_map dmis.simulation.model = m_start # Create a block diagonal regularization wires = maps.Wires(("amp", nC), ("theta", nC), ("phi", nC)) # Create a Combo Regularization # Regularize the amplitude of the vectors reg_a = regularization.Sparse( mesh, gradient_type="total", active_cells=actv, mapping=wires.amp, norms=[0.0, 1.0, 1.0, 1.0], # Only norm on gradients used, reference_model=np.zeros(3 * nC), ) # Regularize the vertical angle of the vectors reg_t = regularization.Sparse( mesh, gradient_type="total", active_cells=actv, mapping=wires.theta, alpha_s=0.0, # No reference angle, norms=[0.0, 1.0, 1.0, 1.0], # Only norm on gradients used, ) reg_t.units = "radian" # Regularize the horizontal angle of the vectors reg_p = regularization.Sparse( mesh, gradient_type="total", active_cells=actv, mapping=wires.phi, alpha_s=0.0, # No reference angle, norms=[0.0, 1.0, 1.0, 1.0], # Only norm on gradients used, ) reg_p.units = "radian" reg = reg_a + reg_t + reg_p reg.reference_model = np.zeros(3 * nC) lower_bound = np.kron(np.asarray([0, -np.inf, -np.inf]), np.ones(nC)) upper_bound = np.kron(np.asarray([10, np.inf, np.inf]), np.ones(nC)) # Add directives to the inversion opt = optimization.ProjectedGNCG( maxIter=20, lower=lower_bound, upper=upper_bound, maxIterLS=20, maxIterCG=30, tolCG=1e-3, stepOffBoundsFact=1e-3, ) opt.approxHinv = None invProb = inverse_problem.BaseInvProblem(dmis, reg, opt, beta=beta) # Here is where the norms are applied irls = directives.Update_IRLS( f_min_change=1e-4, max_irls_iterations=20, minGNiter=1, beta_tol=0.5, coolingRate=1, coolEps_q=True, sphericalDomain=True, ) # Special directive specific to the mag amplitude problem. The sensitivity # weights are updated between each iteration. spherical_projection = directives.ProjectSphericalBounds() sensitivity_weights = directives.UpdateSensitivityWeights() update_Jacobi = directives.UpdatePreconditioner() inv = inversion.BaseInversion( invProb, directiveList=[spherical_projection, irls, sensitivity_weights, update_Jacobi], ) mrec_MVI_S = inv.run(m_start) .. 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 is setting bfgsH0 to the inverse of the eval2Deriv. ***Done using the default solver Pardiso and no solver_opts.*** model has any nan: 0 =============================== Projected GNCG =============================== # beta phi_d phi_m f |proj(x-g)-x| LS Comment ----------------------------------------------------------------------------- x0 has any nan: 0 0 1.61e+04 6.48e+02 1.29e-01 2.72e+03 1.22e+03 0 1 8.04e+03 7.52e+02 1.15e-01 1.68e+03 1.11e+03 1 2 4.02e+03 7.40e+02 1.14e-01 1.20e+03 1.13e+03 5 3 2.01e+03 6.90e+02 1.14e-01 9.19e+02 1.09e+03 3 Skip BFGS 4 1.01e+03 4.90e+02 1.62e-01 6.53e+02 8.83e+02 0 5 5.03e+02 4.56e+02 1.51e-01 5.32e+02 8.92e+02 3 Reached starting chifact with l2-norm regularization: Start IRLS steps... irls_threshold 0.006389732504246946 irls_threshold 1.8403750808867727 irls_threshold 6.276081775408219 6 2.51e+02 3.10e+02 2.67e-01 3.77e+02 8.79e+02 1 Skip BFGS 7 5.26e+02 1.83e+02 3.62e-01 3.73e+02 8.70e+02 1 8 1.27e+03 1.41e+02 4.72e-01 7.40e+02 7.91e+02 1 9 2.75e+03 1.72e+02 5.47e-01 1.68e+03 6.95e+02 2 10 2.75e+03 4.06e+02 6.23e-01 2.12e+03 7.69e+02 1 11 2.75e+03 4.31e+02 6.27e-01 2.16e+03 7.72e+02 3 12 2.75e+03 4.36e+02 7.44e-01 2.48e+03 8.15e+02 1 Skip BFGS 13 1.83e+03 6.87e+02 7.23e-01 2.01e+03 1.00e+03 1 14 1.24e+03 6.66e+02 7.24e-01 1.56e+03 1.12e+03 2 15 1.24e+03 5.21e+02 8.46e-01 1.57e+03 9.02e+02 1 16 1.24e+03 5.67e+02 8.41e-01 1.61e+03 1.05e+03 1 17 1.24e+03 5.06e+02 8.18e-01 1.52e+03 9.05e+02 0 18 1.24e+03 4.99e+02 7.48e-01 1.42e+03 8.99e+02 2 19 1.24e+03 3.82e+02 7.28e-01 1.28e+03 7.78e+02 0 Skip BFGS 20 1.24e+03 4.14e+02 6.82e-01 1.26e+03 7.84e+02 1 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 2.4544e+01 <= tolF*(1+|f0|) = 2.7183e+02 0 : |xc-x_last| = 5.4398e+01 <= tolX*(1+|x0|) = 3.7014e+01 0 : |proj(x-g)-x| = 7.8438e+02 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 7.8438e+02 <= 1e3*eps = 1.0000e-02 1 : maxIter = 20 <= iter = 20 ------------------------- DONE! ------------------------- .. GENERATED FROM PYTHON SOURCE LINES 371-378 Final Plot ---------- Let's compare the smooth and compact model .. GENERATED FROM PYTHON SOURCE LINES 378-441 .. code-block:: Python plt.figure(figsize=(8, 8)) ax = plt.subplot(2, 1, 1) mesh.plot_slice( actv_plot * mrec_MVIC.reshape((nC, 3), order="F"), v_type="CCv", view="vec", ax=ax, normal="Y", ind=66, grid=True, quiver_opts={ "pivot": "mid", "scale": 5 * np.abs(mrec_MVIC).max(), "scale_units": "inches", }, ) ax.set_xlim([-200, 200]) ax.set_ylim([-100, 75]) ax.set_title("Smooth model (Cartesian)") ax.set_xlabel("x") ax.set_ylabel("y") plt.gca().set_aspect("equal", adjustable="box") ax = plt.subplot(2, 1, 2) vec_xyz = utils.mat_utils.spherical2cartesian( mrec_MVI_S.reshape((nC, 3), order="F") ).reshape((nC, 3), order="F") mesh.plot_slice( actv_plot * vec_xyz, v_type="CCv", view="vec", ax=ax, normal="Y", ind=66, grid=True, quiver_opts={ "pivot": "mid", "scale": 5 * np.abs(vec_xyz).max(), "scale_units": "inches", }, ) ax.set_xlim([-200, 200]) ax.set_ylim([-100, 75]) ax.set_title("Sparse model (L0L2)") ax.set_xlabel("x") ax.set_ylabel("y") plt.gca().set_aspect("equal", adjustable="box") plt.show() # 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") .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_003.png :alt: Smooth model (Cartesian), Sparse model (L0L2) :srcset: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_003.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_004.png :alt: Predicted data., Data residual. :srcset: /content/examples/03-magnetics/images/sphx_glr_plot_inv_mag_MVI_Sparse_TreeMesh_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (1 minutes 47.537 seconds) **Estimated memory usage:** 654 MB .. _sphx_glr_download_content_examples_03-magnetics_plot_inv_mag_MVI_Sparse_TreeMesh.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_Sparse_TreeMesh.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_inv_mag_MVI_Sparse_TreeMesh.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_