.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/comparisons/magnetic_source_duality.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_duality.py: 8. Magnetic source using duality ================================ Computing the :math:E and :math:H fields from a magnetic source using the duality principle. 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: - **using the duality principle**, which is what we do in this example, or - creating an electric loop source, see :ref:sphx_glr_gallery_comparisons_magnetic_source_el_loop.py. emg3d solves the following equation, .. math:: :label: eq-maxwell \eta \mathbf{\hat{E}} - \nabla \times \zeta^{-1} \nabla \times \mathbf{\hat{E}} = -\mathbf{\hat{J}}^e_s , where :math:\eta = \sigma - \mathrm{i}\omega\varepsilon, :math:\zeta = \mathrm{i}\omega\mu, :math:\sigma is conductivity (S/m), :math:\omega=2\pi f is the angular frequency (Hz), :math:\mu=\mu_0\mu_\mathrm{r} is magnetic permeability (H/m), :math:\varepsilon=\varepsilon_0\varepsilon_\mathrm{r} is electric permittivity (F/m), :math:\mathbf{\hat{E}} the electric field in the frequency domain (V/m), and :math:\mathbf{\hat{J}}^e_s source current. This is the electric field due to an electric source. One can obtain the magnetic field due to a magnetic field by substituting - :math:\eta \leftrightarrow -\zeta , - :math:\mathbf{\hat{E}} \leftrightarrow -\mathbf{\hat{H}} , - :math:\mathbf{\hat{J}}^e_s \leftrightarrow \mathbf{\hat{J}}^m_s , which is called the **duality principle**. Carrying out the substitution yields .. math:: :label: dualdip \zeta \mathbf{\hat{H}} - \nabla \times \eta^{-1} \nabla \times \mathbf{\hat{H}} = -\mathbf{\hat{J}}^m_s , which is for a magnetic dipole. Changing it for a loop source adds a term :math:\mathrm{i}\omega\mu to the source term, resulting in .. math:: :label: dualloop \zeta \mathbf{\hat{H}} - \nabla \times \eta^{-1} \nabla \times \mathbf{\hat{H}} = -\mathrm{i}\omega\mu\mathbf{\hat{J}}^m_s ; see Dipoles and Loops _ for more information. emg3d is not ideal for the duality principle. Magnetic permeability is implemented isotropically, and discontinuities in magnetic permeabilities can lead to first-order errors in contrary to second-order errors for discontinuities in conductivity. However, we can still abuse the code and use it with the duality principle, at least for isotropic media. The actual implemented equation in emg3d is a slightly modified version of Equation :eq:eq-maxwell, using the diffusive approximation :math:\varepsilon=0, .. math:: :label: dualdiff \mathrm{i}\omega \mu_0 \sigma \mathbf{\hat{E}} - \nabla \times \mu_r^{-1} \nabla \times \mathbf{\hat{E}} = -\mathrm{i}\omega\mu_0\mathbf{\hat{J}}_s . We therefore only need to interchange :math:\sigma with :math:\mu_\mathrm{r}^{-1} or :math:\rho with :math:\mu_\mathrm{r} to get from :eq:dualdiff to :eq:dualloop. This is what we do in this example, for an arbitrarily rotated loop in a homogeneous, isotropic fullspace. We compare the result 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 91-97 .. code-block:: default import emg3d import empymod import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 99-104 Full-space model for a rotated magnetic loop -------------------------------------------- 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 104-107 .. code-block:: default coarse_model = True .. GENERATED FROM PYTHON SOURCE LINES 108-110 Survey and model parameters  .. GENERATED FROM PYTHON SOURCE LINES 110-134 .. 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.TxElectricDipole( coordinates=[0, 0, -300, 10, 70], # [x, y, z, azimuth, elevation] strength=np.pi, # A ) frequency = 0.77 # Hz # Model parameters h_res = 2. # Horizontal resistivity .. GENERATED FROM PYTHON SOURCE LINES 135-139 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 139-162 .. code-block:: default # Collect common input for empymod. inp = { 'src': np.r_[source.coordinates[:2], -source.coordinates[2], source.coordinates[3], -source.coordinates[4]], 'depth': [], 'res': h_res, '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] : 2 aniso [-] : 1 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 dipole(s) > x [m] : 0 > y [m] : 0 > z [m] : 300 > azimuth [°] : 10 > dip [°] : -70 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.167794 :: 8 kernel call(s) :: empymod END; runtime = 0:00:00.181631 :: 9 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 163-165 emg3d  .. GENERATED FROM PYTHON SOURCE LINES 165-187 .. 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
 MESH EXTENT CELL WIDTH FACTOR dir nC TensorMesh 884,736 cells x 96 -5,353.03 5,617.20 40.00 264.17 1.05 y 96 -5,353.03 5,617.20 40.00 264.17 1.05 z 96 -5,835.51 4,995.26 40.00 253.27 1.04

.. GENERATED FROM PYTHON SOURCE LINES 188-207 Abuse the parameters to take advantage of the duality principle --------------------------------------------------------------- See text at the top. We set here :math:\rho=1 and :math:\mu_\mathrm{r} = 1/\rho to get: .. math:: :label: iweta2iwu \mathrm{i}\omega\mu_0(1-\mathrm{i}\omega\varepsilon) = \mathrm{i}\omega\mu_0+\omega^2\mu_0\varepsilon \approx \mathrm{i}\omega\mu (in the diffusive regime), and .. math:: :label: mu2sigma \mu_\mathrm{r} = 1/\rho = \sigma \, . .. GENERATED FROM PYTHON SOURCE LINES 207-216 .. code-block:: default # Define the model => Set property_x = 1 and mu_r = 1./h_res model = emg3d.Model( grid, property_x=1., mu_r=1./h_res, mapping='Resistivity') # Compute the electric field hfield = emg3d.solve_source(model, source, frequency, verb=4, plain=True) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:47:59 :: 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 : 96 x 96 x 96 => 884,736 cells Coarsest grid : 3 x 3 x 3 => 27 cells Coarsest level : 5 ; 5 ; 5 [hh:mm:ss] rel. error [abs. error, last/prev] l s h_ 2h_ \ / 4h_ \ /\ / 8h_ \ /\ / \ / 16h_ \ /\ / \ / \ / 32h_ \/\/ \/ \/ \/ [21:48:04] 3.240e-02 after 1 F-cycles [3.094e-07, 0.032] 0 0 [21:48:09] 2.986e-03 after 2 F-cycles [2.851e-08, 0.092] 0 0 [21:48:13] 3.054e-04 after 3 F-cycles [2.916e-09, 0.102] 0 0 [21:48:18] 3.300e-05 after 4 F-cycles [3.151e-10, 0.108] 0 0 [21:48:23] 3.965e-06 after 5 F-cycles [3.786e-11, 0.120] 0 0 [21:48:27] 8.931e-07 after 6 F-cycles [8.529e-12, 0.225] 0 0 > CONVERGED > MG cycles : 6 > Final rel. error : 8.931e-07 :: emg3d END :: 21:48:27 :: runtime = 0:00:28 .. GENERATED FROM PYTHON SOURCE LINES 217-219 Plot function  .. GENERATED FROM PYTHON SOURCE LINES 219-281 .. 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 282-284 Compare the magnetic field generated from the magnetic source ------------------------------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 284-288 .. code-block:: default e3d_h = hfield.get_receiver((rx, ry, rz, azimuth, elevation)) plot(epm_h, e3d_h, r'Diffusive Fullspace $H$', vmin=-15, vmax=-8) .. image-sg:: /gallery/comparisons/images/sphx_glr_magnetic_source_duality_001.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_duality_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxElectricDipole: 3.1 A; center={0.0; 0.0; -300.0} m; θ=10.0°, φ=70.0°; l=1.0 m - Frequency: 0.77 Hz - Magnetic receivers: z=-400.0 m; θ=25°, φ=10° .. GENERATED FROM PYTHON SOURCE LINES 289-312 Compare the electric field generated from the magnetic source ------------------------------------------------------------- get_magnetic_field gets the :math:H-field from the :math:E-field with Faraday's law, .. math:: :label: faraday2 \nabla \times \mathbf{E} = \rm{i}\omega \mathbf{B} = \rm{i}\omega\mu\mathbf{H}\, . Using the substitutions introduced in the beginning, and using the same function but to get the :math:E-field from the :math:H-field, we have to multiply the result by .. math:: :label: iwu \rm{i}\omega\mu\, . Compute electric field :math:E from the magnetic field  .. GENERATED FROM PYTHON SOURCE LINES 312-320 .. code-block:: default efield = emg3d.get_magnetic_field(model, hfield) efield.field *= efield.smu0 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_duality_002.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_duality_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxElectricDipole: 3.1 A; center={0.0; 0.0; -300.0} m; θ=10.0°, φ=70.0°; l=1.0 m - Frequency: 0.77 Hz - Electric receivers: z=-400.0 m; θ=25°, φ=10° .. GENERATED FROM PYTHON SOURCE LINES 321-322 .. code-block:: default emg3d.Report() .. raw:: html
 Wed Aug 31 21:48:31 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 33.156 seconds) **Estimated memory usage:** 174 MB .. _sphx_glr_download_gallery_comparisons_magnetic_source_duality.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_duality.py  .. container:: sphx-glr-download sphx-glr-download-jupyter :download:Download Jupyter notebook: magnetic_source_duality.ipynb  .. only:: html .. rst-class:: sphx-glr-signature Gallery generated by Sphinx-Gallery _