emg3d

Version: 0.16.1 ~ Date: 09 February 2021

A multigrid solver for 3D electromagnetic diffusion with tri-axial electrical anisotropy. The matrix-free 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- and memory-consuming parts are sped up through jitted numba-functions.

More information

For more information regarding installation, usage, contributing, roadmap, bug reports, and much more, see

• Website: https://emsig.github.io,
• Documentation: https://emg3d.readthedocs.io,
• Source Code: https://github.com/emsig/emg3d,
• Examples: https://emsig.github.io/emg3d-gallery.

Features

• Multigrid solver for 3D electromagnetic (EM) diffusion with regular grids (where source and receiver can be electric or magnetic).
• Compute the 3D EM field in the complex frequency domain or in the real Laplace domain.
• Includes also routines to compute the 3D EM field in the time domain.
• Can be used together with the SimPEG-framework.
• Can be used as a standalone solver or as a pre-conditioner for various Krylov subspace methods implemented in SciPy, e.g., BiCGSTAB (scipy.sparse.linalg.bicgstab) or CGS (scipy.sparse.linalg.cgs).
• Tri-axial electrical anisotropy.
• Isotropic magnetic permeability.
• Semicoarsening and line relaxation.
• Grid-size can be anything.
• As a multigrid method it scales with the number of unknowns N and has therefore optimal complexity O(N).

Installation

You can install emg3d either via conda (preferred):

conda install -c conda-forge emg3d

or via pip:

pip install emg3d

Minimum requirements are Python version 3.7 or higher and the modules scipy and numba. Various other packages are recommended or required for some advanced functionalities (xarray, discretize, matplotlib, h5py, empymod, scooby). Consult the installation notes in the manual for more information regarding installation, requirements, and soft dependencies.

Citation

If you publish results for which you used emg3d, please give credit by citing Werthmüller et al. (2019):

Werthmüller, D., W. A. Mulder, and E. C. Slob, 2019, emg3d: A multigrid solver for 3D electromagnetic diffusion: Journal of Open Source Software, 4(39), 1463; DOI: 10.21105/joss.01463.

All releases have a Zenodo-DOI, which can be found on 10.5281/zenodo.3229006.

See CREDITS for the history of the code.

License information

Copyright 2018-2021 The emg3d Developers.

Licensed under the Apache License, Version 2.0, see the LICENSE-file.

Getting started

The code emg3d ([WeMS19]) is a three-dimensional modeller for electromagnetic (EM) diffusion as used, for instance, in controlled-source EM (CSEM) surveys frequently applied in the search for, amongst other, groundwater, hydrocarbons, and minerals.

The core of the code is primarily based on [Muld06], [Muld07], and [Muld08]. You can read more about the background of the code in the chapter Credits. An introduction to the underlying theory of multigrid methods is given in the chapter Theory, and further literature is provided in the References.

Installation

You can install emg3d either via conda:

conda install -c conda-forge emg3d

or via pip:

pip install emg3d

Minimum requirements are Python version 3.7 or higher and the modules scipy and numba. Various other packages are recommended or required for some advanced functionalities, namely:

• xarray: For the Survey class (many sources and receivers at once).
• discretize: For advanced meshing tools (fancy mesh-representations and plotting utilities).
• matplotlib: To use the plotting utilities within discretize.
• h5py: Save and load data in the HDF5 format.
• empymod: Time-domain modelling (utils.Fourier).
• scooby: For the version and system report (emg3d.Report()).

If you are new to Python we recommend using a Python distribution, which will ensure that all dependencies are met, specifically properly compiled versions of NumPy and SciPy; we recommend using Anaconda. If you install Anaconda you can simply start the Anaconda Navigator, add the channel conda-forge and emg3d will appear in the package list and can be installed with a click.

Using NumPy and SciPy with the Intel Math Kernel Library (mkl) can significantly improve computation time. You can check if mkl is used via conda list: The entries for the BLAS and LAPACK libraries should contain something with mkl, not with openblas. To enforce it you might have to create a file pinned, containing the line libblas[build=*mkl] in the folder path-to-your-conda-env/conda-meta/.

Basic Example

Here we show a very basic example. To see some more realistic models have a look at the gallery. This particular example is also there, with some further explanations and examples to show how to plot the model and the data; see «Minimum working example». It also contains an example without using discretize.

First, we load emg3d and discretize (to create a mesh), along with numpy:

>>> import emg3d
>>> import discretize
>>> import numpy as np

First, we define the mesh (see discretize.TensorMesh for more info). In reality, this task requires some careful considerations. E.g., to avoid edge effects, the mesh should be large enough in order for the fields to dissipate, yet fine enough around source and receiver to accurately model them. This grid is too small, but serves as a minimal example.

>>> grid = discretize.TensorMesh(
>>>         [[(25, 10, -1.04), (25, 28), (25, 10, 1.04)],
>>>          [(50, 8, -1.03), (50, 16), (50, 8, 1.03)],
>>>          [(30, 8, -1.05), (30, 16), (30, 8, 1.05)]],
>>>         x0='CCC')
>>> print(grid)

  TensorMesh: 49,152 cells

                      MESH EXTENT             CELL WIDTH      FACTOR
  dir    nC        min           max         min       max      max
  ---   ---  ---------------------------  ------------------  ------
   x     48       -662.16        662.16     25.00     37.01    1.04
   y     32       -857.96        857.96     50.00     63.34    1.03
   z     32       -540.80        540.80     30.00     44.32    1.05

Next we define a very simple fullspace model with \(\rho_x=1.5\,\Omega\,\text{m}\), \(\rho_y=1.8\,\Omega\,\text{m}\), and \(\rho_z=3.3\,\Omega\,\text{m}\). The source is an x-directed dipole at the origin, with a 10 Hz signal of 1 A.

>>> model = emg3d.models.Model(
>>>     grid, property_x=1.5, property_y=1.8, property_z=3.3)
>>> sfield = emg3d.fields.get_source_field(
>>>     grid, src=[0, 0, 0, 0, 0], freq=10.0)

Now we can compute the electric field with emg3d:

>>> efield = emg3d.solve(grid, model, sfield, verb=4)

:: emg3d START :: 15:24:40 :: v0.13.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      : 3
   Original grid  :  48 x  32 x  32     => 49,152 cells
   Coarsest grid  :   3 x   2 x   2     => 12 cells
   Coarsest level :   4 ;   4 ;   4

   [hh:mm:ss]  rel. error                  [abs. error, last/prev]   l s

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

   [11:18:17]   2.623e-02  after   1 F-cycles   [1.464e-06, 0.026]   0 0
   [11:18:17]   2.253e-03  after   2 F-cycles   [1.258e-07, 0.086]   0 0
   [11:18:17]   3.051e-04  after   3 F-cycles   [1.704e-08, 0.135]   0 0
   [11:18:17]   5.500e-05  after   4 F-cycles   [3.071e-09, 0.180]   0 0
   [11:18:18]   1.170e-05  after   5 F-cycles   [6.531e-10, 0.213]   0 0
   [11:18:18]   2.745e-06  after   6 F-cycles   [1.532e-10, 0.235]   0 0
   [11:18:18]   6.873e-07  after   7 F-cycles   [3.837e-11, 0.250]   0 0

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

:: emg3d END   :: 15:24:42 :: runtime = 0:00:02

So the computation required seven multigrid F-cycles and took just a bit more than 2 seconds. It was able to coarsen in each dimension four times, where the input grid had 49,152 cells, and the coarsest grid had 12 cells.

Related ecosystem

To create advanced meshes it is recommended to use discretize from the SimPEG framework. It also comes with some neat plotting functionalities to plot model parameters and resulting fields. Furthermore, it can serve as a link to use PyVista to create nice 3D plots even within a notebook.

Projects which can be used to compare or validate the results are, e.g., empymod for layered models or SimPEG for 3D models. It is also possible to create a geological model with GemPy and, again via discretize, move it to emg3d to compute CSEM responses for it.

Have a look at the gallery for many examples of how to use emg3d together with the mentioned projects and more!

Tips and Tricks

The function emg3d.solver.solve() is the main entry point, and it takes care whether multigrid is used as a solver or as a preconditioner (or not at all), while the actual multigrid solver is emg3d.solver.multigrid(). Most input parameters for emg3d.solver.solve() are sufficiently described in its docstring. Here a few additional information.

• You can input any three-dimensional tensor mesh into emg3d. However, the implemented multigrid technique works with the existing nodes, meaning there are no new nodes created as coarsening is done by combining adjacent cells. The more times the grid dimension can be divided by two the better it is suited for MG. Ideally, the number should be dividable by two a few times and the dimension of the coarsest grid should be a low prime number \(p\), for which good sizes can then be computed with \(p 2^n\). Good grid sizes (in each direction) up to 1024 are

  • \(2·2^{3, 4, ..., 9}\): 16, 32, 64, 128, 256, 512, 1024,
  • \(3·2^{3, 4, ..., 8}\): 24, 48, 96, 192, 384, 768,
  • \(5·2^{3, 4, ..., 7}\): 40, 80, 160, 320, 640,
  • \(7·2^{3, 4, ..., 7}\): 56, 112, 224, 448, 896,

  and preference decreases from top to bottom row (stick to the first two or three rows if possible). Good grid sizes in sequential order, excluding p=7: 16, 24, 32, 40, 48, 64, 80, 96, 128, 160, 192, 256, 320, 384, 512, 640, 768, 1024. You can get this list via emg3d.meshes.good_mg_cell_nr().

• The multigrid method can be used as a solver or as a preconditioner, for instance for BiCGSTAB. Using multigrid as a preconditioner for BiCGSTAB together with semicoarsening and line relaxation is the most stable version, but expensive, and therefore only recommended on highly stretched grids. Which combination of solver is best (fastest) depends to a large extent on the grid stretching, but also on anisotropy and general model complexity. See «Parameter tests» in the gallery for an example how to run some tests on your particular problem.

Contributing and Roadmap

New contributions, bug reports, or any kind of feedback is always welcomed! Have a look at the Roadmap-project to get an idea of things that could be implemented. The GitHub issues and PR’s are also a good starting point. The best way for interaction is at https://github.com/emsig or by joining the Slack channel «em-x-d» of SimPEG. If you prefer to get in touch outside of GitHub/Slack use the contact form on https://werthmuller.org.

To install emg3d from source, you can download the latest version from GitHub and install it in your python distribution via:

python setup.py install

Please make sure your code follows the pep8-guidelines by using, for instance, the python module flake8, and also that your code is covered with appropriate tests. Just get in touch if you have any doubts.

Tests and benchmarks

The modeller comes with a test suite using pytest. If you want to run the tests, just install pytest and run it within the emg3d-top-directory.

> pytest --cov=emg3d --flake8

It should run all tests successfully. Please let us know if not!

Note that installations of em3gd via conda or pip do not have the test-suite included. To run the test-suite you must download emg3d from GitHub.

There is also a benchmark suite using airspeed velocity, located in the emsig/emg3d-asv-repository. The results of my machine can be found in the emsig/emg3d-bench, its rendered version at emsig.github.io/emg3d-asv.

License

Copyright 2018-2021 The emg3d Developers.

Licensed under the Apache License, Version 2.0 (the “License”); you may not use this file except in compliance with the License. You may obtain a copy of the License at

https://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software distributed under the License is distributed on an “AS IS” BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the License for the specific language governing permissions and limitations under the License.

Multi-what?

