.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/user-guide/examples/20-published/plot_tomo_joint_with_volume.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_20-published_plot_tomo_joint_with_volume.py: Straight Ray with Volume Data Misfit Term ========================================= Based on the SEG abstract Heagy, Cockett and Oldenburg, 2014. Heagy, L. J., Cockett, A. R., & Oldenburg, D. W. (2014, August 5). Parametrized Inversion Framework for Proppant Volume in a Hydraulically Fractured Reservoir. SEG Technical Program Expanded Abstracts 2014. Society of Exploration Geophysicists. doi:10.1190/segam2014-1639.1 This example is a simple joint inversion that consists of a - data misfit for the tomography problem - data misfit for the volume of the inclusions (uses the effective medium theory mapping) - model regularization .. GENERATED FROM PYTHON SOURCE LINES 20-216 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png :alt: Model Transform :srcset: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png :alt: plot tomo joint with volume :srcset: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png :alt: plot tomo joint with volume :srcset: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png :class: sphx-glr-multi-img * .. image-sg:: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_004.png :alt: true, vol: 6.240e-01, recovered(no Volume term), vol: 1.593e-02 , recovered(with Volume term), vol: 2.838e-02 :srcset: /content/user-guide/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_004.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/vsts/work/1/s/simpeg/utils/model_builder.py:188: 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) True Volume: 0.6240000000000001 /home/vsts/work/1/s/examples/20-published/plot_tomo_joint_with_volume.py:157: FutureWarning: The defaults for ProjectedGNCG will change in SimPEG 0.26.0. If you want to maintain the previous behavior, explicitly set 'cg_atol=1E-3' and 'cg_rtol=0.0'. 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 2.50e-01 8.47e+04 0.00e+00 8.47e+04 0 inf inf 1 2.50e-01 4.39e+03 9.96e-01 4.39e+03 5.32e-01 0 5 3.84e-03 8.24e+00 /home/vsts/work/1/s/simpeg/maps/_property_maps.py:1507: UserWarning: Maximum number of iterations reached 2 2.50e-01 4.84e+02 8.10e-01 4.84e+02 3.28e+01 0 5 8.10e-02 5.81e+00 3 2.50e-01 1.71e+02 6.92e-01 1.71e+02 2.57e+01 2 5 2.63e+00 2.39e+02 Skip BFGS 4 2.50e-01 7.60e+01 6.68e-01 7.62e+01 1.60e+01 0 5 2.88e-01 3.35e+01 5 2.50e-01 2.18e+01 6.74e-01 2.19e+01 1.28e+01 0 5 1.96e-01 2.53e+01 6 2.50e-01 1.37e+01 6.76e-01 1.39e+01 9.99e+00 0 5 3.24e-01 1.43e+01 Skip BFGS 7 2.50e-01 4.95e+00 6.78e-01 5.12e+00 9.30e+00 0 5 2.70e-01 1.33e+01 8 2.50e-01 3.18e+00 6.78e-01 3.34e+00 6.66e+00 0 5 4.88e-01 9.19e+00 9 2.50e-01 1.58e+00 6.78e-01 1.75e+00 5.78e+00 0 5 2.51e-01 6.90e+00 10 2.50e-01 9.09e-01 6.78e-01 1.08e+00 3.49e+00 0 5 1.07e+00 7.61e+00 11 2.50e-01 6.07e-01 6.78e-01 7.76e-01 3.32e+00 0 5 5.42e-01 3.77e+00 12 2.50e-01 4.81e-01 6.78e-01 6.50e-01 2.93e+00 0 5 8.02e-01 2.75e+00 13 2.50e-01 3.64e-01 6.78e-01 5.33e-01 2.46e+00 0 5 1.31e+00 3.55e+00 14 2.50e-01 3.16e-01 6.78e-01 4.85e-01 2.33e+00 0 5 5.16e-01 1.94e+00 15 2.50e-01 2.79e-01 6.78e-01 4.48e-01 1.68e+00 0 5 1.93e+00 3.37e+00 Skip BFGS ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 3.7178e-02 <= tolF*(1+|f0|) = 8.4750e+03 1 : |xc-x_last| = 9.6536e-02 <= tolX*(1+|x0|) = 1.0000e-01 0 : |proj(x-g)-x| = 2.1985e+00 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 2.1985e+00 <= 1e3*eps = 1.0000e-02 1 : maxIter = 15 <= iter = 15 ------------------------- DONE! ------------------------- Total recovered volume (no vol misfit term in inversion): 0.01592508882068877 /home/vsts/work/1/s/examples/20-published/plot_tomo_joint_with_volume.py:172: FutureWarning: The defaults for ProjectedGNCG will change in SimPEG 0.26.0. If you want to maintain the previous behavior, explicitly set 'cg_atol=1E-3' and 'cg_rtol=0.0'. 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 2.50e-01 1.20e+05 0.00e+00 1.20e+05 0 inf inf 1 2.50e-01 7.70e+03 8.51e-01 7.70e+03 2.21e+01 0 5 4.88e-02 7.17e+01 2 2.50e-01 7.67e+03 8.36e-01 7.67e+03 2.69e+01 5 5 3.71e-01 1.80e+02 3 2.50e-01 7.66e+03 8.19e-01 7.66e+03 2.64e+01 4 5 5.93e-01 3.12e+02 Skip BFGS 4 2.50e-01 7.64e+03 8.03e-01 7.64e+03 2.61e+01 4 5 4.94e-01 2.77e+02 ------------------------------------------------------------------ 0 : ft = 7.6371e+03 <= alp*descent = 7.6356e+03 1 : maxIterLS = 10 <= iterLS = 10 ------------------------- End Linesearch ------------------------- The linesearch got broken. Boo. Total volume (vol misfit term in inversion): 0.028384850995600085 | .. code-block:: Python import numpy as np import scipy.sparse as sp import matplotlib.pyplot as plt from simpeg.seismic import straight_ray_tomography as tomo import discretize from simpeg import ( maps, utils, regularization, optimization, inverse_problem, inversion, data_misfit, objective_function, ) class Volume(objective_function.BaseObjectiveFunction): r""" A regularization on the volume integral of the model .. math:: \phi_v = || \int_V m dV - \text{knownVolume} ||^2 """ def __init__(self, mesh, knownVolume=0.0, **kwargs): self.mesh = mesh self.knownVolume = knownVolume super().__init__(**kwargs) @property def knownVolume(self): """known volume""" return self._knownVolume @knownVolume.setter def knownVolume(self, value): self._knownVolume = utils.validate_float("knownVolume", value, min_val=0.0) def __call__(self, m): return (self.estVol(m) - self.knownVolume) ** 2 def estVol(self, m): return np.inner(self.mesh.cell_volumes, m) def deriv(self, m): # return (self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, m)) return ( 2 * self.mesh.cell_volumes * (self.knownVolume - np.inner(self.mesh.cell_volumes, m)) ) # factor of 2 from deriv of ||estVol - knownVol||^2 def deriv2(self, m, v=None): if v is not None: return 2 * utils.mkvc( self.mesh.cell_volumes * np.inner(self.mesh.cell_volumes, v) ) else: # TODO: this is inefficent. It is a fully dense matrix return 2 * sp.csc_matrix( np.outer(self.mesh.cell_volumes, self.mesh.cell_volumes) ) def run(plotIt=True): nC = 40 de = 1.0 h = np.ones(nC) * de / nC M = discretize.TensorMesh([h, h]) y = np.linspace(M.cell_centers_y[0], M.cell_centers_x[-1], int(np.floor(nC / 4))) rlocs = np.c_[0 * y + M.cell_centers_x[-1], y] rx = tomo.Rx(rlocs) source_list = [ tomo.Src(location=np.r_[M.cell_centers_x[0], yi], receiver_list=[rx]) for yi in y ] # phi model phi0 = 0 phi1 = 0.65 phitrue = utils.model_builder.create_block_in_wholespace( M.gridCC, [0.4, 0.6], [0.6, 0.4], [phi1, phi0] ) knownVolume = np.sum(phitrue * M.cell_volumes) print("True Volume: {}".format(knownVolume)) # Set up true conductivity model and plot the model transform sigma0 = np.exp(1) sigma1 = 1e4 if plotIt: fig, ax = plt.subplots(1, 1) sigmaMapTest = maps.SelfConsistentEffectiveMedium( nP=1000, sigma0=sigma0, sigma1=sigma1, rel_tol=1e-1, maxIter=150 ) testphis = np.linspace(0.0, 1.0, 1000) sigetest = sigmaMapTest * testphis ax.semilogy(testphis, sigetest) ax.set_title("Model Transform") ax.set_xlabel(r"$\varphi$") ax.set_ylabel(r"$\sigma$") sigmaMap = maps.SelfConsistentEffectiveMedium(M, sigma0=sigma0, sigma1=sigma1) # scale the slowness so it is on a ~linear scale slownessMap = maps.LogMap(M) * sigmaMap # set up the problem and survey survey = tomo.Survey(source_list) problem = tomo.Simulation(M, survey=survey, slownessMap=slownessMap) if plotIt: _, ax = plt.subplots(1, 1) cb = plt.colorbar(M.plot_image(phitrue, ax=ax)[0], ax=ax) survey.plot(ax=ax) cb.set_label(r"$\varphi$") # get observed data data = problem.make_synthetic_data(phitrue, relative_error=0.03, add_noise=True) dpred = problem.dpred(np.zeros(M.nC)) # objective function pieces reg = regularization.WeightedLeastSquares(M) dmis = data_misfit.L2DataMisfit(simulation=problem, data=data) dmisVol = Volume(mesh=M, knownVolume=knownVolume) beta = 0.25 maxIter = 15 # without the volume regularization opt = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0) opt.remember("xc") invProb = inverse_problem.BaseInvProblem(dmis, reg, opt, beta=beta) inv = inversion.BaseInversion(invProb) mopt1 = inv.run(np.zeros(M.nC) + 1e-16) print( "\nTotal recovered volume (no vol misfit term in inversion): " "{}".format(dmisVol(mopt1)) ) # with the volume regularization vol_multiplier = 9e4 reg2 = reg dmis2 = dmis + vol_multiplier * dmisVol opt2 = optimization.ProjectedGNCG(maxIter=maxIter, lower=0.0, upper=1.0) opt2.remember("xc") invProb2 = inverse_problem.BaseInvProblem(dmis2, reg2, opt2, beta=beta) inv2 = inversion.BaseInversion(invProb2) mopt2 = inv2.run(np.zeros(M.nC) + 1e-16) print("\nTotal volume (vol misfit term in inversion): {}".format(dmisVol(mopt2))) # plot results if plotIt: fig, ax = plt.subplots(1, 1) ax.plot(data.dobs) ax.plot(dpred) ax.plot(problem.dpred(mopt1), "o") ax.plot(problem.dpred(mopt2), "s") ax.legend(["dobs", "dpred0", "dpred w/o Vol", "dpred with Vol"]) fig, ax = plt.subplots(1, 3, figsize=(16, 4)) im0 = M.plot_image(phitrue, ax=ax[0])[0] im1 = M.plot_image(mopt1, ax=ax[1])[0] im2 = M.plot_image(mopt2, ax=ax[2])[0] for im in [im0, im1, im2]: im.set_clim([0.0, phi1]) plt.colorbar(im0, ax=ax[0]) plt.colorbar(im1, ax=ax[1]) plt.colorbar(im2, ax=ax[2]) ax[0].set_title("true, vol: {:1.3e}".format(knownVolume)) ax[1].set_title( "recovered(no Volume term), vol: {:1.3e} ".format(dmisVol(mopt1)) ) ax[2].set_title( "recovered(with Volume term), vol: {:1.3e} ".format(dmisVol(mopt2)) ) plt.tight_layout() if __name__ == "__main__": run() plt.show() .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 28.835 seconds) **Estimated memory usage:** 319 MB .. _sphx_glr_download_content_user-guide_examples_20-published_plot_tomo_joint_with_volume.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_tomo_joint_with_volume.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_tomo_joint_with_volume.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: plot_tomo_joint_with_volume.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_