.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/tutorials/total_vs_ps_field.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_tutorials_total_vs_ps_field.py: 05. Total vs primary/secondary field ==================================== We usually use emg3d for total field computations. However, we can also use it in a primary-secondary field formulation, where we compute the primary field with a (semi-)analytical solution. In this example we use emg3d to compute - Total field - Primary field - Secondary field and compare the total field to the primary+secondary field. The primary-field computation could be replaced by a 1D modeller such as empymod. You can play around with the required computation-domain: Using a primary-secondary formulation makes it possible to restrict the required computation domain for the scatterer a lot, therefore speeding up the computation. However, we do not dive into that in this example. Background ---------- The total field is given by .. math:: :label: totalfield s \mu \sigma \mathbf{\hat{E}} + \nabla \times \nabla \times \mathbf{\hat{E}} = -s\mu\mathbf{\hat{J}}_s . We can split it up into a primary field :math:\mathbf{\hat{E}}^p and a secondary field :math:\mathbf{\hat{E}}^s, .. math:: :label: fieldsplit \mathbf{\hat{E}} = \mathbf{\hat{E}}^p + \mathbf{\hat{E}}^s, where we also have to split our conductivity model into .. math:: :label: modelsplit \sigma = \sigma^p + \Delta\sigma. The primary field could be just the direct field, or the direct field plus the air layer, or an entire 1D background, something that can be computed (semi-)analytically. The secondary field is everything that is not included in the primary field. The primary field is then given by .. math:: :label: primaryfield s \mu \sigma^p \mathbf{\hat{E}}^p + \nabla \times \nabla \times \mathbf{\hat{E}}^p = -s\mu\mathbf{\hat{J}}_s , and the secondary field can be computed using the primary field as source, .. math:: :label: secondaryfield s \mu \sigma \mathbf{\hat{E}}^s + \nabla \times \nabla \times \mathbf{\hat{E}}^s = -s\mu\Delta\sigma\mathbf{\hat{E}}^p . .. GENERATED FROM PYTHON SOURCE LINES 74-81 .. code-block:: default import emg3d import numpy as np import matplotlib.pyplot as plt from matplotlib.colors import LogNorm plt.style.use('bmh') .. GENERATED FROM PYTHON SOURCE LINES 83-85 Survey ------ .. GENERATED FROM PYTHON SOURCE LINES 85-94 .. code-block:: default src = [0, 0, -950, 0, 0] # x-dir. source at the origin, 50 m above seafloor off = np.arange(5, 81)*100 # Offsets rec = [off, off*0, -1000] # In-line receivers on the seafloor res = [1e10, 0.3, 1] # 1D resistivities (Ohm.m): [air, water, backgr.] freq = 1.0 # Frequency (Hz) source = emg3d.TxElectricDipole(src) .. GENERATED FROM PYTHON SOURCE LINES 95-104 Mesh ---- We create quite a coarse grid (100 m minimum cell width), to have reasonable fast computation times. Also note that the mesh here includes a large boundary because of the air layer. If you use a semi-analytical solution for the 1D background you could restrict that domain a lot. .. GENERATED FROM PYTHON SOURCE LINES 104-116 .. code-block:: default grid = emg3d.construct_mesh( frequency=freq, properties=[res[1], 100, res[2], 100], center=(src[0], src[1], -1000), min_width_limits=(100, 100, 100), domain=([-100, 8100], [-500, 500], [-2500, 0]), center_on_edge=False, ) grid .. raw:: html
 MESH EXTENT CELL WIDTH FACTOR dir nC TensorMesh 245,760 cells x 128 -33,640.98 48,527.04 100.00 6,886.06 1.20 y 40 -32,849.73 45,222.36 100.00 12,372.63 1.38 z 48 -6,940.24 32,349.73 100.00 8,973.55 1.38

