emg3d

Version: 0.12.0 ~ Date: 25 July 2020

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

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-2020 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, res_x=1.5, res_y=1.8, res_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=3)

:: emg3d START :: 15:24:40 :: v0.9.1

   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.

Tipps 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 grid 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 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^{0, 1, ..., 9}\): 2, 4, 8, 16, 32, 64, 128, 256, 512, 1024,
    • \(3·2^{0, 1, ..., 8}\): 3, 6, 12, 24, 48, 96, 192, 384, 768,
    • \(5·2^{0, 1, ..., 7}\): 5, 10, 20, 40, 80, 160, 320, 640,
    • \(7·2^{0, 1, ..., 7}\): 7, 14, 28, 56, 112, 224, 448, 896,

    and preference decreases from top to bottom row. Good grid sizes in sequential order: 2, 3, 4, 5, 6, 7, 8, 10, 12, 14, 16, 20, 24, 28, 32, 40, 48, 56, 64, 80, 96, 112, 128, 160, 192, 224, 256, 320, 384, 448, 512, 640, 768, 896, 1024.

  • 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. As a rule of thumb:

    • No stretching: Multigrid (MG);
    • Moderate stretching (< 1.04): BiCGSTAB with MG as pre-conditioner;
    • Strong stretching (> 1.04): BicGSTAB with MG as preconditioner and line relaxation/semicoarsening.

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/empymod 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.

The structure of emg3d is:

  • solver: These are the main routines, the flow of the multigrid method;
  • njited: The expensive parts (computation, memory) are here in jitted functions; and
  • utils: Some helper routines.

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 empymod/emg3d-asv-repository. The results of my machine can be found in the empymod/emg3d-bench, its rendered version at empymod.github.io/emg3d-asv.

License

Copyright 2018-2020 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

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

Explanation of smoother

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.

Fine to coarse grid schematic

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.

Schematic of semicoarsening

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

Schematic of line relaxation

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

Staggered grid sketches.

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

Dual volume sketches.

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

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

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

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

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:

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

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

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

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

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

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

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

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:

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 Travis, 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
Automatically deploys if tagged

Main solver routine

emg3d.solver.solve(grid, model, sfield, efield=None, cycle='F', sslsolver=False, semicoarsening=False, linerelaxation=False, verb=2, **kwargs)[source]

Solver for 3D CSEM data with tri-axial electrical anisotropy.

The principal solver of emg3d is using the multigrid method as presented in [Muld06]. Multigrid can be used as a standalone solver, or as a preconditioner for an iterative solver from the scipy.sparse.linalg-library, e.g., scipy.sparse.linalg.bicgstab(). Alternatively, these Krylov subspace solvers can also be used without multigrid at all. See the cycle and sslsolver parameters.

Implemented are the F-, V-, and W-cycle schemes for multigrid (cycle parameter), and the amount of smoothing steps (initial smoothing, pre-smoothing, coarsest-grid smoothing, and post-smoothing) can be set individually (nu_init, nu_pre, nu_coarse, and nu_post, respectively). The maximum level of coarsening can be restricted with the clevel parameter.

Semicoarsening and line relaxation, as presented in [Muld07], are implemented, see the semicoarsening and linerelaxation parameters. Using the BiCGSTAB solver together with multigrid preconditioning with semicoarsening and line relaxation is slow but generally the most robust. Not using BiCGSTAB nor semicoarsening nor line relaxation is fast but may fail on stretched grids.

Parameters:
grid : emg3d.meshes.TensorMesh

The grid. See emg3d.meshes.TensorMesh.

model : emg3d.models.Model

The model. See emg3d.models.Model.

sfield : emg3d.fields.SourceField

The source field. See emg3d.fields.get_source_field().

efield : emg3d.fields.Field, optional

Initial electric field. It is initiated with zeroes if not provided. A provided efield MUST have frequency information (initiated with emg3d.fields.Field(..., freq)).

If an initial efield is provided nothing is returned, but the final efield is directly put into the provided efield.

