.. 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 ` 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: 3.444e-02 , recovered(with Volume term), vol: 3.566e-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 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.27e+03 9.97e-01 4.27e+03 3.72e+01 0 /home/ssoler/simpeg/SimPEG/maps.py:1844: UserWarning: Maximum number of iterations reached 2 2.50e-01 1.03e+03 8.29e-01 1.03e+03 3.28e+01 0 3 2.50e-01 9.21e+02 8.19e-01 9.21e+02 3.26e+01 4 4 2.50e-01 2.81e+02 7.78e-01 2.81e+02 2.47e+01 1 5 2.50e-01 1.50e+02 7.48e-01 1.50e+02 1.68e+01 0 6 2.50e-01 7.16e+01 7.52e-01 7.18e+01 1.45e+01 0 7 2.50e-01 4.52e+01 7.58e-01 4.54e+01 1.60e+01 0 Skip BFGS 8 2.50e-01 3.15e+01 7.59e-01 3.17e+01 1.23e+01 0 9 2.50e-01 2.29e+01 7.67e-01 2.31e+01 9.78e+00 0 Skip BFGS 10 2.50e-01 2.26e+01 7.64e-01 2.28e+01 1.13e+01 0 11 2.50e-01 1.76e+01 7.67e-01 1.78e+01 9.10e+00 0 12 2.50e-01 1.50e+01 7.66e-01 1.52e+01 8.92e+00 1 13 2.50e-01 1.13e+01 7.69e-01 1.15e+01 8.03e+00 0 14 2.50e-01 9.66e+00 7.68e-01 9.85e+00 7.56e+00 0 15 2.50e-01 7.46e+00 7.69e-01 7.65e+00 7.20e+00 0 Skip BFGS ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 2.2026e+00 <= tolF*(1+|f0|) = 8.4760e+03 0 : |xc-x_last| = 5.8258e-01 <= tolX*(1+|x0|) = 1.0000e-01 0 : |proj(x-g)-x| = 7.2049e+00 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 7.2049e+00 <= 1e3*eps = 1.0000e-02 1 : maxIter = 15 <= iter = 15 ------------------------- DONE! ------------------------- Total recovered volume (no vol misfit term in inversion): 0.03444034268723054 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.55e+01 0 1 2.50e-01 7.60e+03 8.46e-01 7.60e+03 2.73e+01 0 2 2.50e-01 7.50e+03 8.19e-01 7.50e+03 2.67e+01 4 3 2.50e-01 7.35e+03 7.68e-01 7.35e+03 2.63e+01 2 Skip BFGS 4 2.50e-01 7.29e+03 7.36e-01 7.29e+03 2.55e+01 4 5 2.50e-01 7.28e+03 7.28e-01 7.28e+03 2.54e+01 6 Skip BFGS 6 2.50e-01 7.26e+03 6.98e-01 7.26e+03 2.46e+01 4 Skip BFGS 7 2.50e-01 7.26e+03 6.92e-01 7.26e+03 2.45e+01 6 8 2.50e-01 7.26e+03 6.86e-01 7.26e+03 2.44e+01 6 Skip BFGS 9 2.50e-01 7.25e+03 6.58e-01 7.25e+03 2.39e+01 4 Skip BFGS ------------------------------------------------------------------ 0 : ft = 7.2537e+03 <= alp*descent = 7.2522e+03 1 : maxIterLS = 10 <= iterLS = 10 ------------------------- End Linesearch ------------------------- The linesearch got broken. Boo. Total volume (vol misfit term in inversion): 0.03566030591476015 | .. 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 40.206 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 ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_tomo_joint_with_volume.py ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_