If you have never heard of the multigrid method before you might ask yourself “multi-what?” The following is an intent to describe the multigrid method without the maths; just some keywords and some figures. It is a heavily simplified intro, using a 2D grid for simplicity. Have a look at the Theory-section for more details. A good, four-page intro with some maths is given by [Muld11]. More in-depth information can be found, e.g., in [BrHM00], [Hack85], and [Wess91].

The multigrid method ([Fedo64])

• is an iterative solver;
• scales almost linearly (CPU & RAM);
• can serve as a pre-conditioner or as a solver on its own.

The main driving motivation to use multigrid is the part about linear scaling.

Matrix-free solver

The implemented multigrid method is a matrix free solver, it never constructs the full matrix. This is how it achieves its relatively low memory consumption. To solve the system, it solves for all fields adjacent to one node, moves then to the next node, and so on until it reaches the last node, see Figure 1, where the red lines indicate the fields which are solved simultaneously per step (the fields on the boundaries are never computed, as they are assumed to be 0).

The multigrid solver solves by default on a node-by-node basis.

Normally, you would have to do this over and over again to achieve a good approximate solution. multigrid typically does it only a few times per grid, typically 2 times (one forward, one backward). This is why it is called smoother, as it only smoothes the error, it does not solve it. The implemented method for this is the Gauss-Seidel method.

Iterative solver which work in this matrix-free manner are typically very fast at solving for the local problem, hence at reducing the high frequency error, but very slow at solving the global problem, hence at reducing the low frequency error. High and low frequency errors are meant relatively to cell-size here.

Moving between different grids

The main thinking behind multigrid is now that we move to coarser grids. This has two advantages:

• Fewer cells means faster computation and less memory.
• Coarser grid size transforms lower frequency error to higher frequency error, relatively to cell size, which means faster convergence.

The implemented multigrid method simply joins two adjacent cells to get from finer to coarser grids, see Figure 2 for an example coarsening starting with a 16 cells by 16 cells grid.

Example of the implemented coarsening scheme.

There are different approaches how to cycle through different grid sizes, see Figures 7 to 9. The downsampling from a finer grid to a coarser grid is often termed restriction, whereas the interpolation from a coarser grid to a finer grid is termed prolongation.

Specialities

The convergence rate of the multigrid method suffers on severely stretched grids or by models with strong anisotropy. Two techniques are implemented, semicoarsening (Figure 3) and line relaxation (Figure 4). Both require more CPU and higher RAM per grid than the standard multigrid, but they can improve the convergence rate, which then in turn improves the overall CPU time.

Example of semicoarsening: The cell size is kept constant in one direction. The direction can be alternated between iterations.

Example of line relaxation: The system is solved for all fields adjacent to a whole line of nodes simultaneously in some direction. The direction can be alternated between iterations.

Theory

The following provides an introduction to the theoretical foundation of the solver emg3d. More specific theory is covered in the docstrings of many functions, have a look at the Other functions-section or follow the links to the corresponding functions here within the theory. If you just want to use the solver, but do not care much about the internal functionality, then the function emg3d.solver.solve() is the only function you will ever need. It is the main entry point, and it takes care whether multigrid is used as a solver or as a preconditioner (or not at all), while the actual multigrid solver is emg3d.solver.multigrid().

Note

This section is not an independent piece of work. Most things are taken from one of the following sources:

• [Muld06], pages 634-639:
  • The Maxwell’s equations and Discretisation sections correspond with some adjustemens and additions to pages 634-636.
  • The start of The Multigrid Method corresponds roughly to page 637.
  • Pages 638 and 639 are in parts reproduced in the code-docstrings of the corresponding functions.
• [BrHM00]: This book is an excellent introduction to multigrid methods. Particularly the Iterative Solvers section is taken to a big extent from the book.

Please consult these original resources for more details, and refer to them for citation purposes and not to this manual. More in-depth information can also be found in, e.g., [Hack85] and [Wess91].

Maxwell’s equations

Maxwell’s equations in the presence of a current source \(\mathbf{J}_\mathrm{s}\) are

()\[\begin{split}\partial_{t} \mathbf{B}(\mathbf{x},t) + \nabla\times\mathbf{E}(\mathbf{x},t) &= 0 , \\ \nabla \times \mathbf{H}(\mathbf{x}, t) - \partial_{t} \mathbf{D}(\mathbf{x}, t) &= \mathbf{J}_{\mathrm{c}}(\mathbf{x}, t) + \mathbf{J}_\mathrm{s}(\mathbf{x}, t) ,\end{split}\]

where the conduction current \(\mathbf{J}_\mathrm{c}\) obeys Ohm’s law,

()\[\mathbf{J}_{\mathrm{c}}(\mathbf{x},t) = \sigma(\mathbf{x})\mathbf{E}(\mathbf{x},t) .\]

Here, \(\sigma(\mathbf{x})\) is the conductivity. \(\mathbf{E}(\mathbf{x}, t)\) is the electric field and \(\mathbf{H}(\mathbf{x}, t)\) is the magnetic field. The electric displacement \(\mathbf{D}(\mathbf{x}, t) = \varepsilon(\mathbf{x})\mathbf{E}(\mathbf{x}, t)\) and the magnetic induction \(\mathbf{B}(\mathbf{x}, t) = \mu(\mathbf{x})\mathbf{H}(\mathbf{x}, t)\). The dielectric constant or permittivity \(\varepsilon\) can be expressed as \(\varepsilon = \varepsilon_r \varepsilon_0\), where \(\varepsilon_r\) is the relative permittivity and \(\varepsilon_0\) is the vacuum value. Similarly, the magnetic permeability \(\mu\) can be written as \(\mu = \mu_r\mu_0\), where \(\mu_r\) is the relative permeability and \(\mu_0\) is the vacuum value.

The magnetic field can be eliminated from Equation (1), yielding the second-order parabolic system of equations,

()\[\varepsilon \partial_{t t} \mathbf{E} + \sigma \partial_{t} \mathbf{E} + \nabla \times \mu^{-1} \nabla \times \mathbf{E} = -\partial_{t} \mathbf{J}_{\mathrm{s}} .\]

To transform from the time domain to the frequency domain, we substitute

