.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "gallery/comparisons/3D_triaxial_SimPEG.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        Click :ref:`here <sphx_glr_download_gallery_comparisons_3D_triaxial_SimPEG.py>`
        to download the full example code

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_gallery_comparisons_3D_triaxial_SimPEG.py:


3. SimPEG: 3D with tri-axial anisotropy
=======================================

`SimPEG <https://simpeg.xyz>`_ is an open source python package for simulation
and gradient based parameter estimation in geophysical applications. Here we
compare ``emg3d`` with ``SimPEG`` using the forward solver ``Pardiso``.

.. GENERATED FROM PYTHON SOURCE LINES 9-21

.. code-block:: default

    import os
    import pooch
    import emg3d
    import numpy as np
    import matplotlib.pyplot as plt
    plt.style.use('bmh')


    # Adjust this path to a folder of your choice.
    data_path = os.path.join('..', 'download', '')









.. GENERATED FROM PYTHON SOURCE LINES 23-25

Model parameters
----------------

.. GENERATED FROM PYTHON SOURCE LINES 25-41

.. code-block:: default


    # Depths (0 is sea-surface);
    # hence a deep sea case where we can ignore the air.
    water_depth = 3000
    target_x = np.r_[-500, 500]
    target_y = target_x
    target_z = -water_depth + np.r_[-400, -100]

    # Resistivities
    res_sea = 0.33
    res_back = [1., 2., 3.]  # Background in x-, y-, and z-directions
    res_target = 100.

    # Acquisition frequency
    frequency = 1.0








.. GENERATED FROM PYTHON SOURCE LINES 42-44

Grid
----

.. GENERATED FROM PYTHON SOURCE LINES 44-56

.. code-block:: default


    vx, vz = np.arange(-20, 21)*100, np.arange(-34, -19)*100
    grid = emg3d.meshes.construct_mesh(
        frequency=frequency,
        properties=1,
        center=(0, 0, 0),
        vector=(vx, vx, vz),
        center_on_edge=True,
    )

    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'>131,072 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;'>64</td>
        <td style='padding: 5px 20px 5px 20px;'>-5,253.29</td>
        <td style='padding: 5px 20px 5px 20px;'>5,253.29</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>515.62</td>
        <td style='padding: 5px 20px 5px 20px;'>1.15</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>y</td>
        <td style='padding: 5px 20px 5px 20px;'>64</td>
        <td style='padding: 5px 20px 5px 20px;'>-5,253.29</td>
        <td style='padding: 5px 20px 5px 20px;'>5,253.29</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>515.62</td>
        <td style='padding: 5px 20px 5px 20px;'>1.15</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>z</td>
        <td style='padding: 5px 20px 5px 20px;'>32</td>
        <td style='padding: 5px 20px 5px 20px;'>-6,583.88</td>
        <td style='padding: 5px 20px 5px 20px;'>1,183.88</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>731.62</td>
        <td style='padding: 5px 20px 5px 20px;'>1.25</td>
      </tr>
    </table>

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

.. GENERATED FROM PYTHON SOURCE LINES 57-59

Survey parameters
-----------------

.. GENERATED FROM PYTHON SOURCE LINES 59-68

.. code-block:: default


    # We take the receiver locations at the actual CCx-locations
    rec_x = grid.cell_centers_x[12:-12]
    rec = (rec_x, 0, -water_depth, 0, 0)
    print(f"Receiver locations:\n{rec_x}\n")

    source = emg3d.TxElectricDipole([-100, 100, 0, 0, -2900, -2900])
    sfield = emg3d.get_source_field(grid, source, frequency)  # Source field





.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Receiver locations:
    [-1950. -1850. -1750. -1650. -1550. -1450. -1350. -1250. -1150. -1050.
      -950.  -850.  -750.  -650.  -550.  -450.  -350.  -250.  -150.   -50.
        50.   150.   250.   350.   450.   550.   650.   750.   850.   950.
      1050.  1150.  1250.  1350.  1450.  1550.  1650.  1750.  1850.  1950.]





.. GENERATED FROM PYTHON SOURCE LINES 69-71

Create model
------------

.. GENERATED FROM PYTHON SOURCE LINES 71-113

