# 9. Magnetic permeability#

The solver emg3d uses the diffusive approximation of Maxwell’s equations; the relative electric permittivity is therefore fixed at $$\varepsilon_\rm{r} = 1$$. The magnetic permeability $$\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.)

import emg3d
import empymod
import numpy as np
import matplotlib.pyplot as plt


## 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.

coarse_model = True


### Survey and model parameters#

# 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


### 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.

# 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))

:: 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)


### emg3d#

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

 MESH EXTENT CELL WIDTH FACTOR dir nC TensorMesh 737,280 cells 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

# 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))

:: 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


### Plot#

# 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()

- 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°

emg3d.Report()

 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

Total running time of the script: ( 0 minutes 30.816 seconds)

Estimated memory usage: 186 MB

Gallery generated by Sphinx-Gallery