.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/comparisons/1D_VTI_empymod.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_gallery_comparisons_1D_VTI_empymod.py: 1. empymod: 1D VTI resistivity ============================== The code ``empymod`` is an open-source code which can model CSEM responses for a layered medium including VTI electrical anisotropy, see `emsig.xyz `_. Content: 1. Full-space VTI model for a finite length, finite strength, rotated bipole. 2. Layered model for a deep water model with a point dipole source. .. GENERATED FROM PYTHON SOURCE LINES 16-22 .. code-block:: default import emg3d import empymod import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 24-29 1. Full-space VTI model for a finite length, finite strength, rotated bipole ---------------------------------------------------------------------------- In order to shorten the build-time of the gallery we use a coarse model. Set ``coarse_model = False`` to obtain a result of higher accuracy. .. GENERATED FROM PYTHON SOURCE LINES 29-32 .. code-block:: default coarse_model = True .. GENERATED FROM PYTHON SOURCE LINES 33-35 Survey and model parameters ``````````````````````````` .. GENERATED FROM PYTHON SOURCE LINES 35-61 .. code-block:: default # Receiver coordinates if coarse_model: x = (np.arange(256))*20-2550 else: x = (np.arange(1025))*5-2560 rx = np.repeat([x, ], np.size(x), axis=0) ry = rx.transpose() frx, fry = rx.ravel(), ry.ravel() rz = -400.0 azimuth = 33 elevation = 18 # Source coordinates, frequency, and strength source = emg3d.TxElectricDipole( coordinates=[-50, 50, -30, 30, -320., -280.], # [x1, x2, y1, y2, z1, z2] strength=np.pi, # A ) frequency = 1.1 # Hz # Model parameters h_res = 1. # Horizontal resistivity aniso = np.sqrt(2.) # Anisotropy v_res = h_res*aniso**2 # Vertical resistivity .. GENERATED FROM PYTHON SOURCE LINES 62-69 1.a Regular VTI case ```````````````````` empymod ``````` Note: The coordinate system of empymod is positive z down, for emg3d it is positive z up. We have to switch therefore src_z, rec_z, and elevation. .. GENERATED FROM PYTHON SOURCE LINES 69-87 .. code-block:: default # Collect common input for empymod. inp = { 'src': np.r_[source.coordinates[:4], -source.coordinates[4:]], 'depth': [], 'res': h_res, 'aniso': aniso, 'strength': source.strength, 'srcpts': 5, 'freqtime': frequency, 'htarg': {'pts_per_dec': -1}, } # Compute epm = empymod.bipole( rec=[frx, fry, -rz, azimuth, -elevation], verb=3, **inp ).reshape(np.shape(rx)) .. rst-class:: sphx-glr-script-out .. code-block:: none :: empymod START :: v2.2.0 depth [m] : res [Ohm.m] : 1 aniso [-] : 1.41421 epermH [-] : 1 epermV [-] : 1 mpermH [-] : 1 mpermV [-] : 1 > MODEL IS A FULLSPACE direct field : Comp. in wavenumber domain frequency [Hz] : 1.1 Hankel : DLF (Fast Hankel Transform) > Filter : Key 201 (2009) > DLF type : Lagged Convolution Loop over : Frequencies Source(s) : 1 bipole(s) > intpts : 5 > length [m] : 123.288 > strength[A] : 3.14159 > x_c [m] : 0 > y_c [m] : 0 > z_c [m] : 300 > azimuth [°] : 30.9638 > dip [°] : -18.9318 Receiver(s) : 65536 dipole(s) > x [m] : -2550 - 2550 : 65536 [min-max; #] > y [m] : -2550 - 2550 : 65536 [min-max; #] > z [m] : 400 > azimuth [°] : 33 > dip [°] : -18 Required ab's : 11 12 13 21 22 23 31 32 33 :: empymod END; runtime = 0:00:00.679417 :: 45 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 88-90 emg3d ````` .. GENERATED FROM PYTHON SOURCE LINES 90-112 .. code-block:: default if coarse_model: min_width_limits = 40 stretching = [1.045, 1.045] else: min_width_limits = 20 stretching = [1.03, 1.045] # Create stretched grid grid = emg3d.construct_mesh( frequency=frequency, properties=h_res, center=source.center, domain=([-2500, 2500], [-2500, 2500], [-2900, 2100]), min_width_limits=min_width_limits, stretching=stretching, lambda_from_center=True, lambda_factor=0.8, center_on_edge=False, ) grid .. raw:: html
TensorMesh 512,000 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 80 -3,668.35 3,851.17 40.00 182.82 1.04
y 80 -3,668.35 3,851.17 40.00 182.82 1.04
z 80 -4,151.17 3,368.35 40.00 182.82 1.04