.. code-block:: default


    # Layered_background
    res_x = res_sea*np.ones(grid.n_cells)
    res_y = res_x.copy()
    res_z = res_x.copy()

    # Tri-axial background.
    res_x[grid.cell_centers[:, 2] <= -water_depth] = res_back[0]
    res_y[grid.cell_centers[:, 2] <= -water_depth] = res_back[1]
    res_z[grid.cell_centers[:, 2] <= -water_depth] = res_back[2]

    res_x_bg = res_x.copy()
    res_y_bg = res_y.copy()
    res_z_bg = res_z.copy()

    # Include the target
    target_inds = (
        (grid.cell_centers[:, 0] >= target_x[0]) &
        (grid.cell_centers[:, 0] <= target_x[1]) &
        (grid.cell_centers[:, 1] >= target_y[0]) &
        (grid.cell_centers[:, 1] <= target_y[1]) &
        (grid.cell_centers[:, 2] >= target_z[0]) &
        (grid.cell_centers[:, 2] <= target_z[1])
    )
    res_x[target_inds] = res_target
    res_y[target_inds] = res_target
    res_z[target_inds] = res_target

    # Create emg3d-models for given frequency
    model = emg3d.Model(
            grid, property_x=res_x, property_y=res_y,
            property_z=res_z, mapping='Resistivity')
    model_bg = emg3d.Model(
            grid, property_x=res_x_bg, property_y=res_y_bg,
            property_z=res_z_bg, mapping='Resistivity')

    # Plot a slice
    grid.plot_3d_slicer(
            model.property_x, zslice=-3200, clim=[0, 2],
            xlim=(-4000, 4000), ylim=(-4000, 4000), zlim=(-4000, -2000)
    )




.. image-sg:: /gallery/comparisons/images/sphx_glr_3D_triaxial_SimPEG_001.png
   :alt: 3D triaxial SimPEG
   :srcset: /gallery/comparisons/images/sphx_glr_3D_triaxial_SimPEG_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 114-116

Compute ``emg3d``
-----------------

.. GENERATED FROM PYTHON SOURCE LINES 116-124

.. code-block:: default


    e3d_ftg = emg3d.solve(model, sfield, verb=1)
    e3d_tg = e3d_ftg.get_receiver(rec)

    e3d_fbg = emg3d.solve(model_bg, sfield, verb=1)
    e3d_bg = e3d_fbg.get_receiver(rec)






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    :: emg3d :: 7.2e-07; 1(4); 0:00:08; CONVERGED
    :: emg3d :: 6.6e-07; 1(4); 0:00:08; CONVERGED




.. GENERATED FROM PYTHON SOURCE LINES 125-127

Fetch and load SimPEG result
----------------------------

.. GENERATED FROM PYTHON SOURCE LINES 127-142

.. code-block:: default


    # Fetch pre-computed data.
    fname = 'simpeg.h5'
    pooch.retrieve(
        'https://raw.github.com/emsig/data/2021-05-21/emg3d/external/'+fname,
        'e0502ccfb6dfec599f4c53d9b8f8a0c79b7d872c7224a9b403cb57f39e729409',
        fname=fname,
        path=data_path,
    )

    # Load pre-computed data.
    spg = emg3d.load(data_path + fname)
    spg_tg, spg_bg = spg['spg_tg'], spg['spg_bg']






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    Data loaded from «/home/dtr/3-GitHub/emsig/emg3d-gallery/examples/download/simpeg.h5»
    [emg3d v0.17.1.dev22+ge23c468 (format 1.0) on 2021-04-15T16:28:20.499498].




.. GENERATED FROM PYTHON SOURCE LINES 143-145

Plot result
-----------

.. GENERATED FROM PYTHON SOURCE LINES 145-187

.. code-block:: default

    def nrmsd(a, b):
        """Return Normalized Root-Mean-Square Difference."""
        return 200 * abs(a - b) / (abs(a) + abs(b))


    fig, axs = plt.subplots(2, 2, figsize=(9, 5), sharex=True, sharey='row')
    ((ax1, ax3), (ax2, ax4)) = axs

    # Real part
    ax1.set_title(r'|Real|')
    ax1.plot(rec_x/1e3, 1e12*np.abs(spg_tg.real), 'C0-', label='SimPEG target')
    ax1.plot(rec_x/1e3, 1e12*np.abs(spg_bg.real), 'C1-', label='SimPEG BG')
    ax1.plot(rec_x/1e3, 1e12*np.abs(e3d_tg.real), 'k:')
    ax1.plot(rec_x/1e3, 1e12*np.abs(e3d_bg.real), 'k--')
    ax1.set_ylabel('$E_x$ (pV/m)')
    ax1.set_yscale('log')
    ax1.legend()

    # Normalized difference real
    ax2.plot(rec_x/1e3, nrmsd(spg_tg.real, e3d_tg.real), 'C0.')
    ax2.plot(rec_x/1e3, nrmsd(spg_bg.real, e3d_bg.real), 'C1.')
    ax2.set_ylabel('Norm. Diff (%)')
    ax2.set_xlabel('Offset (km)')

    # Imaginary part
    ax3.set_title(r'|Imaginary|')
    ax3.plot(rec_x/1e3, 1e12*np.abs(spg_tg.imag), 'C0-')
    ax3.plot(rec_x/1e3, 1e12*np.abs(spg_bg.imag), 'C1-')
    ax3.plot(rec_x/1e3, 1e12*np.abs(e3d_tg.imag), 'k:', label='emg3d target')
    ax3.plot(rec_x/1e3, 1e12*np.abs(e3d_bg.imag), 'k--', label='emg3d BG')
    ax3.legend()

    # Normalized difference imag
    ax4.plot(rec_x/1e3, nrmsd(spg_tg.imag, e3d_tg.imag), 'C0.')
    ax4.plot(rec_x/1e3, nrmsd(spg_bg.imag, e3d_bg.imag), 'C1.')
    ax4.set_xlabel('Offset (km)')
    ax4.set_yscale('log')

    fig.tight_layout()
    fig.show()





