.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "content/user-guide/tutorials/06-ip/plot_inv_3_dcip3d.py"
.. LINE NUMBERS ARE GIVEN BELOW.
.. only:: html
.. note::
:class: sphx-glr-download-link-note
:ref:`Go to the end <sphx_glr_download_content_user-guide_tutorials_06-ip_plot_inv_3_dcip3d.py>`
to download the full example code.
.. rst-class:: sphx-glr-example-title
.. _sphx_glr_content_user-guide_tutorials_06-ip_plot_inv_3_dcip3d.py:
3D Least-Squares Inversion of DC and IP Data
============================================
Here we invert 5 lines of DC and IP data to recover both an electrical
conductivity and a chargeability model. We formulate the corresponding
inverse problems as least-squares optimization problems.
For this tutorial, we focus on the following:
- Generating a mesh based on survey geometry
- Including surface topography
- Defining the inverse problem (data misfit, regularization, directives)
- Applying sensitivity weighting
- Plotting the recovered model and data misfit
The DC data are measured voltages normalized by the source current in V/A
and the IP data are defined as apparent chargeabilities and V/V.
.. GENERATED FROM PYTHON SOURCE LINES 24-27
Import Modules
--------------
.. GENERATED FROM PYTHON SOURCE LINES 27-70
.. code-block:: Python
import os
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import tarfile
from discretize import TreeMesh
from discretize.utils import refine_tree_xyz, active_from_xyz
from simpeg.utils import model_builder
from simpeg.utils.io_utils.io_utils_electromagnetics import read_dcip_xyz
from simpeg import (
maps,
data_misfit,
regularization,
optimization,
inverse_problem,
inversion,
directives,
utils,
)
from simpeg.electromagnetics.static import resistivity as dc
from simpeg.electromagnetics.static import induced_polarization as ip
from simpeg.electromagnetics.static.utils.static_utils import (
apparent_resistivity_from_voltage,
)
# To plot DC/IP data in 3D, the user must have the plotly package
try:
import plotly
from simpeg.electromagnetics.static.utils.static_utils import plot_3d_pseudosection
has_plotly = True
except ImportError:
has_plotly = False
pass
mpl.rcParams.update({"font.size": 16})
# sphinx_gallery_thumbnail_number = 7
.. GENERATED FROM PYTHON SOURCE LINES 71-82
Define File Names
-----------------
Here we provide the file paths to assets we need to run the inversion. The
path to the true model conductivity and chargeability models are also
provided for comparison with the inversion results. These files are stored as a
tar-file on our google cloud bucket:
"https://storage.googleapis.com/simpeg/doc-assets/dcip3d.tar.gz"
.. GENERATED FROM PYTHON SOURCE LINES 82-102
.. code-block:: Python
# storage bucket where we have the data
data_source = "https://storage.googleapis.com/simpeg/doc-assets/dcip3d.tar.gz"
# download the data
downloaded_data = utils.download(data_source, overwrite=True)
# unzip the tarfile
tar = tarfile.open(downloaded_data, "r")
tar.extractall()
tar.close()
# path to the directory containing our data
dir_path = downloaded_data.split(".")[0] + os.path.sep
# files to work with
topo_filename = dir_path + "topo_xyz.txt"
dc_data_filename = dir_path + "dc_data.xyz"
ip_data_filename = dir_path + "ip_data.xyz"
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Downloading https://storage.googleapis.com/simpeg/doc-assets/dcip3d.tar.gz
saved to: /home/vsts/work/1/s/tutorials/06-ip/dcip3d.tar.gz
Download completed!
.. GENERATED FROM PYTHON SOURCE LINES 103-109
Load Data and Topography
------------------------
Here we load the observed data and topography.
.. GENERATED FROM PYTHON SOURCE LINES 109-128
.. code-block:: Python
topo_xyz = np.loadtxt(str(topo_filename))
dc_data = read_dcip_xyz(
dc_data_filename,
"volt",
data_header="V/A",
uncertainties_header="UNCERT",
is_surface_data=False,
)
ip_data = read_dcip_xyz(
ip_data_filename,
"apparent_chargeability",
data_header="APP_CHG",
uncertainties_header="UNCERT",
is_surface_data=False,
)
.. GENERATED FROM PYTHON SOURCE LINES 129-136
Plot Observed DC Data in Pseudosection
--------------------------------------
Here we plot the observed DC data in 3D pseudosection.
To use this utility, you must have Python's *plotly* package.
Here, we represent the DC data as apparent conductivities.
.. GENERATED FROM PYTHON SOURCE LINES 136-166
.. code-block:: Python
# Convert predicted data to apparent conductivities
apparent_conductivity = 1 / apparent_resistivity_from_voltage(
dc_data.survey,
dc_data.dobs,
)
if has_plotly:
# Plot DC Data
fig = plot_3d_pseudosection(
dc_data.survey, apparent_conductivity, scale="log", units="S/m"
)
fig.update_layout(
title_text="Apparent Conductivity",
title_x=0.5,
title_font_size=24,
width=650,
height=500,
scene_camera=dict(
center=dict(x=0, y=0, z=-0.4), eye=dict(x=1.5, y=-1.5, z=1.8)
),
)
plotly.io.show(fig)
else:
print("INSTALL 'PLOTLY' TO VISUALIZE 3D PSEUDOSECTIONS")
.. raw:: html
:file: images/sphx_glr_plot_inv_3_dcip3d_001.html
.. GENERATED FROM PYTHON SOURCE LINES 167-174
Plot Observed IP Data in Pseudosection
--------------------------------------
Here we plot the observed IP data in 3D pseudosection.
To use this utility, you must have Python's *plotly* package.
Here, we represent the IP data as apparent chargeabilities.
.. GENERATED FROM PYTHON SOURCE LINES 174-203
.. code-block:: Python
if has_plotly:
# Plot IP Data
fig = plot_3d_pseudosection(
ip_data.survey,
ip_data.dobs,
scale="linear",
units="V/V",
vlim=[0, np.max(ip_data.dobs)],
marker_opts={"colorscale": "plasma"},
)
fig.update_layout(
title_text="Apparent Chargeability",
title_x=0.5,
title_font_size=24,
width=650,
height=500,
scene_camera=dict(
center=dict(x=0, y=0, z=-0.4), eye=dict(x=1.5, y=-1.5, z=1.8)
),
)
plotly.io.show(fig)
else:
print("INSTALL 'PLOTLY' TO VISUALIZE 3D PSEUDOSECTIONS")
.. raw:: html
:file: images/sphx_glr_plot_inv_3_dcip3d_002.html
.. GENERATED FROM PYTHON SOURCE LINES 204-213
Assign Uncertainties
--------------------
Inversion with SimPEG requires that we define the uncertainties on our data.
This represents our estimate of the standard deviation of the
noise in our data. For DC data, the uncertainties are 10% of the absolute value.
For IP data, the uncertainties are 5e-3 V/V.
.. GENERATED FROM PYTHON SOURCE LINES 213-218
.. code-block:: Python
dc_data.standard_deviation = 0.1 * np.abs(dc_data.dobs)
ip_data.standard_deviation = 5e-3 * np.ones_like(ip_data.dobs)
.. GENERATED FROM PYTHON SOURCE LINES 219-225
Create Tree Mesh
----------------
Here, we create the Tree mesh that will be used to invert both DC
resistivity and IP data.
.. GENERATED FROM PYTHON SOURCE LINES 225-262
.. code-block:: Python
dh = 25.0 # base cell width
dom_width_x = 6000.0 # domain width x
dom_width_y = 6000.0 # domain width y
dom_width_z = 4000.0 # domain width z
nbcx = 2 ** int(np.round(np.log(dom_width_x / dh) / np.log(2.0))) # num. base cells x
nbcy = 2 ** int(np.round(np.log(dom_width_y / dh) / np.log(2.0))) # num. base cells y
nbcz = 2 ** int(np.round(np.log(dom_width_z / dh) / np.log(2.0))) # num. base cells z
# Define the base mesh
hx = [(dh, nbcx)]
hy = [(dh, nbcy)]
hz = [(dh, nbcz)]
mesh = TreeMesh([hx, hy, hz], x0="CCN")
# Mesh refinement based on topography
k = np.sqrt(np.sum(topo_xyz[:, 0:2] ** 2, axis=1)) < 1200
mesh = refine_tree_xyz(
mesh, topo_xyz[k, :], octree_levels=[0, 6, 8], method="surface", finalize=False
)
# Mesh refinement near sources and receivers.
electrode_locations = np.r_[
dc_data.survey.locations_a,
dc_data.survey.locations_b,
dc_data.survey.locations_m,
dc_data.survey.locations_n,
]
unique_locations = np.unique(electrode_locations, axis=0)
mesh = refine_tree_xyz(
mesh, unique_locations, octree_levels=[4, 6, 4], method="radial", finalize=False
)
# Finalize the mesh
mesh.finalize()
.. rst-class:: sphx-glr-script-out
.. code-block:: none
/home/vsts/work/1/s/tutorials/06-ip/plot_inv_3_dcip3d.py:239: 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/tutorials/06-ip/plot_inv_3_dcip3d.py:243: 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.
/home/vsts/work/1/s/tutorials/06-ip/plot_inv_3_dcip3d.py:255: DeprecationWarning:
The radial option is deprecated as of `0.9.0` please update your code to use the `TreeMesh.refine_points` functionality. It will be removed in a future version of discretize.
.. GENERATED FROM PYTHON SOURCE LINES 263-271
Project Electrodes to Discretized Topography
--------------------------------------------
It is important that electrodes are not modeled as being in the air. Even if the
electrodes are properly located along surface topography, they may lie above
the discretized topography. This step is carried out to ensure all electrodes
lie on the discretized surface.
.. GENERATED FROM PYTHON SOURCE LINES 271-287
.. code-block:: Python
# Find cells that lie below surface topography
ind_active = active_from_xyz(mesh, topo_xyz)
# Extract survey from data object
dc_survey = dc_data.survey
ip_survey = ip_data.survey
# Shift electrodes to the surface of discretized topography
dc_survey.drape_electrodes_on_topography(mesh, ind_active, option="top")
ip_survey.drape_electrodes_on_topography(mesh, ind_active, option="top")
# Reset survey in data object
dc_data.survey = dc_survey
ip_data.survey = ip_survey
.. GENERATED FROM PYTHON SOURCE LINES 288-297
Starting/Reference Model and Mapping on OcTree Mesh
---------------------------------------------------
Here, we create starting and/or reference models for the DC inversion as
well as the mapping from the model space to the active cells. Starting and
reference models can be a constant background value or contain a-priori
structures. Here, the starting model is the natural log of 0.01 S/m.
.. GENERATED FROM PYTHON SOURCE LINES 297-312
.. code-block:: Python
# Define conductivity model in S/m (or resistivity model in Ohm m)
air_conductivity = np.log(1e-8)
background_conductivity = np.log(1e-2)
# Define the mapping from active cells to the entire domain
active_map = maps.InjectActiveCells(mesh, ind_active, np.exp(air_conductivity))
nC = int(ind_active.sum())
# Define the mapping from the model to the conductivity of the entire domain
conductivity_map = active_map * maps.ExpMap()
# Define starting model
starting_conductivity_model = background_conductivity * np.ones(nC)
.. GENERATED FROM PYTHON SOURCE LINES 313-319
Define the Physics of the DC Simulation
---------------------------------------
Here, we define the physics of the DC resistivity simulation.
.. GENERATED FROM PYTHON SOURCE LINES 319-324
.. code-block:: Python
dc_simulation = dc.Simulation3DNodal(
mesh, survey=dc_survey, sigmaMap=conductivity_map, storeJ=True
)
.. GENERATED FROM PYTHON SOURCE LINES 325-335
Define DC Inverse Problem
-------------------------
The inverse problem is defined by 3 things:
1) Data Misfit: a measure of how well our recovered model explains the field data
2) Regularization: constraints placed on the recovered model and a priori information
3) Optimization: the numerical approach used to solve the inverse problem
.. GENERATED FROM PYTHON SOURCE LINES 335-362
.. code-block:: Python
# Define the data misfit. Here the data misfit is the L2 norm of the weighted
# residual between the observed data and the data predicted for a given model.
# Within the data misfit, the residual between predicted and observed data are
# normalized by the data's standard deviation.
dc_data_misfit = data_misfit.L2DataMisfit(data=dc_data, simulation=dc_simulation)
# Define the regularization (model objective function)
dc_regularization = regularization.WeightedLeastSquares(
mesh,
active_cells=ind_active,
reference_model=starting_conductivity_model,
)
dc_regularization.reference_model_in_smooth = (
True # Include reference model in smoothness
)
# Define how the optimization problem is solved.
dc_optimization = optimization.InexactGaussNewton(maxIter=15, maxIterCG=30, tolCG=1e-2)
# Here we define the inverse problem that is to be solved
dc_inverse_problem = inverse_problem.BaseInvProblem(
dc_data_misfit, dc_regularization, dc_optimization
)
.. GENERATED FROM PYTHON SOURCE LINES 363-371
Define DC Inversion Directives
------------------------------
Here we define any directives that are carried out during the inversion. This
includes the cooling schedule for the trade-off parameter (beta), stopping
criteria for the inversion and saving inversion results at each iteration.
.. GENERATED FROM PYTHON SOURCE LINES 371-402
.. code-block:: Python
# Apply and update sensitivity weighting as the model updates
update_sensitivity_weighting = directives.UpdateSensitivityWeights()
# Defining a starting value for the trade-off parameter (beta) between the data
# misfit and the regularization.
starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=1e1)
# Set the rate of reduction in trade-off parameter (beta) each time the
# the inverse problem is solved. And set the number of Gauss-Newton iterations
# for each trade-off paramter value.
beta_schedule = directives.BetaSchedule(coolingFactor=2.5, coolingRate=2)
# Options for outputting recovered models and predicted data for each beta.
save_iteration = directives.SaveOutputEveryIteration(save_txt=False)
# Setting a stopping criteria for the inversion.
target_misfit = directives.TargetMisfit(chifact=1)
# Apply and update preconditioner as the model updates
update_jacobi = directives.UpdatePreconditioner()
directives_list = [
update_sensitivity_weighting,
starting_beta,
beta_schedule,
save_iteration,
target_misfit,
update_jacobi,
]
.. GENERATED FROM PYTHON SOURCE LINES 403-409
Running the DC Inversion
------------------------
To define the inversion object, we need to define the inversion problem and
the set of directives. We can then run the inversion.
.. GENERATED FROM PYTHON SOURCE LINES 409-418
.. code-block:: Python
# Here we combine the inverse problem and the set of directives
dc_inversion = inversion.BaseInversion(
dc_inverse_problem, directiveList=directives_list
)
# Run inversion
recovered_conductivity_model = dc_inversion.run(starting_conductivity_model)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Running inversion with SimPEG v0.23.1.dev28+g33c4221d8
/home/vsts/work/1/s/simpeg/base/pde_simulation.py:490: 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.
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver, and solver_opts as the Simulation3DNodal problem***
/usr/share/miniconda/envs/simpeg-test/lib/python3.10/site-packages/pymatsolver/solvers.py:415: FutureWarning:
In Future pymatsolver v0.4.0, passing a vector of shape (n, 1) to the solve method will return an array with shape (n, 1), instead of always returning a flattened array. This is to be consistent with numpy.linalg.solve broadcasting.
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.18e-02 4.29e+04 0.00e+00 4.29e+04 4.82e+03 0
1 3.18e-02 1.07e+04 2.68e+05 1.93e+04 2.48e+02 0
2 1.27e-02 1.14e+04 2.43e+05 1.45e+04 9.75e+02 0
3 1.27e-02 5.47e+03 5.51e+05 1.25e+04 1.42e+02 0
4 5.09e-03 6.10e+03 4.96e+05 8.63e+03 6.01e+02 0
5 5.09e-03 2.55e+03 9.41e+05 7.34e+03 7.61e+01 0
6 2.04e-03 2.80e+03 8.85e+05 4.61e+03 3.41e+02 0
7 2.04e-03 1.05e+03 1.42e+06 3.94e+03 3.78e+01 0
8 8.15e-04 1.14e+03 1.37e+06 2.25e+03 1.82e+02 0
9 8.15e-04 4.04e+02 1.92e+06 1.97e+03 1.79e+01 0
10 3.26e-04 4.30e+02 1.88e+06 1.04e+03 9.17e+01 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 4.2936e+03
1 : |xc-x_last| = 2.1365e+01 <= tolX*(1+|x0|) = 8.6952e+01
0 : |proj(x-g)-x| = 9.1726e+01 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.1726e+01 <= 1e3*eps = 1.0000e-02
0 : maxIter = 15 <= iter = 11
------------------------- DONE! -------------------------
.. GENERATED FROM PYTHON SOURCE LINES 419-422
Recreate True Conductivity Model
--------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 422-439
.. code-block:: Python
background_value = 1e-2
conductor_value = 1e-1
resistor_value = 1e-3
true_conductivity_model = background_value * np.ones(nC)
ind_conductor = model_builder.get_indices_sphere(
np.r_[-350.0, 0.0, -300.0], 160.0, mesh.cell_centers[ind_active, :]
)
true_conductivity_model[ind_conductor] = conductor_value
ind_resistor = model_builder.get_indices_sphere(
np.r_[350.0, 0.0, -300.0], 160.0, mesh.cell_centers[ind_active, :]
)
true_conductivity_model[ind_resistor] = resistor_value
true_conductivity_model_log10 = np.log10(true_conductivity_model)
.. GENERATED FROM PYTHON SOURCE LINES 440-443
Plotting True and Recovered Conductivity Model
----------------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 443-505
.. code-block:: Python
# Plot True Model
fig = plt.figure(figsize=(10, 4))
plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)
ax1 = fig.add_axes([0.15, 0.15, 0.67, 0.75])
mesh.plot_slice(
plotting_map * true_conductivity_model_log10,
ax=ax1,
normal="Y",
ind=int(len(mesh.h[1]) / 2),
grid=False,
clim=(true_conductivity_model_log10.min(), true_conductivity_model_log10.max()),
pcolor_opts={"cmap": mpl.cm.viridis},
)
ax1.set_title("True Conductivity Model")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax1.set_xlim([-1000, 1000])
ax1.set_ylim([-1000, 0])
ax2 = fig.add_axes([0.84, 0.15, 0.03, 0.75])
norm = mpl.colors.Normalize(
vmin=true_conductivity_model_log10.min(), vmax=true_conductivity_model_log10.max()
)
cbar = mpl.colorbar.ColorbarBase(
ax2, cmap=mpl.cm.viridis, norm=norm, orientation="vertical", format="$10^{%.1f}$"
)
cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12)
# Plot recovered model
recovered_conductivity_model_log10 = np.log10(np.exp(recovered_conductivity_model))
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_axes([0.15, 0.15, 0.67, 0.75])
mesh.plot_slice(
plotting_map * recovered_conductivity_model_log10,
ax=ax1,
normal="Y",
ind=int(len(mesh.h[1]) / 2),
grid=False,
clim=(true_conductivity_model_log10.min(), true_conductivity_model_log10.max()),
pcolor_opts={"cmap": mpl.cm.viridis},
)
ax1.set_title("Recovered Conductivity Model")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax1.set_xlim([-1000, 1000])
ax1.set_ylim([-1000, 0])
ax2 = fig.add_axes([0.84, 0.15, 0.03, 0.75])
norm = mpl.colors.Normalize(
vmin=true_conductivity_model_log10.min(), vmax=true_conductivity_model_log10.max()
)
cbar = mpl.colorbar.ColorbarBase(
ax2, cmap=mpl.cm.viridis, norm=norm, orientation="vertical", format="$10^{%.1f}$"
)
cbar.set_label("Conductivity [S/m]", rotation=270, labelpad=15, size=12)
plt.show()
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_003.png
:alt: True Conductivity Model
:srcset: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_003.png
:class: sphx-glr-multi-img
*
.. image-sg:: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_004.png
:alt: Recovered Conductivity Model
:srcset: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_004.png
:class: sphx-glr-multi-img
.. GENERATED FROM PYTHON SOURCE LINES 506-513
Plotting Normalized Data Misfit or Predicted DC Data
----------------------------------------------------
To see how well the recovered model reproduces the observed data,
it is a good idea to compare the predicted and observed data.
Here, we accomplish this by plotting the normalized misfit.
.. GENERATED FROM PYTHON SOURCE LINES 513-548
.. code-block:: Python
# Predicted data from recovered model
dpred_dc = dc_inverse_problem.dpred
# Compute the normalized data misfit
dc_normalized_misfit = (dc_data.dobs - dpred_dc) / dc_data.standard_deviation
if has_plotly:
# Plot IP Data
fig = plot_3d_pseudosection(
dc_data.survey,
dc_normalized_misfit,
scale="linear",
units="",
vlim=[-2, 2],
plane_distance=15,
)
fig.update_layout(
title_text="Normalized Data Misfit",
title_x=0.5,
title_font_size=24,
width=650,
height=500,
scene_camera=dict(
center=dict(x=0, y=0, z=-0.4), eye=dict(x=1.5, y=-1.5, z=1.8)
),
)
plotly.io.show(fig)
else:
print("INSTALL 'PLOTLY' TO VISUALIZE 3D PSEUDOSECTIONS")
.. raw:: html
:file: images/sphx_glr_plot_inv_3_dcip3d_005.html
.. GENERATED FROM PYTHON SOURCE LINES 549-558
Starting/Reference Model for IP Inversion
-----------------------------------------
Here, we would create starting and/or reference models for the IP inversion as
well as the mapping from the model space to the active cells. Starting and
reference models can be a constant background value or contain a-priori
structures. Here, the starting model is the 1e-6 V/V.
.. GENERATED FROM PYTHON SOURCE LINES 558-573
.. code-block:: Python
# Define chargeability model in V/V
air_chargeability = 0.0
background_chargeability = 1e-6
active_map = maps.InjectActiveCells(mesh, ind_active, air_chargeability)
nC = int(ind_active.sum())
chargeability_map = active_map
# Define starting model
starting_chargeability_model = background_chargeability * np.ones(nC)
.. GENERATED FROM PYTHON SOURCE LINES 574-583
Define the Physics of the IP Simulation
---------------------------------------
Here, we define the physics of the IP problem. For the chargeability, we
require a mapping from the model space to the entire mesh. For the background
conductivity/resistivity, we require the conductivity/resistivity on the
entire mesh.
.. GENERATED FROM PYTHON SOURCE LINES 583-592
.. code-block:: Python
ip_simulation = ip.Simulation3DNodal(
mesh,
survey=ip_survey,
etaMap=chargeability_map,
sigma=conductivity_map * recovered_conductivity_model,
storeJ=True,
)
.. GENERATED FROM PYTHON SOURCE LINES 593-598
Define IP Inverse Problem
-------------------------
Here we define the inverse problem in the same manner as the DC inverse problem.
.. GENERATED FROM PYTHON SOURCE LINES 598-623
.. code-block:: Python
# Define the data misfit (Here we use weighted L2-norm)
ip_data_misfit = data_misfit.L2DataMisfit(data=ip_data, simulation=ip_simulation)
# Define the regularization (model objective function)
ip_regularization = regularization.WeightedLeastSquares(
mesh,
active_cells=ind_active,
mapping=maps.IdentityMap(nP=nC),
alpha_s=0.01,
alpha_x=1,
alpha_y=1,
alpha_z=1,
)
# Define how the optimization problem is solved.
ip_optimization = optimization.ProjectedGNCG(
maxIter=15, lower=0.0, upper=10, maxIterCG=30, tolCG=1e-2
)
# Here we define the inverse problem that is to be solved
ip_inverse_problem = inverse_problem.BaseInvProblem(
ip_data_misfit, ip_regularization, ip_optimization
)
.. GENERATED FROM PYTHON SOURCE LINES 624-629
Define IP Inversion Directives
------------------------------
Here we define the directives in the same manner as the DC inverse problem.
.. GENERATED FROM PYTHON SOURCE LINES 629-647
.. code-block:: Python
update_sensitivity_weighting = directives.UpdateSensitivityWeights(threshold_value=1e-3)
starting_beta = directives.BetaEstimate_ByEig(beta0_ratio=1e2)
beta_schedule = directives.BetaSchedule(coolingFactor=2.5, coolingRate=1)
save_iteration = directives.SaveOutputEveryIteration(save_txt=False)
target_misfit = directives.TargetMisfit(chifact=1.0)
update_jacobi = directives.UpdatePreconditioner()
directives_list = [
update_sensitivity_weighting,
starting_beta,
beta_schedule,
save_iteration,
target_misfit,
update_jacobi,
]
.. GENERATED FROM PYTHON SOURCE LINES 648-651
Running the IP Inversion
------------------------
.. GENERATED FROM PYTHON SOURCE LINES 651-660
.. code-block:: Python
# Here we combine the inverse problem and the set of directives
ip_inversion = inversion.BaseInversion(
ip_inverse_problem, directiveList=directives_list
)
# Run inversion
recovered_chargeability_model = ip_inversion.run(starting_chargeability_model)
.. rst-class:: sphx-glr-script-out
.. code-block:: none
Running inversion with SimPEG v0.23.1.dev28+g33c4221d8
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 will set Regularization.reference_model to m0.
simpeg.InvProblem will set Regularization.reference_model to m0.
/home/vsts/work/1/s/simpeg/base/pde_simulation.py:490: 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.
simpeg.InvProblem is setting bfgsH0 to the inverse of the eval2Deriv.
***Done using same Solver, and solver_opts as the Simulation3DNodal problem***
/usr/share/miniconda/envs/simpeg-test/lib/python3.10/site-packages/pymatsolver/solvers.py:415: FutureWarning:
In Future pymatsolver v0.4.0, passing a vector of shape (n, 1) to the solve method will return an array with shape (n, 1), instead of always returning a flattened array. This is to be consistent with numpy.linalg.solve broadcasting.
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 2.23e+01 7.47e+03 0.00e+00 7.47e+03 9.88e+02 0
------------------------- STOP! -------------------------
1 : |fc-fOld| = 0.0000e+00 <= tolF*(1+|f0|) = 7.4694e+02
0 : |xc-x_last| = 1.4540e+00 <= tolX*(1+|x0|) = 1.0002e-01
0 : |proj(x-g)-x| = 9.8696e+02 <= tolG = 1.0000e-01
0 : |proj(x-g)-x| = 9.8696e+02 <= 1e3*eps = 1.0000e-02
0 : maxIter = 15 <= iter = 1
------------------------- DONE! -------------------------
.. GENERATED FROM PYTHON SOURCE LINES 661-664
Recreate True Chargeability Model
---------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 664-675
.. code-block:: Python
background_value = 1e-6
chargeable_value = 1e-1
true_chargeability_model = background_value * np.ones(nC)
ind_chargeable = model_builder.get_indices_sphere(
np.r_[-350.0, 0.0, -300.0], 160.0, mesh.cell_centers[ind_active, :]
)
true_chargeability_model[ind_chargeable] = chargeable_value
.. GENERATED FROM PYTHON SOURCE LINES 676-679
Plot True and Recovered Chargeability Model
--------------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 679-739
.. code-block:: Python
# Plot True Model
fig = plt.figure(figsize=(10, 4))
plotting_map = maps.InjectActiveCells(mesh, ind_active, np.nan)
ax1 = fig.add_axes([0.15, 0.15, 0.67, 0.75])
mesh.plot_slice(
plotting_map * true_chargeability_model,
ax=ax1,
normal="Y",
ind=int(len(mesh.h[1]) / 2),
grid=False,
clim=(true_chargeability_model.min(), true_chargeability_model.max()),
pcolor_opts={"cmap": mpl.cm.plasma},
)
ax1.set_title("True Chargeability Model")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax1.set_xlim([-1000, 1000])
ax1.set_ylim([-1000, 0])
ax2 = fig.add_axes([0.84, 0.15, 0.03, 0.75])
norm = mpl.colors.Normalize(
vmin=true_chargeability_model.min(), vmax=true_chargeability_model.max()
)
cbar = mpl.colorbar.ColorbarBase(
ax2, cmap=mpl.cm.plasma, norm=norm, orientation="vertical", format="%.2f"
)
cbar.set_label("Intrinsic Chargeability [V/V]", rotation=270, labelpad=15, size=12)
# Plot Recovered Model
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_axes([0.15, 0.15, 0.67, 0.75])
mesh.plot_slice(
plotting_map * recovered_chargeability_model,
ax=ax1,
normal="Y",
ind=int(len(mesh.h[1]) / 2),
grid=False,
clim=(true_chargeability_model.min(), true_chargeability_model.max()),
pcolor_opts={"cmap": mpl.cm.plasma},
)
ax1.set_title("Recovered Chargeability Model")
ax1.set_xlabel("x (m)")
ax1.set_ylabel("z (m)")
ax1.set_xlim([-1000, 1000])
ax1.set_ylim([-1000, 0])
ax2 = fig.add_axes([0.84, 0.15, 0.03, 0.75])
norm = mpl.colors.Normalize(
vmin=true_chargeability_model.min(), vmax=true_chargeability_model.max()
)
cbar = mpl.colorbar.ColorbarBase(
ax2, cmap=mpl.cm.plasma, norm=norm, orientation="vertical", format="%.2f"
)
cbar.set_label("Intrinsic Chargeability [V/V]", rotation=270, labelpad=15, size=12)
plt.show()
.. rst-class:: sphx-glr-horizontal
*
.. image-sg:: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_006.png
:alt: True Chargeability Model
:srcset: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_006.png
:class: sphx-glr-multi-img
*
.. image-sg:: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_007.png
:alt: Recovered Chargeability Model
:srcset: /content/user-guide/tutorials/06-ip/images/sphx_glr_plot_inv_3_dcip3d_007.png
:class: sphx-glr-multi-img
.. GENERATED FROM PYTHON SOURCE LINES 740-743
Plotting Normalized Data Misfit or Predicted IP Data
----------------------------------------------------
.. GENERATED FROM PYTHON SOURCE LINES 743-776
.. code-block:: Python
# Predicted data from recovered model
dpred_ip = ip_inverse_problem.dpred
# Normalized misfit
ip_normalized_misfit = (ip_data.dobs - dpred_ip) / ip_data.standard_deviation
if has_plotly:
fig = plot_3d_pseudosection(
ip_data.survey,
ip_normalized_misfit,
scale="linear",
units="",
vlim=[-2, 2],
plane_distance=15,
marker_opts={"colorscale": "plasma"},
)
fig.update_layout(
title_text="Normalized Data Misfit",
title_x=0.5,
title_font_size=24,
width=650,
height=500,
scene_camera=dict(
center=dict(x=0, y=0, z=-0.4), eye=dict(x=1.5, y=-1.5, z=1.8)
),
)
plotly.io.show(fig)
else:
print("INSTALL 'PLOTLY' TO VISUALIZE 3D PSEUDOSECTIONS")
.. raw:: html
:file: images/sphx_glr_plot_inv_3_dcip3d_008.html
.. rst-class:: sphx-glr-timing
**Total running time of the script:** (5 minutes 31.445 seconds)
**Estimated memory usage:** 522 MB
.. _sphx_glr_download_content_user-guide_tutorials_06-ip_plot_inv_3_dcip3d.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_3_dcip3d.ipynb <plot_inv_3_dcip3d.ipynb>`
.. container:: sphx-glr-download sphx-glr-download-python
:download:`Download Python source code: plot_inv_3_dcip3d.py <plot_inv_3_dcip3d.py>`
.. container:: sphx-glr-download sphx-glr-download-zip
:download:`Download zipped: plot_inv_3_dcip3d.zip <plot_inv_3_dcip3d.zip>`
.. only:: html
.. rst-class:: sphx-glr-signature
`Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_