.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "gallery/comparisons/magnetic_permeability.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_permeability.py: 9. Magnetic permeability ======================== The solver ``emg3d`` uses the diffusive approximation of Maxwell's equations; the relative electric permittivity is therefore fixed at :math:`\varepsilon_\rm{r} = 1`. The magnetic permeability :math:`\mu_\rm{r}`, however, is implemented in ``emg3d``, albeit only isotropically. In this example we run the same model as in the example mentioned above: A rotated finite length bipole in a homogeneous VTI fullspace, but here with isotropic magnetic permeability, 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 18-24 .. code-block:: default import emg3d import empymod import numpy as np import matplotlib.pyplot as plt .. GENERATED FROM PYTHON SOURCE LINES 26-31 Full-space model for a finite length, finite strength, rotated bipole --------------------------------------------------------------------- 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 31-34 .. code-block:: default coarse_model = True .. GENERATED FROM PYTHON SOURCE LINES 35-37 Survey and model parameters ``````````````````````````` .. GENERATED FROM PYTHON SOURCE LINES 37-64 .. 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 = 22 elevation = 13 # Source coordinates, frequency, and strength source = emg3d.TxElectricDipole( coordinates=[-50, 50, -30, 30, -320., -280.], # [x1, x2, y1, y2, z1, z2] strength=2.8, # A ) frequency = 0.7 # Hz # Model parameters h_res = 1. # Horizontal resistivity aniso = np.sqrt(2.) # Anisotropy v_res = h_res*aniso**2 # Vertical resistivity mperm = 2.5 # Magnetic permeability .. GENERATED FROM PYTHON SOURCE LINES 65-69 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 69-85 .. code-block:: default # Compute epm = empymod.bipole( src=np.r_[source.coordinates[:4], -source.coordinates[4:]], rec=[frx, fry, -rz, azimuth, -elevation], depth=[], res=h_res, aniso=aniso, strength=source.strength, srcpts=5, freqtime=frequency, mpermH=mperm, htarg={'pts_per_dec': -1}, verb=3, ).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 [-] : 2.5 mpermV [-] : 2.5 > MODEL IS A FULLSPACE direct field : Comp. in wavenumber domain frequency [Hz] : 0.7 Hankel : DLF (Fast Hankel Transform) > Filter : Key 201 (2009) > DLF type : Lagged Convolution Loop over : Frequencies Source(s) : 1 bipole(s) > intpts : 5 > length [m] : 123.288 > strength[A] : 2.8 > x_c [m] : 0 > y_c [m] : 0 > z_c [m] : 300 > azimuth [°] : 30.9638 > dip [°] : -18.9318 Receiver(s) : 65536 dipole(s) > x [m] : -2550 - 2550 : 65536 [min-max; #] > y [m] : -2550 - 2550 : 65536 [min-max; #] > z [m] : 400 > azimuth [°] : 22 > dip [°] : -13 Required ab's : 11 12 13 21 22 23 31 32 33 :: empymod END; runtime = 0:00:00.795417 :: 45 kernel call(s) .. GENERATED FROM PYTHON SOURCE LINES 86-88 emg3d ````` .. GENERATED FROM PYTHON SOURCE LINES 88-110 .. 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 737,280 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 96 -4,274.56 4,457.47 40.00 182.91 1.04
y 96 -4,274.56 4,457.47 40.00 182.91 1.04
z 80 -4,793.87 3,961.21 40.00 232.65 1.05


.. GENERATED FROM PYTHON SOURCE LINES 111-125 .. code-block:: default # Define the model, with magnetic permeability model = emg3d.Model( grid, property_x=h_res, property_z=v_res, mu_r=mperm, mapping='Resistivity' ) # Compute the electric field efield = emg3d.solve_source(model, source, frequency, verb=4, plain=True) # Get responses at receivers e3d = efield.get_receiver((rx, ry, rz, azimuth, elevation)) .. rst-class:: sphx-glr-script-out .. code-block:: none :: emg3d START :: 21:48:33 :: 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 80 => 737,280 cells Coarsest grid : 3 x 3 x 5 => 45 cells Coarsest level : 5 ; 5 ; 4 [hh:mm:ss] rel. error [abs. error, last/prev] l s h_ 2h_ \ / 4h_ \ /\ / 8h_ \ /\ / \ / 16h_ \ /\ / \ / \ / 32h_ \/\/ \/ \/ \/ [21:48:37] 3.578e-02 after 1 F-cycles [2.608e-05, 0.036] 0 0 [21:48:41] 3.810e-03 after 2 F-cycles [2.777e-06, 0.106] 0 0 [21:48:45] 5.278e-04 after 3 F-cycles [3.847e-07, 0.139] 0 0 [21:48:49] 8.297e-05 after 4 F-cycles [6.048e-08, 0.157] 0 0 [21:48:53] 1.414e-05 after 5 F-cycles [1.031e-08, 0.170] 0 0 [21:48:57] 2.595e-06 after 6 F-cycles [1.892e-09, 0.184] 0 0 [21:49:00] 5.524e-07 after 7 F-cycles [4.027e-10, 0.213] 0 0 > CONVERGED > MG cycles : 7 > Final rel. error : 5.524e-07 :: emg3d END :: 21:49:00 :: runtime = 0:00:27 .. GENERATED FROM PYTHON SOURCE LINES 126-128 Plot ```` .. GENERATED FROM PYTHON SOURCE LINES 128-185 .. code-block:: default # Start figure. a_kwargs = {'cmap': "viridis", 'vmin': -12, 'vmax': -6, '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 fig.colorbar(cf0, ax=axs[0, :], label=r"$\log_{10}$ Amplitude (V/m)") 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(r'Diffusive Fullspace, $E$-field', y=1, fontsize=20) print(f"- Source: {source}") print(f"- Frequency: {frequency} Hz") print(f"- Electric receivers: z={rz} m; θ={azimuth}°, φ={elevation}°") fig.show() .. image-sg:: /gallery/comparisons/images/sphx_glr_magnetic_permeability_001.png :alt: Diffusive Fullspace, $E$-field, (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_permeability_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none - Source: TxElectricDipole: 2.8 A; e1={-50.0; -30.0; -320.0} m; e2={50.0; 30.0; -280.0} m - Frequency: 0.7 Hz - Electric receivers: z=-400.0 m; θ=22°, φ=13° .. GENERATED FROM PYTHON SOURCE LINES 186-187 .. code-block:: default emg3d.Report() .. raw:: html
Wed Aug 31 21:49:02 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 30.816 seconds) **Estimated memory usage:** 186 MB .. _sphx_glr_download_gallery_comparisons_magnetic_permeability.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_permeability.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: magnetic_permeability.ipynb ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_