.. GENERATED FROM PYTHON SOURCE LINES 113-122 .. code-block:: default # Define the model model = emg3d.Model( grid, property_x=h_res, property_z=v_res, mapping='Resistivity') # Compute the electric field efield = emg3d.solve_source(model, source, frequency, verb=4, plain=True) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:39:38 :: v1.8.0 MG-cycle : 'F' sslsolver : False semicoarsening : False [0] tol : 1e-06 linerelaxation : False [0] maxit : 50 nu_{i,1,c,2} : 0, 2, 1, 2 verb : 4 Original grid : 80 x 80 x 80 => 512,000 cells Coarsest grid : 5 x 5 x 5 => 125 cells Coarsest level : 4 ; 4 ; 4 [hh:mm:ss] rel. error [abs. error, last/prev] l s h_ 2h_ \ / 4h_ \ /\ / 8h_ \ /\ / \ / 16h_ \/\/ \/ \/ [21:39:41] 3.578e-02 after 1 F-cycles [4.601e-05, 0.036] 0 0 [21:39:43] 3.734e-03 after 2 F-cycles [4.801e-06, 0.104] 0 0 [21:39:46] 4.985e-04 after 3 F-cycles [6.411e-07, 0.134] 0 0 [21:39:48] 7.463e-05 after 4 F-cycles [9.597e-08, 0.150] 0 0 [21:39:51] 1.198e-05 after 5 F-cycles [1.541e-08, 0.161] 0 0 [21:39:53] 2.044e-06 after 6 F-cycles [2.629e-09, 0.171] 0 0 [21:39:56] 3.953e-07 after 7 F-cycles [5.083e-10, 0.193] 0 0 > CONVERGED > MG cycles : 7 > Final rel. error : 3.953e-07 :: emg3d END :: 21:39:56 :: runtime = 0:00:17 .. GENERATED FROM PYTHON SOURCE LINES 123-125 Plot function ````````````` .. GENERATED FROM PYTHON SOURCE LINES 125-187 .. code-block:: default def plot(epm, e3d, title, vmin, vmax): # Start figure. a_kwargs = {'cmap': "viridis", 'vmin': vmin, 'vmax': vmax, 'shading': 'nearest'} e_kwargs = {'cmap': plt.cm.get_cmap("RdBu_r", 8), 'vmin': -2, 'vmax': 2, 'shading': 'nearest'} fig, axs = plt.subplots(2, 3, figsize=(10, 5.5), sharex=True, sharey=True, subplot_kw={'box_aspect': 1}) ((ax1, ax2, ax3), (ax4, ax5, ax6)) = axs x3 = x/1000 # km # Plot Re(data) ax1.set_title(r"(a) |Re(empymod)|") cf0 = ax1.pcolormesh(x3, x3, np.log10(epm.real.amp()), **a_kwargs) ax2.set_title(r"(b) |Re(emg3d)|") ax2.pcolormesh(x3, x3, np.log10(e3d.real.amp()), **a_kwargs) ax3.set_title(r"(c) Error real part") rel_error = 100*np.abs((epm.real - e3d.real) / epm.real) cf2 = ax3.pcolormesh(x3, x3, np.log10(rel_error), **e_kwargs) # Plot Im(data) ax4.set_title(r"(d) |Im(empymod)|") ax4.pcolormesh(x3, x3, np.log10(epm.imag.amp()), **a_kwargs) ax5.set_title(r"(e) |Im(emg3d)|") ax5.pcolormesh(x3, x3, np.log10(e3d.imag.amp()), **a_kwargs) ax6.set_title(r"(f) Error imaginary part") rel_error = 100*np.abs((epm.imag - e3d.imag) / epm.imag) ax6.pcolormesh(x3, x3, np.log10(rel_error), **e_kwargs) # Colorbars unit = "(V/m)" if "E" in title else "(A/m)" fig.colorbar(cf0, ax=axs[0, :], label=r"$\log_{10}$ Amplitude "+unit) cbar = fig.colorbar(cf2, ax=axs[1, :], label=r"Relative Error") cbar.set_ticks([-2, -1, 0, 1, 2]) cbar.ax.set_yticklabels([r"$0.01\,\%$", r"$0.1\,\%$", r"$1\,\%$", r"$10\,\%$", r"$100\,\%$"]) ax1.set_xlim(min(x3), max(x3)) ax1.set_ylim(min(x3), max(x3)) # Axis label fig.text(0.4, 0.05, "Inline Offset (km)", fontsize=14) fig.text(0.05, 0.3, "Crossline Offset (km)", rotation=90, fontsize=14) fig.suptitle(title, y=1, fontsize=20) print(f"- Source: {source}") print(f"- Frequency: {frequency} Hz") rtype = "Electric" if "E" in title else "Magnetic" print(f"- {rtype} receivers: z={rz} m; θ={azimuth}°, φ={elevation}°") fig.show() .. GENERATED FROM PYTHON SOURCE LINES 188-190 Plot ```` .. GENERATED FROM PYTHON SOURCE LINES 190-195 .. code-block:: default e3d = efield.get_receiver((rx, ry, rz, azimuth, elevation)) plot(epm, e3d, r'Diffusive Fullspace $E$', vmin=-12, vmax=-6) .. image-sg:: /gallery/comparisons/images/sphx_glr_1D_VTI_empymod_001.png :alt: Diffusive Fullspace $E$, (a) |Re(empymod)|, (b) |Re(emg3d)|, (c) Error real part, (d) |Im(empymod)|, (e) |Im(emg3d)|, (f) Error imaginary part :srcset: /gallery/comparisons/images/sphx_glr_1D_VTI_empymod_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxElectricDipole: 3.1 A; e1={-50.0; -30.0; -320.0} m; e2={50.0; 30.0; -280.0} m - Frequency: 1.1 Hz - Electric receivers: z=-400.0 m; θ=33°, φ=18° .. GENERATED FROM PYTHON SOURCE LINES 196-201 2. Layered model for a deep water model with a point dipole source ------------------------------------------------------------------ Survey and model parameters ``````````````````````````` .. GENERATED FROM PYTHON SOURCE LINES 201-225 .. code-block:: default # Receiver coordinates if coarse_model: x = (np.arange(256))*20-2550 else: x = (np.arange(1025))*5-2560 rx = np.repeat([x, ], np.size(x), axis=0) ry = rx.transpose() frx, fry = rx.ravel(), ry.ravel() rz = -950.0 azimuth = 30 elevation = 5 # Source coordinates and frequency source = emg3d.TxElectricDipole(coordinates=[0, 0, -900, 0, 0]) frequency = 1.0 # Hz # Model parameters h_res = [1, 50, 1, 0.3, 1e12] # Horizontal resistivity aniso = np.sqrt([2, 2, 2, 1, 1]) # Anisotropy v_res = h_res*aniso**2 # Vertical resistivity depth = np.array([-2200, -2000, -1000, 0]) # Layer boundaries .. GENERATED FROM PYTHON SOURCE LINES 226-228 empymod ``````` .. GENERATED FROM PYTHON SOURCE LINES 228-241 .. code-block:: default epm_d = empymod.bipole( src=source.coordinates, depth=depth, res=h_res, aniso=aniso, freqtime=frequency, htarg={'pts_per_dec': -1}, rec=[frx, fry, rz, azimuth, elevation], verb=3, ).reshape(np.shape(rx)) .. rst-class:: sphx-glr-script-out .. code-block:: none :: empymod START :: v2.2.0 depth [m] : -2200 -2000 -1000 0 res [Ohm.m] : 1 50 1 0.3 1E+12 aniso [-] : 1.41421 1.41421 1.41421 1 1 epermH [-] : 1 1 1 1 1 epermV [-] : 1 1 1 1 1 mpermH [-] : 1 1 1 1 1 mpermV [-] : 1 1 1 1 1 direct field : Comp. in wavenumber domain frequency [Hz] : 1 Hankel : DLF (Fast Hankel Transform) > Filter : Key 201 (2009) > DLF type : Lagged Convolution Loop over : Frequencies Source(s) : 1 dipole(s) > x [m] : 0 > y [m] : 0 > z [m] : -900 > azimuth [°] : 0 > dip [°] : 0 Receiver(s) : 65536 dipole(s) > x [m] : -2550 - 2550 : 65536 [min-max; #] > y [m] : -2550 - 2550 : 65536 [min-max; #] > z [m] : -950 > azimuth [°] : 30 > dip [°] : 5 Required ab's : 11 21 31 :: empymod END; runtime = 0:00:00.074129 :: 3 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 242-244 emg3d ````` .. GENERATED FROM PYTHON SOURCE LINES 244-263 .. code-block:: default if coarse_model: min_width_limits = 40 else: min_width_limits = 20 # Create stretched grid grid = emg3d.construct_mesh( frequency=frequency, properties=[h_res[3], h_res[0]], center=source.center, domain=([-2500, 2500], [-2500, 2500], None), vector=(None, None, -2200 + np.arange(111)*20), min_width_limits=min_width_limits, stretching=[1.1, 1.5], center_on_edge=False, ) grid .. raw:: html
TensorMesh 655,360 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 64 -5,697.64 7,410.84 40.00 1,713.20 1.46
y 64 -5,697.64 7,410.84 40.00 1,713.20 1.46
z 160 -5,508.08 3,308.08 20.00 390.86 1.13


