.. 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 Click :ref:`here ` 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-224 .. 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: 1.154e-02 , recovered(with Volume term), vol: 1.924e-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 4.25e+04 0.00e+00 4.25e+04 3.82e+01 0 1 2.50e-01 1.96e+03 4.98e-01 1.96e+03 3.48e+01 0 /home/vsts/work/1/s/SimPEG/maps.py:1721: UserWarning: Maximum number of iterations reached 2 2.50e-01 2.88e+02 4.20e-01 2.88e+02 2.78e+01 0 3 2.50e-01 2.76e+02 4.14e-01 2.77e+02 2.76e+01 8 4 2.50e-01 1.34e+02 3.51e-01 1.34e+02 1.34e+01 0 5 2.50e-01 2.94e+01 3.50e-01 2.95e+01 9.85e+00 0 6 2.50e-01 1.61e+01 3.52e-01 1.61e+01 7.60e+00 0 Skip BFGS 7 2.50e-01 1.33e+01 3.53e-01 1.34e+01 6.63e+00 0 8 2.50e-01 1.20e+01 3.51e-01 1.21e+01 7.40e+00 0 9 2.50e-01 7.62e+00 3.53e-01 7.71e+00 5.80e+00 0 10 2.50e-01 6.02e+00 3.53e-01 6.10e+00 6.16e+00 0 11 2.50e-01 4.67e+00 3.53e-01 4.76e+00 4.45e+00 0 12 2.50e-01 3.93e+00 3.53e-01 4.02e+00 5.02e+00 0 13 2.50e-01 3.09e+00 3.53e-01 3.18e+00 4.15e+00 0 14 2.50e-01 2.56e+00 3.54e-01 2.65e+00 4.81e+00 0 15 2.50e-01 1.90e+00 3.53e-01 1.98e+00 3.57e+00 0 ------------------------- STOP! ------------------------- 1 : |fc-fOld| = 6.6088e-01 <= tolF*(1+|f0|) = 4.2470e+03 0 : |xc-x_last| = 1.3538e-01 <= tolX*(1+|x0|) = 1.0000e-01 0 : |proj(x-g)-x| = 3.5695e+00 <= tolG = 1.0000e-01 0 : |proj(x-g)-x| = 3.5695e+00 <= 1e3*eps = 1.0000e-02 1 : maxIter = 15 <= iter = 15 ------------------------- DONE! ------------------------- Total recovered volume (no vol misfit term in inversion): 0.011542721350848042 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 6.00e+04 0.00e+00 6.00e+04 1.48e+01 0 1 2.50e-01 3.90e+03 4.24e-01 3.90e+03 2.72e+01 0 2 2.50e-01 3.86e+03 4.09e-01 3.86e+03 2.66e+01 4 3 2.50e-01 3.76e+03 3.89e-01 3.76e+03 2.62e+01 3 Skip BFGS 4 2.50e-01 3.76e+03 3.86e-01 3.76e+03 2.61e+01 6 5 2.50e-01 3.74e+03 3.64e-01 3.74e+03 2.52e+01 3 Skip BFGS 6 2.50e-01 3.73e+03 3.26e-01 3.73e+03 2.43e+01 2 7 2.50e-01 3.72e+03 3.12e-01 3.72e+03 2.35e+01 4 ------------------------------------------------------------------ 0 : ft = 3.7181e+03 <= alp*descent = 3.7166e+03 1 : maxIterLS = 10 <= iterLS = 10 ------------------------- End Linesearch ------------------------- The linesearch got broken. Boo. Total volume (vol misfit term in inversion): 0.019238463032283746 | .. code-block:: default 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, directives, objective_function, ) class Volume(objective_function.BaseObjectiveFunction): """ A regularization on the volume integral of the model .. math:: \phi_v = \frac{1}{2}|| \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 0.5 * (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 self.mesh.cell_volumes * ( self.knownVolume - np.inner(self.mesh.cell_volumes, m) ) def deriv2(self, m, v=None): if v is not None: return 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 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.defineBlock( 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("$\\varphi$") ax.set_ylabel("$\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 true sig model and log model dobs sigtrue = sigmaMap * phitrue # modt = Model.BaseModel(M); slownesstrue = slownessMap * phitrue # true model (m = log(sigma)) # set up the problem and survey survey = tomo.Survey(source_list) problem = tomo.Simulation(M, survey=survey, slownessMap=slownessMap) if plotIt: fig, ax = plt.subplots(1, 1) cb = plt.colorbar(M.plot_image(phitrue, ax=ax)[0], ax=ax) survey.plot(ax=ax) cb.set_label("$\\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.323 seconds) **Estimated memory usage:** 18 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-python :download:`Download Python source code: plot_tomo_joint_with_volume.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_tomo_joint_with_volume.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_