()\[\mathbf{E} (\mathbf{x},t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \mathbf{\hat{E}} (\mathbf{x},\omega) e^{-\mathrm{i}\omega t}\, d\omega,\]

and use a similar representation for \(\mathbf{H}(\mathbf{x}, t)\). The resulting system of equations is

()\[-s \mu_0(\sigma + s\varepsilon) \mathbf{\hat{E}} - \nabla \times \mu_r^{-1} \nabla \times \mathbf{\hat{E}} = s\mu_0\mathbf{\hat{J}}_s ,\]

where \(s = -\mathrm{i}\omega\). The multigrid method converges in the case of the diffusive approximation (with its smoothing and approximation properties), but not in the high-frequency range (at least not in the implemented form of the multigrid method in emg3d). The code emg3d assumes therefore the diffusive approximation, hence only low frequencies are considered that obey \(|\omega\varepsilon| \ll \sigma\). In this case we can set \(\varepsilon=0\), and Equation (5) simplifies to

()\[-s \mu_0 \sigma \mathbf{\hat{E}} - \nabla \times \mu_r^{-1} \nabla \times \mathbf{\hat{E}} = s\mu_0\mathbf{\hat{J}}_s ,\]

From here on, the hats are omitted. We use the perfectly electrically conducting boundary

()\[\mathbf{n}\times\mathbf{E} = 0 \quad \text{and} \quad \mathbf{n}\cdot\mathbf{H} = 0 , \label{eq:sample}\]

where \(\mathbf{n}\) is the outward normal on the boundary of the domain.

The Maxwell’s equations and Ohm’s law are solved in the frequency domain. The time-domain solution can be obtained by taking the inverse Fourier transform.

Note

[Muld06] uses the time convention \(e^{-\mathrm{i}\omega t}\), see Equation (4), with \(s=-\mathrm{i}\omega\). However, the code emg3d uses the convention \(e^{\mathrm{i}\omega t}\), hence \(s=\mathrm{i}\omega\). This is the same convention as used in empymod, and commonly in CSEM.

Laplace domain

It is also possible to solve the problem in the Laplace domain, by using a real value for \(s\) in Equation (6), instead of the complex value \(-\mathrm{i}\omega\). This simplifies the problem from complex numbers to real numbers, which accelerates the computation. It also improves the convergence rate, as the solution is a smoother function. The solver emg3d.solver.solve() is agnostic to the data type of the provided source field, and can solve for real and complex problems, hence frequency and Laplace domain. See the documentation of the functions emg3d.fields.get_source_field() and emg3d.models.Model() to see how you can use emg3d for Laplace-domain computations.

Discretisation

Equation (6) can be discretised by the finite-integration technique ([Weil77], [ClWe01]). This scheme can be viewed as a finite-volume generalization of [Yee66]’s scheme for tensor-product Cartesian grids with variable grid spacings. An error analysis for the constant-coefficient case ([MoSu94]) showed that both the electric and magnetic field components have second-order accuracy.

Consider a tensor-product Cartesian grid with nodes at positions \((x_k, y_l, z_m)\), where \(k=0, \dots, N_x, l=0, \dots, N_y\) and \(m=0, \dots, N_z\). There are \(N_x\times N_y\times N_z\) cells having these nodes as vertices. The cell centres are located at

()\[\begin{split}x_{k+1/2} &= {\textstyle \frac{1}{2}}\left(x_k + x_{k+1}\right) , \\ y_{l+1/2} &= {\textstyle \frac{1}{2}}\left(y_l + y_{l+1}\right) , \\ z_{m+1/2} &= {\textstyle \frac{1}{2}}\left(z_m + z_{m+1}\right) .\end{split}\]

The material properties, \(\sigma\) and \(\mu_\mathrm{r}\), are assumed to be given as cell-averaged values. The electric field components are positioned at the edges of the cells, as shown in Figure 5, in a manner similar to Yee’s scheme. The first component of the electric field \(E_{1, k+1/2, l, m}\) should approximate the average of \(E_1(x, y_l, z_m)\) over the edge from \(x_k\) to \(x_{k+1}\) at given \(y_l\) and \(z_m\). Here, the average is defined as the line integral divided by the length of the integration interval. The other components, \(E_{2, k, l+1/2, m}\) and \(E_{3, k, l, m+1/2}\), are defined in a similar way. Note that these averages may also be interpreted as point values at the midpoint of edges:

()\[\begin{split}E_{1, k+1/2, l, m} \simeq E_1\left(x_{k+1/2}, y_{l}, z_{m}\right) , \\ E_{2, k, l+1/2, m} \simeq E_2\left(x_{k}, y_{l+1/2}, z_{m}\right) , \\ E_{3, k, l, m+1/2} \simeq E_3\left(x_{k}, y_{l}, z_{m+1/2}\right) .\end{split}\]

The averages and point-values are the same within second-order accuracy.

(a) A grid cell with grid nodes and edge-averaged components of the electric field. (b) The face-averaged magnetic field components that are obtained by taking the curl of the electric field.

For the discretisation of the term \(-s\mu_0\sigma\mathbf{E}\) related to Ohm’s law, dual volumes related to edges are introduced. For a given edge, the dual volume is a quarter of the total volume of the four adjacent cells. An example for \(E_1\) is shown in Figure 6(b). The vertices of the dual cell are located at the midpoints of the cell faces.

The first electric field component \(E_{1,k,l,m}\) is located at the intersection of the four cells shown in (a). Four faces of its dual volume are sketched in (b). The first component of the curl of the magnetic field should coincide with the edge on which \(E_1\) is located. The four vectors that contribute to this curl are shown in (a). They are defined as normals to the four faces in (a). Before computing their curl, these vectors are interpreted as tangential components at the faces of the dual volume shown in (b). The curl is evaluated by taking the path integral over a rectangle of the dual volume that is obtained for constant x and by averaging over the interval \([x_k,x_{k+1}]\).

The volume of a normal cell is defined as

()\[V_{k+1/2, l+1/2, m+1/2} = h_{k+1/2}^x h_{l+1/2}^y h_{m+1/2}^z ,\]

where

()\[\begin{split}h_{k+1/2}^x &= x_{k+1}-x_k , \\ h_{l+1/2}^y &= y_{l+1}-y_l , \\ h_{m+1/2}^z &= z_{m+1}-z_m .\end{split}\]

For an edge parallel to the x-axis on which \(E_{1, k+1/2, l, m}\) is located, the dual volume is

()\[V_{k+1/2, l, m} = {\textstyle \frac{1}{4}} h_{k+1/2}^x \sum_{m_2=0}^1 \sum_{m_3=0}^1 h_{l-1/2+m_2}^y h_{m-1/2+m_3}^z .\]

With the definitions,

()\[\begin{split}d_k^x &= x_{k+1/2} - x_{k-1/2} , \\ d_l^y &= y_{l+1/2} - y_{l-1/2} , \\ d_m^z &= z_{m+1/2} - z_{m-1/2} ,\end{split}\]

we obtain

()\[\begin{split}V_{k+1/2, l, m} &= h_{k+1/2}^x d_l^y d_m^z , \\ V_{k, l+1/2, m} &= d_k^x h_{l+1/2}^y d_m^z , \\ V_{k, l, m+1/2} &= d_k^x d_l^y h_{m+1/2}^z .\end{split}\]

Note that Equation (13) does not define \(d_k^x\), etc., at the boundaries. We may simply take \(d^x_0 = h^x_{1/2}\) at \(k = 0\), \(d^x_{N_x} = h^x_{N_x-1/2}\) at \(k = N_x\) and so on, or use half of these values as was done by [MoSu94].

The discrete form of the term \(-s\mu_0\sigma\mathbf{E}\) in Equation (6), with each component multiplied by the corresponding dual volume, becomes \(\mathcal{S}_{k+1/2, l, m}\ E_{1, k+1/2, l, m}\), \(\mathcal{S}_{k, l+1/2, m}\ E_{2, k, l+1/2, m}\) and \(\mathcal{S}_{k, l, m+1/2}\ E_{3, k, l, m+1/2}\) for the first, second and third components, respectively. Here \(\mathcal{S} = -s\mu_0\sigma V\) is defined in terms of cell-averages. At the edges parallel to the x-axis, an averaging procedure similar to (12) gives

()\[\begin{split}\mathcal{S}_{k+1/2, l, m} = &{\textstyle\frac{1}{4}}\left( \mathcal{S}_{k+1/2, l-1/2, m-1/2} + \mathcal{S}_{k+1/2, l+1/2, m-1/2} \right. \\ &+ \left. \mathcal{S}_{k+1/2, l-1/2, m+1/2} + \mathcal{S}_{k+1/2, l+1/2, m+1/2} \right) .\end{split}\]

\(\mathcal{S}_{k, l+1/2, m}\) and \(\mathcal{S}_{k, l, m+1/2}\) are defined in a similar way.

The curl of \(\mathbf{E}\) follows from path integrals around the edges that bound a face of a cell, drawn in Figure 5(a). After division by the area of the faces, the result is a face-averaged value that can be positioned at the centre of the face, as sketched in Figure 5(b). If this result is divided by \(\mathrm{i}\omega\mu\), the component of the magnetic field that is normal to the face is obtained. In order to find the curl of the magnetic field, the magnetic field components that are normal to faces are interpreted as tangential components at the faces of the dual volumes. For \(E_1\), this is shown in Figure 6. For the first component of Equation (6) on the edge \((k+1/2, l, m)\) connecting \((x_k, y_l, z_m)\) and \((x_{k+1}, y_l, z_m)\), the corresponding dual volume comprises the set \([x_k, x_{k+1}] \times [y_{l-1/2}, y_{l+1/2}] \times [z_{m-1/2}, z_{m+1/2}]\) having volume \(V_{k+1/2,l,m}\).

The scaling by \(\mu_r^{-1}\) at the face requires another averaging step because the material properties are assumed to be given as cell-averaged values. We define \(\mathcal{M} = V\mu_r^{-1}\), so

()\[\mathcal{M}_{k+1/2, l+1/2, m+1/2} = h_{k+1/2}^x h_{l+1/2}^y h_{m+1/2}^z \mu_{r, k+1/2, l+1/2, m+1/2}^{-1}\]

for a given cell \((k+1/2, l+1/2, m+1/2)\). An averaging step in, for instance, the z-direction gives

()\[\mathcal{M}_{k+1/2, l+1/2, m} = {\textstyle \frac{1}{2}} \left(\mathcal{M}_{k+1/2, l+1/2, m-1/2} + \mathcal{M}_{k+1/2, l+1/2, m+1/2} \right)\]

at the face \((k+1/2, l+1/2, m)\) between the cells \((k+1/2, l+1/2, m-1/2)\) and \((k+1/2, l+1/2, m+1/2)\).

Starting with \(\mathbf{v}=\nabla \times \mathbf{E}\), we have

()\[\begin{split}v_{1, k, l+1/2, m+1/2} &= e_{l+1/2}^y\left(E_{3, k, l+1, m+1/2} - E_{3, k, l, m+1/2}\right) \\ &-e_{m+1/2}^z\left(E_{2, k, l+1/2, m+1} - E_{2, k, l+1/2, m}\right) , \\ v_{2, k+1/2, l, m+1/2} &= e_{m+1/2}^z\left(E_{1, k+1/2, l, m+1} - E_{1, k+1/2, l, m}\right) \\ &-e_{k+1/2}^x\left(E_{3, k+1, l, m+1/2} - E_{3, k, l, m+1/2}\right) , \\ v_{3, k+1/2, l+1/2, m} &= e_{k+1/2}^x\left(E_{2, k+1/2, l+1, m} - E_{1, k+1/2, l, m}\right) \\ &-e_{l+1/2}^y\left(E_{1, k+1/2, l+1, m} - E_{1, k+1/2, l, m}\right) .\end{split}\]

Here,

()\[e_{k+1/2}^x = 1/h_{k+1/2}^x, \quad e_{l+1/2}^y = 1/h_{l+1/2}^y, \quad e_{m+1/2}^z = 1/h_{m+1/2}^z .\]

Next, we let

()\[\begin{split}u_{1,k,l+1/2,m+1/2} &= \mathcal{M}_{k,l+1/2,m+1/2} v_{1,k,l+1/2,m+1/2} , \\ u_{2,k+1/2,l,m+1/2} &= \mathcal{M}_{k+1/2,l,m+1/2} v_{2,k+1/2,l+1/2,m} , \\ u_{3,k+1/2,l+1/2,m} &= \mathcal{M}_{k+1/2,l+1/2,m} v_{3,k+1/2,l+1/2,m} .\end{split}\]

Note that these components are related to the magnetic field components by

()\[\begin{split}u_{1,k,l+1/2,m+1/2} &= \mathrm{i}\omega\mu_0 V_{k,l+1/2,m+1/2} H_{1,k+1/2,l,m+1/2} , \\ u_{2,k+1/2,l,m+1/2} &= \mathrm{i}\omega\mu_0 V_{k+1/2,l,m+1/2} H_{2,k+1/2,l,m+1/2} , \\ u_{3,k+1/2,l+1/2,m} &= \mathrm{i}\omega\mu_0 V_{k+1/2,l+1/2,m} H_{3,k+1/2,l+1/2,m} ,\end{split}\]

where

()\[\begin{split}V_{k,l+1/2,m+1/2} &= d_k^x h_{l+1/2}^y h_{m+1/2}^z , \\ V_{k+1/2,l,m+1/2} &= h_{k+1/2}^x d_l^y h_{m+1/2}^z , \\ V_{k+1/2,l+1/2,m} &= h_{k+1/2}^x h_{l+1/2}^y d_m^z .\end{split}\]

The discrete representation of the source term \(\mathrm{i}\omega\mu_0\mathbf{J}_\mathrm{s}\), multiplied by the appropriate dual volume, is

()\[\begin{split}s_{1,k+1/2,l,m} &= \mathrm{i}\omega\mu_0 V_{k+1/2,l,m} J_{1,k+1/2,l,m} , \\ s_{2,k,l+1/2,m} &= \mathrm{i}\omega\mu_0 V_{k,l+1/2,m} J_{2,k,l+1/2,m} , \\ s_{3,k,l,m+1/2} &= \mathrm{i}\omega\mu_0 V_{k,l,m+1/2} J_{3,k,l,m+1/2} .\end{split}\]

Let the residual for an arbitrary electric field that is not necessarily a solution to the problem be defined as

()\[\mathbf{r} = V \left(\mathrm{i} \omega \mu_0 \mathbf{J}_\mathrm{s} + -s\mu_0\sigma \mathbf{E} - \nabla \times \mu^{-1}_\mathrm{r} \nabla \times \mathbf{E}\right) .\]

Its discretisation is

()\[\begin{split}r_{1,k+1/2,l,m} = ~&s_{1,k+1/2,l,m} + \mathcal{S}_{k+1/2,l,m} E_{1,k+1/2,l,m} \\ &-\left[e_{l+1/2}^y u_{3,k+1/2,l+1/2,m} - e_{l-1/2}^y u_{3,k+1/2,l-1/2,m]}\right.\\ &+\left[e_{m+1/2}^z u_{2,k+1/2,l,m+1/2} - e_{m-1/2}^z u_{2,k+1/2,l,m-1/2}\right] , \\ % r_{2,k,l+1/2,m} = ~&s_{2,k,l+1/2,m} + \mathcal{S}_{k,l+1/2,m} E_{2,k,l+1/2,m} \\ &-\left[e_{m+1/2}^z u_{1,k,l+1/2, m+1/2} - e_{m-1/2}^z u_{1,k,l+1/2,m-1/2]} \right. \\ &+\left[e_{k+1/2}^x u_{3,k+1/2,l+1/2,m} - e_{k-1/2}^x u_{3,k-1/2,l+1/2,m]}\right] , \\ % r_{3,k,l,m+1/2} = ~&s_{3,k,l,m+1/2} + \mathcal{S}_{k,l,m+1/2} E_{3,k,l,m+1/2} \\ &-\left[e_{k+1/2}^x u_{2,k+1/2,l,m+1/2} - e_{k-1/2}^x u_{2,k-1/2,m+1/2]}\right.\\ &+\left[e_{l+1/2}^y u_{1,k,l+1/2,m+1/2} - e_{l-1/2}^y u_{1,k,l-1/2,m+1/2}\right] .\end{split}\]

The weighting of the differences in \(u_1\), etc., may appear strange. The reason is that the differences have been multiplied by the local dual volume. As already mentioned, the dual volume for \(E_{1,k,l,m}\) is shown in Figure 6(b).

For further details of the discretisation see [Muld06] or [Yee66]. The actual meshing is done using discretize (part of the SimPEG-framework). The coordinate system of discretize uses a coordinate system were positive z is upwards.

