1. empymod: 1D VTI resistivity

The code empymod is an open-source code which can model CSEM responses for a layered medium including VTI electrical anisotropy, see emsig.xyz.

Content:

  1. Full-space VTI model for a finite length, finite strength, rotated bipole.

  2. Layered model for a deep water model with a point dipole source.

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

1. Full-space VTI 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 = 33
elevation = 18

# Source coordinates, frequency, and strength
source = emg3d.TxElectricDipole(
    coordinates=[-50, 50, -30, 30, -320., -280.],  # [x1, x2, y1, y2, z1, z2]
    strength=np.pi,  # A
)
frequency = 1.1  # Hz

# Model parameters
h_res = 1.              # Horizontal resistivity
aniso = np.sqrt(2.)     # Anisotropy
v_res = h_res*aniso**2  # Vertical resistivity

1.a Regular VTI case

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.

# Collect common input for empymod.
inp = {
    'src': np.r_[source.coordinates[:4], -source.coordinates[4:]],
    'depth': [],
    'res': h_res,
    'aniso': aniso,
    'strength': source.strength,
    'srcpts': 5,
    'freqtime': frequency,
    'htarg': {'pts_per_dec': -1},
}

# Compute
epm = empymod.bipole(
    rec=[frx, fry, -rz, azimuth, -elevation], verb=3, **inp
).reshape(np.shape(rx))

Out:

:: empymod START  ::  v2.1.2

   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] :  1.1
   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] :  3.14159
     > 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 [°] :  33
     > dip     [°] :  -18
   Required ab's   :  11 12 13 21 22 23 31 32 33

:: empymod END; runtime = 0:00:01.077291 :: 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,
)
grid
TensorMesh 512,000 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 80 -3,688.35 3,688.35 40.00 176.00 1.04
y 80 -3,688.35 3,688.35 40.00 176.00 1.04
z 80 -4,080.24 3,480.24 40.00 183.05 1.04


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

Out:

:: emg3d START :: 18:07:38 :: v1.1.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_    \/\/  \/    \/

   [18:07:42]   3.453e-02  after   1 F-cycles   [4.302e-05, 0.035]   0 0
   [18:07:46]   3.622e-03  after   2 F-cycles   [4.512e-06, 0.105]   0 0
   [18:07:50]   4.902e-04  after   3 F-cycles   [6.107e-07, 0.135]   0 0
   [18:07:54]   7.451e-05  after   4 F-cycles   [9.282e-08, 0.152]   0 0
   [18:07:58]   1.210e-05  after   5 F-cycles   [1.507e-08, 0.162]   0 0
   [18:08:02]   2.080e-06  after   6 F-cycles   [2.591e-09, 0.172]   0 0
   [18:08:06]   4.151e-07  after   7 F-cycles   [5.171e-10, 0.200]   0 0

   > CONVERGED
   > MG cycles        : 7
   > Final rel. error : 4.151e-07

:: emg3d END   :: 18:08:06 :: runtime = 0:00:28

Plot function

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

Plot

e3d = efield.get_receiver((rx, ry, rz, azimuth, elevation))
plot(epm, e3d, r'Diffusive Fullspace $E$', vmin=-12, vmax=-6)
Diffusive Fullspace $E$, (a) |Re(empymod)|, (b) |Re(emg3d)|, (c) Error real part, (d) |Im(empymod)|, (e) |Im(emg3d)|, (f) Error imaginary part

Out:

- Source: TxElectricDipole: 3.1 A;
    e1={-50.0; -30.0; -320.0} m; e2={50.0; 30.0; -280.0} m
- Frequency: 1.1 Hz
- Electric receivers: z=-400.0 m; θ=33°, φ=18°

2. Layered model for a deep water model with a point dipole source

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 = -950.0
azimuth = 30
elevation = 5

# Source coordinates and frequency
source = emg3d.TxElectricDipole(coordinates=[0, 0, -900, 0, 0])
frequency = 1.0  # Hz

# Model parameters
h_res = [1, 50, 1, 0.3, 1e12]     # Horizontal resistivity
aniso = np.sqrt([2, 2, 2, 1, 1])  # Anisotropy
v_res = h_res*aniso**2            # Vertical resistivity
depth = np.array([-2200, -2000, -1000, 0])  # Layer boundaries

empymod

epm_d = empymod.bipole(
    src=source.coordinates,
    depth=depth,
    res=h_res,
    aniso=aniso,
    freqtime=frequency,
    htarg={'pts_per_dec': -1},
    rec=[frx, fry, rz, azimuth, elevation],
    verb=3,
).reshape(np.shape(rx))

Out:

:: empymod START  ::  v2.1.2

   depth       [m] :  -2200 -2000 -1000 0
   res     [Ohm.m] :  1 50 1 0.3 1E+12
   aniso       [-] :  1.41421 1.41421 1.41421 1 1
   epermH      [-] :  1 1 1 1 1
   epermV      [-] :  1 1 1 1 1
   mpermH      [-] :  1 1 1 1 1
   mpermV      [-] :  1 1 1 1 1
   direct field    :  Comp. in wavenumber domain
   frequency  [Hz] :  1
   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] :  -900
     > azimuth [°] :  0
     > dip     [°] :  0
   Receiver(s)     :  65536 dipole(s)
     > x       [m] :  -2550 - 2550 : 65536  [min-max; #]
     > y       [m] :  -2550 - 2550 : 65536  [min-max; #]
     > z       [m] :  -950
     > azimuth [°] :  30
     > dip     [°] :  5
   Required ab's   :  11 21 31

:: empymod END; runtime = 0:00:00.103962 :: 3 kernel call(s)

emg3d

if coarse_model:
    min_width_limits = 40
else:
    min_width_limits = 20

# Create stretched grid
grid = emg3d.construct_mesh(
    frequency=frequency,
    properties=[h_res[3], h_res[0]],
    center=source.center,
    domain=([-2500, 2500], [-2500, 2500], None),
    vector=(None, None, -2200 + np.arange(111)*20),
    min_width_limits=min_width_limits,
    stretching=[1.1, 1.5],
)
grid
TensorMesh 655,360 cells
MESH EXTENT CELL WIDTH FACTOR
dir nC min max min max max
x 64 -5,681.71 5,681.71 40.00 1,159.16 1.45
y 64 -5,681.71 5,681.71 40.00 1,159.16 1.45
z 160 -5,508.08 3,308.08 20.00 390.86 1.13


# Create the model: horizontal resistivity
res_x_full = h_res[0]*np.ones(grid.n_cells)  # Background
res_x_full[grid.cell_centers[:, 2] >= depth[0]] = h_res[1]  # Target
res_x_full[grid.cell_centers[:, 2] >= depth[1]] = h_res[2]  # Overburden
res_x_full[grid.cell_centers[:, 2] >= depth[2]] = h_res[3]  # Water
res_x_full[grid.cell_centers[:, 2] >= depth[3]] = h_res[4]  # Air

# Create the model: vertical resistivity
res_z_full = v_res[0]*np.ones(grid.n_cells)  # Background
res_z_full[grid.cell_centers[:, 2] >= depth[0]] = v_res[1]
res_z_full[grid.cell_centers[:, 2] >= depth[1]] = v_res[2]
res_z_full[grid.cell_centers[:, 2] >= depth[2]] = v_res[3]
res_z_full[grid.cell_centers[:, 2] >= depth[3]] = v_res[4]

# Get the model
model = emg3d.Model(
        grid, property_x=res_x_full, property_z=res_z_full,
        mapping='Resistivity')
# Compute the electric field
efield = emg3d.solve_source(model, source, frequency, verb=4)

Out:

:: emg3d START :: 18:08:11 :: v1.1.0

   MG-cycle       : 'F'                 sslsolver : 'bicgstab'
   semicoarsening : True [1 2 3]        tol       : 1e-06
   linerelaxation : True [4 5 6]        maxit     : 50 (3)
   nu_{i,1,c,2}   : 0, 2, 1, 2          verb      : 4
   Original grid  :  64 x  64 x 160     => 655,360 cells
   Coarsest grid  :   2 x   2 x   5     => 20 cells
   Coarsest level :   5 ;   5 ;   5

   [hh:mm:ss]  rel. error            solver              MG          l s

       h_
      2h_ \                            /
      4h_  \                  /\      /
      8h_   \          /\    /  \    /
     16h_    \    /\  /  \  /    \  /
     32h_     \/\/  \/    \/      \/

   [18:08:30]   2.394e-03  after                       1 F-cycles    4 1
   [18:08:47]   5.999e-05  after                       2 F-cycles    5 2
   [18:09:02]   1.189e-04  after                       3 F-cycles    6 3
   [18:09:21]   4.446e-06  after                       4 F-cycles    4 1
   [18:09:38]   1.899e-06  after                       5 F-cycles    5 2
   [18:09:54]   1.063e-05  after                       6 F-cycles    6 3
   [18:09:54]   2.611e-06  after   1 bicgstab-cycles
   [18:10:13]   7.479e-08  after                       7 F-cycles    4 1
   [18:10:31]   5.132e-07  after                       8 F-cycles    5 2
   [18:10:31]   5.113e-07  after   2 bicgstab-cycles

   > CONVERGED
   > Solver steps     : 2
   > MG prec. steps   : 8
   > Final rel. error : 5.113e-07

:: emg3d END   :: 18:10:31 :: runtime = 0:02:20

Plot

e3d_d = efield.get_receiver((rx, ry, rz, azimuth, elevation))
plot(epm_d, e3d_d, r'Deep water point dipole $E$', vmin=-14, vmax=-8)
Deep water point dipole $E$, (a) |Re(empymod)|, (b) |Re(emg3d)|, (c) Error real part, (d) |Im(empymod)|, (e) |Im(emg3d)|, (f) Error imaginary part

Out:

- Source: TxElectricDipole: 1.0 A;
    center={0.0; 0.0; -900.0} m; θ=0.0°, φ=0.0°; l=1.0 m
- Frequency: 1.0 Hz
- Electric receivers: z=-950.0 m; θ=30°, φ=5°
Sat Jul 17 18:10:34 2021 CEST
OS Linux CPU(s) 4 Machine x86_64
Architecture 64bit RAM 15.5 GiB Environment Python
Python 3.9.4 | packaged by conda-forge | (default, May 10 2021, 22:13:33) [GCC 9.3.0]
numpy 1.21.0 scipy 1.7.0 numba 0.53.1
emg3d 1.1.0 empymod 2.1.2 xarray 0.18.2
discretize 0.7.0 h5py 3.3.0 matplotlib 3.4.2
tqdm 4.61.2
Intel(R) oneAPI Math Kernel Library Version 2021.2-Product Build 20210312 for Intel(R) 64 architecture applications


Total running time of the script: ( 2 minutes 59.485 seconds)

Estimated memory usage: 471 MB

Gallery generated by Sphinx-Gallery