.. GENERATED FROM PYTHON SOURCE LINES 117-119 Create model ------------ .. GENERATED FROM PYTHON SOURCE LINES 119-145 .. code-block:: default # Layered_background res_x = np.ones(grid.n_cells)*res[0] # Air resistivity res_x[grid.cell_centers[:, 2] < 0] = res[1] # Water resistivity res_x[grid.cell_centers[:, 2] < -1000] = res[2] # Background resistivity # Background model model_pf = emg3d.Model(grid, property_x=res_x.copy(), mapping='Resistivity') # Include the target xx = (grid.cell_centers[:, 0] >= 0) & (grid.cell_centers[:, 0] <= 6000) yy = abs(grid.cell_centers[:, 1]) <= 500 zz = (grid.cell_centers[:, 2] > -2500)*(grid.cell_centers[:, 2] < -2000) res_x[xx*yy*zz] = 100. # Target resistivity # Create target model model = emg3d.Model(grid, property_x=res_x, mapping='Resistivity') # Plot a slice grid.plot_3d_slicer( model.property_x, zslice=-2250, xlim=(-1000, 8000), ylim=(-4000, 4000), zlim=(-3000, 500), pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=200)} ) .. image-sg:: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_001.png :alt: total vs ps field :srcset: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 146-148 Compute total field with emg3d ---------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 148-152 .. code-block:: default em3_tf = emg3d.solve_source(model, source, freq, verb=1) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d :: 7.2e-07; 1(4); 0:00:14; CONVERGED .. GENERATED FROM PYTHON SOURCE LINES 153-158 Compute primary field (1D background) with emg3d ---------------------------------------------------- Here we use emg3d to compute the primary field. This could be replaced by a (semi-)analytical solution. .. GENERATED FROM PYTHON SOURCE LINES 158-161 .. code-block:: default em3_pf = emg3d.solve_source(model_pf, source, freq, verb=1) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d :: 7.0e-07; 1(4); 0:00:15; CONVERGED .. GENERATED FROM PYTHON SOURCE LINES 162-167 Compute secondary field (scatterer) with emg3d -------------------------------------------------- Define the secondary source  .. GENERATED FROM PYTHON SOURCE LINES 167-188 .. code-block:: default # Get the difference of conductivity as volume-average values diff = 1/model.property_x-1/model_pf.property_x dsigma = grid.cell_volumes.reshape(grid.shape_cells, order='F')*diff # Here we use the primary field computed with emg3d. This could be done # with a 1D modeller such as empymod instead. sfield_sf = em3_pf.copy() # Average delta sigma to the corresponding edges sfield_sf.fx[:, 1:-1, 1:-1] *= 0.25*(dsigma[:, :-1, :-1] + dsigma[:, 1:, :-1] + dsigma[:, :-1, 1:] + dsigma[:, 1:, 1:]) sfield_sf.fy[1:-1, :, 1:-1] *= 0.25*(dsigma[:-1, :, :-1] + dsigma[1:, :, :-1] + dsigma[:-1, :, 1:] + dsigma[1:, :, 1:]) sfield_sf.fz[1:-1, 1:-1, :] *= 0.25*(dsigma[:-1, :-1, :] + dsigma[1:, :-1, :] + dsigma[:-1, 1:, :] + dsigma[1:, 1:, :]) # Create field instance -iwu dsigma E sfield_sf = emg3d.Field( sfield_sf.grid, -sfield_sf.field*sfield_sf.smu0, frequency=freq) .. GENERATED FROM PYTHON SOURCE LINES 189-198 Plot the secondary source  Our secondary source is the entire target, the scatterer. Here we look at the :math:E_x secondary source field. But note that the secondary source has all three components :math:E_x, :math:E_y, and :math:E_z, even though our primary source was purely :math:x-directed. (Change fx to fy or fz in the command below, and simultaneously Ex to Ey or Ez, to show the other source fields.) .. GENERATED FROM PYTHON SOURCE LINES 198-205 .. code-block:: default grid.plot_3d_slicer( sfield_sf.fx.ravel('F'), view='abs', v_type='Ex', zslice=-2250, xlim=(-1000, 8000), ylim=(-4000, 4000), zlim=(-3000, 500), pcolor_opts={'norm': LogNorm(vmin=1e-17, vmax=1e-9)} ) .. image-sg:: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_002.png :alt: total vs ps field :srcset: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 206-208 Compute the secondary source  .. GENERATED FROM PYTHON SOURCE LINES 208-211 .. code-block:: default em3_sf = emg3d.solve(model, sfield_sf, verb=1) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d :: 7.7e-07; 2(7); 0:00:25; CONVERGED .. GENERATED FROM PYTHON SOURCE LINES 212-214 Plot result ----------- .. GENERATED FROM PYTHON SOURCE LINES 214-226 .. code-block:: default # E = E^p + E^s em3_ps = em3_pf.copy() em3_ps.field += em3_sf.field # Get the responses at receiver locations rectuple = (rec[0], rec[1], rec[2], 0, 0) em3_pf_rec = em3_pf.get_receiver(rectuple) em3_tf_rec = em3_tf.get_receiver(rectuple) em3_sf_rec = em3_sf.get_receiver(rectuple) em3_ps_rec = em3_ps.get_receiver(rectuple) .. GENERATED FROM PYTHON SOURCE LINES 227-259 .. code-block:: default fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(9, 5), sharex=True) ax1.set_title('|Real part|') ax1.plot(off/1e3, abs(em3_pf_rec.real), 'k', label='Primary Field (1D Background)') ax1.plot(off/1e3, abs(em3_sf_rec.real), '.4', ls='--', label='Secondary Field (Scatterer)') ax1.plot(off/1e3, abs(em3_ps_rec.real)) ax1.plot(off[::2]/1e3, abs(em3_tf_rec[::2].real), '.') ax1.plot(off/1e3, abs(em3_ps_rec.real-em3_tf_rec.real)) ax1.set_xlabel('Offset (km)') ax1.set_ylabel('$E_x$ (V/m)') ax1.set_yscale('log') ax1.legend() ax2.set_title('|Imaginary part|') ax2.plot(off/1e3, abs(em3_pf_rec.imag), 'k') ax2.plot(off/1e3, abs(em3_sf_rec.imag), '.4', ls='--') ax2.plot(off/1e3, abs(em3_ps_rec.imag), label='P/S Field') ax2.plot(off[::2]/1e3, abs(em3_tf_rec[::2].imag), '.', label='Total Field') ax2.plot(off/1e3, abs(em3_ps_rec.imag-em3_tf_rec.imag), label=r'$\Delta$|P/S-Total|') ax2.set_xlabel('Offset (km)') ax2.set_ylabel('$E_x$ (V/m)') ax2.set_yscale('log') ax2.yaxis.tick_right() ax2.yaxis.set_label_position("right") plt.legend() fig.tight_layout() fig.show() .. image-sg:: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_003.png :alt: |Real part|, |Imaginary part| :srcset: /gallery/tutorials/images/sphx_glr_total_vs_ps_field_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 260-262 .. code-block:: default emg3d.Report() .. raw:: html
 Wed Aug 31 21:32:35 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:** ( 0 minutes 58.127 seconds) **Estimated memory usage:** 236 MB .. _sphx_glr_download_gallery_tutorials_total_vs_ps_field.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:Download Python source code: total_vs_ps_field.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: total_vs_ps_field.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _