Simulation#
- class emg3d.simulations.Simulation(survey, model, max_workers=4, gridding='single', **kwargs)[source]#
Bases:
object
Create a simulation for a given survey on a given model.
A simulation can be used to compute responses for an entire survey, hence for an arbitrary amount of sources, receivers, and frequencies. The responses are computed in parallel over sources and frequencies. It can also be used to compute the misfit with the data and to compute the gradient of the misfit function.
The computational grid(s) can either be provided, or automatic gridding can be used; see the description of the parameters
gridding
andgridding_opts
for more details.Note
The automatic gridding does its best to generate meshes that are suitable for the provided model and survey. However, CSEM spans a wide range of acquisition layouts, and both frequencies and conductivities or resistivities span many orders of magnitude. This makes it hard to have a function that fits all purposes. Check the meshes with your expert knowledge. Also, the automatic gridding is conservative in its estimate, trying to be on the save side (correct results over speed). This means, however, that often smaller grids could be used by providing the appropriate options in
gridding_opts
or directly providing your own computational grids.Note
The package
xarray
has to be installed in order to useSimulation
:pip install xarray
orconda install -c conda-forge xarray
.- Parameters
- surveySurvey
The survey; a
emg3d.surveys.Survey
instance. The survey contains sources, receivers, frequencies, and optionally data.The survey-data will be modified in place. Provide
survey.copy()
if you want to avoid this.- modelModel
The model; a
emg3d.models.Model
instance.- max_workersint, default: 4
The maximum number of processes that can be used to execute the given calls.
- griddingstr, default: ‘single’
Method to create the computational grids.
The different methods are:
'same'
: Same grid as for the input model.'single'
: A single grid for all sources and frequencies.'frequency'
: Frequency-dependent grids.'source'
: Source-dependent grids.'both'
: Frequency- and source-dependent grids.'input'
: Same as'single'
, but the grid has to be provided ingridding_opts
instead of being automatically created.'dict'
: Same as'both'
, but the grids have to be provided ingridding_opts
in the form ofdict[source][frequency]
instead of being automatically created.
See the parameter
gridding_opts
for more details.If
gridding
is'same'
or'input'
, the input grid is checked if it is a sensible grid for emg3d; if not, it throws a warning. In the other cases the grids are created by emg3d itself, they will be fine. (If'dict'
we assume the user knows how to provide good grids.)- gridding_opts{dict, TensorMesh}, default: {}
Input format depends on
gridding
:'same'
: Nothing,gridding_opts
is not permitted.'single'
,'frequency'
,'source'
,'both'
: Described below.'input'
: Aemg3d.meshes.TensorMesh
.'dict'
: Dictionary of the formatdict[source][frequency]
containing aemg3d.meshes.TensorMesh
for each source-frequency pair.
The dict in the case of
'single'
,'frequency'
,'source'
,'both
’ is passed toemg3d.meshes.construct_mesh
; consult the corresponding documentation for more information. Parameters that are not provided are estimated from the provided model, grid, and survey usingemg3d.meshes.estimate_gridding_opts
, which documentation contains more information too.There are two notably differences to the parameters described in
emg3d.meshes.construct_mesh
:vector
: besides the normal possibility it can also be a string containing one or several of'x'
,'y'
, and'z'
. In these cases the corresponding dimension of the input mesh is provided as vector. Seeemg3d.meshes.estimate_gridding_opts
.expand
: in the format of[property_sea, property_air]
; if provided, the input model is expanded up to the seasurface with sea water, and an air layer is added. The actual height of the seasurface can be defined with the keyseasurface
. Seeemg3d.models.expand_grid_model
. NOTE:expand
is deprecated in v1.7.0, and will be removed in v1.9.0. A property-complete model has to be provided.
- solver_optsdict, default: {‘verb’: 1’}
Passed through to
emg3d.solver.solve
. The dict can contain any parameter that is accepted by theemg3d.solver.solve
except formodel
,sfield
,efield
,return_info
, andlog
. Default verbosity isverb=2
.- verbint, default: 0
Level of verbosity. Possible options:
-1: Errors.
0: Warnings.
1: Info.
- namestr, default: None
Name of the simulation.
- infostr, default: None
Simulation info or any other info (e.g., what was the purpose of this simulation).
- receiver_interpolationstr, default: ‘cubic’:
Interpolation method to obtain the response at receiver location; ‘cubic’ or ‘linear’. Cubic is more precise. However, if you are interested in the gradient, you need to choose ‘linear’ at the moment, as there are only linearly interpolated source functions. To be the proper adjoint for the gradient the receiver has to be interpolated linearly too.
- file_dirstr, default: None
Absolute or relative path, where temporary files should be stored. By default, everything is done in memory. However, for large models with many sources and frequencies this can become memory consuming. Providing a
file_dir
can help with this. Fields and models are then stored to disk, and each process accesses the files it needs. There is only a gain if there are more source-frequency pairs than concurrent running processes.Note that the directory is created if it does not exist. However, the parent directory must exist.
Also note that the files are stored as .h5-files, and you need to have
h5py
installed to use this feature.- tqdm_opts{bool, dict}, default: True
Boolean if a progress bar should be shown (only if
tqdm
is installed). If a dict is provided it will be forwarded totqdm
.- layeredbool, default: False
If True, the responses are computed with approximated layered (1D) models for each source-receiver pair using empymod.
The computation happens in parallel for each source location. Each source-receiver pair is done separately, but all frequencies at once.
If layered is set, the gradient is computed using the finite-difference method, by perturbing each layer slightly.
Current limitations:
Only point and dipole sources.
Only isotropic and VTI models.
Setting this to True also means:
There are no {e;h}fields, only the fields at receiver locations.
gridding
, most ofgridding_opts
,solver_opts
,receiver_interpolation
, andfile_dir
have no effect.The attribute
emg3d.simulations.Simulation.jvec`
is not implemented.
- layered_optsdict, default: {}
Options passed to
emg3d.models.Model.extract_1d
, defining how the layered model is obtained from the 3D model.p0
andp1
are taken to be the source and receiver locations. Consult that function for more information. Below are only things described which differ.The possible methods are: cylinder, prism, midpoint, source, receiver. The last two are the same as
midpoint
, where just both points are set either to the source or receiver location, respectively.The default method is
'cylinder'
. If the ellipse parameters are not given for the methods cylinder and prism, they are set as follows:factor: 1.2.
minor: 0.8.
radius: one skin depth using the lowest frequency of the survey and the downwards property from the gridding options (or the minimum conductivity value in the lowest vertical layer).
Attributes Summary
Shortcut to survey.data.
Compute the discrete gradient using the adjoint-state method.
If True, use layered computations (empymod).
Misfit or cost function.
Methods Summary
clean
([what])Clean part of the data base.
compute
([observed])Compute efields asynchronously for all sources and frequencies.
copy
([what])Return a copy of the Simulation.
from_dict
(inp)Convert dict into
emg3d.simulations.Simulation
instance.from_file
(fname[, name])Load Simulation from a file.
get_efield
(source, frequency)Return electric field for given source and frequency.
get_efield_info
(source, frequency)Return the solver information of the corresponding computation.
get_grid
(source, frequency)Return computational grid of the given source and frequency.
get_hfield
(source, frequency)Return magnetic field for given source and frequency.
get_model
(source, frequency)Return model on the grid of the given source and frequency.
jtvec
(vector)Compute the sensitivity transpose times a vector.
jvec
(vector)Compute the sensitivity times a vector.
print_grid_info
([verb, return_info])Print info for all generated grids.
print_solver_info
([field, verb, return_info])Print solver info.
to_dict
([what, copy])Store the necessary information of the Simulation in a dict.
to_file
(fname[, what, name])Store Simulation to a file.
Attributes Documentation
- data#
Shortcut to survey.data.
- gradient#
Compute the discrete gradient using the adjoint-state method.
The discrete adjoint-state gradient for a single source at a single frequency is given by Equation (10) in [PlMu08],
\[\begin{split}\nabla_p \phi(\textbf{p}) = -&\sum_{k,l,m}\mathbf{\bar{\lambda}}_{x; k+\frac{1}{2}, l, m} \frac{\partial S_{k+\frac{1}{2}, l, m}}{\partial \textbf{p}} \textbf{E}_{x; k+\frac{1}{2}, l, m}\\ -&\sum_{k,l,m}\mathbf{\bar{\lambda}}_{y; k, l+\frac{1}{2}, m} \frac{\partial S_{k, l+\frac{1}{2}, m}}{\partial \textbf{p}} \textbf{E}_{y; k, l+\frac{1}{2}, m}\\ -&\sum_{k,l,m}\mathbf{\bar{\lambda}}_{z; k, l, m+\frac{1}{2}} \frac{\partial S_{k, l, m+\frac{1}{2}}}{\partial \textbf{p}} \textbf{E}_{z; k, l, m+\frac{1}{2}}\, ,\end{split}\]where \(\textbf{E}\) is the electric (forward) field and \(\mathbf{\lambda}\) is the back-propagated residual field (from electric and magnetic receivers); \(\bar{~}\) denotes conjugate. The \(\partial S\)-part takes care of the volume-averaged model parameters.
Warning
To obtain the proper adjoint-state gradient you have to choose linear interpolation for the receiver responses:
emg3d.Simulation(..., receiver_interpolation='linear')
. The reason is that the point-source is the adjoint of a tri-linear interpolation, so the residual should come from a linear interpolation.Note
The currently implemented gradient does only work for models without relative electric permittivity nor relative magnetic permeability.
Note
If
layered=True
, the gradient is computed using finite differences.- Returns
- gradndarray
Adjoint-state gradient. Shape depends on the anisotropy type:
isotropic: (nx, ny, nz)
HTI/VTI: (2, nx, ny, nz)
triaxial: (3, nx, ny, nz)
- layered#
If True, use layered computations (empymod).
- misfit#
Misfit or cost function.
The data misfit or weighted least-squares functional using an \(l_2\) norm is given by
(37)#\[ \phi = \frac{1}{2} \sum_s\sum_r\sum_f \left\lVert W_{s,r,f} \left( \textbf{d}_{s,r,f}^\text{pred} -\textbf{d}_{s,r,f}^\text{obs} \right) \right\rVert^2 \, ,\]where \(s, r, f\) stand for source, receiver, and frequency, respectively; \(\textbf{d}^\text{obs}\) are the observed electric and magnetic data, and \(\textbf{d}^\text{pred}\) are the synthetic electric and magnetic data. As of now the misfit does not include any regularization term.
The data weight of observation \(d_i\) is given by \(W_i = \varsigma^{-1}_i\), where \(\varsigma_i\) is the standard deviation of the observation, see
emg3d.surveys.Survey.standard_deviation
.Note
You can easily implement your own misfit function (to include, e.g., a regularization term) by monkey patching this misfit function with your own:
@property # misfit is a property def my_misfit_function(self): '''Returns the misfit as a float.''' if self._misfit is None: self.compute() # Ensures fields are computed. # Computing your misfit... self._misfit = your misfit return self._misfit # Monkey patch simulation.misfit: emg3d.simulation.Simulation.misfit = my_misfit_function # And now all the regular stuff, initiate a Simulation etc simulation = emg3d.Simulation(survey, grid, model) simulation.misfit # => will return your misfit # (will also be used for the adjoint-state gradient).
- Returns
- misfitfloat
Value of the misfit function.
Methods Documentation
- clean(what='computed')[source]#
Clean part of the data base.
- Parameters
- whatstr, default: ‘computed’
What to clean. Possibilities:
'computed'
: Removes all computed properties: electric and magnetic fields and responses at receiver locations.'keepresults'
: Removes everything except for the responses at receiver locations.'all'
: Removes everything (leaves it plain as initiated).
- compute(observed=False, **kwargs)[source]#
Compute efields asynchronously for all sources and frequencies.
- Parameters
- observedbool, default: False
If True, it stores the current synthetic responses also as observed responses.
- add_noisebool, default: True
Boolean if to add noise to observed data (if
observed=True
). All remainingkwargs
are forwarded toemg3d.surveys.Survey.add_noise
.
- copy(what='computed')[source]#
Return a copy of the Simulation.
See
to_file
for more information regardingwhat
.
- classmethod from_dict(inp)[source]#
Convert dict into
emg3d.simulations.Simulation
instance.- Parameters
- inpdict
Dictionary as obtained from
emg3d.simulations.Simulation.to_dict
.
- Returns
- simulationSimulation
A
emg3d.simulations.Simulation
instance.
- classmethod from_file(fname, name='simulation', **kwargs)[source]#
Load Simulation from a file.
- Parameters
- fnamestr
Absolute or relative file name including extension.
- namestr, default: ‘simulation’
Name under which the simulation is stored within the file.
- kwargsKeyword arguments, optional
Passed through to
io.load
.
- Returns
- simulationSimulation
A
emg3d.simulations.Simulation
instance.- infostr, returned if verb<0
Info-string.
- get_efield_info(source, frequency)[source]#
Return the solver information of the corresponding computation.
- jtvec(vector)[source]#
Compute the sensitivity transpose times a vector.
If
vector=residual*weights
,jtvec
corresponds to thegradient
.(38)#\[J^H v = G^H A^{-H} P^H v \ ,\]where \(v\) has the shape of the data.
- Parameters
- vectorndarray, DataArray
An array with the shape of the data, or directly a DataArray as stored in the
emg3d.surveys.Survey
.
- Returns
- jtvecndarray
Adjoint-state gradient for the provided vector. Shape depends on the anisotropy type:
isotropic: (nx, ny, nz)
HTI/VTI: (2, nx, ny, nz)
triaxial: (3, nx, ny, nz)
- jvec(vector)[source]#
Compute the sensitivity times a vector.
(39)#\[J v = P A^{-1} G v \ ,\]where \(v\) has size of the model.
- Parameters
- vectorndarray
Vector applied to J. Shape depends on the anisotropy type:
isotropic: (nx, ny, nz) or (1, nx, ny, nz)
HTI/VTI: (2, nx, ny, nz)
triaxial: (3, nx, ny, nz)
- Returns
- jvecndarray
Shape of the data.
- to_dict(what='computed', copy=False)[source]#
Store the necessary information of the Simulation in a dict.
See to_file for more information regarding what.
- to_file(fname, what='computed', name='simulation', **kwargs)[source]#
Store Simulation to a file.
- Parameters
- fnamestr
Absolute or relative file name including ending, which defines the used data format. See
emg3d.io.save
for the options.- whatstr, default: ‘computed’
What to store. Possibilities:
'computed'
,'all'
: Stores all computed properties: electric fields and responses at receiver locations.‘
results'
: Stores only the response at receiver locations.'plain'
: Only stores the plain Simulation (as initiated).
Note that if
file_dir
is set, those files will remain there.- namestr, default: ‘simulation’
Name with which the simulation is stored in the file.
- kwargsKeyword arguments, optional
Passed through to
emg3d.io.save
.