Getting started =============== There are different usages of ``emg3d``. As an end user you might be interested in the :ref:`high-level-usage`: define your survey and your model, and simulate electromagnetic data for them. As a developer you might be interested in the :ref:`solver-level-usage` or :ref:`cli-level-usage`, to implement emg3d for instance as a forward modelling kernel in you inversion code. And then there is also :ref:`time-level-usage`. Here we show first an introductory example of the :ref:`high-level-usage`. Further down you find an overview and more detailed explanation of the different :ref:`usages`. 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. .. ipython:: 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. Grid ~~~~ We first build a simple four-quadrant grid, centred around the origin. .. ipython:: In [1]: hx = np.ones(2)*5000 ...: grid = emg3d.TensorMesh( ...: [[5000, 5000], [10000], [5000, 5000]], origin='CCC') ...: grid # QC Out[1]: ...: ...: 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 Model ~~~~~ We use the constructed grid to create a simple resistivity model: .. ipython:: In [1]: model = emg3d.Model( ...: grid=grid, ...: property_x=[1, 10, 0.3, 0.3], ...: mapping='Resistivity' ...: ) ...: model # QC Out[1]: Model: resistivity; isotropic; 2 x 2 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: .. ipython:: @savefig basic_model.png width=4in In [1]: 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. Survey ~~~~~~ Now that we have a model we need to define our survey. Currently there are three source types implemented, - :class:`emg3d.electrodes.TxElectricDipole`; - :class:`emg3d.electrodes.TxMagneticDipole`; - :class:`emg3d.electrodes.TxElectricWire`; and two receiver types, - :class:`emg3d.electrodes.RxElectricPoint`; - :class:`emg3d.electrodes.RxElectricPoint`. We are going to define a simple survey with an electric dipole source and a line of electric point receivers. .. ipython:: In [1]: 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 Simulation ~~~~~~~~~~ Now that we have a model and a survey we can define our simulation: .. ipython:: :okwarning: In [1]: sim = emg3d.Simulation( ...: survey=survey, ...: model=model, ...: ) ...: ...: sim # QC 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 :class:`emg3d.simulations.Simulation`. Particularly important are the gridding parameters, and as such it is highly recommended to also read :func:`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: .. ipython:: In [1]: comp_model = sim.get_model(source='TxED-1', frequency=1.0) ...: comp_model.grid Out[1]: ...: ...: 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, .. ipython:: In [1]: sim.compute() Results ~~~~~~~ Let's plot the electric field at receiver locations (the responses): .. ipython:: @savefig basic_receivers.png width=4in In [1]: responses = sim.data.synthetic.data.squeeze() ...: 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)') We can also get the entire electric field in a similar way as we got the computational model, and plot it for QC: .. ipython:: @savefig basic_efield.png width=4in In [1]: 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. .. _usages: Usage levels ------------ .. _high-level-usage: Simulations / High-level usage ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. figure:: ../_static/levels1.svg :align: center :alt: High-level usage :name: high-level 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 ~~~~~~~~~~~~~~~~~~ .. figure:: ../_static/levels2.svg :align: center :alt: Solver-level usage :name: solver-level 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 :func:`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 :numref:`Figure %s `. .. figure:: ../_static/levels4.svg :align: center :alt: Solver-source level usage :name: solver-source-level 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). .. _cli-level-usage: Command-line interface (CLI-level) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. figure:: ../_static/levels3.svg :align: center :alt: CLI-level usage :name: cli-level 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 :ref:`io-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-level-usage: 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 :class:`emg3d.time.Fourier`, together with the Solver-level usage mentioned above. Have a look at the repo https://github.com/emsig/article-TDEM. *Note:* In addition to ``emg3d`` this requires the soft dependency ``empymod`` (``discretize`` is recommended).