.. 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 <sphx_glr_download_gallery_tutorials_total_vs_ps_field.py>`
        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

    <div class="output_subarea output_html rendered_html output_result">
    <table>
      <tr>
        <td style='font-weight: bold; font-size: 1.2em; text-align: center;' colspan='3'>TensorMesh</td>
        <td style='font-size: 1.2em; text-align: center;'colspan='4'>245,760 cells</td>
      </tr>
      <tr>
        <th></th>
        <th></th>
        <th colspan='2' style='padding: 5px 20px 5px 20px;'>MESH EXTENT</th>
        <th colspan='2' style='padding: 5px 20px 5px 20px;'>CELL WIDTH</th>
        <th style='padding: 5px 20px 5px 20px;'>FACTOR</th>
      </tr>
      <tr>
        <th style='padding: 5px 20px 5px 20px;'>dir</th>
        <th style='padding: 5px 20px 5px 20px;'>nC</th>
        <th style='padding: 5px 20px 5px 20px;'>min</th>
        <th style='padding: 5px 20px 5px 20px;'>max</th>
        <th style='padding: 5px 20px 5px 20px;'>min</th>
        <th style='padding: 5px 20px 5px 20px;'>max</th>
        <th style='padding: 5px 20px 5px 20px;'>max</th>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>x</td>
        <td style='padding: 5px 20px 5px 20px;'>128</td>
        <td style='padding: 5px 20px 5px 20px;'>-33,640.98</td>
        <td style='padding: 5px 20px 5px 20px;'>48,527.04</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>6,886.06</td>
        <td style='padding: 5px 20px 5px 20px;'>1.20</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>y</td>
        <td style='padding: 5px 20px 5px 20px;'>40</td>
        <td style='padding: 5px 20px 5px 20px;'>-32,849.73</td>
        <td style='padding: 5px 20px 5px 20px;'>45,222.36</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>12,372.63</td>
        <td style='padding: 5px 20px 5px 20px;'>1.38</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>z</td>
        <td style='padding: 5px 20px 5px 20px;'>48</td>
        <td style='padding: 5px 20px 5px 20px;'>-6,940.24</td>
        <td style='padding: 5px 20px 5px 20px;'>32,349.73</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>8,973.55</td>
        <td style='padding: 5px 20px 5px 20px;'>1.38</td>
      </tr>
    </table>

    </div>
    <br />
    <br />

.. 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

    <div class="output_subarea output_html rendered_html output_result">
    <table style='border: 1.5px solid;'>
      <tr>
         <td style='text-align: center; font-weight: bold; font-size: 1.2em; border: 1px solid;' colspan='6'>Wed Aug 31 21:32:35 2022 CEST</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>OS</td>
        <td style='text-align: left; border: 1px solid;'>Linux</td>
        <td style='text-align: right; border: 1px solid;'>CPU(s)</td>
        <td style='text-align: left; border: 1px solid;'>4</td>
        <td style='text-align: right; border: 1px solid;'>Machine</td>
        <td style='text-align: left; border: 1px solid;'>x86_64</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>Architecture</td>
        <td style='text-align: left; border: 1px solid;'>64bit</td>
        <td style='text-align: right; border: 1px solid;'>RAM</td>
        <td style='text-align: left; border: 1px solid;'>15.5 GiB</td>
        <td style='text-align: right; border: 1px solid;'>Environment</td>
        <td style='text-align: left; border: 1px solid;'>Python</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>File system</td>
        <td style='text-align: left; border: 1px solid;'>ext4</td>
      </tr>
      <tr>
         <td style='text-align: center; border: 1px solid;' colspan='6'>Python 3.9.12 | packaged by conda-forge | (main, Mar 24 2022, 23:22:55) 
    [GCC 10.3.0]</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>numpy</td>
        <td style='text-align: left; border: 1px solid;'>1.22.4</td>
        <td style='text-align: right; border: 1px solid;'>scipy</td>
        <td style='text-align: left; border: 1px solid;'>1.9.0</td>
        <td style='text-align: right; border: 1px solid;'>numba</td>
        <td style='text-align: left; border: 1px solid;'>0.55.2</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>emg3d</td>
        <td style='text-align: left; border: 1px solid;'>1.8.0</td>
        <td style='text-align: right; border: 1px solid;'>empymod</td>
        <td style='text-align: left; border: 1px solid;'>2.2.0</td>
        <td style='text-align: right; border: 1px solid;'>xarray</td>
        <td style='text-align: left; border: 1px solid;'>2022.6.0</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>discretize</td>
        <td style='text-align: left; border: 1px solid;'>0.8.2</td>
        <td style='text-align: right; border: 1px solid;'>h5py</td>
        <td style='text-align: left; border: 1px solid;'>3.7.0</td>
        <td style='text-align: right; border: 1px solid;'>matplotlib</td>
        <td style='text-align: left; border: 1px solid;'>3.4.3</td>
      </tr>
      <tr>
        <td style='text-align: right; border: 1px solid;'>tqdm</td>
        <td style='text-align: left; border: 1px solid;'>4.64.0</td>
        <td style='text-align: right; border: 1px solid;'>IPython</td>
        <td style='text-align: left; border: 1px solid;'>8.4.0</td>
        <td style= border: 1px solid;'></td>
        <td style= border: 1px solid;'></td>
      </tr>
      <tr>
         <td style='text-align: center; border: 1px solid;' colspan='6'>Intel(R) oneAPI Math Kernel Library Version 2022.0-Product Build 20211112 for Intel(R) 64 architecture applications</td>
      </tr>
    </table>
    </div>
    <br />
    <br />


.. 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 <total_vs_ps_field.py>`

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: total_vs_ps_field.ipynb <total_vs_ps_field.ipynb>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_