Getting started#

There are different usages of emg3d. As an end user you might be interested in the Simulations / High-level usage: define your survey and your model, and simulate electromagnetic data for them. As a developer you might be interested in the Solver-level usage or Command-line interface (CLI-level), to implement emg3d for instance as a forward modelling kernel in you inversion code. And then there is also Time-domain modelling.

Here we show first an introductory example of the Simulations / High-level usage. Further down you find an overview and more detailed explanation of the different Usage levels.

Basic Example#

The following is a very basic example. To see some more realistic models have a look at the gallery. First, we load emg3d, numpy, and matplotlib. Note that this example requires that you have installed discretize and xarray as well.

In [1]: import emg3d
   ...: import numpy as np
   ...: import matplotlib as mpl

One of the elementary ingredients for modelling is a model, in our case a resistivity or conductivity model. emg3d is not a model builder. We just construct here a very simple dummy model. You can see more realistic models in the models section of the gallery.


We first build a simple four-quadrant grid, centred around the origin.

In [2]: hx = np.ones(2)*5000
   ...: grid = emg3d.TensorMesh(
   ...:         [[5000, 5000], [10000], [5000, 5000]], origin='CCC')
   ...: grid  # QC

  TensorMesh: 4 cells

                      MESH EXTENT             CELL WIDTH      FACTOR
  dir    nC        min           max         min       max      max
  ---   ---  ---------------------------  ------------------  ------
   x      2     -5,000.00      5,000.00  5,000.00  5,000.00    1.00
   y      1     -5,000.00      5,000.00 10,000.00 10,000.00    1.00
   z      2     -5,000.00      5,000.00  5,000.00  5,000.00    1.00


We use the constructed grid to create a simple resistivity model:

In [3]: model = emg3d.Model(
   ...:     grid=grid,
   ...:     property_x=[1, 10, 0.3, 0.3],
   ...:     mapping='Resistivity'
   ...: )
   ...: model  # QC
Out[3]: Model: resistivity; isotropic; 2 x 1 x 2 (4)

We defined an isotropic model in terms of resistivities, but through the mapping parameter one can also define a model in terms of conductivities or the logarithms thereof.

Quick check how the model looks like:

In [4]: fig, ax = plt.subplots()
   ...: cf = ax.pcolormesh(
   ...:     grid.cell_centers_x/1e3,
   ...:     grid.cell_centers_z/1e3,
   ...:     model.property_x[:, 0, :].T,
   ...:     shading='nearest',
   ...:     norm=mpl.colors.LogNorm(),
   ...: )
   ...: fig.colorbar(cf)
   ...: ax.set_title(r'Depth slice ($\Omega$ m)');
   ...: ax.set_xlabel('Easting (km)');
   ...: ax.set_ylabel('Depth (km)');

So we have an upper halfspace of 0.3 Ohm.m, a lower-left quadrant of 1 Ohm.m, and a lower-right quadrant of 10 Ohm.m.


Now that we have a model we need to define our survey. Currently there are three source types implemented,

and two receiver types,

We are going to define a simple survey with an electric dipole source and a line of electric point receivers.

In [5]: source = emg3d.TxElectricDipole(
   ...:     coordinates=(-3000, 0, 0, 0, 0)  # x, y, z, azimuth, elevation
   ...: )
   ...: offsets = np.linspace(-2000, 3000, 21)
   ...: receivers = emg3d.surveys.txrx_coordinates_to_dict(
   ...:     emg3d.RxElectricPoint,
   ...:     coordinates=(offsets, 0, 0, 0, 0),  # x, y, z, azimuth, elevation
   ...: )
   ...: survey = emg3d.Survey(
   ...:     sources=source,
   ...:     receivers=receivers,
   ...:     frequencies=1.0,       # Hz
   ...: )
   ...: survey  # QC
:: Survey ::

Dimensions:   (src: 1, rec: 21, freq: 1)
  * src       (src) <U6 'TxED-1'
  * rec       (rec) <U7 'RxEP-01' 'RxEP-02' 'RxEP-03' ... 'RxEP-20' 'RxEP-21'
  * freq      (freq) <U3 'f-1'
Data variables:
    observed  (src, rec, freq) complex128 (nan+nanj) (nan+nanj) ... (nan+nanj)
    noise_floor:     None
    relative_error:  None


Now that we have a model and a survey we can define our simulation:

In [6]: sim = emg3d.Simulation(
   ...:     survey=survey,
   ...:     model=model,
   ...: )
   ...: sim  # QC
:: Simulation ::

- Survey: 1 sources; 21 receivers; 1 frequencies
- Model: resistivity; isotropic; 2 x 1 x 2 (4)
- Gridding: Single grid for all sources and frequencies; 96 x 48 x 64 (294,912)

From the output we see that we defined a survey with one source, 21 receivers, and one frequency. We see that we have a model, defined as resistivity, which consists of four cells. And we see that the simulation created a computational grid of 96x48x64 cells.

