Code¶

Electromagnetic modeller in the diffusive limit (low frequencies) for 3D media with tri-axial electrical anisotropy. The matrix-free multigrid solver can be used as main solver or as preconditioner for one of the Krylov subspace methods implemented in scipy.sparse.linalg, and the governing equations are discretized on a staggered Yee grid. The code is written completely in Python using the numpy/scipy-stack, where the most time-consuming parts are sped-up through jitted numba-functions.

`solver` – Multigrid solver¶

The actual solver routines. The most computationally intensive parts, however, are in the emg3d.core as numba-jitted functions.

emg3d.solver.multigrid(grid, model, sfield, efield, var, **kwargs)[source]¶

Multigrid solver for 3D controlled-source electromagnetic (CSEM) data.

Multigrid solver as presented in [Muld06], including semicoarsening and line relaxation as presented in and [Muld07].

The electric field is stored in-place in efield.
The number of multigrid cycles is stored in var.it.
The current error (l2-norm) is stored in var.l2.
The reference error (l2-norm of sfield) is stored in var.l2_refe.

This function is called by solve().

Parameters:

grid : emg3d.meshes.TensorMesh: The grid. See emg3d.meshes.TensorMesh.
model : emg3d.models.VolumeModel: The Model. See emg3d.models.VolumeModel.
sfield : emg3d.fields.SourceField: The source field. See emg3d.fields.get_source_field().
efield : emg3d.fields.Field: The electric field. See emg3d.fields.Field.
var : MGParameters instance: As returned by multigrid().
**kwargs : Recursion parameters.: Do not use; only used internally by recursion; level (current coarsening level) and new_cycmax (new maximum of MG cycles, takes care of V/W/F-cycling).

emg3d.solver.smoothing(grid, model, sfield, efield, nu, lr_dir)[source]¶

Reducing high-frequency error by smoothing.

Solves the linear equation system \(A x = b\) iteratively using the Gauss-Seidel method. This acts as smoother or, on the coarsest grid, as a direct solver.

This is a simple wrapper for the jitted calculation in emg3d.core.gauss_seidel(), emg3d.core.gauss_seidel_x(), emg3d.core.gauss_seidel_y(), and emg3d.core.gauss_seidel_z() (@njit can not [yet] access class attributes). See these functions for more details and corresponding theory.

The electric fields are updated in-place.

This function is called by multigrid().

Parameters:

grid : emg3d.meshes.TensorMesh: Input grid.
model : emg3d.models.VolumeModel: Input model.
sfield : emg3d.fields.SourceField: Input source field.
efield : emg3d.fields.Field: Input electric field.
nu : int: Number of Gauss-Seidel steps; odd numbers are forward, even numbers are reversed. E.g., nu=2 is one symmetric Gauss-Seidel iteration, with a forward and a backward step.
lr_dir : int: Direction of line relaxation {0, 1, 2, 3, 4, 5, 6, 7}.

emg3d.solver.restriction(grid, model, sfield, residual, sc_dir)[source]¶

Downsampling of grid, model, and fields to a coarser grid.

The restriction of the residual is used as source term for the coarse grid.

Corresponds to Equations 8 and 9 in [Muld06] and surrounding text. In the case of the restriction of the residual, this function is a wrapper for the jitted functions emg3d.core.restrict_weights() and emg3d.core.restrict() (@njit can not [yet] access class attributes). See these functions for more details and corresponding theory.

This function is called by multigrid().

Parameters:	grid : `emg3d.meshes.TensorMesh` Input grid. model : `emg3d.models.VolumeModel` Input model. sfield : `emg3d.fields.SourceField` Input source field. sc_dir : int Direction of semicoarsening (0, 1, 2, or 3).
Returns:	cgrid : `emg3d.meshes.TensorMesh` Coarse grid. cmodel : `emg3d.models.VolumeModel` Coarse model. csfield : `emg3d.fields.SourceField` Coarse source field. Corresponds to restriction of fine-grid residual. cefield : `emg3d.fields.Field` Coarse electric field, complex zeroes.

emg3d.solver.prolongation(grid, efield, cgrid, cefield, sc_dir)[source]¶

Interpolating the electric field from coarse grid to fine grid.

The prolongation from a coarser to a finer grid is the inverse process of the restriction (restriction()) from a finer to a coarser grid. The interpolated values of the coarse grid electric field are added to the fine grid electric field, in-place. Piecewise constant interpolation is used in the direction of the field, and bilinear interpolation in the other two directions.

See Equation 10 in [Muld06] and surrounding text.

This function is called by multigrid().

Parameters:	grid, cgrid : `emg3d.meshes.TensorMesh` Fine and coarse grids. efield, cefield : `emg3d.fields.Field` Fine and coarse grid electric fields. sc_dir : int Direction of semicoarsening (0, 1, 2, or 3).

emg3d.solver.residual(grid, model, sfield, efield, norm=False)[source]¶

Calculating the residual.

Returns the complete residual as given in [Muld06], page 636, middle of the right column:

\[\mathbf{r} = V \left( \mathrm{i}\omega\mu_0\mathbf{J_s} + \mathrm{i}\omega\mu_0 \tilde{\sigma} \mathbf{E} - \nabla \times \mu_\mathrm{r}^{-1} \nabla \times \mathbf{E} \right) .\]

This is a simple wrapper for the jitted calculation in emg3d.core.amat_x() (@njit can not [yet] access class attributes). See emg3d.core.amat_x() for more details and corresponding theory.

This function is called by multigrid().

Parameters:	grid : `emg3d.meshes.TensorMesh` Input grid. model : `emg3d.models.VolumeModel` Input model. sfield : `emg3d.fields.SourceField` Input source field. efield : `emg3d.fields.Field` Input electric field. norm : bool If True, the error (l2-norm) of the residual is returned, not the residual.
Returns:	residual : Field Returned if `norm=False`. The residual field; `emg3d.fields.Field` instance. norm : float Returned if `norm=True`. The error (l2-norm) of the residual

emg3d.solver.krylov(grid, model, sfield, efield, var)[source]¶

Krylov Subspace iterative solver for 3D CSEM data.

Using a Krylov subspace iterative solver (defined in var.sslsolver) implemented in SciPy with or without multigrid as a pre-conditioner ([Muld06]).

The electric field is stored in-place in efield.
The current error (l2-norm) is stored in var.l2.
The reference error (l2-norm of sfield) is stored in var.l2_refe.

This function is called by solve().

Parameters:

grid : emg3d.meshes.TensorMesh: The grid. See emg3d.meshes.TensorMesh.
model : emg3d.models.VolumeModel: The Model. See emg3d.models.VolumeModel.
sfield : emg3d.fields.SourceField: The source field. See emg3d.fields.get_source_field().
efield : emg3d.fields.Field: The electric field. See emg3d.fields.Field.
var : MGParameters instance: As returned by multigrid().

class emg3d.solver.MGParameters(verb: int, cycle: str, sslsolver: str, linerelaxation: int, semicoarsening: int, vnC: tuple, tol: float = 1e-06, maxit: int = 50, nu_init: int = 0, nu_pre: int = 2, nu_coarse: int = 1, nu_post: int = 2, clevel: int = -1, return_info: bool = False)[source]¶

Collect multigrid solver settings.

This dataclass is used by the main solve()-routine. See solve() for a description of the mandatory and optional input parameters and more information .

Returns:	var : class:MGParameters As required by `multigrid()`.

cprint(self, info, verbosity, **kwargs)[source]¶

Conditional printing.

Prints info if self.verb > verbosity.

Parameters:	info : str String to be printed. verbosity : int Verbosity of info. kwargs : optional Arguments passed to print.

max_level¶

Sets dimension-dependent level variable clevel.

Requires at least two cells in each direction (for nCx, nCy, and nCz).

one_liner(self, l2_last, last=False)[source]¶

Print continuously updated one-liner.

Parameters:	l2_last : float Current error. last : bool If True, adds exit_message and finishes line.

class emg3d.solver.RegularGridProlongator(x, y, cxy)[source]¶

Prolongate field from coarse to fine grid.

This is a heavily modified and adapted version of scipy.interpolate.RegularGridInterpolator.

The main difference (besides the pre-sets) is that this version allows to initiate an instance with the coarse and fine grids. This initialize will calculate the required weights, and it has therefore only to be done once.

After this, interpolating values from the coarse to the fine grid can be carried out much faster.

Simplifications in comparison to scipy.interpolate.RegularGridInterpolator:

No sanity checks what-so-ever.
Only 2D data;
method='linear';
bounds_error=False;
fill_value=None.

It results in a speed-up factor of about 2, independent of grid size, for this particular case. The prolongation is the second-most expensive part of multigrid after the smoothing. Trying to improve this further might therefore be useful.

Parameters:	x, y : ndarray The x/y-coordinates defining the coarse grid. cxy : ndarray of shape (…, 2) The ([[x], [y]]).T-coordinates defining the fine grid.

`core` – Number crunching¶

The core functionalities, the most computationally demanding parts, of the emg3d.solver as just-in-time (jit) compiled functions using numba.

emg3d.core.amat_x(rx, ry, rz, ex, ey, ez, eta_x, eta_y, eta_z, zeta, hx, hy, hz)[source]¶

Residual without or with source term.

Calculate the residual as given in [Muld06] in middle of the right column on page 636, but without the source term:

\[\mathbf{r} = V \left( \mathrm{i}\omega\mu_0 \tilde{\sigma} \mathbf{E} - \nabla \times \mu_\mathrm{r}^{-1} \nabla \times \mathbf{E} \right) .\]

The calculation is carried out in a matrix-free manner; on said page 636 (or in the Theory) are the various steps laid out to discretise the different parts, for instance involved curls. This can also be understood as the left-hand-side of \(A x = b\), as given in Equation 2 in [Muld06] (here without the cell volumes V),

\[\mathrm{i}\omega\mu_0 \tilde{\sigma} \mathrm{E} - \nabla \times \zeta^{-1} \nabla \times \mathrm{E} = - \mathrm{i} \omega \mu_0 \mathrm{J_s} .\]

It can therefore be used as matvec to create a LinearOperator, which can be passed to a solver.

It is assumed that ex, ey, and ez have PEC boundaries; otherwise the output will not have PEC boundaries.

The residuals are subtracted in-place from rx, ry, and rz. That means that if rx, ry, and rz contain the source field, they will contain the total residual afterwards; if they are empty fields, they will contain the negative partial residual afterwards.

Parameters:

rx, ry, rz : ndarray: Source field or pre-allocated zero residual field in x-, y-, and z-directions.
ex, ey, ez : ndarray: Electric fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
eta_x, eta_y, eta_z, zeta : ndarray: VolumeModel parameters (multiplied by volumes) as obtained from emg3d.models.VolumeModel().
hx, hy, hz : ndarray: Cell widths in x-, y-, and z-directions.

emg3d.core.blocks_to_amat(amat, bvec, middle, left, rhs, im, nC)[source]¶

Insert middle, left, and rhs into main arrays amat and bvec.