If an initial field is provided and a sslsolver is used, then it first carries out one multigrid cycle without semicoarsening nor line relaxation. The sslsolver is at times unstable with an initial guess, carrying out one MG cycle helps to stabilize it.

cycle : str; optional.

Type of multigrid cycle. Default is ‘F’.

  • ‘V’: V-cycle, simplest version;
  • ‘W’: W-cycle, most expensive version;
  • ‘F’: F-cycle, sort of a compromise between ‘V’ and ‘W’;
  • None: Does not use multigrid, only sslsolver.

If None, sslsolver must be provided, and the sslsolver will be used without multigrid pre-conditioning.

Comparison of V (left), F (middle), and W (right) cycles for the case of four grids (three relaxation and prolongation steps):

 h_
2h_   \    /   \          /   \            /
4h_    \  /     \    /\  /     \    /\    /
8h_     \/       \/\/  \/       \/\/  \/\/
sslsolver : str, optional

A scipy.sparse.linalg-solver, to use with MG as pre-conditioner or on its own (if cycle=None). Default is False.

Current possibilities:

It does currently not work with ‘cg’, ‘bicg’, ‘qmr’, and ‘minres’ for various reasons (e.g., some require rmatvec in addition to matvec).

semicoarsening : int; optional

Semicoarsening. Default is False.

  • True: Cycling over 1, 2, 3.
  • 0 or False: No semicoarsening.
  • 1: Semicoarsening in x direction.
  • 2: Semicoarsening in y direction.
  • 3: Semicoarsening in z direction.
  • Multi-digit number containing digits from 0 to 3. Multigrid will cycle over these values, e.g., semicoarsening=1213 will cycle over [1, 2, 1, 3].
linerelaxation : int; optional

Line relaxation. Default is False.

This parameter is not respected on the coarsest grid, except if it is set to 0. If it is bigger than zero line relaxation on the coarsest grid is carried out along all dimensions which have more than 2 cells.

  • True: Cycling over [4, 5, 6].
  • 0 or False: No line relaxation.
  • 1: line relaxation in x direction.
  • 2: line relaxation in y direction.
  • 3: line relaxation in z direction.
  • 4: line relaxation in y and z directions.
  • 5: line relaxation in x and z directions.
  • 6: line relaxation in x and y directions.
  • 7: line relaxation in x, y, and z directions.
  • Multi-digit number containing digits from 0 to 7. Multigrid will cycle over these values, e.g., linerelaxation=1213 will cycle over [1, 2, 1, 3].

Note: Smoothing is generally done in lexicographical order, except for line relaxation in y direction; the reason is speed (memory access).

verb : int; optional

Level of verbosity (the higher the more verbose). Default is 2.

  • 0: Print nothing.
  • 1: Print warnings.
  • 2: Print runtime and information about the method.
  • 3: Print additional information for each MG-cycle.
  • 4: Print everything (slower due to additional error computations).
  • -1: Print one-liner (dynamically updated).
**kwargs : Optional solver options:

  • tol : float

    Convergence tolerance. Default is 1e-6.

    Iterations stop as soon as the norm of the residual has decreased by this factor, relative to the residual norm obtained for a zero electric field.

  • maxit : int

    Maximum number of multigrid iterations. Default is 50.

    If sslsolver is used, this applies to the sslsolver.

    In the case that multigrid is used as a pre-conditioner for the sslsolver, the maximum iteration for multigrid is defined by the maximum length of the linerelaxation and semicoarsening-cycles.

  • nu_init : int

    Number of initial smoothing steps, before MG cycle. Default is 0.

  • nu_pre : int

    Number of pre-smoothing steps. Default is 2.

  • nu_coarse : int

    Number of smoothing steps on coarsest grid. Default is 1.

  • nu_post : int

    Number of post-smoothing steps. Default is 2.

  • clevel : int

    The maximum coarsening level can be different for each dimension and is, by default, automatically determined (clevel=-1). The parameter clevel can be used to restrict the maximum coarsening level in any direction by its value. Default is -1.

  • return_info : bool

    If True, a dictionary is returned with runtime info (final norm and number of iterations of MG and the sslsolver).
Returns:
efield : emg3d.fields.Field