A simulation takes many input parameters, and you should really read the API reference of emg3d.simulations.Simulation. Particularly important are the gridding parameters, and as such it is highly recommended to also read emg3d.meshes.construct_mesh. The simulation will try its best to create suitable computational grids. However, these are not necessarily the best ones, and the user has to pay attention to the grid creation. Also, the default grids will most likely be quite big, to be on the safe side. Providing suitable user inputs can significantly reduce the grid sizes and therefore computational time!

So lets check in more detail the computational grid it created:

In [7]: comp_model = sim.get_model(source='TxED-1', frequency=1.0)
   ...: comp_model.grid

  TensorMesh: 294,912 cells

                      MESH EXTENT             CELL WIDTH      FACTOR
  dir    nC        min           max         min       max      max
  ---   ---  ---------------------------  ------------------  ------
   x     96     -6,939.37     13,615.81     91.89  3,045.59    1.42
   y     48    -11,286.30     11,286.30     91.89  3,045.59    1.42
   z     64    -13,651.53      2,245.90     91.89  1,866.31    1.21

We can see that the grid extends further in the directions where there are higher resistivities. For instance, in positive z we have 0.3 Ohm.m, so the grid does only extend roughly to +2.2 km. In positive x and z as well as in positive and negative y we have 10 Ohm.m, so the grid goes to over 10 km.

Computing the fields is now a simple command,

In [8]: sim.compute()


Let’s plot the electric field at receiver locations (the responses):

In [9]: responses =
   ...: fig, ax = plt.subplots()
   ...: ax.semilogy(offsets/1e3, abs(responses.real), 'C0o-', label='|Real|')
   ...: ax.semilogy(offsets/1e3, abs(responses.imag), 'C1o-', label='|Imag|')
   ...: ax.legend()
   ...: ax.set_title('Electric field at receivers')
   ...: ax.set_xlabel('Easting (km)')
   ...: ax.set_ylabel('E-field (V/m)')
Out[9]: Text(0, 0.5, 'E-field (V/m)')

We can also get the entire electric field in a similar way as we got the computational model, and plot it for QC:

In [10]: efield = sim.get_efield(source='TxED-1', frequency=1.0)
   ....: fig, ax = plt.subplots()
   ....: cf = ax.pcolormesh(
   ....:     efield.grid.cell_centers_x/1e3,
   ....:     efield.grid.nodes_z/1e3,
   ....:     abs(efield.fx[:, 24, :].T),
   ....:     shading='gouraud',
   ....:     norm=mpl.colors.LogNorm(vmin=1e-16, vmax=1e-8),
   ....: )
   ....: fig.colorbar(cf)
   ....: ax.set_xlim([-3.5, 4])
   ....: ax.set_ylim([-3, 1])
   ....: ax.set_title(r'|$E_x$| Field (V/m)');
   ....: ax.set_xlabel('Easting (km)');
   ....: ax.set_ylabel('Depth (km)');

This is frankly a very fast rundown, and many things are only scratched at the surface or not explained at all. However, you should find much more information and explanation in the rest of the manual, and many examples in the gallery.

Usage levels#

Simulations / High-level usage#

High-level usage

Fig. 1 Workflow for the high-level usage: A Simulation needs a Model and a Survey. A survey contains all acquisition parameters such as sources, receivers, frequencies, and data, if available. A model contains the subsurface properties such as conductivities or resistivities, and the grid information.#

Simulate responses for electric and magnetic receivers due to electric and magnetic sources, in parallel. If data is provided it can also compute the misfit and the gradient of the misfit function. It includes automatic, source and frequency dependent gridding.

Note: In addition to emg3d this requires the soft dependency xarray (tqdm and discretize are recommended).

Solver-level usage#

Solver-level usage

Fig. 2 Workflow for the solver-level usage: The solve function requires a Model A and a Source-Field b. It then solves Ax=b and returns x, the electric field, corresponding to the provided subsurface model and source field.#

The solver level is the core of emg3d: It solves Maxwell’s equations for the provided subsurface model and the provided source field using the multigrid method, returning the resulting electric field.

The function emg3d.solver.solve_source simplifies the solver scheme. It takes a model, a source, and a frequency, avoiding the need to generate the source field manually, as shown in Figure 3.

Solver-source level usage

Fig. 3 Simplified solver-level workflow: The solve_source function requires a Model, a Source, and a frequency. It generates the source field internally, and returns x, the electric field, corresponding to the provided input.#

Note: This requires only emg3d (discretize is recommended).

Command-line interface (CLI-level)#

CLI-level usage

Fig. 4 CLI-level usage: file-driven command-line usage of the high-level (Simulation) functionality of emg3d.#

The command-line interface is a terminal utility for the high-level (Simulation) usage of emg3d. The model and the survey have to be provided as files (HDF5, npz, or json), various settings can be defined in the config file, and the output will be written to the output file (see also I/O & Persistence).

Note: In addition to emg3d this requires the soft dependency xarray (tqdm and discretize are recommended), and h5py if the provided files are in the HDF5 format.

Time-domain modelling#

Time-domain modelling with emg3d is possible, but it is not implemented in the high-level class Simulation. It has to be carried out by using emg3d.time.Fourier, together with the Solver-level usage mentioned above. Have a look at the repo emsig/article-TDEM.

Note: In addition to emg3d this requires the soft dependency empymod (discretize is recommended).