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

.. only:: html

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

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

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

.. _sphx_glr_gallery_comparisons_2D_triaxial_MARE2DEM.py:


2. MARE2DEM: 2D with tri-axial anisotropy
=========================================

``MARE2DEM`` is an open-source, finite element 2.5D code for controlled-source
electromagnetic (CSEM) and magnetotelluric (MT) data, see
`mare2dem.bitbucket.io <https://mare2dem.bitbucket.io>`_.


.. note::

    The ``MARE2DEM`` results are pre-computed. All input files to reproduce the
    results are available on
    https://github.com/emsig/data/tree/main/emg3d/external/MARE2DEM .

.. GENERATED FROM PYTHON SOURCE LINES 17-30

.. code-block:: default

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


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









.. GENERATED FROM PYTHON SOURCE LINES 32-34

Fetch and load MARE2DEM result
------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 34-57

.. code-block:: default


    url = 'https://raw.github.com/emsig/data/2021-05-21/emg3d/external/MARE2DEM/'
    fname1 = 'triaxial.0.resp'
    pooch.retrieve(
        url + fname1,
        '29ec8e3dbfc615bcb430df5cbd89fea6302bb3867d90ae969907314013dc871b',
        fname=fname1,
        path=data_path,
    )
    mar_tg = np.loadtxt(data_path + fname1, skiprows=93, usecols=6)
    mar_tg = mar_tg[::2] + 1j*mar_tg[1::2]

    fname2 = 'triaxial-BG.0.resp'
    pooch.retrieve(
        url + fname2,
        '036f72e30b7794304c45ef73403cdd8318ca0fc5c2fdbe7d05a33731cf3f2cf6',
        fname=fname2,
        path=data_path,
    )
    mar_bg = np.loadtxt(data_path + fname2, skiprows=93, usecols=6)
    mar_bg = mar_bg[::2] + 1j*mar_bg[1::2]









.. GENERATED FROM PYTHON SOURCE LINES 58-63

emg3d
-----

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 63-70

.. code-block:: default

    coarse_model = True

    # Source location [x, y, z, azimuth, elevation]
    source = emg3d.TxElectricDipole((50, 0, -1950, 0, 0))
    rec = (np.arange(80)*100+2050, 0, -1999.9, 0, 0)
    frequency = 0.5             # Frequency (Hz)








.. GENERATED FROM PYTHON SOURCE LINES 71-90