.. GENERATED FROM PYTHON SOURCE LINES 264-284 .. code-block:: default # Create the model: horizontal resistivity res_x_full = h_res[0]*np.ones(grid.n_cells) # Background res_x_full[grid.cell_centers[:, 2] >= depth[0]] = h_res[1] # Target res_x_full[grid.cell_centers[:, 2] >= depth[1]] = h_res[2] # Overburden res_x_full[grid.cell_centers[:, 2] >= depth[2]] = h_res[3] # Water res_x_full[grid.cell_centers[:, 2] >= depth[3]] = h_res[4] # Air # Create the model: vertical resistivity res_z_full = v_res[0]*np.ones(grid.n_cells) # Background res_z_full[grid.cell_centers[:, 2] >= depth[0]] = v_res[1] res_z_full[grid.cell_centers[:, 2] >= depth[1]] = v_res[2] res_z_full[grid.cell_centers[:, 2] >= depth[2]] = v_res[3] res_z_full[grid.cell_centers[:, 2] >= depth[3]] = v_res[4] # Get the model model = emg3d.Model( grid, property_x=res_x_full, property_z=res_z_full, mapping='Resistivity') .. GENERATED FROM PYTHON SOURCE LINES 285-290 .. code-block:: default # Compute the electric field efield = emg3d.solve_source(model, source, frequency, verb=4) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:39:59 :: v1.8.0 MG-cycle : 'F' sslsolver : 'bicgstab' semicoarsening : True [1 2 3] tol : 1e-06 linerelaxation : True [4 5 6] maxit : 50 (3) nu_{i,1,c,2} : 0, 2, 1, 2 verb : 4 Original grid : 64 x 64 x 160 => 655,360 cells Coarsest grid : 2 x 2 x 5 => 20 cells Coarsest level : 5 ; 5 ; 5 [hh:mm:ss] rel. error solver MG l s h_ 2h_ \ / 4h_ \ /\ / 8h_ \ /\ / \ / 16h_ \ /\ / \ / \ / 32h_ \/\/ \/ \/ \/ [21:40:11] 2.609e-03 after 1 F-cycles 4 1 [21:40:22] 6.118e-05 after 2 F-cycles 5 2 [21:40:32] 1.494e-04 after 3 F-cycles 6 3 [21:40:44] 5.038e-06 after 4 F-cycles 4 1 [21:40:56] 1.933e-06 after 5 F-cycles 5 2 [21:41:06] 1.344e-05 after 6 F-cycles 6 3 [21:41:06] 3.922e-06 after 1 bicgstab-cycles [21:41:17] 1.480e-07 after 7 F-cycles 4 1 [21:41:29] 6.793e-07 after 8 F-cycles 5 2 [21:41:30] 6.771e-07 after 2 bicgstab-cycles > CONVERGED > Solver steps : 2 > MG prec. steps : 8 > Final rel. error : 6.771e-07 :: emg3d END :: 21:41:30 :: runtime = 0:01:31 .. GENERATED FROM PYTHON SOURCE LINES 291-293 Plot ```` .. GENERATED FROM PYTHON SOURCE LINES 293-297 .. code-block:: default e3d_d = efield.get_receiver((rx, ry, rz, azimuth, elevation)) plot(epm_d, e3d_d, r'Deep water point dipole $E$', vmin=-14, vmax=-8) .. image-sg:: /gallery/comparisons/images/sphx_glr_1D_VTI_empymod_002.png :alt: Deep water point dipole $E$, (a) |Re(empymod)|, (b) |Re(emg3d)|, (c) Error real part, (d) |Im(empymod)|, (e) |Im(emg3d)|, (f) Error imaginary part :srcset: /gallery/comparisons/images/sphx_glr_1D_VTI_empymod_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxElectricDipole: 1.0 A; center={0.0; 0.0; -900.0} m; θ=0.0°, φ=0.0°; l=1.0 m - Frequency: 1.0 Hz - Electric receivers: z=-950.0 m; θ=30°, φ=5° .. GENERATED FROM PYTHON SOURCE LINES 298-300 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 21:41:31 2022 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
File system ext4
Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) [GCC 10.3.0]
numpy 1.22.4 scipy 1.9.0 numba 0.55.2
emg3d 1.8.0 empymod 2.2.0 xarray 2022.6.0
discretize 0.8.2 h5py 3.7.0 matplotlib 3.4.3
tqdm 4.64.0 IPython 8.4.0
Intel(R) oneAPI Math Kernel Library Version 2022.0-Product Build 20211112 for Intel(R) 64 architecture applications


.. rst-class:: sphx-glr-timing **Total running time of the script:** ( 1 minutes 54.798 seconds) **Estimated memory usage:** 456 MB .. _sphx_glr_download_gallery_comparisons_1D_VTI_empymod.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 1D_VTI_empymod.py <1D_VTI_empymod.py>` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 1D_VTI_empymod.ipynb <1D_VTI_empymod.ipynb>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_