.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/comparisons/magnetic_source_el_loop.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_magnetic_source_el_loop.py: 7. Magnetic source using an electric loop ========================================= Computing the :math:`E` and :math:`H` fields generated by a magnetic source We know that we can get the magnetic fields from the electric fields using Faraday's law, see :ref:`sphx_glr_gallery_comparisons_magnetic_field.py`. However, what about computing the fields generated by a magnetic source? There are two ways we can achieve that: - **creating an electric loop source**, which is what we do in this example, or - using the duality principle, see :ref:`sphx_glr_gallery_comparisons_magnetic_source_duality.py`. We create a "magnetic dipole" through an electric loop perpendicular to the defined dipole in a homogeneous VTI fullspace, and compare it to the semi-analytical solution of ``empymod``. (The code ``empymod`` is an open-source code which can model CSEM responses for a layered medium including VTI electrical anisotropy, see `emsig.xyz `_.) .. GENERATED FROM PYTHON SOURCE LINES 25-31 .. code-block:: default import emg3d import empymod import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 33-38 Full-space model for a loop dipole ---------------------------------- 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 38-41 .. code-block:: default coarse_model = True .. GENERATED FROM PYTHON SOURCE LINES 42-47 Survey and model parameters ``````````````````````````` ``emg3d.TxMagneticDipole`` creates an electric square loop perpendicular to the defined dipole, where the area of the square loop corresponds to the length of the dipole. .. GENERATED FROM PYTHON SOURCE LINES 47-73 .. 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 = 25 elevation = 10 # Source coordinates, frequency, and strength source = emg3d.TxMagneticDipole( coordinates=[-0.5, 0.5, -0.3, 0.3, -300.5, -299.5], # [x1,x2,y1,y2,z1,z2] strength=np.pi, # A ) frequency = 0.77 # 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 74-78 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 78-101 .. 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, 'freqtime': frequency, 'htarg': {'pts_per_dec': -1}, } # Compute e-field epm_e = -empymod.loop( rec=[frx, fry, -rz, azimuth, -elevation], mrec=False, verb=3, **inp ).reshape(np.shape(rx)) # Compute h-field epm_h = -empymod.loop( rec=[frx, fry, -rz, azimuth, -elevation], **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] : 0.77 Hankel : DLF (Fast Hankel Transform) > Filter : Key 201 (2009) > DLF type : Lagged Convolution Loop over : Frequencies Source(s) : 1 bipole(s) > intpts : 1 (as dipole) > length [m] : 1.53623 > strength[A] : 3.14159 > x_c [m] : 0 > y_c [m] : 0 > z_c [m] : 300 > azimuth [°] : 30.9638 > dip [°] : -40.6129 Receiver(s) : 65536 dipole(s) > x [m] : -2550 - 2550 : 65536 [min-max; #] > y [m] : -2550 - 2550 : 65536 [min-max; #] > z [m] : 400 > azimuth [°] : 25 > dip [°] : -10 Required ab's : 14 15 16 24 25 26 34 35 36 :: empymod END; runtime = 0:00:00.180658 :: 8 kernel call(s) :: empymod END; runtime = 0:00:00.177053 :: 9 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 102-104 emg3d ````` .. GENERATED FROM PYTHON SOURCE LINES 104-126 .. 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 -4,155.20 4,378.72 40.00 223.51 1.04
y 80 -4,155.20 4,378.72 40.00 223.51 1.04
z 80 -4,678.72 3,855.20 40.00 223.51 1.04


.. GENERATED FROM PYTHON SOURCE LINES 127-136 .. 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:47:42 :: 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:47:45] 4.451e-03 after 1 F-cycles [4.617e-09, 0.004] 0 0 [21:47:48] 1.643e-04 after 2 F-cycles [1.704e-10, 0.037] 0 0 [21:47:51] 8.138e-06 after 3 F-cycles [8.443e-12, 0.050] 0 0 [21:47:54] 4.803e-07 after 4 F-cycles [4.983e-13, 0.059] 0 0 > CONVERGED > MG cycles : 4 > Final rel. error : 4.803e-07 :: emg3d END :: 21:47:54 :: runtime = 0:00:12 .. GENERATED FROM PYTHON SOURCE LINES 137-139 Plot function ````````````` .. GENERATED FROM PYTHON SOURCE LINES 139-201 .. 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 202-204 Compare the electric field generated from the magnetic source ------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 204-208 .. code-block:: default e3d_e = efield.get_receiver((rx, ry, rz, azimuth, elevation)) plot(epm_e, e3d_e, r'Diffusive Fullspace $E$', vmin=-17, vmax=-10) .. image-sg:: /gallery/comparisons/images/sphx_glr_magnetic_source_el_loop_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_magnetic_source_el_loop_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxMagneticDipole: 3.1 A; e1={-0.5; 0.8; -300.0} m; e2={-0.5; -0.3; -299.3} m - Frequency: 0.77 Hz - Electric receivers: z=-400.0 m; θ=25°, φ=10° .. GENERATED FROM PYTHON SOURCE LINES 209-211 Compare the magnetic field generated from the magnetic source ------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 211-218 .. code-block:: default # Get the magnetic field :math:`H` from the electric field hfield = emg3d.get_magnetic_field(model, efield) e3d_h = hfield.get_receiver((rx, ry, rz, azimuth, elevation)) plot(epm_h, e3d_h, r'Diffusive Fullspace $H$', vmin=-13, vmax=-8) .. image-sg:: /gallery/comparisons/images/sphx_glr_magnetic_source_el_loop_002.png :alt: Diffusive Fullspace $H$, (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_magnetic_source_el_loop_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxMagneticDipole: 3.1 A; e1={-0.5; 0.8; -300.0} m; e2={-0.5; -0.3; -299.3} m - Frequency: 0.77 Hz - Magnetic receivers: z=-400.0 m; θ=25°, φ=10° .. GENERATED FROM PYTHON SOURCE LINES 219-220 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 21:47:57 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 16.924 seconds) **Estimated memory usage:** 99 MB .. _sphx_glr_download_gallery_comparisons_magnetic_source_el_loop.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: magnetic_source_el_loop.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: magnetic_source_el_loop.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_