The banded matrix amat contains the main diagonal and the first five lower off-diagonals. They are stored one column after the other, in a 6*n ndarray.

The complete main matrix amat and the middle and left blocks are given by:

.-0
|X|\   0
0-.-0       left:  middle:  right:
 \|X|\                      (not used)
  0-.-0      0-     .-      0
   \|X|\      \     |X      |\
    0-.-0
 0   \|X|
      0-.

. 1*1, - 4*1, | 1*4, X 4*4, \ 4*4 upper or lower

Both, middle and left, are 5x5 matrices. The corresponding right-hand-side rhs is filled into bvec. The matrices left and middle provided in a single call are horizontally aligned (not vertically). The sorting of amat (banded matrix) and bvec are given by:

 amat (66,)             example: n = 11                   bvec (11,)
 --------------                                                 --
|01            |                    FIRST CALL                  01
|02 07         |                    Only `middle` and `rhs`     02
|03 08 13      |                    are used, not `left`.       03
|04 09 14 19   |                                                04
|05 10 15 20 25|                                                05
 -------------- --------------                                  --
| 0 11 16 21 26|31            |     SUBSEQUENT CALLS            06
|   12 17 22 27|32 37         |     (normal case)               07
|      18 23 28|33 38 43      |     Complete `left`,            08
|         24 29|34 39 44 49   |     `middle` and `rhs`          09
|            30|35 40 45 50 55|     are used.                   10
 -------------- -------------- ---                              --
               | 0 41 46 51 56|61   LAST CALL                   11
               |    0  0  0  0| 0   Only top row of `left`
               |       0  0  0| 0   and the first elements
               |          0  0| 0   of `middle` and `rhs`
               |             0| 0   are used.
                -------------- ---
                              | 0

Single zeros (0) show elements in amat which are 0, hence not used.
Their location in amat can be deduced from their neighbours.

Parameters:

amat : ndarray: Main banded matrix (stored as array) of length 6*n.
bvec : ndarray: Main right-hand-side of length n.
middle : ndarray: Middle block of size 5x5, as ndarray of length 25. Only the diagonal and the lower triangular part are used.
left : ndarray: Left block of size 5x5, as ndarray of length 25. Only the diagonal and the first row are used.
rhs : ndarray: Corresponding right-hand-side of length 5.
im : int: Current minus-index of direction of line relaxation, {ixm, iym, izm}.
nC : int: Total number of cells in direction of line relaxation, {nCx, nCy, nCz}.

emg3d.core.gauss_seidel(ex, ey, ez, sx, sy, sz, eta_x, eta_y, eta_z, zeta, hx, hy, hz, nu)[source]¶

Gauss-Seidel method.

Solves the linear equation system \(A x = b\) iteratively using the following method:

\[\mathbf{x}^{(k+1)} = L_*^{-1} \left(\mathbf{b} - U \mathbf{x}^{(k)} \right) \ ,\]

where \(L_*\) is the lower triangular component, and \(U\) the strictly upper triangular component, \(A = L_* + U\):

\[\begin{split}L_* = \left[ \begin{array} {cccc} a_{11} & 0 & \cdots & 0 \\ a_{21} & a_{22} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{array} \right] \ , \quad U = \left[ \begin{array} {cccc} 0 & a_{12} & \cdots & a_{1n} \\ 0 & 0 & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{array} \right] \ .\end{split}\]

On the coarsest grid it acts as direct solver, whereas on the fine grid it acts as a smoother with only few iterations, defined by \(\nu\) (nu). Odd numbers of nu use forward ordering, even numbers use backwards ordering. nu=2 is therefore one symmetric Gauss-Seidel iteration, one forward ordered iteration followed by one backward ordered iteration.

From [Muld06]: The method proposed by [ArFW00] is chosen as a smoother. It selects one node of the grid and simultaneously solves for the six degrees of freedom on the six edges attached to the node. If node \((x_k, y_l, z_m)\) is selected, the six equations, \(r_{x;k\pm1/2,l,m} = 0\), \(r_{y;k,l\pm1/2,m} = 0\) and \(r_{z;k,l,m\pm1/2} = 0\), are solved for \(e_{x;k\pm1/2,l,m}\), \(e_{y;k,l\pm1/2,m}\) and \(e_{z;k,l,m\pm1/2}\). Here, this smoother is applied in a symmetric Gauss-Seidel fashion, following the lexicographical ordering of the nodes \((x_k, y_l, z_m)\), with fastest index \(k=1, \dots, N_x-1\), intermediate index \(l=1, \dots, N_y-1\), and slowest index \(m=1, \ldots, N_z-1\).

To actually solve the system of six equations a non-standard Cholesky factorisation is used, solve().

Tangential components at the boundaries are assumed to be zero (PEC boundaries).

The result is stored in the provided electric fields ex, ey, and ez.

Parameters:

ex, ey, ez : ndarray: Electric fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
sx, sy, sz :: Source fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
eta_x, eta_y, eta_z, zeta :: VolumeModel parameters (multiplied by volumes) as obtained from emg3d.models.VolumeModel().
hx, hy, hz : ndarray: Cell widths in x-, y-, and z-directions.
nu : int: Number of Gauss-Seidel iterations.

emg3d.core.gauss_seidel_x(ex, ey, ez, sx, sy, sz, eta_x, eta_y, eta_z, zeta, hx, hy, hz, nu)[source]¶

Gauss-Seidel method with line relaxation in x-direction.

This is the equivalent to gauss_seidel(), but with line relaxation in the x-direction. See gauss_seidel() for more details.

The resulting system A x = b to solve consists of n unknowns (x-vector), and the corresponding matrix A is a banded matrix with the main diagonal and five upper and lower diagonals:

.-0
|X|\   0
0-.-0       left:  middle:  right:
 \|X|\                      (not used)
  0-.-0      0-     .-      0
   \|X|\      \     |X      |\
    0-.-0
 0   \|X|
      0-.

. 1*1, - 4*1, | 1*4, X 4*4, \ 4*4 upper or lower

The matrix A is complex and symmetric (A = A^T), and therefore only the main diagonal and the lower five off-diagonals are required.

The right-hand-side b has length 5*nCx-4 (nCx even).
The matrix A has length of b and 1+2*5 diagonals; we use for it an array of length 6*len(b).

The values are calculated in rows of 5 lines, with the indicated middle and left matrices as indicated in the above scheme. These blocks are filled into the main matrix A and vector b, and subsequently solved with a non-standard Cholesky factorisation, solve().

Tangential components at the boundaries are assumed to be 0 (PEC boundaries).

The result is stored in the provided electric fields ex, ey, and ez.

Parameters:

ex, ey, ez : ndarray: Electric fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
sx, sy, sz :: Source fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
eta_x, eta_y, eta_z, zeta :: VolumeModel parameters (multiplied by volumes) as obtained from emg3d.models.VolumeModel().
hx, hy, hz : ndarray: Cell widths in x-, y-, and z-directions.
nu : int: Number of Gauss-Seidel iterations.

emg3d.core.gauss_seidel_y(ex, ey, ez, sx, sy, sz, eta_x, eta_y, eta_z, zeta, hx, hy, hz, nu)[source]¶

Gauss-Seidel method with line relaxation in y-direction.

This is the equivalent to gauss_seidel(), but with line relaxation in the y-direction. See gauss_seidel() for more details.

The resulting system A x = b to solve consists of n unknowns (x-vector), and the corresponding matrix A is a banded matrix with the main diagonal and five upper and lower diagonals:

.-0
|X|\   0
0-.-0       left:  middle:  right:
 \|X|\                      (not used)
  0-.-0      0-     .-      0
   \|X|\      \     |X      |\
    0-.-0
 0   \|X|
      0-.

. 1*1, - 4*1, | 1*4, X 4*4, \ 4*4 upper or lower

The matrix A is complex and symmetric (A = A^T), and therefore only the main diagonal and the lower five off-diagonals are required.

The right-hand-side b has length 5*nCy-4 (nCy even).
The matrix A has length of b and 1+2*5 diagonals; we use for it an array of length 6*len(b).

The values are calculated in rows of 5 lines, with the indicated middle and left matrices as indicated in the above scheme. These blocks are filled into the main matrix A and vector b, and subsequently solved with a non-standard Cholesky factorisation, solve().

Note: The smoothing with linerelaxation in y-direction is carried out in reversed lexicographical order, in order to improve speed (memory access). All other smoothers (gauss_seidel(), gauss_seidel_x(), and gauss_seidel_z()) use lexicographical order.

Tangential components at the boundaries are assumed to be 0 (PEC boundaries).

The result is stored in the provided electric fields ex, ey, and ez.

Parameters:

ex, ey, ez : ndarray: Electric fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
sx, sy, sz :: Source fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
eta_x, eta_y, eta_z, zeta :: VolumeModel parameters (multiplied by volumes) as obtained from emg3d.models.VolumeModel().
hx, hy, hz : ndarray: Cell widths in x-, y-, and z-directions.
nu : int: Number of Gauss-Seidel iterations.

emg3d.core.gauss_seidel_z(ex, ey, ez, sx, sy, sz, eta_x, eta_y, eta_z, zeta, hx, hy, hz, nu)[source]¶

Gauss-Seidel method with line relaxation in z-direction.

This is the equivalent to gauss_seidel(), but with line relaxation in the z-direction. See gauss_seidel() for more details.

The resulting system A x = b to solve consists of n unknowns (x-vector), and the corresponding matrix A is a banded matrix with the main diagonal and five upper and lower diagonals:

.-0
|X|\   0
0-.-0       left:  middle:  right:
 \|X|\                      (not used)
  0-.-0      0-     .-      0
   \|X|\      \     |X      |\
    0-.-0
 0   \|X|
      0-.

. 1*1, - 4*1, | 1*4, X 4*4, \ 4*4 upper or lower

The matrix A is complex and symmetric (A = A^T), and therefore only the main diagonal and the lower five off-diagonals are required.

The right-hand-side b has length 5*nCz-4 (nCz even).
The matrix A has length of b and 1+2*5 diagonals; we use for it an array of length 6*len(b).

The values are calculated in rows of 5 lines, with the indicated middle and left matrices as indicated in the above scheme. These blocks are filled into the main matrix A and vector b, and subsequently solved with a non-standard Cholesky factorisation, solve().

Tangential components at the boundaries are assumed to be 0 (PEC boundaries).

The result is stored in the provided electric fields ex, ey, and ez.

Parameters:

ex, ey, ez : ndarray: Electric fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
sx, sy, sz :: Source fields in x-, y-, and z-directions, as obtained from emg3d.fields.Field.
eta_x, eta_y, eta_z, zeta :: VolumeModel parameters (multiplied by volumes) as obtained from emg3d.models.VolumeModel().
hx, hy, hz : ndarray: Cell widths in x-, y-, and z-directions.
nu : int: Number of Gauss-Seidel iterations.

emg3d.core.restrict(crx, cry, crz, rx, ry, rz, wx, wy, wz, sc_dir)[source]¶

Restriction of residual from fine to coarse grid.

Corresponds to Equation 8 in [Muld06]. The equation for the x-direction, using the notation \(\{x,y,z\}\) instead of \(\{1,2,3\}\), is given by

\[\begin{split}r_{x,K+1/2,L,M}^{2h} = &\sum_{j_y=-1}^1\sum_{j_z=-1}^1 w_{L,j_y}^y w_{M,j_z}^z \\ &\times \left(r_{x,k+1/2,l+j_y,m+j_z}^h+r_{x,k+3/2,l+j_y,m+j_z}^h\right) .\end{split}\]

The superscripts \(h, 2h\) indicate quantities defined on the coarse grid and on the fine grid, respectively. The indices \(\{K, L, M\}\) on the coarse grid correspond to \(\{k, l, m\} = 2\{K, L, M\}\) on the fine grid. The weights \(w\) are obtained from restrict_weights().

The restrictions of rx, ry, and rz are stored directly in crx, cry, and crz.

Parameters:	crx, cry, crz : ndarray Coarse grid {x,y,z}-directed residual (pre-allocated empty arrays). rx, ry, rz : ndarray Fine grid {x,y,z}-directed residual. wx, wy, wz: tuple Tuples containing the weights (wl, w0, wr) as returned from `restrict_weights()` for the x-, y-, and z-directions. sc_dir : int Direction of semicoarsening; 0 for no semicoarsening.

emg3d.core.restrict_weights(vectorN, vectorCC, h, cvectorN, cvectorCC, ch)[source]¶

Restriction weights for the coarse-grid correction operator.

Corresponds to Equation 9 in [Muld06]. A generalized version of that equation is given by

\[\begin{split}w_{Q,-1}^v &= \left(v_{q-1/2}^h-v_{Q-1/2}^{2h}\right)/d_{q-1}^v ,\\ w_{Q,0}^v &= 1 ,\\ w_{Q,1}^v &= \left(v_{Q+1/2}^{2h}-v_{q+1/2}^h \right)/d_{q+1}^v ,\end{split}\]

where \(d\) are the dual grid cell widths, \(v\) is one of \(\{x, y, z\}\), and \(Q, q\) the corresponding entries of \(\{K, L, M\}, \{k, l, m\}\). The superscripts \(h, 2h\) indicate quantities defined on the coarse grid and on the fine grid, respectively. The indices \(\{K, L, M\}\) on the coarse grid correspond to \(\{k, l, m\} = 2\{K, L, M\}\) on the fine grid.

For the dual volume cell widths at the boundaries the scheme of [MoSu94] is applied, where \(d_0^x = h_{1/2}^x/2\) at \(k = 0\), \(d_{N_x}^x = h_{N_x-1/2}^x\) at \(k = N_x\), and so on.

The following parameters must all be in the same direction, hence, all must be either for the x, the y, or the z direction. The returned weights are for this direction.

Parameters:	vectorN, cvectorN : ndarray Cell edges of the fine (vectorN) and coarse (cvectorN) grids. vectorCC, cvectorCC : ndarray Cell centers of the fine (vectorCC) and coarse (cvectorCC) grids. h, ch : ndarray Cell widths of the fine (h) and coarse (ch) grids.
Returns:	wl, w0, wr : ndarray Left, central, and right weights in the direction provided in the input.

emg3d.core.solve(amat, bvec)[source]¶

Solve A x = b using a non-standard Cholesky factorisation.

Solve the system A x = b using a non-standard Cholesky factorisation without pivoting for a symmetric, complex matrix A tailored to the problem of the multigrid solver. The matrix A (amat) is an array of length 6*n, containing the main diagonal and the first five lower off-diagonals (ordered so that the first element of the main diagonal is followed by the first elements of the off diagonals, then the second elements and so on). The vector bvec has length b.

The solution is placed in b (bvec), and A (amat) is replaced by its decomposition.

Non-standard Cholesky factorisation.

From [Muld07]: We use a non-standard Cholesky factorisation. The standard factorisation factors a hermitian matrix A into L L^H, where L is a lower triangular matrix and L^H its complex conjugate transpose. In our case, the discretisation is based on the Finite Integration Technique ([Weil77]) and provides a matrix A that is complex-valued and symmetric: A = A^T, where the superscript T denotes the transpose. The line relaxation scheme takes a matrix B that is a subset of A along the line. B is a complex symmetric band matrix with eleven diagonals. The non-standard Cholesky factorisation factors the matrix B into L L^T. Because of the symmetry, only the main diagonal and five lower diagonal elements of B need to be computed. The Cholesky factorisation replaces this matrix by L, containing six diagonals, after which the line relaxation can be carried out by simple back-substitution.

\(A = L D L^T\) factorisation without pivoting:

\[\begin{split}D(j) &= A(j,j)-\sum_{k=1}^{j-1} L(j,k)^2 D(k),\ j=1,..,n ;\\ L(i,j) &= \frac{1}{D(j)} \left[A(i,j)-\sum_{k=1}^{j-1} L(i,k)L(j,k)D(k)\right], \ i=j+1,..,n .\end{split}\]

A and L are in this case arrays, where \(A(i, j) \rightarrow A(i+5j)\).
Solve A x = b.

Solve A x = b, given L which is the result from the factorisation in the first step (and stored in A), hence, solve L x = b, where x is stored in b:

\[b(j) = b(j) - \sum_{k=1}^{j-1} L(j,k) x(k), j = 2,..,n .\]

The result is equivalent with simply using numpy.linalg.solve(), but faster for the particular use-case of this code.

Note that in this custom solver there is no pivoting, and the diagonals of the matrix cannot be zero.

Parameters:	amat : ndarray Banded matrix A provided as a vector of length 6n, containing main diagonal plus first five lower diagonals. bvec* : ndarray Right-hand-side vector b of length n.

`utils` – Utilities¶

Utility functions for the multigrid solver.

class emg3d.utils.Fourier(time, fmin, fmax, signal=0, ft='dlf', ftarg=None, **kwargs)[source]¶

Time-domain CSEM calculation.

Class to carry out time-domain modelling with the frequency-domain code emg3d. Instances of the class take care of calculating the required frequencies, the interpolation from coarse, limited-band frequencies to the required frequencies, and carrying out the actual transform.

Everything related to the Fourier transform is done by utilising the capabilities of the 1D modeller empymod. The input parameters time, signal, ft, and ftarg are passed to the function empymod.utils.check_time() to obtain the required frequencies. The actual transform is subsequently carried out by calling empymod.model.tem(). See these functions for more details about the exact implementations of the Fourier transforms and its parameters. Note that also the verb-argument follows the definition in empymod.

The mapping from calculated frequencies to the frequencies required for the Fourier transform is done in three steps:

Data for \(f>f_\mathrm{max}\) is set to 0+0j.
Data for \(f<f_\mathrm{min}\) is interpolated by adding an additional data point at a frequency of 1e-100 Hz. The data for this point is data.real[0]+0j, hence the real part of the lowest calculated frequency and zero imaginary part. Interpolation is carried out using PCHIP scipy.interpolate.pchip_interpolate().
Data for \(f_\mathrm{min}\le f \le f_\mathrm{max}\) is calculated with cubic spline interpolation (on a log-scale) scipy.interpolate.InterpolatedUnivariateSpline.

Note that fmin and fmax should be chosen wide enough such that the mapping for \(f>f_\mathrm{max}\) \(f<f_\mathrm{min}\) does not matter that much.

Parameters:

time : ndarray

Desired times (s).

fmin, fmax : float

Minimum and maximum frequencies (Hz) to calculate:

Data for freq > fmax is set to 0+0j.

Data for freq < fmin is interpolated, using an extra data-point at f = 1e-100 Hz, with value data.real[0]+0j. (Hence zero imaginary part, and the lowest calculated real value.)

signal : {0, 1, -1}, optional

Source signal, default is 0:

None: Frequency-domain response
-1 : Switch-off time-domain response
0 : Impulse time-domain response
+1 : Switch-on time-domain response

ft : {‘sin’, ‘cos’, ‘fftlog’}, optional

Flag to choose either the Digital Linear Filter method (Sine- or Cosine-Filter) or the FFTLog for the Fourier transform. Defaults to ‘sin’.

ftarg : dict, optional

Depends on the value for ft:

If ft=’dlf’:

dlf: string of filter name in empymod.filters or the filter method itself. (Default: empymod.filters.key_201_CosSin_2012())

pts_per_dec: points per decade; (default: -1)

If 0: Standard DLF.

If < 0: Lagged Convolution DLF.

If > 0: Splined DLF

If ft=’fftlog’:

pts_per_dec: sampels per decade (default: 10)

add_dec: additional decades [left, right] (default: [-2, 1])

q: exponent of power law bias (default: 0); -1 <= q <= 1

freq_inp : array

Frequencies to use for calculation. Mutually exclusive with every_x_freq.

every_x_freq : int

Every every_x_freq-th frequency of the required frequency-range is used for calculation. Mutually exclusive with freq_calc.

every_x_freq¶: If set, freq_coarse is every_x_freq-frequency of freq_req.

fmax¶: Maximum frequency (Hz) to calculate.

fmin¶: Minimum frequency (Hz) to calculate.

fourier_arguments(self, ft, ftarg)[source]¶: Set Fourier type and its arguments.

freq2time(self, fdata, off)[source]¶

Calculate corresponding time-domain signal.

Carry out the actual Fourier transform.

Parameters:	fdata : ndarray Frequency-domain data corresponding to freq_calc. off : float Corresponding offset (m).
Returns:	tdata : ndarray Time-domain data corresponding to Fourier.time.

freq_calc¶: Frequencies at which the model has to be calculated.

freq_calc_i¶: Indices of freq_coarse which have to be calculated.

freq_coarse¶: Coarse frequency range, can be different from freq_req.

freq_extrapolate¶

These are the frequencies to extrapolate.

In fact, it is dow via interpolation, using an extra data-point at f = 1e-100 Hz, with value data.real[0]+0j. (Hence zero imaginary part, and the lowest calculated real value.)

freq_extrapolate_i¶: Indices of the frequencies to extrapolate.

freq_inp¶: If set, freq_coarse is set to freq_inp.

freq_interpolate¶

These are the frequencies to interpolate.

If freq_req is equal freq_coarse, then this is eual to freq_calc.

freq_interpolate_i¶

Indices of the frequencies to interpolate.

If freq_req is equal freq_coarse, then this is eual to freq_calc_i.

freq_req¶: Frequencies required to carry out the Fourier transform.

ft¶: Type of Fourier transform. Set via fourier_arguments(ft, ftarg).

ftarg¶: Fourier transform arguments. Set via fourier_arguments(ft, ftarg).

interpolate(self, fdata)[source]¶

Interpolate from calculated data to required data.

Parameters:	fdata : ndarray Frequency-domain data corresponding to freq_calc.
Returns:	full_data : ndarray Frequency-domain data corresponding to freq_req.

signal¶: Signal in time domain {0, 1, -1}.

time¶: Desired times (s).

class emg3d.utils.Time[source]¶

Class for timing (now; runtime).

elapsed¶: Return runtime in seconds since time zero.

now¶: Return string of current time.

runtime¶: Return string of runtime since time zero.

t0¶: Return time zero of this class instance.

class emg3d.utils.Report(add_pckg=None, ncol=3, text_width=80, sort=False)[source]¶

Print date, time, and version information.

Use scooby to print date, time, and package version information in any environment (Jupyter notebook, IPython console, Python console, QT console), either as html-table (notebook) or as plain text (anywhere).

Always shown are the OS, number of CPU(s), numpy, scipy, emg3d, numba, sys.version, and time/date.

Additionally shown are, if they can be imported, IPython and matplotlib. It also shows MKL information, if available.

All modules provided in add_pckg are also shown.

Note

The package scooby has to be installed in order to use Report: pip install scooby.

Parameters:	add_pckg : packages, optional Package or list of packages to add to output information (must be imported beforehand). ncol : int, optional Number of package-columns in html table (no effect in text-version); Defaults to 3. text_width : int, optional The text width for non-HTML display modes sort : bool, optional Sort the packages when the report is shown

Examples

>>> import pytest
>>> import dateutil
>>> from emg3d import Report
>>> Report()                            # Default values
>>> Report(pytest)                      # Provide additional package
>>> Report([pytest, dateutil], ncol=5)  # Set nr of columns

class emg3d.utils.Field[source]¶

Create a Field instance with x-, y-, and z-views of the field.

A Field is an ndarray with additional views of the x-, y-, and z-directed fields as attributes, stored as fx, fy, and fz. The default array contains the whole field, which can be the electric field, the source field, or the residual field, in a 1D array. A Field instance has additionally the property ensure_pec which, if called, ensures Perfect Electric Conductor (PEC) boundary condition. It also has the two attributes amp and pha for the amplitude and phase, as common in frequency-domain CSEM.

A Field can be initiated in three ways:

Field(grid, dtype=complex): Calling it with a TensorMesh instance returns a Field instance of correct dimensions initiated with zeroes of data type dtype.
Field(grid, field): Calling it with a TensorMesh instance and an ndarray returns a Field instance of the provided ndarray, of same data type.
Field(fx, fy, fz): Calling it with three ndarray’s which represent the field in x-, y-, and z-direction returns a Field instance with these views, of same data type.

Sort-order is ‘F’.

Parameters:

fx_or_grid : TensorMesh or ndarray

Either a TensorMesh instance or an ndarray of shape grid.nEx or grid.vnEx. See explanations above. Only mandatory parameter; if the only one provided, it will initiate a zero-field of dtype.

fy_or_field : Field or ndarray, optional

Either a Field instance or an ndarray of shape grid.nEy or grid.vnEy. See explanations above.

fz : ndarray, optional

An ndarray of shape grid.nEz or grid.vnEz. See explanations above.

dtype : dtype, optional

Only used if fy_or_field=None and fz=None; the initiated zero-field for the provided TensorMesh has data type dtype. Default: complex.

freq : float, optional

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

Just added as info if provided.

amp(self)[source]¶: Amplitude of the electromagnetic field.

copy(self)[source]¶: Return a copy of the Field.

ensure_pec¶: Set Perfect Electric Conductor (PEC) boundary condition.

field¶: Entire field, 1D [fx, fy, fz].

freq¶: Return frequency.

classmethod from_dict(inp)[source]¶

Convert dictionary into Field instance.

Parameters:	inp : dict Dictionary as obtained from `Field.to_dict()`. The dictionary needs the keys field, freq, vnEx, vnEy, and vnEz.
Returns:	obj : `Field` instance

fx¶: View of the x-directed field in the x-direction (nCx, nNy, nNz).

fy¶: View of the field in the y-direction (nNx, nCy, nNz).

fz¶: View of the field in the z-direction (nNx, nNy, nCz).

is_electric¶: Returns True if Field is electric, False if it is magnetic.

pha(self, deg=False, unwrap=True, lag=True)[source]¶

Phase of the electromagnetic field.

Parameters:	deg : bool If True the returned phase is in degrees, else in radians. Default is False (radians). unwrap : bool If True the returned phase is unwrapped. Default is True (unwrapped). lag : bool If True the returned phase is lag, else lead defined. Default is True (lag defined).

smu0¶: Return s*mu_0; mu_0 = Magn. permeability of free space [H/m].

sval¶: Return s; s=iw in frequency domain; s=freq in Laplace domain.

to_dict(self, copy=False)[source]¶: Store the necessary information of the Field in a dict.

class emg3d.utils.SourceField[source]¶

Create a Source-Field instance with x-, y-, and z-views of the field.

A subclass of Field. Additional properties are the real-valued source vector (vector, vx, vy, vz), which sum is always one. For a SourceField frequency is a mandatory parameter, unlike for a Field (recommended also for Field though),

Parameters:

fx_or_grid : TensorMesh or ndarray

Either a TensorMesh instance or an ndarray of shape grid.nEx or grid.vnEx. See explanations above. Only mandatory parameter; if the only one provided, it will initiate a zero-field of dtype.

fy_or_field : Field or ndarray, optional

Either a Field instance or an ndarray of shape grid.nEy or grid.vnEy. See explanations above.

fz : ndarray, optional

An ndarray of shape grid.nEz or grid.vnEz. See explanations above.

dtype : dtype, optional

Only used if fy_or_field=None and fz=None; the initiated zero-field for the provided TensorMesh has data type dtype. Default: complex.

freq : float

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

In difference to Field, the frequency has to be provided for a SourceField.

copy(self)[source]¶: Return a copy of the SourceField.

classmethod from_dict(inp)[source]¶

Convert dictionary into SourceField instance.

Parameters:	inp : dict Dictionary as obtained from `SourceField.to_dict()`. The dictionary needs the keys field, freq, vnEx, vnEy, and vnEz.
Returns:	obj : `SourceField` instance

vector¶: Entire vector, 1D [vx, vy, vz].

vx¶: View of the x-directed vector in the x-direction (nCx, nNy, nNz).

vy¶: View of the vector in the y-direction (nNx, nCy, nNz).

vz¶: View of the vector in the z-direction (nNx, nNy, nCz).

emg3d.utils.get_source_field(grid, src, freq, strength=0)[source]¶

Return the source field.

The source field is given in Equation 2 in [Muld06],

\[s \mu_0 \mathbf{J}_\mathrm{s} ,\]

where \(s = \mathrm{i} \omega\). Either finite length dipoles or infinitesimal small point dipoles can be defined, whereas the return source field corresponds to a normalized (1 Am) source distributed within the cell(s) it resides (can be changed with the strength-parameter).

The adjoint of the trilinear interpolation is used to distribute the point(s) to the grid edges, which corresponds to the discretization of a Dirac ([PlDM07]).

Parameters:

grid : TensorMesh

Model grid; a TensorMesh instance.

src : list of floats

Source coordinates (m). There are two formats:

Finite length dipole: [x0, x1, y0, y1, z0, z1].

Point dipole: [x, y, z, azimuth, dip].

freq : float

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

strength : float or complex, optional

Source strength (A):

If 0, output is normalized to a source of 1 m length, and source strength of 1 A.

If != 0, output is returned for given source length and strength.

Default is 0.

Returns:

sfield : SourceField() instance: Source field, normalized to 1 A m.

emg3d.utils.get_receiver(grid, values, coordinates, method='cubic', extrapolate=False)[source]¶

Return values corresponding to grid at coordinates.

Works for electric fields as well as magnetic fields obtained with get_h_field(), and for model parameters.

Parameters:

grid : TensorMesh

Model grid; a TensorMesh instance.

values : ndarray

Field instance, or a particular field (e.g. field.fx); Model parameters.

coordinates : tuple (x, y, z)

Coordinates (x, y, z) where to interpolate values; e.g. receiver locations.

method : str, optional

The method of interpolation to perform, ‘linear’ or ‘cubic’. Default is ‘cubic’ (forced to ‘linear’ if there are less than 3 points in any direction).

extrapolate : bool

If True, points on new_grid which are outside of grid are filled by the nearest value (if method='cubic') or by extrapolation (if method='linear'). If False, points outside are set to zero.

Default is False.

Returns:

new_values : ndarray or empymod.utils.EMArray

Values at coordinates.

If input was a field it returns an EMArray, which is a subclassed ndarray with .pha and .amp attributes.

If input was an entire Field instance, output is a tuple (fx, fy, fz).

See also

grid2grid: Interpolation of model parameters or fields to a new grid.

emg3d.utils.get_h_field(grid, model, field)[source]¶

Return magnetic field corresponding to provided electric field.

Retrieve the magnetic field \(\mathbf{H}\) from the electric field \(\mathbf{E}\) using Farady’s law, given by

\[\nabla \times \mathbf{E} = \rm{i}\omega\mu\mathbf{H} .\]

Note that the magnetic field in x-direction is defined in the center of the face defined by the electric field in y- and z-directions, and similar for the other field directions. This means that the provided electric field and the returned magnetic field have different dimensions:

E-field:  x: [grid.vectorCCx,  grid.vectorNy,  grid.vectorNz]
          y: [ grid.vectorNx, grid.vectorCCy,  grid.vectorNz]
          z: [ grid.vectorNx,  grid.vectorNy, grid.vectorCCz]

H-field:  x: [ grid.vectorNx, grid.vectorCCy, grid.vectorCCz]
          y: [grid.vectorCCx,  grid.vectorNy, grid.vectorCCz]
          z: [grid.vectorCCx, grid.vectorCCy,  grid.vectorNz]

Parameters:	grid : TensorMesh Model grid; `TensorMesh` instance. model : Model Model; `Model` instance. field : Field Electric field; `Field` instance.
Returns:	hfield : Field Magnetic field; `Field` instance.

class emg3d.utils.Model(grid, res_x=1.0, res_y=None, res_z=None, mu_r=None, epsilon_r=None)[source]¶

Create a model instance.

Class to provide model parameters (x-, y-, and z-directed resistivities, electric permittivity and magnetic permeability) to the solver. Relative magnetic permeability \(\mu_\mathrm{r}\) is by default set to one and electric permittivity \(\varepsilon_\mathrm{r}\) is by default set to zero, but they can also be provided (isotropically). Keep in mind that the multigrid method as implemented in emg3d only works for the diffusive approximation. As soon as the displacement-part in the Maxwell’s equations becomes too dominant it will fail (high frequencies or very high electric permittivity).

Parameters:

grid : TensorMesh: Grid on which to apply model.
res_x, res_y, res_z : float or ndarray; default to 1.: Resistivity in x-, y-, and z-directions. If ndarray, they must have the shape of grid.vnC (F-ordered) or grid.nC. Resistivities have to be bigger than zero and smaller than infinity.
mu_r : None, float, or ndarray: Relative magnetic permeability (isotropic). If ndarray it must have the shape of grid.vnC (F-ordered) or grid.nC. Default is None, which corresponds to 1., but avoids the calculation of zeta. Magnetic permeability has to be bigger than zero and smaller than infinity.
epsilon_r : None, float, or ndarray: Relative electric permittivity (isotropic). If ndarray it must have the shape of grid.vnC (F-ordered) or grid.nC. The displacement part is completely neglected (diffusive approximation) if set to None, which is the default. Electric permittivity has to be bigger than zero and smaller than infinity.

copy(self)[source]¶: Return a copy of the Model.

epsilon_r¶: Electric permittivity.

classmethod from_dict(inp)[source]¶

Convert the dictionary into a Model instance.

Parameters:	inp : dict Dictionary as obtained from `Model.to_dict()`. The dictionary needs the keys res_x, res_y, res_z, mu_r, epsilon_r, and vnC.
Returns:	obj : `Model` instance

mu_r¶: Magnetic permeability.

res_x¶: Resistivity in x-direction.

res_y¶: Resistivity in y-direction.

res_z¶: Resistivity in z-direction.

to_dict(self, copy=False)[source]¶: Store the necessary information of the Model in a dict.

class emg3d.utils.VolumeModel(grid, model, sfield)[source]¶

Return a volume-averaged version of provided model.

Takes a Model instance and returns the volume averaged values. This is used by the solver internally.

\[\eta_{\{x,y,z\}} = -V\mathrm{i}\omega\mu_0 \left(\rho^{-1}_{\{x,y,z\}} + \mathrm{i}\omega\varepsilon\right)\]

\[\zeta = V\mu_\mathrm{r}^{-1}\]

Parameters:	grid : TensorMesh Grid on which to apply model. model : Model Model to transform to volume-averaged values. sfield : SourceField A VolumeModel is frequency-dependent. The frequency-information is taken from the provided source filed.

static calculate_eta(name, grid, model, field)[source]¶: eta: volume divided by resistivity.

static calculate_zeta(name, grid, model)[source]¶: zeta: volume divided by mu_r.

eta_x¶: eta in x-direction.

eta_y¶: eta in y-direction.

eta_z¶: eta in z-direction.

zeta¶: zeta.

emg3d.utils.grid2grid(grid, values, new_grid, method='linear', extrapolate=True, log=False)[source]¶

Interpolate values located on grid to new_grid.

Note 1: The default method is ‘linear’, because it works with fields and model parameters. However, recommended are ‘volume’ for model parameters and ‘cubic’ for fields.

Note 2: For model parameters with method=’volume’ the result is quite different if you provide resistivity, conductivity, or the logarithm of any of the two. The recommended way is to provide the logarithm of resistivity or conductivity, in which case the output of one is indeed the inverse of the output of the other.

Parameters:

grid, new_grid : TensorMesh

Input and output model grids; TensorMesh instances.

values : ndarray

Model parameters; emg3d.fields.Field instance, or a particular field (e.g. field.fx). For fields the method cannot be ‘volume’.

method : {<’linear’>, ‘volume’, ‘cubic’}, optional

The method of interpolation to perform. The volume averaging method ensures that the total sum of the property stays constant.

Volume averaging is only implemented for model parameters, not for fields. The method ‘cubic’ requires at least three points in any direction, otherwise it will fall back to ‘linear’.

Default is ‘linear’, because it works with fields and model parameters. However, recommended are ‘volume’ for model parameters and ‘cubic’ for fields.

extrapolate : bool

If True, points on new_grid which are outside of grid are filled by the nearest value (if method='cubic') or by extrapolation (if method='linear'). If False, points outside are set to zero.

For method='volume' it always uses the nearest value for points outside of grid.

Default is True.

log : bool

If True, the interpolation is carried out on a log10-scale; hence the same as 10**grid2grid(grid, np.log10(values), ...). Default is False.

Returns:

new_values : ndarray: Values corresponding to new_grid.

See also

get_receiver: Interpolation of model parameters or fields to (x, y, z).

emg3d.utils.interp3d(points, values, new_points, method, fill_value, mode)[source]¶

Interpolate values in 3D either linearly or with a cubic spline.

Return values corresponding to a regular 3D grid defined by points on new_points.

This is a modified version of scipy.interpolate.interpn(), using scipy.interpolate.RegularGridInterpolator if method='linear' and a custom-wrapped version of scipy.ndimage.map_coordinates() if method='cubic'. If speed is important then choose ‘linear’, as it can be significantly faster.

Parameters:

points : tuple of ndarray of float, with shapes ((nx, ), (ny, ) (nz, )): The points defining the regular grid in three dimensions.
values : array_like, shape (nx, ny, nz): The data on the regular grid in three dimensions.
new_points : tuple (rec_x, rec_y, rec_z): Coordinates (x, y, z) of new points.
method : {‘cubic’, ‘linear’}, optional: The method of interpolation to perform, ‘linear’ or ‘cubic’. Default is ‘cubic’ (forced to ‘linear’ if there are less than 3 points in any direction).
fill_value : float or None: Passed to scipy.interpolate.RegularGridInterpolator if method='linear': The value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated.
mode : {‘constant’, ‘nearest’, ‘mirror’, ‘reflect’, ‘wrap’}: Passed to scipy.ndimage.map_coordinates() if method='cubic': Determines how the input array is extended beyond its boundaries.

Returns:

new_values : ndarray: Values corresponding to new_points.

class emg3d.utils.TensorMesh(h, x0)[source]¶

Rudimentary mesh for multigrid calculation.

The tensor-mesh discretize.TensorMesh is a powerful tool, including sophisticated mesh-generation possibilities in 1D, 2D, and 3D, plotting routines, and much more. However, in the multigrid solver we have to generate a mesh at each level, many times over and over again, and we only need a very limited set of attributes. This tensor-mesh class provides all required attributes. All attributes here are the same as their counterparts in discretize.TensorMesh (both in name and value).

Warning

This is a slimmed-down version of discretize.TensorMesh, meant principally for internal use by the multigrid modeller. It is highly recommended to use discretize.TensorMesh to create the input meshes instead of this class. There are no input-checks carried out here, and there is only one accepted input format for h and x0.

Parameters:	h : list of three ndarrays Cell widths in [x, y, z] directions. x0 : ndarray of dimension (3, ) Origin (x, y, z).

copy(self)[source]¶: Return a copy of the TensorMesh.

classmethod from_dict(inp)[source]¶

Convert dictionary into TensorMesh instance.

Parameters:	inp : dict Dictionary as obtained from `TensorMesh.to_dict()`. The dictionary needs the keys hx, hy, hz, and x0.
Returns:	obj : `TensorMesh` instance

to_dict(self, copy=False)[source]¶: Store the necessary information of the TensorMesh in a dict.

vol¶: Construct cell volumes of the 3D model as 1D array.

emg3d.utils.get_hx_h0(freq, res, domain, fixed=0.0, possible_nx=None, min_width=None, pps=3, alpha=None, max_domain=100000.0, raise_error=True, verb=1, return_info=False)[source]¶

Return cell widths and origin for given parameters.

Returns cell widths for the provided frequency, resistivity, domain extent, and other parameters using a flexible amount of cells. See input parameters for more details. A maximum of three hard/fixed boundaries can be provided (one of which is the grid center).

The minimum cell width is calculated through \(\delta/\rm{pps}\), where the skin depth is given by \(\delta = 503.3 \sqrt{\rho/f}\), and the parameter pps stands for ‘points-per-skindepth’. The minimum cell width can be restricted with the parameter min_width.

The actual calculation domain adds a buffer zone around the (survey) domain. The thickness of the buffer is six times the skin depth. The field is basically zero after two wavelengths. A wavelength is \(2\pi\delta\), hence roughly 6 times the skin depth. Taking a factor 6 gives therefore almost two wavelengths, as the field travels to the boundary and back. The actual buffer thickness can be steered with the res parameter.

One has to take into account that the air is very resistive, which has to be considered not just in the vertical direction, but also in the horizontal directions, as the airwave will bounce back from the sides otherwise. In the marine case this issue reduces with increasing water depth.

Parameters:

freq : float

Frequency (Hz) to calculate the skin depth. The skin depth is a concept defined in the frequency domain. If a negative frequency is provided, it is assumed that the calculation is carried out in the Laplace domain. To calculate the skin depth, the value of freq is then multiplied by \(-2\pi\), to simulate the closest frequency-equivalent.

res : float or list

Resistivity (Ohm m) to calculate the skin depth. The skin depth is used to calculate the minimum cell width and the boundary thicknesses. Up to three resistivities can be provided:

float: Same resistivity for everything;
[min_width, boundaries];
[min_width, left boundary, right boundary].

domain : list

Contains the survey-domain limits [min, max]. The actual calculation domain consists of this domain plus a buffer zone around it, which depends on frequency and resistivity.

fixed : list, optional

Fixed boundaries, one, two, or maximum three values. The grid is centered around the first value. Hence it is the center location with the smallest cell. Two more fixed boundaries can be added, at most one on each side of the first one. Default is 0.

possible_nx : list, optional

List of possible numbers of cells. See get_cell_numbers(). Default is get_cell_numbers(500, 5, 3), which corresponds to [16, 24, 32, 40, 48, 64, 80, 96, 128, 160, 192, 256, 320, 384].

min_width : float, list or None, optional

Minimum cell width restriction:

None : No restriction;
float : Fixed to this value, ignoring skin depth and pps.
list [min, max] : Lower and upper bounds.

Default is None.

pps : int, optional

Points per skindepth; minimum cell width is calculated via dmin = skindepth/pps. Default = 3.

alpha : list, optional

Maximum alpha and step size to find a good alpha. The first value is the maximum alpha of the survey domain, the second value is the maximum alpha for the buffer zone, and the third value is the step size. Default = [1, 1.5, .01], hence no stretching within the survey domain and a maximum stretching of 1.5 in the buffer zone; step size is 0.01.

max_domain : float, optional

Maximum calculation domain from fixed[0] (usually source position). Default is 100,000.

raise_error : bool, optional

If True, an error is raised if no suitable grid is found. Otherwise it just prints a message and returns None’s. Default is True.

verb : int, optional

Verbosity, 0 or 1. Default = 1.

return_info : bool

If True, a dictionary is returned with some grid info (min and max cell width and alpha).

Returns:

hx : ndarray

Cell widths of mesh.

x0 : float

Origin of the mesh.

info : dict

Dictionary with mesh info; only if return_info=True.

Keys:

dmin: Minimum cell width;
dmax: Maximum cell width;
amin: Minimum alpha;
amax: Maximum alpha.

See also

get_stretched_h: Get hx for a fixed number nx and within a fixed domain.

emg3d.utils.get_cell_numbers(max_nr, max_prime=5, min_div=3)[source]¶

Returns ‘good’ cell numbers for the multigrid method.

‘Good’ cell numbers are numbers which can be divided by 2 as many times as possible. At the end there will be a low prime number.

The function adds all numbers \(p 2^n \leq M\) for \(p={2, 3, ..., p_\text{max}}\) and \(n={n_\text{min}, n_\text{min}+1, ..., \infty}\); \(M, p_\text{max}, n_\text{min}\) correspond to max_nr, max_prime, and min_div, respectively.

Parameters:	max_nr : int Maximum number of cells. max_prime : int Highest permitted prime number p for p2^n. {2, 3, 5, 7} are good upper limits in order to avoid too big lowest grids in the multigrid method. Default is 5. min_div* : int Minimum times the number can be divided by two. Default is 3.
Returns:	numbers : array Array containing all possible cell numbers from lowest to highest.

emg3d.utils.get_stretched_h(min_width, domain, nx, x0=0, x1=None, resp_domain=False)[source]¶

Return cell widths for a stretched grid within the domain.

Returns nx cell widths within domain, where the minimum cell width is min_width. The cells are not stretched within x0 and x1, and outside uses a power-law stretching. The actual stretching factor and the number of cells left and right of x0 and x1 are find in a minimization process.

The domain is not completely respected. The starting point of the domain is, but the endpoint of the domain might slightly shift (this is more likely the case for small nx, for big nx the shift should be small). The new endpoint can be obtained with domain[0]+np.sum(hx). If you want the domain to be respected absolutely, set resp_domain=True. However, be aware that this will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Parameters:

min_width : float: Minimum cell width. If x1 is provided, the actual minimum cell width might be smaller than min_width.
domain : list: [start, end] of model domain.
nx : int: Number of cells.
x0 : float: Center of the grid. x0 is restricted to domain. Default is 0.
x1 : float: If provided, then no stretching is applied between x0 and x1. The non-stretched part starts at x0 and stops at the first possible location at or after x1. x1 is restricted to domain. This will min_width so that an integer number of cells fit within x0 and x1.
resp_domain : bool: If False (default), then the domain-end might shift slightly to assure that the same stretching factor is applied throughout. If set to True, however, the domain is respected absolutely. This will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Returns:

hx : ndarray: Cell widths of mesh.

See also

get_hx_x0: Get hx and x0 for a flexible number of nx with given bounds.

emg3d.utils.get_domain(x0=0, freq=1, res=0.3, limits=None, min_width=None, fact_min=0.2, fact_neg=5, fact_pos=None)[source]¶

Get domain extent and minimum cell width as a function of skin depth.

Returns the extent of the calculation domain and the minimum cell width as a multiple of the skin depth, with possible user restrictions on minimum calculation domain and range of possible minimum cell widths.

\[\begin{split}\delta &= 503.3 \sqrt{\frac{\rho}{f}} , \\ x_\text{start} &= x_0-k_\text{neg}\delta , \\ x_\text{end} &= x_0+k_\text{pos}\delta , \\ h_\text{min} &= k_\text{min} \delta .\end{split}\]

Parameters:

x0 : float

Center of the calculation domain. Normally the source location. Default is 0.

freq : float

Frequency (Hz) to calculate the skin depth. The skin depth is a concept defined in the frequency domain. If a negative frequency is provided, it is assumed that the calculation is carried out in the Laplace domain. To calculate the skin depth, the value of freq is then multiplied by \(-2\pi\), to simulate the closest frequency-equivalent.

Default is 1 Hz.

res : float, optional

Resistivity (Ohm m) to calculate skin depth. Default is 0.3 Ohm m (sea water).

limits : None or list

[start, end] of model domain. This extent represents the minimum extent of the domain. The domain is therefore only adjusted if it has to reach outside of [start, end]. Default is None.

min_width : None, float, or list of two floats

Minimum cell width is calculated as a function of skin depth: fact_min*sd. If min_width is a float, this is used. If a list of two values [min, max] are provided, they are used to restrain min_width. Default is None.

fact_min, fact_neg, fact_pos : floats

The skin depth is multiplied with these factors to estimate:

Minimum cell width (fact_min, default 0.2)

Domain-start (fact_neg, default 5), and

Domain-end (fact_pos, defaults to fact_neg).

Returns:

h_min : float: Minimum cell width.
domain : list: Start- and end-points of calculation domain.

emg3d.utils.get_hx(alpha, domain, nx, x0, resp_domain=True)[source]¶

Return cell widths for given input.

Find the number of cells left and right of x0, nl and nr respectively, for the provided alpha. For this, we solve

\[\frac{x_\text{max}-x_0}{x_0-x_\text{min}} = \frac{a^{nr}-1}{a^{nl}-1}\]

where \(a = 1+\alpha\).

Parameters:

alpha : float: Stretching factor a is given by a=1+alpha.
domain : list: [start, end] of model domain.
nx : int: Number of cells.
x0 : float: Center of the grid. x0 is restricted to domain.
resp_domain : bool: If False (default), then the domain-end might shift slightly to assure that the same stretching factor is applied throughout. If set to True, however, the domain is respected absolutely. This will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Returns:

hx : ndarray: Cell widths of mesh.

`meshes` – Discretization¶

Everything related to meshes appropriate for the multigrid solver.

class emg3d.meshes.TensorMesh(h, x0)[source]¶

Rudimentary mesh for multigrid calculation.

The tensor-mesh discretize.TensorMesh is a powerful tool, including sophisticated mesh-generation possibilities in 1D, 2D, and 3D, plotting routines, and much more. However, in the multigrid solver we have to generate a mesh at each level, many times over and over again, and we only need a very limited set of attributes. This tensor-mesh class provides all required attributes. All attributes here are the same as their counterparts in discretize.TensorMesh (both in name and value).

Warning

This is a slimmed-down version of discretize.TensorMesh, meant principally for internal use by the multigrid modeller. It is highly recommended to use discretize.TensorMesh to create the input meshes instead of this class. There are no input-checks carried out here, and there is only one accepted input format for h and x0.

Parameters:	h : list of three ndarrays Cell widths in [x, y, z] directions. x0 : ndarray of dimension (3, ) Origin (x, y, z).

copy(self)[source]¶: Return a copy of the TensorMesh.

classmethod from_dict(inp)[source]¶

Convert dictionary into TensorMesh instance.

Parameters:	inp : dict Dictionary as obtained from `TensorMesh.to_dict()`. The dictionary needs the keys hx, hy, hz, and x0.
Returns:	obj : `TensorMesh` instance

to_dict(self, copy=False)[source]¶: Store the necessary information of the TensorMesh in a dict.

vol¶: Construct cell volumes of the 3D model as 1D array.

emg3d.meshes.get_hx_h0(freq, res, domain, fixed=0.0, possible_nx=None, min_width=None, pps=3, alpha=None, max_domain=100000.0, raise_error=True, verb=1, return_info=False)[source]¶

Return cell widths and origin for given parameters.

Returns cell widths for the provided frequency, resistivity, domain extent, and other parameters using a flexible amount of cells. See input parameters for more details. A maximum of three hard/fixed boundaries can be provided (one of which is the grid center).

The minimum cell width is calculated through \(\delta/\rm{pps}\), where the skin depth is given by \(\delta = 503.3 \sqrt{\rho/f}\), and the parameter pps stands for ‘points-per-skindepth’. The minimum cell width can be restricted with the parameter min_width.

The actual calculation domain adds a buffer zone around the (survey) domain. The thickness of the buffer is six times the skin depth. The field is basically zero after two wavelengths. A wavelength is \(2\pi\delta\), hence roughly 6 times the skin depth. Taking a factor 6 gives therefore almost two wavelengths, as the field travels to the boundary and back. The actual buffer thickness can be steered with the res parameter.

One has to take into account that the air is very resistive, which has to be considered not just in the vertical direction, but also in the horizontal directions, as the airwave will bounce back from the sides otherwise. In the marine case this issue reduces with increasing water depth.

Parameters:

freq : float

Frequency (Hz) to calculate the skin depth. The skin depth is a concept defined in the frequency domain. If a negative frequency is provided, it is assumed that the calculation is carried out in the Laplace domain. To calculate the skin depth, the value of freq is then multiplied by \(-2\pi\), to simulate the closest frequency-equivalent.

res : float or list

Resistivity (Ohm m) to calculate the skin depth. The skin depth is used to calculate the minimum cell width and the boundary thicknesses. Up to three resistivities can be provided:

float: Same resistivity for everything;
[min_width, boundaries];
[min_width, left boundary, right boundary].

domain : list

Contains the survey-domain limits [min, max]. The actual calculation domain consists of this domain plus a buffer zone around it, which depends on frequency and resistivity.

fixed : list, optional

Fixed boundaries, one, two, or maximum three values. The grid is centered around the first value. Hence it is the center location with the smallest cell. Two more fixed boundaries can be added, at most one on each side of the first one. Default is 0.

possible_nx : list, optional

List of possible numbers of cells. See get_cell_numbers(). Default is get_cell_numbers(500, 5, 3), which corresponds to [16, 24, 32, 40, 48, 64, 80, 96, 128, 160, 192, 256, 320, 384].

min_width : float, list or None, optional

Minimum cell width restriction:

None : No restriction;
float : Fixed to this value, ignoring skin depth and pps.
list [min, max] : Lower and upper bounds.

Default is None.

pps : int, optional

Points per skindepth; minimum cell width is calculated via dmin = skindepth/pps. Default = 3.

alpha : list, optional

Maximum alpha and step size to find a good alpha. The first value is the maximum alpha of the survey domain, the second value is the maximum alpha for the buffer zone, and the third value is the step size. Default = [1, 1.5, .01], hence no stretching within the survey domain and a maximum stretching of 1.5 in the buffer zone; step size is 0.01.

max_domain : float, optional

Maximum calculation domain from fixed[0] (usually source position). Default is 100,000.

raise_error : bool, optional

If True, an error is raised if no suitable grid is found. Otherwise it just prints a message and returns None’s. Default is True.

verb : int, optional

Verbosity, 0 or 1. Default = 1.

return_info : bool

If True, a dictionary is returned with some grid info (min and max cell width and alpha).

Returns:

hx : ndarray

Cell widths of mesh.

x0 : float

Origin of the mesh.

info : dict

Dictionary with mesh info; only if return_info=True.

Keys:

dmin: Minimum cell width;
dmax: Maximum cell width;
amin: Minimum alpha;
amax: Maximum alpha.

See also

get_stretched_h: Get hx for a fixed number nx and within a fixed domain.

emg3d.meshes.get_cell_numbers(max_nr, max_prime=5, min_div=3)[source]¶

Returns ‘good’ cell numbers for the multigrid method.

‘Good’ cell numbers are numbers which can be divided by 2 as many times as possible. At the end there will be a low prime number.

The function adds all numbers \(p 2^n \leq M\) for \(p={2, 3, ..., p_\text{max}}\) and \(n={n_\text{min}, n_\text{min}+1, ..., \infty}\); \(M, p_\text{max}, n_\text{min}\) correspond to max_nr, max_prime, and min_div, respectively.

Parameters:	max_nr : int Maximum number of cells. max_prime : int Highest permitted prime number p for p2^n. {2, 3, 5, 7} are good upper limits in order to avoid too big lowest grids in the multigrid method. Default is 5. min_div* : int Minimum times the number can be divided by two. Default is 3.
Returns:	numbers : array Array containing all possible cell numbers from lowest to highest.

emg3d.meshes.get_stretched_h(min_width, domain, nx, x0=0, x1=None, resp_domain=False)[source]¶

Return cell widths for a stretched grid within the domain.