The method is implemented in a matrix-free manner: the large sparse linear matrix that describes the discretised problem is never explicitly formed, only its action is evaluated on the latest estimate of the solution, thereby reducing storage requirements.

Iterative Solvers

The multigrid method is an iterative (or relaxation) method and shares as such the underlying idea of iterative solvers. We want to solve the linear equation system

()\[A \mathbf{x} = \mathbf{b} ,\]

where \(A\) is the \(n\times n\) system matrix and \(x\) the unknown. If \(v\) is an approximation to \(x\), then we can define two important measures. One is the error \(e\)

()\[\mathbf{e} = \mathbf{x} - \mathbf{v} ,\]

which magnitude can be measured by any standard vector norm, for instance the maximum norm and the Euclidean or 2-norm defined respectively, by

\[\|\mathbf{e}\|_\infty = \max_{1\leq j \leq n}|e_j| \quad \text{and} \quad \|\mathbf{e}\|_{2} = \sqrt{\sum_{j=1}^{n} e_{j}^{2}} .\]

However, as the solution is not known the error cannot be computed either. The second important measure, however, is a computable measure, the residual \(r\) (computed in emg3d.solver.residual())

()\[\mathbf{r} = \mathbf{b} - A\mathbf{v} .\]

Using Equation (27) we can rewrite Equation (26) as

\[A\mathbf{e} = \mathbf{b} - A\mathbf{v} ,\]

from which we obtain with Equation (28) the Residual Equation

()\[A\mathbf{e} = \mathbf{r} .\]

The Residual Correction is given by

()\[\mathbf{x} = \mathbf{v}+\mathbf{e} .\]

The Multigrid Method

Note

If you have never heard of multigrid methods before you might want to read through the Multi-what?-section.

Multigrid is a numerical technique for solving large, often sparse, systems of equations, using several grids at the same time. An elementary introduction can be found in [BrHM00]. The motivation for this approach follows from the observation that it is fairly easy to determine the local, short-range behaviour of the solution, but more difficult to find its global, long-range components. The local behaviour is characterized by oscillatory or rough components of the solution. The slowly varying smooth components can be accurately represented on a coarser grid with fewer points. On coarser grids, some of the smooth components become oscillatory and again can be easily determined.