.. code-block:: default

    if coarse_model:
        min_width = 100
        stretching = ([1.02, 1.5], [1.05, 1.5], [1, 1.5])
    else:
        min_width = 50
        stretching = [1, 1.5]

    # Create grid.
    grid = emg3d.construct_mesh(
        frequency=frequency,
        properties=[0.3, 1, 100],
        center=(0, 0, -2000),
        domain=([-100, 10100], [-1000, 1000], [-4200, 0]),
        stretching=stretching,
        min_width_limits=min_width,
        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'>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;'>80</td>
        <td style='padding: 5px 20px 5px 20px;'>-50,416.82</td>
        <td style='padding: 5px 20px 5px 20px;'>55,590.88</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>16,432.88</td>
        <td style='padding: 5px 20px 5px 20px;'>1.48</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>y</td>
        <td style='padding: 5px 20px 5px 20px;'>48</td>
        <td style='padding: 5px 20px 5px 20px;'>-46,181.63</td>
        <td style='padding: 5px 20px 5px 20px;'>46,181.63</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>13,445.46</td>
        <td style='padding: 5px 20px 5px 20px;'>1.42</td>
      </tr>
      <tr>
        <td style='padding: 5px 20px 5px 20px;'>z</td>
        <td style='padding: 5px 20px 5px 20px;'>64</td>
        <td style='padding: 5px 20px 5px 20px;'>-9,432.76</td>
        <td style='padding: 5px 20px 5px 20px;'>45,181.63</td>
        <td style='padding: 5px 20px 5px 20px;'>100.00</td>
        <td style='padding: 5px 20px 5px 20px;'>13,445.46</td>
        <td style='padding: 5px 20px 5px 20px;'>1.42</td>
      </tr>
    </table>

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

.. GENERATED FROM PYTHON SOURCE LINES 91-132

.. code-block:: default


    xx = (grid.cell_centers[:, 0] > 0)*(grid.cell_centers[:, 0] <= 6000)
    zz = (grid.cell_centers[:, 2] > -4200)*(grid.cell_centers[:, 2] < -4000)

    # Background
    res_x_full = 2*np.ones(grid.n_cells)
    res_y_full = 1*np.ones(grid.n_cells)
    res_z_full = 3*np.ones(grid.n_cells)

    # Water - isotropic
    res_x_full[grid.cell_centers[:, 2] >= -2000] = 0.3
    res_y_full[grid.cell_centers[:, 2] >= -2000] = 0.3
    res_z_full[grid.cell_centers[:, 2] >= -2000] = 0.3

    # Air - isotropic
    res_x_full[grid.cell_centers[:, 2] >= 0] = 1e10
    res_y_full[grid.cell_centers[:, 2] >= 0] = 1e10
    res_z_full[grid.cell_centers[:, 2] >= 0] = 1e10

    # Target
    res_x_full_tg = res_x_full.copy()
    res_y_full_tg = res_y_full.copy()
    res_z_full_tg = res_z_full.copy()
    res_x_full_tg[xx*zz] = 200
    res_y_full_tg[xx*zz] = 100
    res_z_full_tg[xx*zz] = 300

    # Collect models
    model_bg = emg3d.Model(
            grid, property_x=res_x_full, property_y=res_y_full,
            property_z=res_z_full, mapping='Resistivity')
    model_tg = emg3d.Model(
            grid, property_x=res_x_full_tg, property_y=res_y_full_tg,
            property_z=res_z_full_tg, mapping='Resistivity')

    # QC model
    grid.plot_3d_slicer(
            model_tg.property_x, zlim=[-6000, 500],
            pcolor_opts={'norm': LogNorm(vmin=0.3, vmax=300)})





.. image-sg:: /gallery/comparisons/images/sphx_glr_2D_triaxial_MARE2DEM_001.png
   :alt: 2D triaxial MARE2DEM
   :srcset: /gallery/comparisons/images/sphx_glr_2D_triaxial_MARE2DEM_001.png
   :class: sphx-glr-single-img





.. GENERATED FROM PYTHON SOURCE LINES 133-135

Model background
````````````````

.. GENERATED FROM PYTHON SOURCE LINES 135-140

.. code-block:: default


    efield_bg = emg3d.solve_source(model_bg, source, frequency)
    e3d_bg = efield_bg.get_receiver(rec)









.. GENERATED FROM PYTHON SOURCE LINES 141-143

Model target
````````````

.. GENERATED FROM PYTHON SOURCE LINES 143-148

.. code-block:: default


    efield_tg = emg3d.solve_source(model_tg, source, frequency)
    e3d_tg = efield_tg.get_receiver(rec)









.. GENERATED FROM PYTHON SOURCE LINES 149-151

Plot
----

.. GENERATED FROM PYTHON SOURCE LINES 151-196

.. 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[0]/1e3, 1e12*np.abs(mar_tg.real), 'C0-', label='MARE2DEM target')
    ax1.plot(rec[0]/1e3, 1e12*np.abs(mar_bg.real), 'C1-', label='MARE2DEM BG')
    ax1.plot(rec[0]/1e3, 1e12*np.abs(e3d_tg.real), 'k--')
    ax1.plot(rec[0]/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[0]/1e3, nrmsd(mar_tg.real, e3d_tg.real), 'C0.')
    ax2.plot(rec[0]/1e3, nrmsd(mar_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[0]/1e3, 1e12*np.abs(mar_tg.imag), 'C0-')
    ax3.plot(rec[0]/1e3, 1e12*np.abs(mar_bg.imag), 'C1-')
    ax3.plot(rec[0]/1e3, 1e12*np.abs(e3d_tg.imag), 'k--', label='emg3d target')
    ax3.plot(rec[0]/1e3, 1e12*np.abs(e3d_bg.imag), 'k-.', label='emg3d BG')
    ax3.legend()

    # Normalized difference imag
    ax4.plot(rec[0]/1e3, nrmsd(mar_tg.imag, e3d_tg.imag), 'C0.')
    ax4.plot(rec[0]/1e3, nrmsd(mar_bg.imag, e3d_bg.imag), 'C1.')
    ax4.set_xlabel('Offset (km)')
    ax4.set_yscale('log')
    ax4.set_ylim([8e-3, 120])
    ax4.set_yticks([0.01, 0.1, 1, 10, 100])
    ax4.set_yticklabels(('0.01', '0.1', '1', '10', '100'))

    fig.tight_layout()
    fig.show()




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





.. GENERATED FROM PYTHON SOURCE LINES 197-199

.. 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:23 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  50.938 seconds)

**Estimated memory usage:**  214 MB


.. _sphx_glr_download_gallery_comparisons_2D_triaxial_MARE2DEM.py:

.. only:: html

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


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

      :download:`Download Python source code: 2D_triaxial_MARE2DEM.py <2D_triaxial_MARE2DEM.py>`

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

      :download:`Download Jupyter notebook: 2D_triaxial_MARE2DEM.ipynb <2D_triaxial_MARE2DEM.ipynb>`


.. only:: html

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

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