Returns nx cell widths within domain, where the minimum cell width is min_width. The cells are not stretched within x0 and x1, and outside uses a power-law stretching. The actual stretching factor and the number of cells left and right of x0 and x1 are find in a minimization process.

The domain is not completely respected. The starting point of the domain is, but the endpoint of the domain might slightly shift (this is more likely the case for small nx, for big nx the shift should be small). The new endpoint can be obtained with domain[0]+np.sum(hx). If you want the domain to be respected absolutely, set resp_domain=True. However, be aware that this will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Parameters:

min_width : float: Minimum cell width. If x1 is provided, the actual minimum cell width might be smaller than min_width.
domain : list: [start, end] of model domain.
nx : int: Number of cells.
x0 : float: Center of the grid. x0 is restricted to domain. Default is 0.
x1 : float: If provided, then no stretching is applied between x0 and x1. The non-stretched part starts at x0 and stops at the first possible location at or after x1. x1 is restricted to domain. This will min_width so that an integer number of cells fit within x0 and x1.
resp_domain : bool: If False (default), then the domain-end might shift slightly to assure that the same stretching factor is applied throughout. If set to True, however, the domain is respected absolutely. This will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Returns:

hx : ndarray: Cell widths of mesh.

See also

get_hx_x0: Get hx and x0 for a flexible number of nx with given bounds.

emg3d.meshes.get_domain(x0=0, freq=1, res=0.3, limits=None, min_width=None, fact_min=0.2, fact_neg=5, fact_pos=None)[source]¶

Get domain extent and minimum cell width as a function of skin depth.

Returns the extent of the calculation domain and the minimum cell width as a multiple of the skin depth, with possible user restrictions on minimum calculation domain and range of possible minimum cell widths.

\[\begin{split}\delta &= 503.3 \sqrt{\frac{\rho}{f}} , \\ x_\text{start} &= x_0-k_\text{neg}\delta , \\ x_\text{end} &= x_0+k_\text{pos}\delta , \\ h_\text{min} &= k_\text{min} \delta .\end{split}\]

Parameters:

x0 : float

Center of the calculation domain. Normally the source location. Default is 0.

freq : float

Frequency (Hz) to calculate the skin depth. The skin depth is a concept defined in the frequency domain. If a negative frequency is provided, it is assumed that the calculation is carried out in the Laplace domain. To calculate the skin depth, the value of freq is then multiplied by \(-2\pi\), to simulate the closest frequency-equivalent.

Default is 1 Hz.

res : float, optional

Resistivity (Ohm m) to calculate skin depth. Default is 0.3 Ohm m (sea water).

limits : None or list

[start, end] of model domain. This extent represents the minimum extent of the domain. The domain is therefore only adjusted if it has to reach outside of [start, end]. Default is None.

min_width : None, float, or list of two floats

Minimum cell width is calculated as a function of skin depth: fact_min*sd. If min_width is a float, this is used. If a list of two values [min, max] are provided, they are used to restrain min_width. Default is None.

fact_min, fact_neg, fact_pos : floats

The skin depth is multiplied with these factors to estimate:

Minimum cell width (fact_min, default 0.2)

Domain-start (fact_neg, default 5), and

Domain-end (fact_pos, defaults to fact_neg).

Returns:

h_min : float: Minimum cell width.
domain : list: Start- and end-points of calculation domain.

emg3d.meshes.get_hx(alpha, domain, nx, x0, resp_domain=True)[source]¶

Return cell widths for given input.

Find the number of cells left and right of x0, nl and nr respectively, for the provided alpha. For this, we solve

\[\frac{x_\text{max}-x_0}{x_0-x_\text{min}} = \frac{a^{nr}-1}{a^{nl}-1}\]

where \(a = 1+\alpha\).

Parameters:

alpha : float: Stretching factor a is given by a=1+alpha.
domain : list: [start, end] of model domain.
nx : int: Number of cells.
x0 : float: Center of the grid. x0 is restricted to domain.
resp_domain : bool: If False (default), then the domain-end might shift slightly to assure that the same stretching factor is applied throughout. If set to True, however, the domain is respected absolutely. This will introduce one stretch-factor which is different from the other stretch factors, to accommodate the restriction. This one-off factor is between the left- and right-side of x0, or, if x1 is provided, just after x1.

Returns:

hx : ndarray: Cell widths of mesh.

`models` – Earth properties¶

Everything to create model-properties for the multigrid solver.

class emg3d.models.Model(grid, res_x=1.0, res_y=None, res_z=None, mu_r=None, epsilon_r=None)[source]¶

Create a model instance.

Class to provide model parameters (x-, y-, and z-directed resistivities, electric permittivity and magnetic permeability) to the solver. Relative magnetic permeability \(\mu_\mathrm{r}\) is by default set to one and electric permittivity \(\varepsilon_\mathrm{r}\) is by default set to zero, but they can also be provided (isotropically). Keep in mind that the multigrid method as implemented in emg3d only works for the diffusive approximation. As soon as the displacement-part in the Maxwell’s equations becomes too dominant it will fail (high frequencies or very high electric permittivity).

Parameters:

grid : TensorMesh: Grid on which to apply model.
res_x, res_y, res_z : float or ndarray; default to 1.: Resistivity in x-, y-, and z-directions. If ndarray, they must have the shape of grid.vnC (F-ordered) or grid.nC. Resistivities have to be bigger than zero and smaller than infinity.
mu_r : None, float, or ndarray: Relative magnetic permeability (isotropic). If ndarray it must have the shape of grid.vnC (F-ordered) or grid.nC. Default is None, which corresponds to 1., but avoids the calculation of zeta. Magnetic permeability has to be bigger than zero and smaller than infinity.
epsilon_r : None, float, or ndarray: Relative electric permittivity (isotropic). If ndarray it must have the shape of grid.vnC (F-ordered) or grid.nC. The displacement part is completely neglected (diffusive approximation) if set to None, which is the default. Electric permittivity has to be bigger than zero and smaller than infinity.

copy(self)[source]¶: Return a copy of the Model.

epsilon_r¶: Electric permittivity.

classmethod from_dict(inp)[source]¶

Convert the dictionary into a Model instance.

Parameters:	inp : dict Dictionary as obtained from `Model.to_dict()`. The dictionary needs the keys res_x, res_y, res_z, mu_r, epsilon_r, and vnC.
Returns:	obj : `Model` instance

mu_r¶: Magnetic permeability.

res_x¶: Resistivity in x-direction.

res_y¶: Resistivity in y-direction.

res_z¶: Resistivity in z-direction.

to_dict(self, copy=False)[source]¶: Store the necessary information of the Model in a dict.

class emg3d.models.VolumeModel(grid, model, sfield)[source]¶

Return a volume-averaged version of provided model.

Takes a Model instance and returns the volume averaged values. This is used by the solver internally.

\[\eta_{\{x,y,z\}} = -V\mathrm{i}\omega\mu_0 \left(\rho^{-1}_{\{x,y,z\}} + \mathrm{i}\omega\varepsilon\right)\]

\[\zeta = V\mu_\mathrm{r}^{-1}\]

Parameters:	grid : TensorMesh Grid on which to apply model. model : Model Model to transform to volume-averaged values. sfield : SourceField A VolumeModel is frequency-dependent. The frequency-information is taken from the provided source filed.

static calculate_eta(name, grid, model, field)[source]¶: eta: volume divided by resistivity.

static calculate_zeta(name, grid, model)[source]¶: zeta: volume divided by mu_r.

eta_x¶: eta in x-direction.

eta_y¶: eta in y-direction.

eta_z¶: eta in z-direction.

zeta¶: zeta.

`maps` – Interpolation routines¶

Interpolation routines mapping grids to grids, grids to fields, and fields to grids.

emg3d.maps.grid2grid(grid, values, new_grid, method='linear', extrapolate=True, log=False)[source]¶

Interpolate values located on grid to new_grid.

Note 1: The default method is ‘linear’, because it works with fields and model parameters. However, recommended are ‘volume’ for model parameters and ‘cubic’ for fields.

Note 2: For model parameters with method=’volume’ the result is quite different if you provide resistivity, conductivity, or the logarithm of any of the two. The recommended way is to provide the logarithm of resistivity or conductivity, in which case the output of one is indeed the inverse of the output of the other.

Parameters:

grid, new_grid : TensorMesh

Input and output model grids; TensorMesh instances.

values : ndarray

Model parameters; emg3d.fields.Field instance, or a particular field (e.g. field.fx). For fields the method cannot be ‘volume’.

method : {<’linear’>, ‘volume’, ‘cubic’}, optional

The method of interpolation to perform. The volume averaging method ensures that the total sum of the property stays constant.

Volume averaging is only implemented for model parameters, not for fields. The method ‘cubic’ requires at least three points in any direction, otherwise it will fall back to ‘linear’.

Default is ‘linear’, because it works with fields and model parameters. However, recommended are ‘volume’ for model parameters and ‘cubic’ for fields.

extrapolate : bool

If True, points on new_grid which are outside of grid are filled by the nearest value (if method='cubic') or by extrapolation (if method='linear'). If False, points outside are set to zero.

For method='volume' it always uses the nearest value for points outside of grid.

Default is True.

log : bool

If True, the interpolation is carried out on a log10-scale; hence the same as 10**grid2grid(grid, np.log10(values), ...). Default is False.

Returns:

new_values : ndarray: Values corresponding to new_grid.

See also

get_receiver: Interpolation of model parameters or fields to (x, y, z).

emg3d.maps.interp3d(points, values, new_points, method, fill_value, mode)[source]¶

Interpolate values in 3D either linearly or with a cubic spline.

Return values corresponding to a regular 3D grid defined by points on new_points.

This is a modified version of scipy.interpolate.interpn(), using scipy.interpolate.RegularGridInterpolator if method='linear' and a custom-wrapped version of scipy.ndimage.map_coordinates() if method='cubic'. If speed is important then choose ‘linear’, as it can be significantly faster.

Parameters:

points : tuple of ndarray of float, with shapes ((nx, ), (ny, ) (nz, )): The points defining the regular grid in three dimensions.
values : array_like, shape (nx, ny, nz): The data on the regular grid in three dimensions.
new_points : tuple (rec_x, rec_y, rec_z): Coordinates (x, y, z) of new points.
method : {‘cubic’, ‘linear’}, optional: The method of interpolation to perform, ‘linear’ or ‘cubic’. Default is ‘cubic’ (forced to ‘linear’ if there are less than 3 points in any direction).
fill_value : float or None: Passed to scipy.interpolate.RegularGridInterpolator if method='linear': The value to use for points outside of the interpolation domain. If None, values outside the domain are extrapolated.
mode : {‘constant’, ‘nearest’, ‘mirror’, ‘reflect’, ‘wrap’}: Passed to scipy.ndimage.map_coordinates() if method='cubic': Determines how the input array is extended beyond its boundaries.

Returns:

new_values : ndarray: Values corresponding to new_points.

`fields` – Electric and magnetic fields¶

Everything related to the multigrid solver that is a field: source field, electric and magnetic fields, and fields at receivers.