.. image-sg:: /gallery/comparisons/images/sphx_glr_3D_triaxial_SimPEG_002.png
   :alt: |Real|, |Imaginary|
   :srcset: /gallery/comparisons/images/sphx_glr_3D_triaxial_SimPEG_002.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 188-248

Reproduce ``SimPEG`` result
---------------------------

In order to reduce (a) the number of dependencies to generate the gallery
and, more importantly, (b) the runtime and memory requirements of the gallery
the SimPEG result is pre-computed.

.. note::

    The following cell needs to be carried out to compute the SimPEG results
    from scratch. For this you  have to install ``simpeg`` and
    ``pymatsolver``. The code example and the ``simpeg.h5``-file used above
    were created on 2021-04-14 with ``simpeg=0.14.3``, ``pymatsolver=0.1.1``,
    and ``discretize=0.6.3``.


.. code-block:: python

    # Note, in order to use the ``Pardiso``-solver ``pymatsolver`` has to be
    # installed via ``conda``, not via ``pip``!
    import SimPEG
    import discretize
    import pymatsolver
    import SimPEG.electromagnetics.frequency_domain as FDEM

    # Set up the receivers
    rx_locs = discretize.utils.ndgrid([rec_x, np.r_[0], np.r_[-water_depth]])
    rx_list = [
        FDEM.receivers.PointElectricField(
            orientation='x', component="real", locations=rx_locs),
        FDEM.receivers.PointElectricField(
            orientation='x', component="imag", locations=rx_locs)
    ]

    # We use the emg3d-source-vector, to ensure we use the same in both cases
    svector = np.real(sfield.field/-sfield.smu0)
    src_sp = FDEM.sources.RawVec_e(rx_list, s_e=svector, frequency=frequency)
    src_list = [src_sp]
    survey = FDEM.Survey(src_list)

    # Define the Simulation
    mesh = discretize.TensorMesh(grid.h, grid.origin)
    sim = FDEM.simulation.Simulation3DElectricField(
            mesh,
            survey=survey,
            sigmaMap=SimPEG.maps.IdentityMap(mesh),
            solver=pymatsolver.Pardiso,
    )

    spg_tg_dobs = sim.dpred(np.vstack([1./res_x, 1./res_y, 1./res_z]).T)
    spg_ftg = SimPEG.survey.Data(survey, dobs=spg_tg_dobs)

    spg_bg_dobs = sim.dpred(
            np.vstack([1./res_x_bg, 1./res_y_bg, 1./res_z_bg]).T)
    spg_fbg = SimPEG.survey.Data(survey, dobs=spg_bg_dobs)

    spg_tg = spg_ftg[src_sp, rx_list[0]] + 1j*spg_ftg[src_sp, rx_list[1]]
    spg_bg = spg_fbg[src_sp, rx_list[0]] + 1j*spg_fbg[src_sp, rx_list[1]]

    # emg3d.save('simpeg.h5', spg_tg=spg_tg, spg_bg=spg_bg)

.. GENERATED FROM PYTHON SOURCE LINES 251-252

.. 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:42:43 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  19.683 seconds)

**Estimated memory usage:**  101 MB


.. _sphx_glr_download_gallery_comparisons_3D_triaxial_SimPEG.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example


    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: 3D_triaxial_SimPEG.py <3D_triaxial_SimPEG.py>`

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

      :download:`Download Jupyter notebook: 3D_triaxial_SimPEG.ipynb <3D_triaxial_SimPEG.ipynb>`


.. only:: html

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

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