.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "content/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 <sphx_glr_download_content_examples_20-published_plot_tomo_joint_with_volume.py>` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_content_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/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png :alt: Model Transform :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_001.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png :alt: plot tomo joint with volume :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_002.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png :alt: plot tomo joint with volume :srcset: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_003.png :class: sphx-glr-multi-img * .. image-sg:: /content/examples/20-published/images/sphx_glr_plot_tomo_joint_with_volume_004.png :alt: true, vol: 6.240e-01, recovered(no Volume term), vol: 2.652e-02 , recovered(with Volume term), vol: 3.795e-02 :srcset: /content/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 True Volume: 0.6240000000000001 Running inversion with SimPEG v0.22.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 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 2.50e-01 8.48e+04 0.00e+00 8.48e+04 3.83e+01 0 1 2.50e-01 4.24e+03 9.97e-01 4.24e+03 3.67e+01 0 /home/vsts/work/1/s/simpeg/maps/_property_maps.py:1456: UserWarning: Maximum number of iterations reached 2 2.50e-01 6.98e+02 8.59e-01 6.98e+02 3.24e+01 0 3 2.50e-01 6.67e+02 8.35e-01 6.67e+02 3.22e+01 4 4 2.50e-01 3.43e+02 7.28e-01 3.44e+02 2.57e+01 0 5 2.50e-01 8.63e+01 7.20e-01 8.64e+01 1.57e+01 0 6 2.50e-01 5.73e+01 7.23e-01 5.74e+01 1.51e+01 0 7 2.50e-01 4.99e+01 7.23e-01 5.00e+01 1.29e+01 0 8 2.50e-01 3.15e+01 7.26e-01 3.17e+01 1.10e+01 0 9 2.50e-01 2.31e+01 7.28e-01 2.32e+01 1.03e+01 0 10 2.50e-01 1.57e+01 7.32e-01 1.59e+01 8.89e+00 0 11 2.50e-01 1.26e+01 7.33e-01 1.28e+01 8.47e+00 0 12 2.50e-01 9.58e+00 7.36e-01 9.76e+00 7.66e+00 0 13 2.50e-01 6.81e+00 7.38e-01 6.99e+00 7.18e+00 0 14 2.50e-01 4.86e+00 7.38e-01 5.04e+00 6.23e+00 0 15 2.50e-01 3.60e+00 7.38e-01 3.78e+00 5.31e+00 0 Skip BFGS ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 1.2589e+00 <= tolF*(1+|f0|) = 8.4765e+03 0 : |xc-x_last| = 4.5379e-01 <= tolX*(1+|x0|) = 1.0000e-01 0 : |proj(x-g)-x| = 5.3051e+00 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 5.3051e+00 <= 1e3*eps = 1.0000e-02 1 : maxIter = 15 <= iter = 15 ------------------------- DONE! ------------------------- Total recovered volume (no vol misfit term in inversion): 0.026515074507646754 Running inversion with SimPEG v0.22.0 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 2.50e-01 1.20e+05 0.00e+00 1.20e+05 1.56e+01 0 1 2.50e-01 7.79e+03 8.43e-01 7.79e+03 2.71e+01 0 2 2.50e-01 7.78e+03 8.15e-01 7.78e+03 2.65e+01 4 3 2.50e-01 7.77e+03 7.91e-01 7.78e+03 2.61e+01 4 Skip BFGS 4 2.50e-01 7.74e+03 7.57e-01 7.74e+03 2.55e+01 4 5 2.50e-01 7.64e+03 7.06e-01 7.64e+03 2.48e+01 3 6 2.50e-01 7.49e+03 6.73e-01 7.49e+03 2.39e+01 3 ------------------------------------------------------------------ 0 : ft = 7.4923e+03 <= alp*descent = 7.4873e+03 1 : maxIterLS = 10 <= iterLS = 10 ------------------------- End Linesearch ------------------------- The linesearch got broken. Boo. Total volume (vol misfit term in inversion): 0.03795357151511156 | .. 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 35.356 seconds) **Estimated memory usage:** 9 MB .. _sphx_glr_download_content_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 <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 <plot_tomo_joint_with_volume.py>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_