From 3b9a3c989be3bcb0a23524d6b35ae90b73cb1772 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Dieter=20Werthm=C3=BCller?= Date: Wed, 31 Aug 2022 16:46:49 +0200 Subject: [PATCH] Layered computation using empymod. (#302) The simulation class takes new the parameters ``layered`` and ``layered_opts``, where the default values are False and None, respectively. If ``layered=True``, there will be no 3D computations. Instead, it will create a local layered (1D) model for each source-receiver pair, and compute the response using the semi-analytical code ``empymod``. In this case an eventual gradient is computed using the finite-difference method, not the adjoint-state method, perturbing each layer slightly. Layered computation is also possible through the CLI, through the new flag ``-l`` or ``--layered``, and a new section ``[layered]`` in the config file. --- CHANGELOG.rst | 20 +- docs/manual/cli.rst | 16 +- emg3d/_multiprocessing.py | 316 ++++++++++++++++++++++++- emg3d/cli/main.py | 8 + emg3d/cli/parser.py | 56 ++++- emg3d/cli/run.py | 6 + emg3d/simulations.py | 427 ++++++++++++++++++++++++++-------- tests/test_cli.py | 49 +++- tests/test_multiprocessing.py | 180 ++++++++++++++ tests/test_simulations.py | 161 ++++++++++++- 10 files changed, 1138 insertions(+), 101 deletions(-) diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 5656c40b..f3383c65 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -6,8 +6,24 @@ Changelog """""""""" -latest ------- +v1.8.0 : Layered modelling +-------------------------- + +**2022-08-31** + +The simulation class takes new the parameters ``layered`` and ``layered_opts``, +where the default values are False and None, respectively. If ``layered=True``, +there will be no 3D computations. Instead, it will create a local layered (1D) +model for each source-receiver pair, and compute the response using the +semi-analytical code ``empymod`` (which needs to be installed manually, as it +is a soft dependency). In this case an eventual gradient is computed using the +finite-difference method, not the adjoint-state method, perturbing each layer +slightly. The main purpose of these layered computations is for quick checks, +QC, verifications, etc. Layered computation is also possible through the CLI, +through the new flag ``-l`` or ``--layered``, and a new section ``[layered]`` +in the config file. + +Other changes (many of them related to the above): - Model instances have a new attribute ``exctract_1d``, which returns a layered (1D) model, extracted from the 3D model according the provided parameters; diff --git a/docs/manual/cli.rst b/docs/manual/cli.rst index 5267fd99..915f347c 100644 --- a/docs/manual/cli.rst +++ b/docs/manual/cli.rst @@ -59,8 +59,9 @@ remove the comment signs to use them. # max_workers = 4 # Also via `-n` or `--nproc`. # gridding = single # name = MyTestSimulation - # file_dir = None # For file-based comp; absolute or relative path. + # file_dir = None # For file-based comp; absolute or relative path. # receiver_interpolation = cubic # Set it to for the gradient. + # layered = False # Also via `-l` or `--layered`. # Solver options # -------------- @@ -102,7 +103,6 @@ remove the comment signs to use them. # center = # list, e.g.: 0, 0, 0 # cell_number = # list, e.g.: 8, 16, 32, 64, 128 # min_width_pps = # list, e.g.: 5, 3, 3 - # expand = # list, e.g.: 0.3, 1e8 # domain = # list of lists, e.g.: -10000, 10000; None; None # distance = # list of lists, e.g., None; None; -10000, 10000 # stretching = # list of lists, e.g.: None; None; 1.05, 1.5 @@ -136,3 +136,15 @@ remove the comment signs to use them. # receivers = RxEP-01, RxMP-10 # frequencies = f-1, f-3 # remove_empty = False # CLI uses False by default. + + # Layered computation + # ------------------- + # The following parameters are only used if `-l`/`--layered` is set or the + # simulation section has set `layered` to True. + [layered] + # method = # str + # radius = # float + # factor = # float + # minor = # float + # merge = # bool + # check_foci = # bool diff --git a/emg3d/_multiprocessing.py b/emg3d/_multiprocessing.py index 15f7ac45..2dd43f94 100644 --- a/emg3d/_multiprocessing.py +++ b/emg3d/_multiprocessing.py @@ -16,15 +16,18 @@ # WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. See the # License for the specific language governing permissions and limitations under # the License. +from copy import deepcopy from concurrent.futures import ProcessPoolExecutor +import numpy as np + try: import tqdm import tqdm.contrib.concurrent except ImportError: tqdm = None -from emg3d import io, solver +from emg3d import io, solver, utils def process_map(fn, *iterables, max_workers, **kwargs): @@ -143,3 +146,314 @@ def solve(inp): return fname, fname else: return efield, info + + +@utils._requires('empymod') +def layered(inp): + """Returns response or gradient using layered models; for a `process_map`. + + Used within a Simulation to call empymod in parallel for layered models. + Depending on the input it returns either the forward responses or the + finite-difference gradient. + + The parameters section describes the content of the input dict. + + Parameters + ---------- + model : Model + The model; a :class:`emg3d.models.Model` instance. Must be isotropic or + VTI. + + src : Tx* + Any dipole or point source of the available sources, e.g., + :class:`emg3d.electrodes.TxElectricDipole`. + + receivers : dict of Rx* + Receiver dict (:attr:`emg3d.surveys.Survey.receivers`). + + frequencies : dict + Frequency dict (:attr:`emg3d.surveys.Survey.frequencies`). + + empymod_opts : dict + Options passed to empymod ({src;rec}pts, {h;f}t, {h;f}targ, xdirect, + loop). + + observed : DataArray + Observed data of this source. + + layered_opts : dict + Options passed to :attr:`emg3d.models.Model.extract_1d`. + + gradient : bool + If False, the electromagnetic responses are returned; if True, the + gradient is returned. + + If True, the following things _have_ to be provided: ``observed``, + ``weights``, and ``residual``; otherwise a zero gradient is returned. + + weights : DataArray, optional + Data weights corresponding to the data; only required if + ``gradient=True``. + + residual : DataArray + Residuals using the current model; only required if ``gradient=True``. + + + Returns + ------- + out : ndarray + If ``gradient=False``, the output are the electromagnetic responses + (synthetic data) of shape (nrec, nfreq). + + If ``gradient=True``, the output is the finite-difference gradient of + shape (3, nx, ny, nz). + + """ + + # Extract input. + model = inp['model'] + src = inp['src'] + receivers = inp['receivers'] + frequencies = np.array([f for f in inp['frequencies'].values()]) + empymod_opts = inp['empymod_opts'] + observed = inp['observed'] + lopts = deepcopy(inp['layered_opts']) + gradient = inp['gradient'] + + # Get method and set to return_imat. + method = lopts.pop('method') + lopts['return_imat'] = True + + # Collect rec-independent empymod options. + empymod_opts = { + # User input ({src;rec}pts, {h;f}t, {h;f}targ, xdirect, loop). + # Contains also verb from simulation class. + **empymod_opts, + # + # Source properties, same for all receivers. + 'src': src.coordinates, + 'msrc': src.xtype != 'electric', + 'strength': src.strength, + # + # Enforced properties (not implemented). + 'signal': None, + 'epermV': None, + 'mpermV': None, + 'squeeze': True, + } + + # Create some flags. + epsilon_r = model.epsilon_r is not None + mu_r = model.mu_r is not None + vti = model.case == 'VTI' + + # Pre-allocate output array. + if gradient: + # Gradient. + out = np.zeros((3, *model.shape)) + + # Get weights and residual if the gradient is wanted. + weights = inp.get('weights', None) + residual = inp.get('residual', None) + if weights is None or residual is None or observed is None: + return out + + else: + # Responses. + out = np.full((len(receivers), frequencies.size), np.nan+1j*np.nan) + + # Loop over receivers. + for i, (rkey, rec) in enumerate(receivers.items()): + + # Check observed data, limit to finite values if provided. + if observed is not None: + fi = np.isfinite(observed.loc[rkey, :].data) + if fi.sum() == 0: + continue + freqs = frequencies[fi] + + # Generating obs data for all. + else: + fi = np.ones(frequencies.size, dtype=bool) + freqs = frequencies + + # Get 1D model. + # Note: if method='source', this would be faster outside the loop. + oned, imat = model.extract_1d(**_get_points(method, src, rec), **lopts) + + # Collect input. + empymod_inp = { + **empymod_opts, + 'rec': rec.coordinates, + 'mrec': rec.xtype != 'electric', + 'depth': oned.grid.nodes_z[1:-1], + 'freqtime': freqs, + 'epermH': None if not epsilon_r else oned.epsilon_r[0, 0, :], + 'mpermH': None if not mu_r else oned.mu_r[0, 0, :], + } + + # Get horizontal and vertical conductivities. + map2cond = oned.map.backward + cond_h = map2cond(oned.property_x[0, 0, :]) + cond_v = None if not vti else map2cond(oned.property_z[0, 0, :]) + + # Compute gradient. + if gradient: + + # Get misfit of this src-rec pair. + obs = observed.loc[rkey, :].data[fi] + wgt = weights.loc[rkey, :].data[fi] + res = residual.loc[rkey, :].data[fi] + misfit = np.sum(wgt*(res.conj()*res)).real/2 + + # Get horizontal gradient. + out[0, ...] += _fd_gradient(cond_h, cond_v, obs, wgt, misfit, + empymod_inp, imat, vertical=False) + + # Get vertical gradient if VTI. + if vti: + out[2, ...] += _fd_gradient(cond_h, cond_v, obs, wgt, misfit, + empymod_inp, imat, vertical=True) + + # Compute response. + else: + out[i, fi] = _empymod_fwd(cond_h, cond_v, empymod_inp) + + return out + + +@utils._requires('empymod') +def _empymod_fwd(cond_h, cond_v, empymod_inp): + """Thin wrapper for empymod.bipole(). + + Parameters + ---------- + cond_h, cond_v : ndarray + Horizontal and vertical conductivities (S/m). ``cond_v`` can be None, + in which case an isotropic medium is assumed. + + empymod_inp : dict + Passed through to :func:`empymod.model.bipole`. Any parameter except + for ``res`` and ``aniso``. + + + Returns + ------- + resp : EMArray + Electromagnetic field as returned from :func:`empymod.model.bipole`. + + + """ + from empymod import bipole + aniso = None if cond_v is None else np.sqrt(cond_h/cond_v) + return bipole(res=1/cond_h, aniso=aniso, **empymod_inp) + + +def _get_points(method, src, rec): + """Returns correct method and points for model.extract_1d. + + Parameters + ---------- + method : str + All methods accepted by :attr:`emg3d.models.Model.extract_1d` plus + ``'source'``, ``'receiver'``. + + src, rec : {Tx*, Rx*) + Any of the available point and dipole sources or receivers, e.g., + :class:`emg3d.electrodes.TxElectricDipole`. + + Returns + ------- + out : dict + Can be passed directly to :attr:`emg3d.models.Model.extract_1d` for the + parameters ``method``, ``p0``, and ``p1``. + + """ + + # Get default points. + p0 = src.center[:2] + p1 = rec.center[:2] + + # If source or receiver, we re-set one point and rename the method + if method == 'source': + p1 = p0 + method = 'midpoint' + + elif method == 'receiver': + p0 = p1 + method = 'midpoint' + + return {'method': method, 'p0': p0, 'p1': p1} + + +def _fd_gradient(cond_h, cond_v, data, weight, misfit, empymod_inp, imat, + vertical): + """Computes the finite-difference gradient using empymod. + + The finite difference is obtained by adding a relative difference of 0.01 % + to the layer (currently hard-coded). + + Parameters + ---------- + cond_h, cond_v : ndarray + Horizontal and vertical conductivities (S/m). ``cond_v`` can be None, + in which case an isotropic medium is assumed. + + data : ndarray + Observed data. + + weight : ndarray + Weights corresponding to these data. + + misfit : float + Misfit using the current model. + + empymod_inp : dict + Passed through to :func:`empymod.model.bipole`. Any parameter except + for ``res`` and ``aniso``. + + imat : ndarray + Interpolation matrix as returned by + :attr:`emg3d.models.Model.extract_1d`. + + vertical : bool + If True, the gradient for the vertical conductivities is assumed, + otherwise the gradient for the horizontal conductivities. + If ``vertical=True``, ``cond_v`` cannot be None (not checked, will fail + with an AttributeError). + + + Returns + ------- + gradient : ndarray of shape (nx, ny, nz) + Gradient. + + + """ + # Relative difference fixed to 0.01 %; could be made an input parameter. + rel_diff = 0.0001 + + # Loop over layers and compute FD gradient for each. + grad = np.zeros(cond_h.size) + for iz in range(cond_h.size): + + # Get 1D model. + cond_p = cond_h.copy() if not vertical else cond_v.copy() + + # Add relative difference to the layer. + delta = cond_p[iz] * rel_diff + cond_p[iz] += delta + + # Call empymod. + if vertical: + response = _empymod_fwd(cond_h, cond_p, empymod_inp) + else: + response = _empymod_fwd(cond_p, cond_v, empymod_inp) + + # Calculate gradient and add it. + residual = response - data + fd_misfit = np.sum(weight*(residual.conj()*residual)).real/2 + grad[iz] = (fd_misfit - misfit)/delta + + # Bring back to full grid and return. + return imat[..., None] * grad[None, :] diff --git a/emg3d/cli/main.py b/emg3d/cli/main.py index 294d2291..04e91937 100644 --- a/emg3d/cli/main.py +++ b/emg3d/cli/main.py @@ -145,6 +145,14 @@ def main(args=None): help="replace model and all computed data of loaded simulation" ) + # arg: Use layered models with empymod. + parser.add_argument( + "-l", "--layered", + action="store_true", + default=None, + help="use a layered model with empymod for each src-rec pair" + ) + # arg: Run without emg3d-computation. parser.add_argument( "-d", "--dry-run", diff --git a/emg3d/cli/parser.py b/emg3d/cli/parser.py index 73c9ddbf..7b1a41d4 100644 --- a/emg3d/cli/parser.py +++ b/emg3d/cli/parser.py @@ -66,7 +66,7 @@ def parse_config_file(args_dict): term['config_file'] = configfile # Get terminal input. - for key in ['verbosity', 'nproc', 'dry_run', 'clean']: + for key in ['verbosity', 'nproc', 'dry_run', 'clean', 'layered']: term[key] = args_dict.pop(key) for key in ['forward', 'misfit', 'gradient']: @@ -171,6 +171,15 @@ def parse_config_file(args_dict): simulation[key] = cfg.getint('simulation', key) del term['nproc'] + # Check layered. + key = 'layered' + _ = all_sim.pop(key, None) + if term[key] is not None: + simulation[key] = term[key] + elif cfg.has_option('simulation', key): + simulation[key] = cfg.getboolean('simulation', key) + del term[key] + # Check gridding. key = 'gridding' if cfg.has_option('simulation', key): @@ -252,6 +261,51 @@ def parse_config_file(args_dict): f"{list(all_noise.keys())}." ) + # # Layered parameters # # + + # Check if parameter-file has a layered-section, add it otherwise. + layered_opts = {} + if 'layered' in cfg.sections(): + + all_layered = dict(cfg.items('layered')) + + key = 'method' + if cfg.has_option('layered', key): + _ = all_layered.pop(key) + layered_opts[key] = cfg.get('layered', key) + + key = 'merge' + if cfg.has_option('layered', key): + _ = all_layered.pop(key) + layered_opts[key] = cfg.getboolean('layered', key) + + # Collect Ellipse-Options + ellipse = {} + + for key in ['radius', 'minor', 'factor']: + if cfg.has_option('layered', key): + _ = all_layered.pop(key) + ellipse[key] = float(cfg.get('layered', key)) + + key = 'check_foci' + if cfg.has_option('layered', key): + _ = all_layered.pop(key) + ellipse[key] = cfg.getboolean('layered', key) + + if ellipse: + layered_opts['ellipse'] = ellipse + + # Ensure no keys are left. + if all_layered: + raise TypeError( + f"Unexpected parameter in [layered]: " + f"{list(all_layered.keys())}." + ) + + # Add to simulation dict if not empty. + if layered_opts: + simulation['layered_opts'] = layered_opts + # # Solver parameters # # # Check if parameter-file has a solver-section, add it otherwise. diff --git a/emg3d/cli/run.py b/emg3d/cli/run.py index 6656be64..3814a895 100644 --- a/emg3d/cli/run.py +++ b/emg3d/cli/run.py @@ -103,6 +103,12 @@ def simulation(args_dict): sim.model = models.expand_grid_model( sim.model, expand, interface) + # Set layered according to user input if it differs. + layered = cfg['simulation_options'].get('layered', False) + if sim.layered != layered: + logger.info(f"Change «layered» of simulation to {layered}.") + sim.layered = layered + else: # Load input. diff --git a/emg3d/simulations.py b/emg3d/simulations.py index 1e00398f..379ad748 100644 --- a/emg3d/simulations.py +++ b/emg3d/simulations.py @@ -191,6 +191,49 @@ class Simulation: Boolean if a progress bar should be shown (only if ``tqdm`` is installed). If a dict is provided it will be forwarded to ``tqdm``. + layered : bool, 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 of ``gridding_opts``, ``solver_opts``, + ``receiver_interpolation``, and ``file_dir`` have no effect. + - The attribute :attr:`emg3d.simulations.Simulation.jvec`` is not + implemented. + + layered_opts : dict, default: {} + Options passed to :attr:`emg3d.models.Model.extract_1d`, defining how + the layered model is obtained from the 3D model. ``p0`` and ``p1`` + 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). + """ # Gridding descriptions (for repr's). @@ -232,6 +275,7 @@ def __init__(self, survey, model, max_workers=4, gridding='single', self._dict_efield_info = self._dict_initiate self._gradient = None self._misfit = None + self._computed = False # Initiate file_dir self.file_dir = kwargs.pop('file_dir', None) @@ -243,6 +287,8 @@ def __init__(self, survey, model, max_workers=4, gridding='single', # Get model taking gridding_opts into account. # Sets self.model and self.gridding_opts. self._set_model(model, kwargs) + self._set_layered_opts(kwargs.pop('layered', False), + kwargs.pop('layered_opts', {})) # Initiate synthetic data with NaN's if they don't exist. if 'synthetic' not in self.survey.data.keys(): @@ -279,8 +325,7 @@ def __repr__(self): f"{self.survey.shape[1]} receivers; " f"{self.survey.shape[2]} frequencies\n" f"- {self.model.__repr__()}\n" - f"- Gridding: {self._gridding_descr[self.gridding]}; " - f"{self._info_grids}") + f"- {self._info_grids}") def _repr_html_(self): """HTML representation.""" @@ -294,8 +339,7 @@ def _repr_html_(self): f" {self.survey.shape[1]} receivers;" f" {self.survey.shape[2]} frequencies" f"
  • {self.model.__repr__()}
    • " - f"
  • Gridding: {self._gridding_descr[self.gridding]}; " - f" {self._info_grids}
    • " + f"
  • {self._info_grids}
    • " f"") def clean(self, what='computed'): @@ -348,6 +392,7 @@ def clean(self, what='computed'): # Clean data. if what in ['computed', 'all']: + self._computed = False for key in ['residual', 'weights']: if key in self.data.keys(): del self.data[key] @@ -393,6 +438,8 @@ def to_dict(self, what='computed', copy=False): 'name': self.name, 'info': self.info, 'tqdm_opts': self._tqdm_opts, + 'layered': self.layered, + 'layered_opts': self.layered_opts, 'receiver_interpolation': self.receiver_interpolation, 'file_dir': self.file_dir, '_input_sc2': self._input_sc2, @@ -416,6 +463,7 @@ def to_dict(self, what='computed', copy=False): if what in ['computed', 'results', 'all']: out['gradient'] = self._gradient out['misfit'] = self._misfit + out['computed'] = self._computed if copy: return deepcopy(out) @@ -455,6 +503,8 @@ def from_dict(cls, inp): 'receiver_interpolation', 'cubic') cls_inp['file_dir'] = inp.pop('file_dir', None) cls_inp['tqdm_opts'] = inp.pop('tqdm_opts', {}) + cls_inp['layered'] = inp.pop('layered', False) + cls_inp['layered_opts'] = inp.pop('layered_opts', {}) # Instantiate the class. out = cls(**cls_inp) @@ -473,7 +523,7 @@ def from_dict(cls, inp): setattr(out, name, values) # Add gradient and misfit. - data = ['gradient', 'misfit'] + data = ['gradient', 'misfit', 'computed'] for name in data: if name in inp.keys(): setattr(out, '_'+name, inp.pop(name)) @@ -743,13 +793,35 @@ def compute(self, observed=False, **kwargs): :meth:`emg3d.surveys.Survey.add_noise`. """ + source = kwargs.pop('source', None) + frequency = kwargs.pop('frequency', None) + if self.layered: + if source or frequency: + raise NotImplementedError("No fields if `layered` is used.") + self._compute_1d() + else: + # Undocumented (internal): + # If the call is from `get_efield`, it will have source/frequency. + # This use is only internal. End users should use `get_efield()`. + self._compute([(source, frequency), ]) + + # If it shall be used as observed data save a copy. + if observed: + + # Copy synthetic to observed. + self.data['observed'] = self.data['synthetic'].copy() + + # Add noise. + if kwargs.pop('add_noise', True): + self.survey.add_noise(**kwargs) + + elif source is None and frequency is None: + self._computed = True + + def _compute(self, srcfreq): + """Compute efields and responses asynchronously using emg3d.solve().""" + from emg3d import _multiprocessing as _mp - # Undocumented (internal): - # If the call is from `get_efield`, it will have source/frequency. - # This use is only internal. End users should use `get_efield()`. - srcfreq = [ - (kwargs.pop('source', None), kwargs.pop('frequency', None)), - ] if not srcfreq[0][0]: # "Normal" case: all source-frequency pairs. srcfreq = self._srcfreq @@ -771,7 +843,12 @@ def collect_efield_inputs(inp): return self._data_or_file('efield', source, freq, data) # Compute fields in parallel. - out = self._compute(collect_efield_inputs, 'Compute efields', srcfreq) + out = _mp.process_map( + _mp.solve, + list(map(collect_efield_inputs, srcfreq)), + max_workers=self.max_workers, + **{'desc': 'Compute efields', **self._tqdm_opts}, + ) # Loop over src-freq combinations to extract and store. for i, (src, freq) in enumerate(srcfreq): @@ -787,26 +864,66 @@ def collect_efield_inputs(inp): # Print solver info. self.print_solver_info('efield', verb=self.verb) - # If it shall be used as observed data save a copy. - if observed: + def _compute_1d(self, gradient=False): + """Compute responses asynchronously using empymod.bipole().""" + from emg3d import _multiprocessing as _mp - # Copy synthetic to observed. - self.data['observed'] = self.data['synthetic'].copy() + # Check if there are observed data. + has_data = np.isfinite(self.data.observed.data).sum() > 0 - # Add noise. - if kwargs.pop('add_noise', True): - self.survey.add_noise(**kwargs) + # Create iterable from src-list for parallel computation. + def collect_empymod_inputs(source): + """Collect inputs.""" - def _compute(self, fn, description, srcfreq=None): - """Use process_map to call solver.solve asynchronously.""" - from emg3d import _multiprocessing as _mp - return _mp.process_map( - _mp.solve, - list(map(fn, self._srcfreq if srcfreq is None else srcfreq)), + data = { + 'model': self.model, + 'src': self.survey.sources[source], + 'receivers': self.survey.receivers, + 'frequencies': self.survey.frequencies, + 'empymod_opts': self._empymod_opts, + 'observed': None, + 'layered_opts': self.layered_opts, + 'gradient': gradient, + } + + # If there is data, add it. + if has_data: + data['observed'] = self.data.observed.loc[source, :, :] + + # For the gradient we also need the residuals and weights. + if gradient: + data['residual'] = self.data.residual.loc[source, :, :] + data['weights'] = self.data.weights.loc[source, :, :] + + return data + + # Compute responses in parallel. + out = _mp.process_map( + _mp.layered, + list(map(collect_empymod_inputs, self.survey.sources.keys())), max_workers=self.max_workers, - **{'desc': description, **self._tqdm_opts}, + **{'desc': 'Compute empymod', **self._tqdm_opts}, ) + # If gradient, return it. + if gradient: + + # Pre-allocate the gradient on the mesh. + grad = np.zeros((3, *self.model.grid.shape_cells), order='F') + + # Sum all sources up and return. + for val in out: + grad += val + + return grad + + # If forward, store responses in synthetic. + else: + + # Loop over sources to extract and store. + for i, src in enumerate(self.survey.sources.keys()): + self.data['synthetic'].loc[src, :, :] = out[i] + # OPTIMIZATION @property def gradient(self): @@ -850,6 +967,11 @@ def gradient(self): without relative electric permittivity nor relative magnetic permeability. + .. note:: + + If ``layered=True``, the gradient is computed using finite + differences. + Returns ------- @@ -863,71 +985,73 @@ def gradient(self): """ if self._gradient is None: - # Warn that cubic is not good for adjoint-state gradient. - if self.receiver_interpolation == 'cubic': - msg = ( - "emg3d: Receiver responses were obtained with cubic " - "interpolation. This will not yield the exact gradient. " - "Change `receiver_interpolation='linear'` in the call to " - "Simulation()." - ) - warnings.warn(msg, UserWarning) - - # Check limitation: No epsilon_r, mu_r. - var = (self.model.epsilon_r, self.model.mu_r) - for v, n in zip(var, ('el. permittivity', 'magn. permeability')): - if v is not None and not np.allclose(v, 1.0): - raise NotImplementedError( - f"Gradient not implemented for {n}." - ) - # Ensure misfit has been computed # (and therefore the electric fields). _ = self.misfit - # Compute back-propagating electric fields. - self._bcompute() + if self.layered: # 1D finite-difference gradient. + gradient = self._compute_1d(gradient=True) - # Inversion grid and its shape. - igrid = self.model.grid - ishape = igrid.shape_cells + else: # 3D adjoint-state gradient - # Pre-allocate the gradient on the mesh. - gradient = np.zeros((3, *ishape), order='F') + # Warn that cubic is not good for adjoint-state gradient. + if self.receiver_interpolation == 'cubic': + msg = ( + "emg3d: Receiver responses were obtained with cubic " + "interpolation. This will not yield the exact " + "gradient. Change `receiver_interpolation='linear'` " + "in the call to Simulation()." + ) + warnings.warn(msg, UserWarning) - # Loop over source-frequency pairs. - for src, freq in self._srcfreq: + # Check limitation: No epsilon_r, mu_r. + var = (self.model.epsilon_r, self.model.mu_r) + nam = ('el. permittivity', 'magn. permeability') + for v, n in zip(var, nam): + if v is not None and not np.allclose(v, 1.0): + raise NotImplementedError( + f"Gradient not implemented for {n}." + ) - efield = self._dict_get('efield', src, freq) - bfield = self._dict_get('bfield', src, freq) + # Compute back-propagating electric fields. + self._bcompute() - # Multiply forward field with backward field; take real part. - gfield = fields.Field( - grid=efield.grid, - data=np.real(bfield.field * efield.smu0 * efield.field), - ) + # Pre-allocate the gradient on the mesh. + gradient = np.zeros((3, *self.model.shape), order='F') + + # Loop over source-frequency pairs. + for src, freq in self._srcfreq: - # Pre-allocate the gradient for the computational grid. - shape = gfield.grid.shape_cells - grad = np.zeros((3, *shape), order='F') - - # Map the field to cell centers times volume. - cell_volumes = gfield.grid.cell_volumes - maps.interp_edges_to_vol_averages( - ex=gfield.fx, ey=gfield.fy, ez=gfield.fz, - volumes=cell_volumes.reshape(shape, order='F'), - ox=grad[0, ...], oy=grad[1, ...], oz=grad[2, ...]) - - # Bring gradient back from computation grid to inversion grid - # and add it to the total gradient. - if igrid != gfield.grid: - # Requires discretize; for this reason wrapped in own fct. - maps._interp_volume_average_adj( - oval=gradient, ogrid=igrid, - nval=grad, ngrid=gfield.grid, + efield = self._dict_get('efield', src, freq) + bfield = self._dict_get('bfield', src, freq) + + # Multiply forward field with backward; take real part. + gfield = fields.Field( + grid=efield.grid, + data=np.real(bfield.field*efield.smu0*efield.field), ) - else: - gradient += grad + + # Pre-allocate the gradient for the computational grid. + shape = gfield.grid.shape_cells + grad = np.zeros((3, *shape), order='F') + + # Map the field to cell centers times volume. + cell_volumes = gfield.grid.cell_volumes + maps.interp_edges_to_vol_averages( + ex=gfield.fx, ey=gfield.fy, ez=gfield.fz, + volumes=cell_volumes.reshape(shape, order='F'), + ox=grad[0, ...], oy=grad[1, ...], oz=grad[2, ...]) + + # Bring gradient back from computation grid to inversion + # grid and add it to the total gradient. + if self.model.grid != gfield.grid: + # Wrapped in own function as it requires discretize. + maps._interp_volume_average_adj( + oval=gradient, ogrid=self.model.grid, + nval=grad, ngrid=gfield.grid, + ) + else: + gradient += grad # Apply derivative-chain of property-map (only relevant if # `mapping` is something else than conductivity) and collect. @@ -1016,10 +1140,12 @@ def my_misfit_function(self): Value of the misfit function. """ + if self._misfit is None: # Ensure efields are computed - self.compute() + if not self._computed: + self.compute() # Check if weights are stored already. (weights are currently # simply 1/std^2; but might change in the future). @@ -1051,6 +1177,7 @@ def my_misfit_function(self): def _bcompute(self): """Compute bfields asynchronously for all sources and frequencies.""" + from emg3d import _multiprocessing as _mp # Initiate back-propagated electric field and info dicts. if not hasattr(self, '_dict_bfield'): @@ -1072,7 +1199,12 @@ def collect_bfield_inputs(inp): return self._data_or_file('bfield', source, freq, data) # Compute fields in parallel. - out = self._compute(collect_bfield_inputs, 'Back-propagate') + out = _mp.process_map( + _mp.solve, + list(map(collect_bfield_inputs, self._srcfreq)), + max_workers=self.max_workers, + **{'desc': 'Back-propagate', **self._tqdm_opts}, + ) # Loop over src-freq combinations to extract and store. for i, (src, freq) in enumerate(self._srcfreq): @@ -1147,6 +1279,12 @@ def jvec(self, vector): Shape of the data. """ + from emg3d import _multiprocessing as _mp + + if self.layered: + msg = "`jvec` is not implemented for `layered`." + raise NotImplementedError(msg) + # Missing for jvec/jtvec # - Refactor `compute/gradient/_bcompute/_get_rfield/jvec/jtvec`. # - Document properly jvec and jtvec. @@ -1221,7 +1359,12 @@ def collect_gfield_inputs(inp, vector=vector): return self._data_or_file('gfield', source, freq, data) # Compute fields (A^-1 * G * vector) in parallel. - out = self._compute(collect_gfield_inputs, 'Compute jvec') + out = _mp.process_map( + _mp.solve, + list(map(collect_gfield_inputs, self._srcfreq)), + max_workers=self.max_workers, + **{'desc': 'Compute jvec', **self._tqdm_opts}, + ) # Initiate jvec data with NaN's if it doesn't exist. if 'jvec' not in self.data.keys(): @@ -1318,6 +1461,21 @@ def _freq_inp2key(self, frequency): def _info_grids(self): """Return a string with "min {- max}" grid size.""" + info = "Gridding: " + + if self.layered: + info += "layered computation using method " + info += f"'{self.layered_opts['method']}'" + + if self.layered_opts['method'] in ['prism', 'cylinder']: + opts = '; '.join( + [f"{k}: {v:.2f}" for k, v in + self.layered_opts['ellipse'].items()] + ) + info += "; "+opts + + return info + # Single grid for all sources and receivers. if self.gridding in ['same', 'single', 'input']: grid = self.get_grid(*self._srcfreq[0]) @@ -1341,7 +1499,8 @@ def _info_grids(self): has_minmax = min_nc != max_nc # Assemble info. - info = f"{min_vc[0]} x {min_vc[1]} x {min_vc[2]} ({min_nc:,})" + info += f"{self._gridding_descr[self.gridding]}; " + info += f"{min_vc[0]} x {min_vc[1]} x {min_vc[2]} ({min_nc:,})" if has_minmax: info += f" - {max_vc[0]} x {max_vc[1]} x {max_vc[2]} ({max_nc:,})" @@ -1350,6 +1509,13 @@ def _info_grids(self): def print_grid_info(self, verb=1, return_info=False): """Print info for all generated grids.""" + # Act depending on gridding: + out = "" + + # If layered, return. + if self.layered: + return out if return_info else None + def get_grid_info(src, freq): """Return grid info for given source and frequency.""" grid = self.get_grid(src, freq) @@ -1359,9 +1525,6 @@ def get_grid_info(src, freq): out += grid.__repr__() return out - # Act depending on gridding: - out = "" - # Frequency-dependent. if self.gridding == 'frequency': for freq in self.survey.frequencies.values(): @@ -1394,13 +1557,13 @@ def get_grid_info(src, freq): def print_solver_info(self, field='efield', verb=1, return_info=False): """Print solver info.""" - # If not verbose, return. - if verb < 0: - return - # Get info dict. out = "" + # If not verbose or layered, return. + if verb < 0 or self.layered: + return out if return_info else None + # Loop over sources and frequencies. for src, freq in self._srcfreq: cinfo = self._dict_get(f"{field}_info", src, freq) @@ -1484,3 +1647,85 @@ def _set_model(self, model, kwargs): self.gridding_opts = gridding_opts self.model = model + + @property + def layered(self): + """If True, use layered computations (empymod).""" + return self._layered + + @layered.setter + def layered(self, layered): + """Update layered and therefore layered_opts.""" + self._set_layered_opts(layered, self.layered_opts) + + def _set_layered_opts(self, layered, layered_opts): + """Set self.layered and self.layered_opts.""" + + # Set layered. + self._layered = layered + + # If not layered, just store layered_opts and return. + if not layered: + self.layered_opts = layered_opts + return + + # Ensure we can handle the sources and receivers. + srlist = list(self.survey.sources.values()) + srlist = srlist + list(self.survey.receivers.values()) + for sr in srlist: + name = sr.__class__.__name__ + if 'Point' not in name and 'Dipole' not in name: + raise ValueError( + "Layered: Only Points and Dipoles supported, " + f"provided: {sr}!" + ) + + # Check limitation: Only isotropic and VTI. + if self.model.case not in ['isotropic', 'VTI']: + raise NotImplementedError( + f"Layered compute not implemented for {self.model.case} case." + ) + + # Make a copy of layered to not overwrite. + layered_opts = deepcopy(layered_opts) + + # Ensure method is defined; default: cylinder + layered_opts['method'] = layered_opts.get('method', 'cylinder') + + # For cylinder/prism, ensure there is ellipse['radius']. + if layered_opts['method'] in ['prism', 'cylinder']: + + # Initiate or get ellipse dict. + ellipse = layered_opts.get('ellipse', {}) + + # Try to estimate radius if not given. + if not ellipse.get('radius'): + + # Check if negz-cond is in gridding_opts. + try: + prop = self.gridding_opts['properties'] + prop = np.atleast_1d(prop) + m = getattr(maps, 'Map' + self.gridding_opts['mapping'])() + # Take the negative z property + ind = -1 if prop.size < 3 else -2 + cond = m.backward(prop[ind]) + + # If not, calculate it. + except (KeyError, TypeError): + zneg = self.model.property_x[:, :, 0] + cond = np.min(self.model.map.backward(zneg)) + + # Lowest frequency. + freq = min(self.survey.frequencies.values()) + + # Set the radius to one skin depth. + ellipse['radius'] = meshes.skin_depth(freq, cond) + + # Set factor/minor and store back. + ellipse['factor'] = ellipse.get('factor', 1.2) + ellipse['minor'] = ellipse.get('minor', 0.8) + layered_opts['ellipse'] = ellipse + + # Store layered options. + self.layered_opts = layered_opts + self._empymod_opts = {'verb': self.verb+1} diff --git a/tests/test_cli.py b/tests/test_cli.py index 8c344518..de620d80 100644 --- a/tests/test_cli.py +++ b/tests/test_cli.py @@ -89,6 +89,7 @@ class TestParser: 'forward': False, 'misfit': False, 'gradient': False, + 'layered': False, 'path': None, 'survey': None, 'model': None, @@ -178,6 +179,7 @@ def test_term_various(self, tmpdir): args_dict['nproc'] = -1 args_dict['verbosity'] = 20 args_dict['dry_run'] = True + args_dict['layered'] = True args_dict['clean'] = True args_dict['gradient'] = True args_dict['path'] = tmpdir @@ -190,6 +192,7 @@ def test_term_various(self, tmpdir): assert term['clean'] is True assert term['function'] == 'gradient' assert cfg['simulation_options']['max_workers'] == 1 + assert cfg['simulation_options']['layered'] is True assert cfg['files']['survey'] == join(tmpdir, 'testit.h5') assert cfg['files']['model'] == join(tmpdir, 'model.json') assert cfg['files']['output'] == join(tmpdir, 'output.npz') @@ -312,6 +315,41 @@ def test_simulation(self, tmpdir): sim_opts = cfg['simulation_options'] assert sim_opts['receiver_interpolation'] == 'cubic' + def test_layered(self, tmpdir): + + # Write a config file. + config = os.path.join(tmpdir, 'emg3d.cfg') + with open(config, 'w') as f: + f.write("[simulation]\n") + f.write("layered=True\n") + f.write("[layered]\n") + f.write("method=prism\n") + f.write("radius=3000\n") + f.write("factor=1.2\n") + f.write("minor=0.8\n") + f.write("check_foci=True\n") + f.write("merge=False") + + args_dict = self.args_dict.copy() + args_dict['config'] = config + args_dict['layered'] = None + cfg, term = cli.parser.parse_config_file(args_dict) + sim_opts = cfg['simulation_options'] + lay_opts = sim_opts['layered_opts'] + assert sim_opts['layered'] is True + assert lay_opts['ellipse']['radius'] == 3000.0 + assert lay_opts['ellipse']['factor'] == 1.2 + assert lay_opts['ellipse']['minor'] == 0.8 + assert lay_opts['ellipse']['check_foci'] is True + assert lay_opts['merge'] is False + + with pytest.raises(TypeError, match="Unexpected parameter in"): + with open(config, 'a') as f: + f.write("\nanother=True") + args_dict = self.args_dict.copy() + args_dict['config'] = config + _ = cli.parser.parse_config_file(args_dict) + def test_noise(self, tmpdir): # Write a config file. @@ -468,6 +506,7 @@ class TestRun: 'forward': False, 'misfit': False, 'gradient': True, + 'layered': False, 'path': None, 'survey': 'survey.npz', 'model': 'model.npz', @@ -505,7 +544,7 @@ class TestRun: model_vti = emg3d.Model(grid, 1., property_z=2.0) model_tri = emg3d.Model(grid, 1., 1.5, 2.0) - def test_basic(self, tmpdir, capsys): + def test_basic(self, tmpdir): # Store survey and model. self.survey.to_file(os.path.join(tmpdir, 'survey.npz'), verb=0) @@ -521,7 +560,6 @@ def test_basic(self, tmpdir, capsys): args_dict['path'] = tmpdir args_dict['config'] = 'bla' args_dict['verbosity'] = 2 - _, _ = capsys.readouterr() with pytest.raises(SystemExit) as e: cli.run.simulation(args_dict) @@ -656,7 +694,7 @@ def test_run(self, tmpdir): assert 'misfit' not in res3 assert 'gradient' not in res3 - # Redo for misfit, loading existing simulation. + # Redo for misfit, loading existing simulation, setting to layered. args_dict = self.args_dict.copy() args_dict['config'] = os.path.join(tmpdir, 'emg3d.cfg') args_dict['path'] = tmpdir @@ -664,9 +702,14 @@ def test_run(self, tmpdir): args_dict['misfit'] = True args_dict['gradient'] = False args_dict['dry_run'] = False + args_dict['layered'] = True # Change to layered! args_dict['load'] = 'mysim.npz' args_dict['output'] = 'output3.npz' cli.run.simulation(args_dict) + with open(os.path.join(tmpdir, 'output3.log'), 'r') as f: + log = f.read() + assert "Change «layered» of simulation to True." in log + assert "Gridding: layered computation using method 'cylinder'" in log res3 = emg3d.load(os.path.join(tmpdir, 'output3.npz')) assert 'misfit' in res3 assert 'gradient' not in res3 diff --git a/tests/test_multiprocessing.py b/tests/test_multiprocessing.py index 2b14d6bd..9ef46442 100644 --- a/tests/test_multiprocessing.py +++ b/tests/test_multiprocessing.py @@ -1,5 +1,7 @@ from os.path import join, dirname +import pytest +import numpy as np from numpy.testing import assert_allclose import emg3d @@ -10,6 +12,11 @@ except ImportError: tqdm = None +try: + import empymod +except ImportError: + empymod = None + # Data generated with tests/create_data/regression.py REGRES = emg3d.load(join(dirname(__file__), 'data', 'regression.npz')) @@ -65,3 +72,176 @@ def test__solve(): 'solver_opts': {'plain': True}} efield, info = _mp.solve(inp) assert_allclose(dat['Fresult'].field, efield.field) + + +@pytest.mark.skipif(empymod is None, reason="empymod not installed.") +def test_layered(): + + src = emg3d.TxElectricDipole((-2000, 0, 200, 20, 5)) + off = np.array([2000, 6000]) + rz = -250 + rec = emg3d.surveys.txrx_coordinates_to_dict( + emg3d.RxElectricPoint, (off, 0, rz, 0, 0)) + freqs = [1.0, 2.0, 3.0] + + grid = emg3d.TensorMesh( + h=[[4000, 4000, 4000], [2000, 2000], [100, 250, 100]], + origin=(-4000, -2000, -350), + ) + res = [2e14, 0.3, 1.0] + model = emg3d.Model(grid, property_x=1.0, property_z=1.0) + model.property_x[:, :, 1] = res[1] + model.property_z[:, :, 1] = res[1] + model.property_x[:, :, 2] = res[0] + model.property_z[:, :, 2] = res[0] + + obs = empymod.bipole( + src=src.coordinates, + rec=(off, off*0, rz, 0, 0), + depth=[0, -250], + res=[res[0], res[1], 2*res[2]], + freqtime=freqs, + verb=1, + ) + syn = empymod.bipole( + src=src.coordinates, + rec=(off, off*0, rz, 0, 0), + depth=[0, -250], + res=res, + freqtime=freqs, + verb=1, + ) + + survey = emg3d.Survey(src, rec, freqs) + survey.data.observed[0, ...] = obs.T + survey.data['synthetic'] = survey.data['observed'].copy() + survey.data.synthetic[0, ...] = obs.T + + inp = { + 'model': model, + 'src': survey.sources['TxED-1'], + 'receivers': survey.receivers, + 'frequencies': survey.frequencies, + 'empymod_opts': {'verb': 1}, + 'observed': None, + 'layered_opts': {'method': 'receiver'}, + 'gradient': False, + } + + # Forward - generate data + out = _mp.layered(inp) + assert_allclose(out, syn.T) + + # Gradient without required data. + inp['gradient'] = True + out = _mp.layered(inp) + assert_allclose(out, 0.0) + + # Gradient but all-NaN obs. + inp['observed'] = survey.data.observed[0, :, :]*np.nan + inp['weights'] = survey.data['observed'].copy()*0+1 + inp['weights'] = inp['weights'][0, :, :] + inp['residual'] = survey.data.synthetic - survey.data.observed + inp['residual'] = inp['residual'][0, :, :] + out = _mp.layered(inp) + assert_allclose(out, 0.0) + + # Gradient + inp['observed'] = survey.data.observed[0, :, :] + out = _mp.layered(inp) + # Only checking locations, not actual gradient + # (that is done in _fd_gradient) + assert_allclose(out[1, ...], 0.0) + assert_allclose(out[::2, 0, :, :], 0.0) + assert_allclose(out[::2, :, 0, :], 0.0) + assert np.all(out[::2, 1:, 1, :] != 0.0) + + +@pytest.mark.skipif(empymod is None, reason="empymod not installed.") +def test_empymod_fwd(): + # Simple check for status quo. + empymod_inp = { + 'src': [0, 0, -150, 40, -20], + 'rec': [4000, 0, -200, 13, 58], + 'depth': [0, -200, -1000], + 'freqtime': [0.01, np.pi, 10], + } + res = np.array([2e14, 0.3, 1, 2]) + aniso = np.array([1., 1., np.sqrt(2), 2]) + resp1 = empymod.bipole(res=res, aniso=aniso, **empymod_inp) + resp2 = _mp._empymod_fwd(1/res, 1/(aniso**2 * res), empymod_inp) + + assert_allclose(resp1, resp2) + + +def test_get_points(): + + class DummySrcRec: + pass + + src = DummySrcRec() + src.center = [1, 2] + rec = DummySrcRec() + rec.center = [10, 20] + + # Method which is not source/receiver: + out = _mp._get_points('cylinder', src, rec) + assert out['method'] == 'cylinder' + assert out['p0'] == src.center + assert out['p1'] == rec.center + + # Method source: + out = _mp._get_points('source', src, rec) + assert out['method'] == 'midpoint' + assert out['p0'] == src.center + assert out['p1'] == src.center + + # Method receiver: + out = _mp._get_points('receiver', src, rec) + assert out['method'] == 'midpoint' + assert out['p0'] == rec.center + assert out['p1'] == rec.center + + +@pytest.mark.skipif(empymod is None, reason="empymod not installed.") +def test_fd_gradient(): + + res = np.array([0.9876, ]) + + empymod_inp = { + 'src': (0, 0, 0, 20, 5), + 'rec': (100, 0, 0, -5, -5), + 'depth': [], + 'freqtime': 1.0, + 'verb': 1, + } + + obs = empymod.bipole(res=res, **empymod_inp) + + d = 1/res*0.0001 + dobs_h = empymod.bipole(res=d+res, **empymod_inp) + dobs_v = empymod.bipole( + res=res, aniso=np.sqrt((d+res)/res), **empymod_inp) + + weight = np.pi + misfit = np.e + inp = { + 'data': obs, + 'weight': weight, + 'misfit': misfit, + 'empymod_inp': empymod_inp, + 'imat': np.array([[0.5, 0.0], [0.0, 0.0], [.1, 0.4]]), + } + + rh = dobs_h-obs + gh = (np.real(weight*(rh.conj()*rh)/2) - misfit)/d + rv = dobs_v-obs + gv = (np.real(weight*(rv.conj()*rv)/2) - misfit)/d + + cond = 1/res + ghc = _mp._fd_gradient(cond_h=cond, cond_v=None, vertical=False, **inp) + gvc = _mp._fd_gradient(cond_h=cond, cond_v=cond, vertical=True, **inp) + + assert_allclose(gh, 2*ghc[0, 0]) + assert_allclose(gv, 10*gvc[2, 0]) + assert_allclose(gv, 2.5*gvc[2, 1]) diff --git a/tests/test_simulations.py b/tests/test_simulations.py index d2a6c030..0ef2254c 100644 --- a/tests/test_simulations.py +++ b/tests/test_simulations.py @@ -6,7 +6,7 @@ import emg3d from emg3d import simulations -from . import alternatives +from . import alternatives, helpers # Soft dependencies @@ -541,6 +541,160 @@ def test_tqdm(self): assert sim._tqdm_opts == tqdm_opts +@pytest.mark.skipif(xarray is None, reason="xarray not installed.") +class TestLayeredSimulation(): + if xarray is not None: + # Create a simple survey + + # Sources: 1 Electric Dipole, 1 Magnetic Dipole. + s1 = emg3d.TxElectricDipole((0, 1000, -950, 0, 0)) + s2 = emg3d.TxMagneticDipole((0, 3000, -950, 90, 0)) + sources = emg3d.surveys.txrx_lists_to_dict([s1, s2]) + + # Receivers: 1 Electric Point, 1 Magnetic Point + e_rec = emg3d.surveys.txrx_coordinates_to_dict( + emg3d.RxElectricPoint, + (np.arange(6)*1000, 0, -1000, 0, 0)) + m_rec = emg3d.surveys.txrx_coordinates_to_dict( + emg3d.RxMagneticPoint, + (np.arange(6)*1000, 0, -1000, 90, 0)) + receivers = emg3d.surveys.txrx_lists_to_dict([e_rec, m_rec]) + + # Frequencies + frequencies = (1.0, 2.0) + + survey = emg3d.Survey( + sources, receivers, frequencies, name='TestSurv', + noise_floor=1e-15, relative_error=0.05) + + # Create a simple grid and model + grid = emg3d.TensorMesh( + [np.ones(32)*250, np.ones(16)*500, np.ones(16)*500], + np.array([-1250, -1250, -2250])) + model = emg3d.Model(grid, 1) + + # Create a simulation, compute all fields. + simulation = simulations.Simulation( + survey, model, name='TestSim', max_workers=1, gridding='single', + layered=True, layered_opts={'ellipse': {'minor': 0.99}} + ) + + # Compute observed. + simulation.compute(observed=True, min_offset=1100) + del simulation.survey.data['synthetic'] + + def test_reprs(self): + test = self.simulation.__repr__() + + assert "Simulation «TestSim»" in test + assert "Survey «TestSurv»: 2 sources; 12 receivers; 2 frequenc" in test + assert "Model: resistivity; isotropic; 32 x 16 x 16 (8,192)" in test + assert "Gridding: layered computation using method 'cylinder'" in test + + test = self.simulation._repr_html_() + assert "Simulation «TestSim»" in test + assert "Survey «TestSurv»: 2 sources; 12 receivers; 2" in test + assert "Model: resistivity; isotropic; 32 x 16 x 16 (8,192)" in test + assert "Gridding: layered computation using method 'cylinder'" in test + + # Also _info_grids + txt = ("Gridding: layered computation using method 'cylinder'; " + "minor: 0.99; radius: 503.29; factor: 1.20") + assert txt == self.simulation._info_grids + + def test_copy(self, tmpdir): + + sim2 = self.simulation.copy() + assert self.simulation.layered == sim2.layered + assert helpers.compare_dicts(self.simulation.layered_opts, + sim2.layered_opts) + assert sim2.layered_opts['ellipse']['minor'] == 0.99 + + def test_layered_opts(self): + sminp = {'survey': self.survey, 'model': self.model} + + # layered=False: layered_opts just stored as is + inp = {'ellipse': {'minor': 0.99}} + simulation = simulations.Simulation(**sminp, layered_opts=inp) + assert helpers.compare_dicts(simulation.layered_opts, inp) + + # layered=True: layered_opts complete + simulation.layered = True + assert simulation.layered_opts['ellipse']['factor'] == 1.2 + + # If method not in 'prism'/'cylinder': layered_opts just stored as is + inp = {'method': 'aoeuaoeu', 'whatever_else': 'yes'} + sim = simulations.Simulation(**sminp, layered=True, layered_opts=inp) + assert sim.layered_opts['method'] == 'aoeuaoeu' + assert sim.layered_opts['whatever_else'] == 'yes' + + # no gridding -> estimate from model + inp = {'method': 'prism'} + sim = simulations.Simulation( + **sminp, gridding='same', layered=True, layered_opts=inp) + assert sim.layered_opts['method'] == 'prism' + assert round(sim.layered_opts['ellipse']['radius']) == 503 + + # Check defaults for cylinder. + inp = {'ellipse': {'check_foci': True, 'merge': True}} + sim = simulations.Simulation(**sminp, layered=True, layered_opts=inp) + assert sim.layered_opts['method'] == 'cylinder' + assert round(sim.layered_opts['ellipse']['radius']) == 503 + # Skindepth 1 Hz; 1 Ohm.m + assert sim.layered_opts['ellipse']['factor'] == 1.2 + assert sim.layered_opts['ellipse']['minor'] == 0.8 + assert sim.layered_opts['ellipse']['merge'] + + # Ensure it doesn't change. + inp = {'mapping': 'Conductivity', + 'properties': [1/0.3, 1.0, 1.0, 1e-8]} + inp = {'ellipse': {'factor': 2.0, 'minor': 0.5}} + sim = simulations.Simulation(**sminp, layered=True, layered_opts=inp) + assert sim.layered_opts['method'] == 'cylinder' + assert round(sim.layered_opts['ellipse']['radius']) == 503 + # Skindepth 1 Hz; 1 Ohm.m + assert sim.layered_opts['ellipse']['factor'] == 2.0 + assert sim.layered_opts['ellipse']['minor'] == 0.5 + + def test_print_info(self): + assert self.simulation.print_solver_info(return_info=True) == "" + assert self.simulation.print_grid_info(return_info=True) == "" + + def test_not_implemented_things(self): + with pytest.raises(NotImplementedError, match="for `layered`"): + self.simulation.jvec('vector') + + with pytest.raises(NotImplementedError, match="No fields if `laye"): + self.simulation.get_efield('TxED-1', 'f-1') + + # Survey with unsupported electrodes. + survey = emg3d.Survey( + sources=emg3d.TxElectricWire([[0, 0, 0], [2, 2, 2]]), + receivers=emg3d.RxElectricPoint((1000, 500, 100, 0, 0)), + frequencies=1.0, + ) + with pytest.raises(ValueError, match='Layered: Only Points and D'): + simulations.Simulation(survey, self.model, layered=True) + + model = emg3d.Model(self.grid, 1.0, 2.0, 3.0) + with pytest.raises(NotImplementedError, match='Layered compute not'): + simulations.Simulation(self.survey, model, layered=True) + + def test_compute_1d(self): + + # Just checking the workflow, not the actual result. + assert not np.isfinite(self.simulation.data.synthetic.data).sum() > 0 + assert not self.simulation._computed + self.simulation.compute() + assert self.simulation._computed + assert np.isfinite(self.simulation.data.synthetic.data).sum() > 0 + + assert self.simulation.misfit > 0.0 + + grad = self.simulation.gradient + assert grad.shape == self.model.shape + + @pytest.mark.skipif(xarray is None, reason="xarray not installed.") def test_misfit(): data = 1 @@ -645,9 +799,13 @@ def test_errors(self): sim_inp = {'survey': survey, 'gridding': 'same', 'receiver_interpolation': 'linear'} + # Dummy misfit, so it doesn't want to compute it. + misfit = np.sum(self.survey.data.observed) + # Model with electric permittivity. simulation = simulations.Simulation( model=emg3d.Model(mesh, epsilon_r=3), **sim_inp) + simulation._misfit = misfit with pytest.raises(NotImplementedError, match='for el. permittivity'): simulation.gradient @@ -655,6 +813,7 @@ def test_errors(self): simulation = simulations.Simulation( model=emg3d.Model(mesh, mu_r=np.ones(mesh.shape_cells)*np.pi), **sim_inp) + simulation._misfit = misfit with pytest.raises(NotImplementedError, match='for magn. permeabili'): simulation.gradient