The following constituents are required to carry out multigrid. First, a sequence of grids is needed. If the finest grid on which the solution is to be found has a constant grid spacing \(h\), then it is natural to define coarser grids with spacings of \(2h\), \(4h\), etc. Let the problem on the finest grid be defined by \(A^h \mathbf{x}^h = \mathbf{b}^h\). The residual is \(\mathbf{r}^h = \mathbf{b}^h - A^h \mathbf{x}^h\) (see the corresponding function emg3d.solver.residual(), and for more details also the function emg3d.core.amat_x()). To find the oscillatory components for this problem, a smoother or relaxation scheme is applied. Such a scheme is usually based on an approximation of \(A^h\) that is easy to invert. After one or more smoothing steps (see the corresponding function emg3d.solver.smoothing()), say \(\nu_1\) in total, convergence will slow down because it is generally difficult to find the smooth, long-range components of the solution. At this point, the problem is mapped to a coarser grid, using a restriction operator \(\tilde{I}^{2h}_h\) (see the corresponding function emg3d.solver.restriction(), and for more details, the functions emg3d.core.restrict_weights() and emg3d.core.restrict(). On the coarse-grid, \(\mathbf{b}^{2h} = \tilde{I}^{2h}_h\mathbf{r}^h\). The problem \(\mathbf{r}^{2h} = \mathbf{b}^{2h} - A^{2h} \mathbf{x}^{2h} = 0\) is now solved for \(\mathbf{x}^{2h}\), either by a direct method if the number of points is sufficiently small or by recursively applying multigrid. The resulting approximate solution needs to be interpolated back to the fine grid and added to the solution. An interpolation operator \(I^h_{2h}\), usually called prolongation in the context of multigrid, is used to update \(\mathbf{x}^h := \mathbf{x}^h + I^h_{2h}\mathbf{x}^{2h}\) (see the corresponding function emg3d.solver.prolongation()). Here \(I^h_{2h}\mathbf{x}^{2h}\) is called the coarse-grid correction. After prolongation, \(\nu_2\) additional smoothing steps can be applied. This constitutes one multigrid iteration.

So far, we have not specified the coarse-grid operator \(A^{2h}\). It can be formed by using the same discretisation scheme as that applied on the fine grid. Another popular choice, \(A^{2h} = \tilde{I}^{2h}_h A^h I^h_{2h}\), has not been considered here. Note that the tilde is used to distinguish restriction of the residual from operations on the solution, because these act on elements of different function spaces.

If multigrid is applied recursively, a strategy is required for moving through the various grids. The simplest approach is the V-cycle shown in Figure 7 for the case of four grids. Here, the same number of pre- and post-smoothing steps is used on each grid, except perhaps on the coarsest. In many cases, the V-cycle does not solve the coarse-grid equations sufficiently well. The W-cycle, shown in Figure 8, will perform better in that case. In a W-cycle, the number of coarse-grid corrections is doubled on subsequent coarser grids, starting with one coarse-grid correction on the finest grid. Because of its cost, it is often replaced by the F-cycle (Figure 9). In the F-cycle, the number of coarse-grid corrections increases by one on each subsequent coarser grid.

V-cycle with \(\nu_1\) pre-smoothing steps and \(\nu_2\) post-smoothing steps. On the coarsest grid, \(\nu_c\) smoothing steps are applied or an exact solver is used. The finest grid has a grid spacing \(h\) and the coarsest \(8h\). A single coarse-grid correction is computed for all grids but the coarsest.

W-cycle with \(\nu_1\) pre-smoothing steps and \(\nu_2\) post-smoothing steps. On each grid except the coarsest, the number of coarse-grid corrections is twice that of the underlying finer grid.

F-cycle with \(\nu_1\) pre-smoothing steps and \(\nu_2\) post-smoothing steps. On each grid except the coarsest, the number of coarse-grid corrections increases by one compared to the underlying finer grid.

One reason why multigrid methods may fail to reach convergence is strong anisotropy in the coefficients of the governing partial differential equation or severely stretched grids (which has the same effect as anisotropy). In that case, more sophisticated smoothers or coarsening strategies may be required. Two strategies are currently implemented, semicoarsening and line relaxation, which can be used on their own or combined. Semicoarsening is when the grid is only coarsened in some directions. Line relaxation is when in some directions the whole gridlines of values are found simultaneously. If slow convergence is caused by just a few components of the solution, a Krylov subspace method can be used to remove them. In this way, multigrid is accelerated by a Krylov method. Alternatively, multigrid might be viewed as a preconditioner for a Krylov method.

Gauss-Seidel

The smoother implemented in emg3d is a Gauss-Seidel smoother. The Gauss-Seidel method solves the linear equation system \(A \mathbf{x} = \mathbf{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\). On the coarsest grid it acts as direct solver, whereas on the finer grid it acts as a smoother with only few iterations.

See the function emg3d.solver.smoothing(), and for more details, the functions emg3d.core.gauss_seidel(), emg3d.core.gauss_seidel_x(), emg3d.core.gauss_seidel_y(), emg3d.core.gauss_seidel_z(), and also emg3d.core.blocks_to_amat().

Choleski factorisation

The actual solver of the system \(A\mathbf{x}=\mathbf{b}\) is a non-standard Cholesky factorisation without pivoting for a symmetric, complex matrix \(A\) tailored to the problem of the multigrid solver, using only the main diagonal and five lower off-diagonals of the banded matrix \(A\). The result is the same as simply using, e.g., numpy.linalg.solve(), but faster for the particular use-case of this code.

See emg3d.core.solve() for more details.

CPU & RAM

The multigrid method is attractive because it shows optimal scaling for both runtime and memory consumption. In the following are a few notes regarding memory and runtime requirements. It also contains information about what has been tried and what still could be tried in order to improve the current code.

Runtime

The gallery contains a script to do some testing with regards to runtime, see the Tools Section. An example output of that script is shown in Figure 10.

Runtime as a function of cell size, which shows nicely the linear scaling of multigrid solvers (using a single thread).

The costliest functions (for big models) are:

• >90 %: emg3d.solver.smoothing() (emg3d.core.gauss_seidel())

• <5 % each, in decreasing importance:

  • emg3d.solver.prolongation() (emg3d.solver.RegularGridProlongator)
  • emg3d.solver.residual() (emg3d.core.amat_x())
  • emg3d.solver.restriction()

Example with 262,144 / 2,097,152 cells (nu_{i,1,c,2}=0,2,1,2; sslsolver=False; semicoarsening=True; linerelaxation=True):

• 93.7 / 95.8 % smoothing
• 3.6 / 2.0 % prolongation
• 1.9 / 1.9 % residual
• 0.6 / 0.4 % restriction

The rest can be ignored. For small models, the percentage of smoothing goes down and of prolongation and restriction go up. But then the modeller is fast anyway.

emg3d.core.gauss_seidel() and emg3d.core.amat_x() are written in numba; jitting emg3d.solver.RegularGridProlongator turned out to not improve things, and many functions used in the restriction are jitted too. The costliest functions (RAM- and CPU-wise) are therefore already written in numba.

Any serious attempt to improve the speed will have to tackle the smoothing itself.

Things which could be tried

• Not much has been tested with the numba-options parallel; prange; and nogil.
• There might be an additional gain by making emg3d.meshes.TensorMesh, emg3d.models.Model, and emg3d.fields.Field instances jitted classes.

Things which have been tried

• One important aspect of the smoothing part is the memory layout. emg3d.core.gauss_seidel() and emg3d.core.gauss_seidel_x() are ideal for F-arrays (loop z-y-x, hence slowest to fastest axis). emg3d.core.gauss_seidel_y() and emg3d.core.gauss_seidel_z(), however, would be optimal for C-arrays. But copying the arrays to C-order and afterwards back is costlier in most cases for both CPU and RAM. The one possible and therefore implemented solution was to swap the loop-order in emg3d.core.gauss_seidel_y().
• Restriction and prolongation information could be saved in a dictionary instead of recomputing it every time. Turns out to be not worth the trouble.
• Rewrite emg3d.RegularGridInterpolator as jitted function, but the iterator approach seems to be better for large grids.

Memory

Most of the memory requirement comes from storing the data itself, mainly the fields (source field, electric field, and residual field) and the model parameters (resistivity, eta, mu). For a big model, they some up; e.g., almost 3 GB for an isotropic model with 256x256x256 cells.

The gallery contains a script to do some testing with regards to the RAM usage, see the Tools Section. An example output of that script is shown in Figure 11.

RAM usage, showing the optimal behaviour of multigrid methods. “Data RAM” is the memory required by the fields (source field, electric field, residual field) and by the model parameters (resistivity; and eta, mu). “MG Base” is for solving one Gauss-Seidel iteration on the original grid. “MG full RAM” is for solving one multigrid F-Cycle.

The theory of multigrid says that in an ideal scenario, multigrid requires 8/7 (a bit over 1.14) the memory requirement of carrying out one Gauss-Seidel step on the finest grid. As can be seen in the figure, for models up to 2 million cells that holds pretty much, afterwards it becomes a bit worse.

However, for this estimation one has to run the model first. Another way to estimate the requirement is by starting from the RAM used to store the fields and parameters. As can be seen in the figure, for big models one is on the save side estimating the required RAM as 1.35 times the storage required for the fields and model parameters.

The figure also shows nicely the linear behaviour of multigrid; for twice the number of cells twice the memory is required (from a certain size onwards).

Attempts at improving memory usage should focus on the difference between the red line (actual usage) and the dashed black line (1.14 x base usage).

CLI interface

Command-line interface for certain specific tasks, such as forward modelling and gradient computation of the misfit function. The command is emg3d, consult the inbuilt help to get started:

emg3d --help

The CLI is driven by command-line parameters and a configuration file. The default configuration file is emg3d.cfg, but another name can be provided as first positional argument to emg3d. Note that arguments provided in the command line overwrite the settings in the configuration file.

Format of the config file

The shown values are the defaults. All values are commented out in this example; remove the comment signs to use them.

# Files
# -----
# If the files are provided without ending the suffix `.h5` will be appended.
# The log has the same name as `output`, but with the suffix `.log`.
[files]
# path = .                   # Path (absolute or relative) to the data
# survey = survey.h5         # Also via `--survey`
# model = model.h5           # Also via `--model`
# output = emg3d_out.h5      # Also via `--output`
# store_simulation = False   # Stores entire simulation in output if True

# Simulation parameters
# ---------------------
# Input parameters for the `Simulation` class, except for `solver_opts`
# (defined in their own section), but including the parameter `min_offset`
# for `compute()`.
[simulation]
# max_workers = 4    # Also via `-n` or `--nproc`.
# gridding = single  # One grid for all sources and frequencies.
# min_offset = 0.0   # Only relevant if `observed=True` (r<r_min set to NaN).

# Solver options
# --------------
# Input parameters for the solver.
# See https://emg3d.readthedocs.io/en/stable/api/emg3d.solver.solve.html
# for a list of all parameters. The only parameters that cannot be provided
# here are grid, model, sfield, efield, and return_info.
#
# Note that currently sslsolver, semicoarsening, and linerelaxation only
# accept True/False through the CLI.
[solver_opts]
# sslsolver = True
# semicoarsening = True
# linerelaxation = True
# verb = 0

# Gridding options
# ----------------
# Input parameters for the automatic gridding.
# See the description of `gridding_opts` and link therein in
# https://emg3d.readthedocs.io/en/stable/api/emg3d.simulations.Simulation.html
# for more details.
#
# List of lists: lists are comma-separated values, lists are separated by
# semi-colons.
#
# One of the limitation of the CLI is that `vector` has to be a string.
[gridding_opts]
# properties =          # list, e.g.: 0.3, 1, 1e5
# center =              # list, e.g.: 0, 0, 0
# cell_number =         # list, e.g.: 8, 16, 32, 64, 128
# min_width_pps =       # list, e.g.: 5, 3, 3
# expand =              # list, e.g.: 0.3, 1e8
# domain =              # list of lists, e.g.: -10000, 10000; None; None
# stretching =          # list of lists, e.g.: None; None; 1.05, 1.5
# min_width_limits =    # list of lists, e.g.: 10, 100; None; 50
# mapping =             # string, e.g.: Resistivity
# vector =              # string, e.g.: xy
# frequency =           # float, e.g.: 1.0
# seasurface =          # float, e.g.: 0.0
# max_buffer =          # float, e.g.: 100000.0
# lambda_factor =       # float, e.g.: 1.0
# verb =                # int, e.g.: 0
# lambda_from_center =  # bool, e.g.: False


# Data
# ----
# Select which sources, receivers, and frequencies of the survey are used. By
# default all data is used. These are comma-separated lists.
[data]
# sources = Tx02, Tx08, Tx14
# receivers = Rx01, Rx10
# frequencies = 0.5, 0.75

Gallery

The gallery with many examples can be found at emsig.github.io/emg3d-gallery.

References

[ArFW00]

Arnold, D. N., R. S. Falk, and R. Winther, 2000, Multigrid in H(div) and H(curl): Numerische Mathematik, 85, 197–217; DOI: 10.1007/PL00005386.

[BrHM00]

Briggs, W., V. Henson, and S. McCormick, 2000, A Multigrid Tutorial, Second Edition: Society for Industrial and Applied Mathematics; DOI: 10.1137/1.9780898719505.

[ClWe01]

Clemens, M., and T. Weiland, 2001, Discrete electromagnetism with the finite integration technique: PIER, 32, 65–87; DOI: 10.2528/PIER00080103.

[Fedo64]

Fedorenko, R. P., 1964, The speed of convergence of one iterative process: USSR Computational Mathematics and Mathematical Physics, 4, 227–235; DOI 10.1016/0041-5553(64)90253-8.

[JoOM06]

Jönsthövel, T. B., C. W. Oosterlee, and W. A. Mulder, 2006, Improving multigrid for 3-D electro-magnetic diffusion on stretched grids: European Conference on Computational Fluid Dynamics; UUID: df65da5c-e43f-47ab-b80d-2f8ee7f35464.

[Hack85]

Hackbusch, W., 1985, Multi-grid methods and applications: Springer, Berlin, Heidelberg, Volume 4 of Springer Series in Computational Mathematics; DOI: 10.1007/978- 3-662-02427-0.

[MoSu94]

Monk, P., and E. Süli, 1994, A convergence analysis of Yee’s scheme on nonuniform grids: SIAM Journal on Numerical Analysis, 31, 393–412; DOI 10.1137/0731021.

[Muld06]

Mulder, W. A., 2006, A multigrid solver for 3D electromagnetic diffusion: Geophysical Prospecting, 54, 633–649; DOI: 10.1111/j.1365-2478.2006.00558.x.

[Muld07]

Mulder, W. A., 2007, A robust solver for CSEM modelling on stretched grids: EAGE Technical Program Expanded Abstracts, D036; DOI 10.3997/2214-4609.201401567.

[Muld08]

Mulder, W. A., 2008, Geophysical modelling of 3D electromagnetic diffusion with multigrid: Computing and Visualization in Science, 11, 29–138; DOI: 10.1007/s00791-007-0064-y.

[Muld11]

Mulder, W. A., 2011, in Numerical Methods, Multigrid: Springer Netherlands, 895–900; DOI 10.1007/978-90-481-8702-7_153.

[MuWS08]

Mulder, W. A., M. Wirianto, and E. C. Slob, 2008, Time-domain modeling of electromagnetic diffusion with a frequency-domain code: Geophysics, 73, F1–F8; DOI: 10.1190/1.2799093.

[PlDM07]

Plessix, R.-É., M. Darnet, and W. A. Mulder, 2007, An approach for 3D multisource, multifrequency CSEM modeling: Geophysics, 72, SM177–SM184; DOI: 10.1190/1.2744234.

[PlMu08]

Plessix, R.-É., and W. A. Mulder, 2008, Resistivity imaging with controlled-source electromagnetic data: depth and data weighting: Inverse Problems, 24, no. 3, 034012; DOI: 10.1088/0266-5611/24/3/034012.

[SlHM10]

Slob, E., J. Hunziker, and W. A. Mulder, 2010, Green’s tensors for the diffusive electric field in a VTI half-space: PIER, 107, 1–20: DOI: 10.2528/PIER10052807.

[Weil77]

Weiland, T., 1977, Eine Methode zur Lösung der Maxwellschen Gleichungen für sechskomponentige Felder auf diskreter Basis: Archiv für Elektronik und Übertragungstechnik, 31, 116–120; pdf: leibniz-publik.de/de/fs1/object/display/bsb00064886_00001.html.

[WeMS19]

Werthmüller, D., W. A. Mulder, and E. C. Slob, 2019, emg3d: A multigrid solver for 3D electromagnetic diffusion: Journal of Open Source Software, 4(39), 1463; DOI: 10.21105/joss.01463.

[Wess91]

Wesseling, P., 1991, An introduction to multigrid methods: John Wiley & Sons. Pure and Applied Mathematics; ISBN: 0-471-93083-0.

[WiMS10]

Wirianto, M., W. A. Mulder, and E. C. Slob, 2010, A feasibility study of land CSEM reservoir monitoring in a complex 3-D model: Geophysical Journal International, 181, 741–755; DOI: 10.1111/j.1365-246X.2010.04544.x.

[WiMS11]

Wirianto, M., W. A. Mulder, and E. C. Slob, 2011, Applying essentially non-oscillatory interpolation to controlled-source electromagnetic modelling: Geophysical Prospecting, 59, 161–175; DOI: 10.1111/j.1365-2478.2010.00899.x.

[Yee66]

Yee, K., 1966, Numerical solution of initial boundary value problems involving maxwell’s equations in isotropic media: IEEE Transactions on Antennas and Propagation, 14, 302–307; DOI: 10.1109/TAP.1966.1138693.

Credits

This project was started by Dieter Werthmüller. Every contributor will be listed here and is considered to be part of «The emg3d Developers»:

• Dieter Werthmüller

Various bits got improved through discussions on Slack at SWUNG and at SimPEG, thanks to both communities. Special thanks to @jokva (general), @banesullivan (visualization), @joferkington (interpolation), and @jcapriot (volume averaging).

Historical credits

The core of emg3d is a complete rewrite and redesign of the multigrid code by Wim A. Mulder ([Muld06], [Muld07], [Muld08], [MuWS08]), developed at Shell and at TU Delft. Various authors contributed to the original code, amongst others, Tom Jönsthövel ([JoOM06]; improved solver for strongly stretched grids), Marwan Wirianto ([WiMS10], [WiMS11]; computation of time-domain data), and Evert C. Slob ([SlHM10]; analytical solutions). The original code was written in Matlab, where the most time- and memory-consuming parts were sped up through mex-files (written in C). It included multigrid with or without BiCGSTAB, VTI resistivity, semicoarsening, and line relaxation; the number of cells had to be powers of two, and coarsening was done only until the first dimension was at two cells. As such it corresponded roughly to emg3d v0.3.0.

See the References in the manual for the full citations and a more extensive list.

Note

This software was initially (till 05/2021) developed at Delft University of Technology (https://www.tudelft.nl) within the Gitaro.JIM project funded through MarTERA as part of Horizon 2020, a funding scheme of the European Research Area (ERA-NET Cofund, https://www.martera.eu).

Changelog

recent versions

v0.16.1: Verbosity & Logging

2021-02-09

• Solve has a new keyword log, which enables to log the solver messages in the returned info dictionary instead of printing them to screen. This is utilized in the CLI and in the Simulation class to log the solver info.
• Survey has a new attribute select, which returns a reduced survey containing the selected sources, receivers, and frequencies.
• CLI:
  • Configuration info is added to output data.
  • Checks now first if all required files and directories exist, and exits gracefully otherwise informing the user. (The default thrown Python errors would be good enough; but user of the CLI interface might not be familiar with Python, so it is better to throw a very simple, clear message.)
  • Log is more verbose with regards to solver (rel. error, time, nr of it.).
• Dipole throws new an error instead of a warning if it received an unknown keyword.
• Various small things with regard to how things are logged or shown on screen.
• Changed all DeprecationWarnings to FutureWarnings, meaning they will be removed in the next release.
• Bug fix with regards to data selection in the CLI; moved to Survey (see above).

v0.16.0: Arbitrarily shaped sources

2021-01-13

• fields.get_source_field:
  • Arbitrarily shaped sources (and therefore also loops) can now be created by providing a src that consists of x-, y-, and z-coordinates of all endpoints of the individual segments.
  • Simple “magnetic dipole” sources can now be created by providing a point dipole ([x, y, z, azm, dip]) and set msrc=True. This will create a square loop of length``x``length m perpendicular to the defined point dipole, hence simulating a magnetic source. Default length is 1 meter.
  • Point dipoles and finite length dipoles were before treated differently. Point dipoles are new converted into finite length dipoles of provided length (default is 1 meter), and treated as finite length dipoles. This is backwards incompatible and means that the source field for point dipoles might not be exactly the same as before. However, in any properly set-up simulation this should have no influence on the result.
  • Bugfix: Fix floating point issue when the smaller coordinate of a finite length dipole source was very close to a node, but not exactly. This is done by rounding the grid locations and source position, and the precision can be controlled via decimals; default is micrometer.
• fields: Values outside the grid in get_receiver and get_receiver_response are new set to NaN’s instead of zeroes. Additionally, the first and last values in each direction of the fields are ignored, to avoid effects form the boundary condition (receivers should not be placed that close to the boundary anyway).
• simulations:
  • Within the automatic gridding the properties are estimated much more conservative now, if not provided: before the log10-average of the last slice in a given direction was used; now it uses the maximum resistivity. This is usually the air value for x/y and positive z. This is very conservative, but avoids that users use too small computational domains in the case of land and shallow marine surveys. The downside is that it heavily over-estimates the required domain in the deep marine case. However, slower but safe is better in this case.
  • New method print_grids, which prints the info of all created grids. This is also used for logging in the CLI interface.
• maps: interp3d takes a new keyword cval, which is passed to map_coordinates.

v0.11.0 - v0.15.3

v0.15.3: Move to EMSiG

2020-12-09

Various small things, mostly related to the automatic meshing.

• New parameter distance for get_origin_widths, as an alternative for domain and vector: distance defines the survey domain as distance from the center. This is then also available in construct_mesh and Simulation, including the CLI.
• Removed precision from skin_depth, wavelength, min_cell_width; all in meshes. It caused problems for high frequencies.
• All data is stored in the Survey, not partly in Survey and partly in Simulation.
• Deprecated collect_classes in io.
• Expanded the what-parameter in the Simulation-class to include properties related to the gradient.
• Moved from github.com/empymod to github.com/emsig.

v0.15.2 : Bugfix deploy II

2020-12-04

• Fixing deploy script with GHA.

v0.15.1 : Bugfix deploy

2020-12-04

Small bugfix release, as v0.15.0 never got deployed.

• Fix CI deploy script.
• Makefile for the most common dev-tasks.

v0.15.0 : discretize restructure

2020-12-04

The package discretize went through a major restructuring with many name changes and consequent deprecations (see below for a list of affected mesh-properties for emg3d). This version updates emg3d to be compatible with discretize>=0.6.0 in the long run. It also means that emg3d will, from emg3d>=0.15.0 onwards, only work with discretize>=0.6.0.

Other notable changes:

• Bug fix re storing/loading synthetics
• Moved from Travis CI to GitHub Actions.

The relevant aliases and deprecations for emg3d are (consult the release notes of discretize for all changes):

Aliases: Aliases (left) remain valid pointers to the new names (right).

• x0 => origin
• nC => n_cells
• vnC => shape_cells
• nN => n_nodes
• vnN => shape_nodes
• nE => n_edges
• nEx => n_edges_x
• nEy => n_edges_y
• nEz => n_edges_z
• vnE => n_edges_per_direction
• vnEx => shape_edges_x
• vnEy => shape_edges_y
• vnEz => shape_edges_z

Deprecations: Deprecated properties (left) raise a deprecation warning and will be removed in the future. Currently, they still work and point to the new names (right).

• hx => h[0]
• hy => h[1]
• hz => h[2]
• nCx => shape_cells[0]
• nCy => shape_cells[1]
• nCz => shape_cells[2]
• nNx => shape_nodes[0]
• nNy => shape_nodes[1]
• nNz => shape_nodes[2]
• vectorNx => nodes_x
• vectorNy => nodes_y
• vectorNz => nodes_z
• vectorCCx => cell_centers_x
• vectorCCy => cell_centers_y
• vectorCCz => cell_centers_z
• vol => cell_volumes

v0.14.3 : Bug fix

2020-11-19

• Bug fix for discretize>=0.6.0.

v0.14.2 : Bug fix

2020-11-18

• Bug fix for Windows affecting good_mg_cell_nr (int32 issue).

v0.14.1 : Bug fix

2020-11-14

• Fix for h5py>=3.0.
• Improved docs re automatic gridding.

v0.14.0 : Automatic gridding

2020-11-07

The simulation class comes new with an automatic gridding functionality, which should make it much easier to compute CSEM data. With that the entire optimization routine was improved too. See the API docs for more info of the relevant implementation.

• simulation:
  • Simulation: New gridding options 'single', 'frequency' 'source', and 'both'; new default is 'single'.
  • compute() takes a new argument, min_offset. If observed=True, it will add Gaussian random noise according to the standard deviation of the data; it will set receivers responses below the minimum offset to NaN.
  • There is no longer a reference model.
  • misfit and gradient can now handle observations with NaN’s.
• survey: A Survey has new attributes standard_error, noise_floor, and relative_error.
• optimize: Completely changed misfit and data-weighting to more sensible functions.
• cli:
  • As a consequence of the changes the data_weight_opts got removed.
  • New sections [data] to select the wanted data and [gridding_opts] for options of the automatic gridding.
  • Section [simulation] has a new parameter min_offset (for creating observed data).
  • Output has a new parameter n_observations if misfit or gradient were called, which is the number of observations that were used to compute the misfit.
• meshes:
  • New functions construct_mesh, get_origin_widths, good_mg_cell_nr and other, smaller helper routines.
  • Deprecated the old meshing routines get_hx_h0, get_cell_numbers, get_stretched_h, get_domain, get_hx; they will be removed in the future.
  • Default of good_mg_cell_nr changed, and the documentation (and verbosity) with regards to «good» number of cells was improved.
• Bug fixes:
  • maps: Fixed the mapping of the gradients (Conductivity is the only mapping that was not affected by this bug).
• Removed deprecated features:
  • models.Model: Removed parameters res_{x;y;z}.
  • io.save: Removed deprecated parameter backend.
  • io.save: Removed default, file extension has to be provided.

v0.13.0 : CLI

2020-09-22

• New Module cli for command-line interaction:

  The command-line interface can currently be used to forward model an entire Simulation, and also to compute the misfit of it with respect to some data and the gradient of the misfit function. See the section “CLI interface” in the documentation for more info.

Note that, while cli (v0.13.0) and optimize (v0.12.0) are implemented, they are still in development and are likely going to change throughout the next two minor releases or so.

• Other changes:
  • solver: Changes in verbosity for emg3d.solve:
    • New default verbosity is 1 (only warnings; before it was 2).
    • Verbosities {-1;0;1} remain unchanged.
    • Verbosities {2;3;4} => {3;4;5}.
    • New verbosity 2: Only shows a one-liner at the end (plus warnings).
  • survey and simulation: to_file and from_file have new a parameter name, to store and load with a particular name instead of the default survey/simulation (useful when storing, e.g., many surveys in one file).
  • survey: stores new also the reference-data; different data (observed, reference) is contained in a data-dict when storing.
  • simulation: takes new a verb parameter.
  • optimize:
    • Gradient now possible for arbitrarily rotated sources and receivers.
    • Falls back to synthetic instead of observed now if reference not found.
  • io: np.bool_ are converted back to bool when loading.
  • Re-arrange, improve, and update documentation.

v0.12.0 : Survey & Simulation

2020-07-25

This is a big release with many new features, and unfortunately not completely backwards compatible. The main new features are the new Survey and Simulation classes, as well as some initial work for optimization (misfit, gradient). Also, a Model can now be a resistivity model, a conductivity model, or the logarithm (natural or base 10) therefore. Receivers can now be arbitrarily rotated, just as the sources. In addition to the existing soft-dependencies empymod, discretize, and h5py there are the new soft-dependencies xarray and tqm; discretize is now much tighter integrated. For the new survey and simulation classes xarray is a required dependency. However, the only hard dependency remain scipy and numba, if you use emg3d purely as a solver. Data reading and writing has new a JSON-backend, in addition to the existing HDF5 and NumPy-backends.

In more detail:

• Modules:
  • surveys (new; requires xarray):
    • Class surveys.Survey, which combines sources, receivers, and data.
    • Class surveys.Dipole, which defines electric or magnetic point dipoles and finite length dipoles.
  • simulations (new; requires xarray; soft-dependency tqdm):
    • Class simulations.Simulation, which combines a survey with a model. A simulation computes the e-field (and h-field) asynchronously using concurrent.futures. This class will include automatic, source- and frequency-dependent gridding in the future. If tqdm is installed it displays a progress bar for the asynchronous computation. Note that the simulation class has still some limitations, consult the class documentation.
  • models:
    • Model instances take new the parameters property_{x;y;z} instead of res_{x;y;z}. The properties can be either resistivity, conductivity, or log_{e;10} thereof. What is actually provided has to be defined with the parameter mapping. By default, it remains resistivity, as it was until now. The keywords res_{x;y;z} are deprecated, but still accepted at the moment. The attributes model.res_{x;y;z} are still available too, but equally deprecated. However, it is no longer possible to assign values to these attributes, which is a backwards incompatible change.
    • A model knows now how to interpolate itself from its grid to another grid (interpolate2grid).
  • maps:
    • New mappings for models.Model instances: The mappings take care of how to transform the investigation variable to conductivity and back, and how it affects its derivative.
    • New interpolation routine edges2cellaverages.
  • fields:
    • Function get_receiver_response (new), which returns the response for arbitrarily rotated receivers.
    • Improvements to Field and SourceField:
      • _sval and _smu0 not stored any longer, derived from _freq.
      • SourceField is now using the copy() and from_dict() from its parents class Field.
  • io:
    • File-format json (new), writes to a hierarchical, plain json file.
    • Deprecated the use of backend, it uses the file extension of fname instead.
    • This means .npz (instead of numpy), .h5 (instead of h5py), and new .json.
    • New parameter collect_classes, which can be used to switch-on collection of the main classes in root-level dictionaries. By default, they are no longer collected (changed).
  • meshes:
    • meshes.TensorMesh new inherits from discretize if installed.
    • Added __eq__ to models.TensorMesh to compare meshes.
  • optimize (new)
    • Functionalities related to inversion (data misfit, gradient, data weighting, and depth weighting). This module is in an early stage, and the API will likely change in the future. Current functions are misfit, gradient (using the adjoint-state method), and data_weighting. These functionalities are best accessed through the Simulation class.
• Dependencies:
  • empymod is now a soft dependency (no longer a hard dependency), only required for utils.Fourier (time-domain modelling).
  • Existing soft dependency discretize is now baked straight into meshes.
  • New soft dependency xarray for the Survey class (and therefore also for the Simulation class and the optimize module).
  • New soft dependency tqdm for nice progress bars in asynchronous computation.
• Deprecations and removals:
  • Removed deprecated functions data_write and data_read.
  • Removed all deprecated functions from utils.
• Miscellaneous:
  • Re-organise API-docs.
  • Much bookkeeping (improve error raising and checking; chaining errors, numpy types, etc).

v0.11.0 : Refactor

2020-05-05

Grand refactor with new internal layout. Mainly splitting-up utils into smaller bits. Most functionalities (old names) are currently retained in utils and it should be mostly backwards compatible for now, but they are deprecated and will eventually be removed. Some previously deprecated functions were removed, however.

• Removed deprecated functions:
  • emg3d.solver.solver (use emg3d.solver.solve instead).
  • Aliases of emg3d.io.data_write and emg3d.io.data_read in emg3d.utils.
• Changes:
  • SourceField has now the same signature as Field (this might break your code if you called SourceField directly, with positional arguments, and not through get_source_field).
  • More functions and classes in the top namespace.
  • Replaced core.l2norm with scipy.linalg.norm, as SciPy 1.4 got the following PR: https://github.com/scipy/scipy/pull/10397 (reason to raise minimum SciPy to 1.4).
  • Increased minimum required versions of dependencies to
    • scipy>=1.4.0 (raised from 1.1, see note above)
    • empymod>=2.0.0 (no min requirement before)
    • numba>=0.45.0 (raised from 0.40)
• New layout
  • njitted -> core.
  • utils split in fields, meshes, models, maps, and utils.
• Bugfixes:
  • Fixed to_dict, from_dict, and copy for the SourceField.
  • Fixed io for SourceField, that was not implemented properly.

v0.8.0 - v0.10.1

v0.10.1 : Zero Source

2020-04-29

• Bug fixes:
  • Checks now if provided source-field is zero, and exists gracefully if so, returning a zero electric field. Until now it failed with a division-by-zero error.
• Improvements:
  • Warnings: If verb=1 it prints a warning in case it did not converge (it finished silently until now).
  • Improvements to docs (figures-scaling; intersphinx).
  • Adjust Fields.pha and Fields.amp in accordance with empymod v2: .pha and .amp are now methods; uses directly empymod.utils.EMArray.
  • Adjust tests for empymod v2 (Fields, Fourier).

v0.10.0 : Data persistence

2020-03-25

• New:
  • New functions emg3d.save and emg3d.load to save and load all sort of emg3d instances. The currently implemented backends are h5py for .h5-files (default, but requires h5py to be installed) and numpy for .npz-files.
  • Classes emg3d.utils.Field, emg3d.utils.Model, and emg3d.utils.TensorMesh have new methods .copy(), .to_dict(), and .from_dict().
  • emg3d.utils.Model: Possible to create new models by adding or subtracting existing models, and comparing two models (+, -, == and !=). New attributes shape and size.
  • emg3d.utils.Model does not store the volume any longer (just vnC).
• Deprecations:
  • Deprecated data_write and data_read.
• Internal and bug fixes:
  • All I/O-related stuff moved to its own file io.py.
  • Change from NUMBA_DISABLE_JIT to use py_func for testing and coverage.
  • Bugfix: emg3d.njitted.restrict did not store the {x;y;z}-field if sc_dir was {4;5;6}, respectively.

v0.9.3 : Sphinx gallery

2020-02-11

• Rename solver.solver to solver.solve; load solve also into the main namespace as emg3d.solve.
• Adjustment to termination criterion for STAGNATION: The current error is now compared to the last error of the same cycle type. Together with this the workaround for sslsolver when called with an initial efield introduced in v0.8.0 was removed.
• Adjustment to utils.get_hx_h0 (this might change your boundaries): The computation domain is now computed so that the distance for the signal travelling from the source to the boundary and back to the most remote receiver is at least two wavelengths away. If this is within the provided domain, then now extra buffer is added around the domain. Additionally, the function has a new parameter max_domain, which is the maximum distance from the center to the boundary; defaults to 100 km.
• New parameter log for utils.grid2grid; if True, then the interpolation is carried out on a log10-scale.
• Change from the notebook-based emg3d-examples-repo to the sphinx-based emg3d-gallery-repo.

v0.9.2 : Complex sources

2019-12-26

• Strength input for get_source_field can now be complex; it also stores now the source location and its strength and moment.
• get_receiver can now take entire Field instances, and returns in that case (fx, fy, fz) at receiver locations.
• Krylov subspace solvers:
  • Solver now finishes in the middle of preconditioning cycles if tolerance is reached.
  • Solver now aborts if solution diverges or stagnates also for the SSL solvers; it fails and returns a zero field.
  • Removed gmres and lgmres from the supported SSL solvers; they do not work nice for this problem. Supported remain bicgstab (default), cgs, and gcrotmk.
• Various small things:
  • New attribute Field.is_electric, so the field knows if it is electric or magnetic.
  • New verb-possibility: verb=-1 is a continuously updated one-liner, ideal to monitor large sets of computations or in inversions.
  • The returned info dictionary contains new keys:
    • runtime_at_cycle: accumulated total runtime at each cycle;
    • error_at_cycle: absolute error at each cycle.
  • Simple __repr__ for TensorMesh, Model, Fourier, Time.
• Bugfixes:
  • Related to get_hx_h0, data_write, printing in Fourier.

v0.9.1 : VolumeModel

2019-11-13

• New class VolumeModel; changes in Model:
  • Model now only contains resistivity, magnetic permeability, and electric permittivity.
  • VolumeModel contains the volume-averaged values eta and zeta; called from within emg3d.solver.solver.
  • Full wave equation is enabled again, via epsilon_r; by default it is set to None, hence diffusive approximation.
  • Model parameters are now internally stored as 1D arrays.
  • An {isotropic, VTI, HTI} initiated model can be changed by providing the missing resistivities.
• Bugfix: Up and till version 0.8.1 there was a bug. If resistivity was set with slices, e.g., model.res[:, :, :5]=1e10, it DID NOT update the corresponding eta. This bug was unintentionally fixed in 0.9.0, but only realised now.
• Various:
  • The log now lists the version of emg3d.
  • PEP8: internal imports now use absolute paths instead of relative ones.
  • Move from conda-channel prisae to conda-forge.
  • Automatic deploy for PyPi and conda-forge.

v0.9.0 : Fourier

2019-11-07

• New routine:
  • emg3d.utils.Fourier, a class to handle Fourier-transform related stuff for time-domain modelling. See the example notebooks for its usage.
• Utilities:
  • Fields and returned receiver-arrays (EMArray) both have amplitude (.amp) and phase (.pha) attributes.
  • Fields have attributes containing frequency-information (freq, smu0).
  • New class SourceField; a subclass of Field, adding vector and v{x,y,z} attributes for the real valued source vectors.
  • The Model is not frequency-dependent any longer and does NOT take a freq-parameter any more (currently it still takes it, but it is deprecated and will be removed in the future).
  • data_write automatically removes _vol from TensorMesh instances and _eta_{x,y,z}, _zeta from Model instances. This makes the archives smaller, and they are not required, as they are simply reconstructed if needed.
• Internal changes:
  • The multigrid method, as implemented, only works for the diffusive approximation. Nevertheless, we always used \sigma-i\omega\epsilon, hence a complex number. This is now changed and \epsilon set to 0, leaving only \sigma.
  • Change time convention from exp(-iwt) to exp(iwt), as used in empymod and commonly in CSEM. Removed the parameter conjugate from the solver, to simplify.
  • Change own private class variables from __ to _.
  • res and mu_r are now checked to ensure they are >0; freq is checked to ensure !=0.
• New dependencies and maintenance:
  • empymod is a new dependency.
  • Travis now checks all the url’s in the documentation, so there should be no broken links down the road. (Check is allowed to fail, it is visual QC.)
• Bugfixes:
  • Fixes to the setuptools_scm-implementation (MANIFEST.in).

v0.8.1 : setuptools_scm

2019-10-22

• Implement setuptools_scm for versioning (adds git hashes for dev-versions).

v0.8.0 : Laplace

2019-10-04

• Laplace-domain computation: By providing a negative freq-value to utils.get_source_field and utils.Model, the computation is carried out in the real Laplace domain s = freq instead of the complex frequency domain s = 2i*pi*freq.
• New meshing helper routines (particularly useful for transient modelling where frequency-dependent/adaptive meshes are inevitable):
  • utils.get_hx_h0 to get cell widths and origin for given parameters including a few fixed interfaces (center plus two, e.g. top anomaly, sea-floor, and sea-surface).
  • utils.get_cell_numbers to get good values of number of cells for given primes.
• Speed-up njitted.volume_average significantly thanks to @jcapriot.
• Bugfixes and other minor things:
  • Abort if l2-norm is NaN (only works for MG).
  • Workaround for the case where a sslsolver is used together with a provided initial efield.
  • Changed parameter rho to res for consistency reasons in utils.get_domain.
  • Changed parameter h_min to min_width for consistency reasons in utils.get_stretched_h.

v0.1.0 - v0.7.1

v0.7.1 : JOSS article

2019-07-17

• Version of the JOSS article, https://doi.org/10.21105/joss.01463 .
• New function utils.grid2grid to move from one grid to another. Both functions (utils.get_receiver and utils.grid2grid) can be used for fields and model parameters (with or without extrapolation). They are very similar, the former taking coordinates (x, y, z) as new points, the latter one another TensorMesh instance.
• New jitted function njitted.volume_average for interpolation using the volume-average technique.
• New parameter conjugate in solver.solver to permit both Fourier transform conventions.
• Added exit_status and exit_message to info_dict.
• Add section Related ecosystem to documentation.

v0.7.0 : H-field

2019-07-05

• New routines:
  • utils.get_h_field: Small routine to compute the magnetic field from the electric field using Faraday’s law.
  • utils.get_receiver: Small wrapper to interpolate a field at receiver positions. Added 3D spline interpolation; is the new default.
• Re-implemented the possibility to define isotropic magnetic permeabilities in utils.Model. Magnetic permeability is not tri-axially included in the solver currently; however, it would not be too difficult to include if there is a need.
• CPU-graph added on top of RAM-graph.
• Expand utils.Field to work with pickle/shelve.
• Jit np.linalg.norm (njitted.l2norm).
• Use scooby (soft dependency) for versioning, rename Version to Report (backwards incompatible).
• Bug fixes:
  • Small bugfix introduced in ebd2c9d5: sc_cycle and lr_cycle was not updated any longer at the end of a cycle (only affected sslsolver=True.
  • Small bugfix in utils.get_hx.

v0.6.2 : CPU & RAM

2019-06-03

Further speed and memory improvements:

• Add CPU & RAM-page to documentation.
• Change loop-order from x-z-y to z-x-y in Gauss-Seidel smoothing with line relaxation in y-direction. Hence reversed lexicographical order. This results in a significant speed-up, as x is the fastest changing axis.
• Move total residual computation from solver.residual into njitted.amat_x.
• Simplifications in utils:
  • Simplify utils.get_source_field.
  • Simplify utils.Model.
  • Removed unused timing-stuff from early development.

v0.6.1 : Memory

2019-05-28

Memory and speed improvements:

• Only compute residual and l2-norm when absolutely necessary.
• Inplace computations for np.conjugate in solver.solver and np.subtract in solver.residual.

v0.6.0 : RegularGridInterpolator

2019-05-26

• Replace scipy.interpolate.RegularGridInterpolator with a custom tailored version of it (class:emg3d.solver.RegularGridProlongator); results in twice as fast prolongation.
• Simplify the fine-grid computation in prolongation without using gridE*; memory friendlier.
• Submission to JOSS.
• Add Multi-what?-page to documentation.
• Some major refactoring, particularly in solver.
• Removed discretize as hard dependency.
• Rename rdir and ldir (and related p*dir; *cycle) to the more descriptive sc_dir and lr_dir.

v0.5.0 : Accept any grid size

2019-05-01

• First open-source version.
• Include RTD, Travis, Coveralls, Codacy, and Zenodo. No benchmarks yet.
• Accepts now any grid size (warns if a bad grid size for MG is provided).
• Coarsens now to the lowest level of each dimension, not only to the coarsest level of the smallest dimension.
• Combined restrict_rx, restrict_ry, and restrict_rz to restrict.
• Improve speed by passing pre-allocated arrays to jitted functions.
• Store res_y, res_z and corresponding eta_y, eta_z only if res_y, res_z were provided in initial call to utils.model.
• Change zeta to v_mu_r.
• Include rudimentary TensorMesh-class in utils; removes hard dependency on discretize.
• Bugfix: Take a provided efield into account; don’t return if provided.

v0.4.0 : Cholesky

2019-03-29

• Use solve_chol for everything, remove solve_zlin.
• Moved mesh.py and some functionalities from solver.py into utils.py.
• New mesh-tools. Should move to discretize eventually.
• Improved source generation tool. Might also move to discretize.
• printversion is now included in utils.
• Many bug fixes.
• Lots of improvements to tests.
• Lots of improvements to documentation. Amongst other, moved docs from __init__.py into the docs rst.

v0.3.0 : Semicoarsening

2019-01-18

• Semicoarsening option.
• Number of cells must still be 2^n, but n can be different in the x-, y-, and z-directions.
• Many other iterative solvers from scipy.sparse.linalg can be used. It seems to work fine with the following methods:
  • scipy.sparse.linalg.bicgstab(): BIConjugate Gradient STABilize;
  • scipy.sparse.linalg.cgs(): Conjugate Gradient Squared;
  • scipy.sparse.linalg.gmres(): Generalized Minimal RESidual;
  • scipy.sparse.linalg.lgmres(): Improvement of GMRES using alternating residual vectors;
  • scipy.sparse.linalg.gcrotmk(): GCROT: Generalized Conjugate Residual with inner Orthogonalization and Outer Truncation.
• The SciPy-solver or MG can be used all in combination or on its own, hence only MG, SciPy-solver with MG preconditioning, only SciPy-solver.

v0.2.0 : Line relaxation

2019-01-14

• Line relaxation option.

v0.1.0 : Initial

2018-12-28

• Standard multigrid with or without BiCGSTAB.
• Tri-axial anisotropy.
• Number of cells must be 2^n, and n has to be the same in the x-, y-, and z-directions.

Maintainers Guide

Making a release

1. Update CHANGELOG.rst.
2. Push it to GitHub, create a release tagging it.
3. Tagging it on GitHub will automatically deploy it to PyPi, which in turn will create a PR for the conda-forge feedstock. Merge that PR.
4. Check that:

• PyPi deployed;
• conda-forge deployed;
• Zenodo minted a DOI;
• emg3d.rtfd.io created a tagged version.

Useful things

• If there were changes to README, check it with:

  python setup.py --long-description | rst2html.py --no-raw > index.html
  

• If unsure, test it first on testpypi (requires ~/.pypirc):

  ~/anaconda3/bin/twine upload dist/* -r testpypi
  

• If unsure, test the test-pypi for conda if the skeleton builds:

  conda skeleton pypi --pypi-url https://test.pypi.io/pypi/ emg3d
  

• If it fails, you might have to install python3-setuptools:

  sudo apt install python3-setuptools
  

CI

Automatic bits

• Testing on Github Actions includes:
  • Tests using pytest
  • Linting / code style with pytest-flake8
  • Ensure all http(s)-links work (sphinx linkcheck)
• Line-coverage with pytest-cov on Coveralls
• Code-quality on Codacy
• Manual on ReadTheDocs
• DOI minting on Zenodo

Manual things

• Benchmarks with Airspeed Velocity (asv)
• Gallery in emg3d-gallery (sphinx-gallery)

Automatically deploys if tagged

• PyPi
• conda -c conda-forge

Solver

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.

emg3d.solver Module

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

Functions

solve(grid, model, sfield[, efield, cycle, …])	Solver for 3D CSEM data with tri-axial electrical anisotropy.
multigrid(grid, model, sfield, efield, var, …)	Multigrid solver for 3D controlled-source electromagnetic (CSEM) data.
smoothing(grid, model, sfield, efield, nu, …)	Reducing high-frequency error by smoothing.
restriction(grid, model, sfield, residual, …)	Downsampling of grid, model, and fields to a coarser grid.
prolongation(grid, efield, cgrid, cefield, …)	Interpolating the electric field from coarse grid to fine grid.
residual(grid, model, sfield, efield[, norm])	Computing the residual.
krylov(grid, model, sfield, efield, var)	Krylov Subspace iterative solver for 3D CSEM data.

Classes

MGParameters(verb, cycle, sslsolver, …)	Collect multigrid solver settings.
RegularGridProlongator(x, y, cxy)	Prolongate field from coarse to fine grid.

emg3d.core Module

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

Functions

amat_x(rx, ry, rz, ex, ey, ez, eta_x, eta_y, …)	Residual without or with source term.
blocks_to_amat(amat, bvec, middle, left, …)	Insert middle, left, and rhs into main arrays amat and bvec.
gauss_seidel(ex, ey, ez, sx, sy, sz, eta_x, …)	Gauss-Seidel method.
gauss_seidel_x(ex, ey, ez, sx, sy, sz, …)	Gauss-Seidel method with line relaxation in x-direction.
gauss_seidel_y(ex, ey, ez, sx, sy, sz, …)	Gauss-Seidel method with line relaxation in y-direction.
gauss_seidel_z(ex, ey, ez, sx, sy, sz, …)	Gauss-Seidel method with line relaxation in z-direction.
restrict(crx, cry, crz, rx, ry, rz, wx, wy, …)	Restriction of residual from fine to coarse grid.
restrict_weights(vectorN, vectorCC, h, …)	Restriction weights for the coarse-grid correction operator.
solve(amat, bvec)	Solve A x = b using a non-standard Cholesky factorisation.

Meshes, Models, and Fields

emg3d.meshes Module

Everything related to meshes appropriate for the multigrid solver.

Functions

construct_mesh(frequency, properties, center)	Return a TensorMesh for given parameters.
get_origin_widths(frequency, properties, center)	Return origin and cell widths for given parameters.
skin_depth(frequency, conductivity[, mu])	Return skin depth for provided frequency and conductivity.
wavelength(skin_depth)	Return the wavelength.
good_mg_cell_nr([max_nr, max_prime, min_div])	Returns ‘good’ cell numbers for the multigrid method.
min_cell_width(skin_depth[, pps, limits])	Return the minimum cell width.
get_hx_h0(freq, res, domain[, fixed, …])	Return cell widths and origin for given parameters.
get_cell_numbers(max_nr[, max_prime, min_div])
get_stretched_h(min_width, domain, nx[, x0, …])	Return cell widths for a stretched grid within the domain.
get_domain([x0, freq, res, limits, …])	Get domain extent and minimum cell width as a function of skin depth.
get_hx(alpha, domain, nx, x0[, resp_domain])	Return cell widths for given input.

Classes

TensorMesh([h, origin])

A slightly modified discretize.TensorMesh.

emg3d.models Module

Everything to create model-properties for the multigrid solver.

Classes

Model(grid[, property_x, property_y, …])	Create a model instance.
VolumeModel(grid, model, sfield)	Return a volume-averaged version of provided model.

emg3d.maps Module

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

Functions

grid2grid(grid, values, new_grid[, method, …])	Interpolate values located on grid to new_grid.
interp3d(points, values, new_points, method, …)	Interpolate values in 3D either linearly or with a cubic spline.
edges2cellaverages(ex, ey, ez, vol, out_x, …)	Interpolate fields defined on edges to volume-averaged cell values.

Classes

MapConductivity()	Maps σ to computational variable σ (conductivity).
MapLgConductivity()	Maps log_10(σ) to computational variable σ (conductivity).
MapLnConductivity()	Maps log_e(σ) to computational variable σ (conductivity).
MapResistivity()	Maps ρ to computational variable σ (conductivity).
MapLgResistivity()	Maps log_10(ρ) to computational variable σ (conductivity).
MapLnResistivity()	Maps log_e(ρ) to computational variable σ (conductivity).

emg3d.fields Module

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

Functions

get_source_field(grid, src, freq[, …])	Return the source field.
get_receiver(grid, values, coordinates[, …])	Return values corresponding to grid at coordinates.
get_receiver_response(grid, field, rec)	Return the field (response) at receiver coordinates.
get_h_field(grid, model, field)	Return magnetic field corresponding to provided electric field.

Classes

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

Surveys and Simulations

emg3d.surveys Module

A survey stores a set of sources, receivers, and the measured data.

Classes

Survey(name, sources, receivers, frequencies)	Create a survey with sources, receivers, and data.
Dipole(name, coordinates[, electric])	Finite length dipole or point dipole.
PointDipole(name, xco, yco, zco, azm, dip, …)	Infinitesimal small electric or magnetic point dipole.

emg3d.simulations Module

A simulation is the computation (modelling) of electromagnetic responses of a resistivity (conductivity) model for a given survey.

In its heart, emg3d is a multigrid solver for 3D electromagnetic diffusion with tri-axial electrical anisotropy. However, it contains most functionalities to also act as a modeller. The simulation module combines all these things by combining surveys with computational meshes and fields and providing high-level, specialised modelling routines.

Functions

expand_grid_model(grid, model, expand, interface)	Expand model and grid according to provided parameters.
estimate_gridding_opts(gridding_opts, grid, …)	Estimate parameters for automatic gridding.

Classes

Simulation(name, survey, grid, model[, …])

Create a simulation for a given survey on a given model.

Optimize

emg3d.optimize Module

Inversion

Functionalities related to optimization (inversion), e.g., misfit function, gradient of the misfit function, or data- and depth-weighting.

Currently it follows the implementation of [PlMu08], using the adjoint-state technique for the gradient.

Functions

gradient(simulation)	Compute the discrete gradient using the adjoint-state method.
misfit(simulation)	Return the misfit function.

I/O and Utils

emg3d.io Module

Utility functions for writing and reading data.

Functions

save(fname, **kwargs)	Save surveys, meshes, models, fields, and more to disk.
load(fname, **kwargs)	Load meshes, models, fields, and other data from disk.

emg3d.utils Module

Utility functions for the multigrid solver.

Classes

Fourier(time, fmin, fmax[, signal, ft, ftarg])	Time-domain CSEM computation.
Time()	Class for timing (now; runtime).
Report([add_pckg, ncol, text_width, sort])	Print date, time, and version information.
EMArray	Create an EM-ndarray: add amplitude <amp> and phase <pha> methods.

Command Line Interface

Functions related to the command-line interface (CLI) of emg3d.

Consult the CLI interface section in the documentation for more information.

emg3d.cli.main Module

Entry point for the command-line interface (CLI).

Functions

main([args])

Parsing command line inputs of CLI interface.

emg3d.cli.parser Module

Parser for the configuration file of the command-line interface.

Functions

parse_config_file(args_dict)

Read and parse the configuration file and set defaults.

emg3d.cli.run Module

Functions that actually call emg3d within the CLI interface.

Functions

check_files(cfg, term)	Ensure all paths and files exist.
initiate_logger(cfg, runtime, verb)	Initiate logger for CLI of emg3d.
simulation(args_dict)	Run emg3d invoked by CLI.