class emg3d.fields.Field[source]¶

Create a Field instance with x-, y-, and z-views of the field.

A Field is an ndarray with additional views of the x-, y-, and z-directed fields as attributes, stored as fx, fy, and fz. The default array contains the whole field, which can be the electric field, the source field, or the residual field, in a 1D array. A Field instance has additionally the property ensure_pec which, if called, ensures Perfect Electric Conductor (PEC) boundary condition. It also has the two attributes amp and pha for the amplitude and phase, as common in frequency-domain CSEM.

A Field can be initiated in three ways:

Field(grid, dtype=complex): Calling it with a TensorMesh instance returns a Field instance of correct dimensions initiated with zeroes of data type dtype.
Field(grid, field): Calling it with a TensorMesh instance and an ndarray returns a Field instance of the provided ndarray, of same data type.
Field(fx, fy, fz): Calling it with three ndarray’s which represent the field in x-, y-, and z-direction returns a Field instance with these views, of same data type.

Sort-order is ‘F’.

Parameters:

fx_or_grid : TensorMesh or ndarray

Either a TensorMesh instance or an ndarray of shape grid.nEx or grid.vnEx. See explanations above. Only mandatory parameter; if the only one provided, it will initiate a zero-field of dtype.

fy_or_field : Field or ndarray, optional

Either a Field instance or an ndarray of shape grid.nEy or grid.vnEy. See explanations above.

fz : ndarray, optional

An ndarray of shape grid.nEz or grid.vnEz. See explanations above.

dtype : dtype, optional

Only used if fy_or_field=None and fz=None; the initiated zero-field for the provided TensorMesh has data type dtype. Default: complex.

freq : float, optional

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

Just added as info if provided.

amp(self)[source]¶: Amplitude of the electromagnetic field.

copy(self)[source]¶: Return a copy of the Field.

ensure_pec¶: Set Perfect Electric Conductor (PEC) boundary condition.

field¶: Entire field, 1D [fx, fy, fz].

freq¶: Return frequency.

classmethod from_dict(inp)[source]¶

Convert dictionary into Field instance.

Parameters:	inp : dict Dictionary as obtained from `Field.to_dict()`. The dictionary needs the keys field, freq, vnEx, vnEy, and vnEz.
Returns:	obj : `Field` instance

fx¶: View of the x-directed field in the x-direction (nCx, nNy, nNz).

fy¶: View of the field in the y-direction (nNx, nCy, nNz).

fz¶: View of the field in the z-direction (nNx, nNy, nCz).

is_electric¶: Returns True if Field is electric, False if it is magnetic.

pha(self, deg=False, unwrap=True, lag=True)[source]¶

Phase of the electromagnetic field.

Parameters:	deg : bool If True the returned phase is in degrees, else in radians. Default is False (radians). unwrap : bool If True the returned phase is unwrapped. Default is True (unwrapped). lag : bool If True the returned phase is lag, else lead defined. Default is True (lag defined).

smu0¶: Return s*mu_0; mu_0 = Magn. permeability of free space [H/m].

sval¶: Return s; s=iw in frequency domain; s=freq in Laplace domain.

to_dict(self, copy=False)[source]¶: Store the necessary information of the Field in a dict.

class emg3d.fields.SourceField[source]¶

Create a Source-Field instance with x-, y-, and z-views of the field.

A subclass of Field. Additional properties are the real-valued source vector (vector, vx, vy, vz), which sum is always one. For a SourceField frequency is a mandatory parameter, unlike for a Field (recommended also for Field though),

Parameters:

fx_or_grid : TensorMesh or ndarray

Either a TensorMesh instance or an ndarray of shape grid.nEx or grid.vnEx. See explanations above. Only mandatory parameter; if the only one provided, it will initiate a zero-field of dtype.

fy_or_field : Field or ndarray, optional

Either a Field instance or an ndarray of shape grid.nEy or grid.vnEy. See explanations above.

fz : ndarray, optional

An ndarray of shape grid.nEz or grid.vnEz. See explanations above.

dtype : dtype, optional

Only used if fy_or_field=None and fz=None; the initiated zero-field for the provided TensorMesh has data type dtype. Default: complex.

freq : float

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

In difference to Field, the frequency has to be provided for a SourceField.

copy(self)[source]¶: Return a copy of the SourceField.

classmethod from_dict(inp)[source]¶

Convert dictionary into SourceField instance.

Parameters:	inp : dict Dictionary as obtained from `SourceField.to_dict()`. The dictionary needs the keys field, freq, vnEx, vnEy, and vnEz.
Returns:	obj : `SourceField` instance

vector¶: Entire vector, 1D [vx, vy, vz].

vx¶: View of the x-directed vector in the x-direction (nCx, nNy, nNz).

vy¶: View of the vector in the y-direction (nNx, nCy, nNz).

vz¶: View of the vector in the z-direction (nNx, nNy, nCz).

emg3d.fields.get_source_field(grid, src, freq, strength=0)[source]¶

Return the source field.

The source field is given in Equation 2 in [Muld06],

\[s \mu_0 \mathbf{J}_\mathrm{s} ,\]

where \(s = \mathrm{i} \omega\). Either finite length dipoles or infinitesimal small point dipoles can be defined, whereas the return source field corresponds to a normalized (1 Am) source distributed within the cell(s) it resides (can be changed with the strength-parameter).

The adjoint of the trilinear interpolation is used to distribute the point(s) to the grid edges, which corresponds to the discretization of a Dirac ([PlDM07]).

Parameters:

grid : TensorMesh

Model grid; a TensorMesh instance.

src : list of floats

Source coordinates (m). There are two formats:

Finite length dipole: [x0, x1, y0, y1, z0, z1].

Point dipole: [x, y, z, azimuth, dip].

freq : float

Source frequency (Hz), used to calculate the Laplace parameter s. Either positive or negative:

freq > 0: Frequency domain, hence \(s = -\mathrm{i}\omega = -2\mathrm{i}\pi f\) (complex);
freq < 0: Laplace domain, hence \(s = f\) (real).

strength : float or complex, optional

Source strength (A):

If 0, output is normalized to a source of 1 m length, and source strength of 1 A.

If != 0, output is returned for given source length and strength.

Default is 0.

Returns:

sfield : SourceField() instance: Source field, normalized to 1 A m.

emg3d.fields.get_receiver(grid, values, coordinates, method='cubic', extrapolate=False)[source]¶

Return values corresponding to grid at coordinates.

Works for electric fields as well as magnetic fields obtained with get_h_field(), and for model parameters.

Parameters:

grid : TensorMesh

Model grid; a TensorMesh instance.

values : ndarray

Field instance, or a particular field (e.g. field.fx); Model parameters.

coordinates : tuple (x, y, z)

Coordinates (x, y, z) where to interpolate values; e.g. receiver locations.

method : str, optional

The method of interpolation to perform, ‘linear’ or ‘cubic’. Default is ‘cubic’ (forced to ‘linear’ if there are less than 3 points in any direction).

extrapolate : bool

If True, points on new_grid which are outside of grid are filled by the nearest value (if method='cubic') or by extrapolation (if method='linear'). If False, points outside are set to zero.

Default is False.

Returns:

new_values : ndarray or empymod.utils.EMArray

Values at coordinates.

If input was a field it returns an EMArray, which is a subclassed ndarray with .pha and .amp attributes.

If input was an entire Field instance, output is a tuple (fx, fy, fz).

See also

grid2grid: Interpolation of model parameters or fields to a new grid.

emg3d.fields.get_h_field(grid, model, field)[source]¶

Return magnetic field corresponding to provided electric field.

Retrieve the magnetic field \(\mathbf{H}\) from the electric field \(\mathbf{E}\) using Farady’s law, given by

\[\nabla \times \mathbf{E} = \rm{i}\omega\mu\mathbf{H} .\]

Note that the magnetic field in x-direction is defined in the center of the face defined by the electric field in y- and z-directions, and similar for the other field directions. This means that the provided electric field and the returned magnetic field have different dimensions:

E-field:  x: [grid.vectorCCx,  grid.vectorNy,  grid.vectorNz]
          y: [ grid.vectorNx, grid.vectorCCy,  grid.vectorNz]
          z: [ grid.vectorNx,  grid.vectorNy, grid.vectorCCz]

H-field:  x: [ grid.vectorNx, grid.vectorCCy, grid.vectorCCz]
          y: [grid.vectorCCx,  grid.vectorNy, grid.vectorCCz]
          z: [grid.vectorCCx, grid.vectorCCy,  grid.vectorNz]

Parameters:	grid : TensorMesh Model grid; `TensorMesh` instance. model : Model Model; `Model` instance. field : Field Electric field; `Field` instance.
Returns:	hfield : Field Magnetic field; `Field` instance.

`io` – I/O utilities¶

Utility functions for writing and reading data.

emg3d.io.save(fname, backend='h5py', compression='gzip', **kwargs)[source]¶

Save meshes, models, fields, and other data to disk.

Serialize and save emg3d.meshes.TensorMesh, emg3d.fields.Field, and emg3d.models.Model instances and add arbitrary other data, where instances of the same type are grouped together.

The serialized instances will be de-serialized if loaded with load().

Parameters:

fname : str

File name.

backend : str, optional

Backend to use. Implemented are currently:

h5py (default): Uses h5py to store inputs to a hierarchical, compressed binary hdf5 file with the extension ‘.h5’. Recommended and default backend, but requires the module h5py. Use numpy if you don’t want to install h5py.
numpy: Uses numpy to store inputs to a flat, compressed binary file with the extension ‘.npz’.

compression : int or str, optional

Passed through to h5py, default is ‘gzip’.

kwargs : Keyword arguments, optional

Data to save using its key as name. The following instances will be properly serialized: emg3d.meshes.TensorMesh, emg3d.fields.Field, and emg3d.models.Model and serialized again if loaded with load(). These instances are collected in their own group if h5py is used.

emg3d.io.load(fname, **kwargs)[source]¶

Load meshes, models, fields, and other data from disk.

Load and de-serialize emg3d.meshes.TensorMesh, emg3d.fields.Field, and emg3d.models.Model instances and add arbitrary other data that were saved with save().

Parameters:	fname : str File name including extension. Used backend depends on the file extensions: ‘.npz’: numpy-binary ‘.h5’: h5py-binary (needs h5py) verb : int If 1 (default) verbose, if 0 silent.
Returns:	out : dict A dictionary containing the data stored in fname; `emg3d.meshes.TensorMesh`, `emg3d.fields.Field`, and `emg3d.models.Model` instances are de-serialized and returned as instances.

Code¶

solver – Multigrid solver¶

core – Number crunching¶

utils – Utilities¶

meshes – Discretization¶

models – Earth properties¶

maps – Interpolation routines¶

fields – Electric and magnetic fields¶

io – I/O utilities¶

`solver` – Multigrid solver¶

`core` – Number crunching¶

`utils` – Utilities¶

`meshes` – Discretization¶

`models` – Earth properties¶

`maps` – Interpolation routines¶

`fields` – Electric and magnetic fields¶

`io` – I/O utilities¶