Resulting electric field. Is not returned but replaced in-place if an initial efield was provided.

info_dict : dict

Dictionary with runtime info; only if return_info=True.

Keys:

  • exit: Exit status, 0=Success, 1=Failure;
  • exit_message: Exit message, check this if exit=1;
  • abs_error: Absolute error;
  • rel_error: Relative error;
  • ref_error: Reference error [norm(sfield)];
  • tol: Tolerance (abs_error<ref_error*tol);
  • it_mg: Number of multigrid iterations;
  • it_ssl: Number of SSL iterations;
  • time: Runtime (s).
  • runtime_at_cycle: Runtime after each cycle (s).
  • error_at_cycle: Absolute error after each cycle.

Examples

>>> import emg3d
>>> import numpy as np
>>> # Create a simple grid, 8 cells of length 1 in each direction,
>>> # starting at the origin.
>>> grid = emg3d.meshes.TensorMesh(
>>>         [np.ones(8), np.ones(8), np.ones(8)],
>>>         x0=np.array([0, 0, 0]))
>>> # The model is a fullspace with tri-axial anisotropy.
>>> model = emg3d.models.Model(grid, res_x=1.5, res_y=1.8, res_z=3.3)
>>> # The source is a x-directed, horizontal dipole at (4, 4, 4)
>>> # with a frequency of 10 Hz.
>>> sfield = emg3d.fields.get_source_field(
>>>         grid, src=[4, 4, 4, 0, 0], freq=10)
>>> # Compute the electric signal.
>>> efield = emg3d.solve(grid, model, sfield, verb=3)
>>> # Get the corresponding magnetic signal.
>>> hfield = emg3d.fields.get_h_field(grid, model, efield)
.
:: emg3d START :: 10:27:25 :: v0.9.1
.
   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  :   8 x   8 x   8     => 512 cells
   Coarsest grid  :   2 x   2 x   2     => 8 cells
   Coarsest level :   2 ;   2 ;   2
.
   [hh:mm:ss]  rel. error                  [abs. error, last/prev]   l s
.
       h_
      2h_ \    /
      4h_  \/\/
.
   [10:27:25]   2.284e-02  after   1 F-cycles   [1.275e-06, 0.023]   0 0
   [10:27:25]   1.565e-03  after   2 F-cycles   [8.739e-08, 0.069]   0 0
   [10:27:25]   1.295e-04  after   3 F-cycles   [7.232e-09, 0.083]   0 0
   [10:27:25]   1.197e-05  after   4 F-cycles   [6.685e-10, 0.092]   0 0
   [10:27:25]   1.233e-06  after   5 F-cycles   [6.886e-11, 0.103]   0 0
   [10:27:25]   1.415e-07  after   6 F-cycles   [7.899e-12, 0.115]   0 0
.
   > CONVERGED
   > MG cycles        : 6
   > Final rel. error : 1.415e-07
.
:: emg3d END   :: 10:27:25 :: runtime = 0:00:00

Code

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

emg3d.meshes Module

meshes – Discretization

Everything related to meshes appropriate for the multigrid solver.

Functions
get_hx_h0(freq, res, domain[, fixed, …]) Return cell widths and origin for given parameters.
get_cell_numbers(max_nr[, max_prime, min_div]) Returns ‘good’ cell numbers for the multigrid method.
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, x0) A slightly modified discretize.TensorMesh.

emg3d.models Module

models – Earth properties

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.fields Module

fields – Electric and magnetic fields

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

Functions
get_source_field(grid, src, freq[, strength]) 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.

emg3d.io Module

io – I/O utilities

Utility functions for writing and reading data.

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

emg3d.surveys Module

survey – Surveys

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

simulation – Model a survey

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.

Classes
Simulation(name, survey, grid, model[, …]) Create a simulation for a given survey on a given model.

emg3d.maps Module

maps – Interpolation routines

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.optimize Module

optimize – 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.
data_weighting(simulation) Return weighted residual.

emg3d.utils Module

utils – Utilities

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.

emg3d.solver Module

solver – Multigrid solver

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

Functions
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

core – Number crunching

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.