Changelog

Various bug fixes and quality-of-life improvements when using a SphericalChunkDomain together with a spherical polygon to further constrain the spatial extent. Most importantly conversions between WGS84 and fully spherical coordinates are now consistent with the rest of Salvus. The in-notebook domain plot now also works and shows the original spherical chunk together with the domain polygon.

Add the sampling interval for point output to the WaveformSimulationConfiguration.

Fix a bug where side-sets involving fluid layers on UTM domains where incorrectly attached.

The polynomial degree of the spectral-element basis is now configurable in the WaveformSimulationConfiguration.

Copy
w = WaveformSimulationConfiguration(spectral_element_order=6)

If not specified, the default remains at 4 as in all previous releases.

A new parameter can optionally be passed to a cartesian SurfaceTopography object that restricts the mesh deformation above a given z-value. This is useful, for instance, when placing a layer (i.e. the atmosphere) above a free surface with topography. Related constructors (such as SurfaceTopography.from_gmrt_file() optionally take this parameter as well.

Bug fix for the gradient widget of the inversion dashboard, where previously some event-dependent gradients were not displayed correctly.

Bug fix in the SimulationConfiguration to ensure the extra_output_configuration passed to launch overwrites additional output values defined in the WaveformSimulationConfiguration. The output options in the WaveformSimulationConfiguration will be deprecated in a future release.

Add a new tab to the inversion dashboard to visualize the progress of model updates and trust-region radius.

Add a toggle button to the misfit widget to switch between absolute and relative misfits - normalized per iteration or event, respectively.

Much faster seismological window picking for some cases. Additionally it is now a bit more conservative and will not pick some sketchy windows it did not pick before. All in all the picked windows should be very similar to results from the previous version of this algorithm but a few less windows will be picked.

The p.waveforms.has_data() method now also works for synthetic data.

The partial copy of the Salvus job array results has been integrated into all parts of SalvusProject. If one or more jobs fail for a forward or an adjoint simulation, it will still raise an exception but copy all successfully completed jobs to the project folder. Any subsequent simulations launched will be aware of this and not submit the already run ones.

This makes it much easier to recover from sporadic cluster failures without having to directly modify the project directory.

Add memory_per_rank_in_MB as optional parameter to SiteConfig. This allows to adjust the memory settings for different machines during an inversion.

Improve memory management of checkpoints in SalvusProject.

This gives more flexibility to tune the strategy regarding checkpoints. By default, checkpoints will be stored for every forward run, and automatically removed after (a) gradients for the simulation configuration have been computed and/or (b) the model got rejected by the trust region algorithm.

This behavior can be modified using optional boolean flags store_checkpoints and cleanup_checkpoints in the job_submission settings.

Copy
job_submission={
  "forward": ...,
  "preconditioner": ...,
  "store_checkpoints": False,
  "cleanup_checkpoints": False,
}

If the former is False checkpoints will only be stored when gradients are requested. If the latter is False checkpoints will never be deleted automatically.

The _plot method of an EventData object now takes an additional optional parameter: _exclude_cbar. If set to true, a colorbar won't be added to the axis if a shotgather plot is requested. This makes it easier to produce shotgathers for multi-component data (i.e. for strain output).

Fix a bug in plotting gradients in the inversion dashboard or using p.viz.nb.gradients.

Enable stopping criteria to limit the maximum number of iterations in an inverse problem. This is mainly useful for automated monolithic inversion workflows.

Copy
p.inversions.set_stopping_criteria(
  inverse_problem_configuration="...",
  criteria={
    "max_iterations_global": 10,
    "max_iterations": 5,
  }
)

The criteria can be overwritten anytime to continue with the inversion. global settings apply to the entire inversion tree, the other settings apply to individual branches.

Remove some path dependencies which previously prevented renaming the project.

Various improvements to make 2-D seismological domains more useful.

The interpolation of volumetric models onto seismological meshes was previously only available for 3-D domains. Support is now extended to 2-D domains, with "longitude" and "radius" OR "depth" being required coordinate axes in the parent xarray Dataset. A new tutorial (teleseismic_2d) has been added, fleshing out this feature.

Add the ability to ignore elements flagged with a certain value during volume model interpolation for seismological models. This is useful, for instance, if one wants to limit interpolation to a specific region (say, the elastic subsurface, and not any water layers). See the description in the seismology GenericModel constructor for more details. Crustal and Mantle models will now ignore fluid elements by default, to ensure that elastic models are not interpolated into oceans. Is this behavior is undesirable, one has more control with the GenericModel interface.

When re-launching simulations for which data already exist (e.g. to compute checkpoints), previously computed outputs are kept in the project file structure and only overwritten once the new jobs have finished. Additionally, a utility to manually clean the simulation store has been added.

The internal workflow for handling checkpoints and adjoint simulations has been improved to avoid duplicated simulations for missing checkpoints on remote sites. This changes the hash values of adjoint simulations, which do not depend on synthetic receiver files anymore.

The waveform plotting widget can now also plot windows without a reference time.

The p.simulations.launch() method can now take extra output configurations. These can be any setting that does not modify the receiver output which currently are surface and volume data settings as well as the memory per rank buffer settings.

Copy
p.simulations.launch(
    simulation_configuration="basic",
    events=p.events.list(),
    site_name="my_computer",
    ranks_per_job=12,
    extra_output_configuration={
        "volume_data": {
            "sampling_interval_in_time_steps": 50,
            "fields": ["displacement"],
        },
        "surface_data": {
            "sampling_interval_in_time_steps": 20,
            "fields": ["acceleration", "velocity"],
            "side_sets": ["x0", "x1"],
        },
        "memory_per_rank_in_MB": 2000.0,
    },
)

The Project.from_domain() constructor now takes an optional argument load_if_exists which defaults to False. If this is set to True, and if the domains are identical, future runs of this constructor will simply load the existing project from disk rather than raising a "project already exists" error. This helps remove some boilerplate code which the user previously had to write.

Add widget to query status of simulations and and improve the widgets of waveforms and inversions.

Minor bug fix in visualizing 2D domains using p.viz.nb.domain that did not display some source types correctly.

Previous to this change, the parameters governing the behavior of the linear solids in viscous simulations would only be properly attached if they came in via a background model. Now these parameters can be supplied to a ModelConfiguration via the LinearSolids class (which can be found in salvus.namespace). Additional validations and checks have also been added to help the user when viscous simulations are either requested or thought to be desired.

Automatically clean-up checkpoint files after an adjoint simulation.

Add utilities to handle remote output data and to manually delete jobs in SalvusProject.

Chaining processing function in SalvusProject now works as expected.

Fix a bug in the horizontal spacing when generating 2D layered meshes from splines and extruding horizontally for absorbing boundaries.

Add a function derive_bm_file() to derive a 1-D background model from a 2- or 3-D Cartesian volume model. This is useful, for instance, when the 3-D model has a significant increase in velocity with depth, and when no appropriate 1-D model exists a-priori. If the velocity profile allows, using the 1-D model output by this function may result in the insertion of doubling layers when a simulation mesh is created, and as a consequence may allow for a reduction of computational cost.

Minor modifications to create_events_from_segy_file to handle 2-D events.

Implement p.viz.nb.domain() method for 2D box domains.

Move job submission settings to the InverseProblemConfiguration and enable different site configurations for simulations and the preconditioner. The previous way of providing site_name, ranks_per_job and wall_time_per_job_in_seconds through the iterate and resume functions still work, but will be deprecated in the future.

Add interactive misfit plots to the inverse problem Jupyter notebook widget.

Better workflow management for misfit and gradient computations. The newly added functionality ensures consistent settings for site name and the number of ranks and automatically resubmits simulations if the checkpoints are no longer available.

With very large models and / or meshes, interpolating volumetric models can take some time. Here a 'verbose' argument is now accepted by toolbox.interpolate_spherical, which will print a progress bar if the verbosity is 1 or greater. The verbose behavior is the default when generating spherical meshes with project.

Better error handling and management of entities using UnstructuredMeshSimulationConfiguration objects.

Adding an already existing inverse problem configuration again with exactly the same parameters will no longer throw.

Projects created with an older Salvus version using the windows component cannot be imported anymore.

The windows component has been removed in favor of a more generic DataSelectionConfiguration. This is for one more in line with the rest of SalvusProject but also more powerful and flexible.

A direct consequence is that any function or object that took a window_set_name argument now takes a data_selection_configuration argument.

A DataSelectionConfiguration is either a function that is applied on the fly to select pieces of data (think generalized windows and weights) or one EventWindowAndWeightSet set per event (or a combination of both).

Projects that used an on-the-fly temporal weighting/windowing function now have to use it by creating a new DataSelectionConfiguration with an attached data_selection_function. The function itself does not have to be changed.

Copy
p.add_to_project(
    sn.DataSelectionConfiguration(
        name="mute_with_direct_phase",
        receiver_field="phi",
        data_selection_function=mute_with_direct_phase,
    )
)

The pick_windows() seismological action also requires changing window_set_name to data_selection_configuration (if that configuration does not yet exist, it will be created). In addition the window_taper_width_in_seconds argument is now exposed and required. It controls the width of the taper for the windows.

Copy
p.actions.seismology.pick_windows(
    data_selection_configuration="custom_windows",
    ...
    window_taper_width_in_seconds=10.0,
)

Salvus contains a lot of convenience functions to pick windows and weight individual receivers but it is also possible to directly operate on these data structures for full control.

Copy
# Get the data selection configuration.
dsc = p.entities.get(
    entity_type="data_selection_configuration", entity_name="custom_windows"
)
# If interval windows have been picked before an event, its
# `EventWindowAndWeightSet` can be retrieved.
ewws = dsc.get_event_window_and_weight_set(event=p.events.get("event_0000"))

This can than be used to exert detailed control on everything. Once this has been done it will be applied in all subsequent steps that utilize that particular DataSelectionConfiguration.

Copy
In  [1]: ewws.receivers
Out [1]:
{'XX.A0000.': {'receiver_weight': 1.5462720296177797,
  'components': {'Z': {'component_weight': 1.0,
    'windows': [{'window_start': 229.67554,
      'window_end': 1389.219377,
      'window_weight': 1.0}]},
...

In [2]: ewws.receivers["XX.A0000."]["receiver_weight"] = 0.7

In [3]: ewws.write(overwrite=True)

Can now add weights to individual receivers and components which will then be applied to each misfit and adjoint source computation. This could have been previously done with a custom windowing function with weights but this is more convenient and easier to control. It is detailed in other changelog items of this release.

It can either happen via window picking and receiver weighting or directly in the data selection function:

Copy
def data_selection_function(st, receiver, sources):

    # Receiver and component weights.
    weights = {
        # This is optional but allows to specify additional weights
        # important for the misfit and adjoint source computations.
        "receiver_and_component_weights": {
            # Pass the receiver weight here.
            "receiver_weight": 0.75,
            # It is also possible to give weights to individual
            # receivers.
            "component_weights": {"X": 1.0, "Y": 0.8", "Z": 0.5},
        }
    }

    for tr in st:
        component = tr.stats.channel[-1]
        temporal_weights = compute_window(...)

        # Individual window weights.
        weights[component] = [
            {"values": temporal_weights, "misfit_weight": 0.8},
            ...
        ]

    return weights

# Add to the project.
p.add_to_project(
    sn.DataSelectionConfiguration(
        name="custom_data_selection",
        receiver_field="phi",
        data_selection_function=data_selection_function,
    ),
    overwrite=True,
)

Added implementations for a few common seismological receiver weighting schemes with a general and easy to extend interface.

Copy
# Custom receiver weighting function.
def my_favorite_station(event, receivers, receiver_name):
    weights = {}
    for name in receivers:
        if name == receiver_name:
            weight = 2.0
        else:
            weight = 1.0
        weights[name] = weight
    return weights


p.actions.seismology.add_receiver_weights_to_windows(
    # Add the weights to existing windows in this data selection
    # configuration.
    data_selection_configuration="custom_windows",
    events=p.events.list(),
    # The weights for each item in the chain will be multiplied together.
    weighting_chain=[
        {
            "weighting_scheme": "ruan_et_al_2019",
            "function_kwargs": {"ref_distance_condition_fraction": 1.0 / 3.0},
        },
        {
            "weighting_scheme": "mute_near_source_receivers",
            "function_kwargs": {
                # All receivers closer than this to the sources
                # will have zero weight.
                "minimum_receiver_distance_in_m": 200e3,
                # All receivers further away than this will not be muted.
                "maximum_receiver_mute_distance_in_m": 1000e3,
                # How to go from weight 0 to weight 1 for receivers in the
                # transition region. Currently only "hanning" is supported.
                "taper_type": "hanning",
            },
        },
        # Custom weighting function.
        {
            "weighting_scheme": my_favorite_station,
            "function_kwargs": {"receiver_name": "XX.A0044."},
        },
    ],
    # Normalize the sum of all receiver weights to the receiver count.
    normalize=True,
)

New seismological receiver weight map plot.

Copy
p.viz.seismology.receiver_weights_map(
    data_selection_configuration="custom_windows"
)

Various part of SalvusProject can serialize functions to disc. This is now more powerful and easier to use as the serialized functions can now access closure values as well as other variables in the global scope of the functions which will also be serialized alongside. Thus external variables and imports will now work with these functions.

Basically it means that it should now just work in more cases without any extra work or rewrites.

The p.viz.seismology.misfit_map() method is now a widget and the event can be selected interactively. Thus passing the event argument is no longer necessary.

New method to visualize and compare histograms of misfits.

Copy
p.viz.misfit_histogram(
    simulation_configuration_a="initial_model_70_120_seconds",
    simulation_configuration_b=p.inversions.get_simulation_name(
            inverse_problem_configuration="inversion", iteration_id=25
    ),
    misfit_configuration="phase_misfit_70_120_seconds",
    events=p.events.list(),
    merge_all_components=False
)

A few new utilities to analyze the statistics of interval windows.

One is the ability to get a pandas.DataFrame object containing statistics about all windows for a list of events. This is very useful to do some custom analysis and figure out which events to use in the final inversion.

Copy
df = p.entities.get(
    entity_type="data_selection_configuration", entity_name="selection_a"
).get_interval_window_statistics_table(events=p.events.get_all())

A styled version of this table with all events in a project can be visualized with:

Copy
p.viz.interval_window_statistics("selection_a")

Warn when receiver fields other than "displacement" or "phi" are used in a MisfitConfiguration object. Otherwise the adjoint sources are currently not fully consistent.

New p.misfits.get_misfit_comparison_table() method yielding a pandas.DataFrame containing detailed per-receiver misfit information that can be used for further custom analysis.

Copy
In [1]: p.misfits.get_misfit_comparison_table(
            reference_data="initial_model",
            other_data=[
                p.inversions.get_simulation_name(
                    inverse_problem_configuration="inv", iteration_id=4
                )
            ],
            misfit_configuration="L2-misfit-to-target-model",
            event="event_0000",
        )

Out [1]:
           initial_model (ref)  inv_it_3__trial   Reduction inv_it_3__trial
XX.A0000.      1.266640e-03     0.000224           0.001043
XX.A0001.      1.438440e-03     0.000361           0.001078
XX.A0002.      1.450685e-03     0.000369           0.001081
XX.A0003.      1.607995e-03     0.000430           0.001178
XX.A0004.      9.473526e-04     0.000238           0.000709
...                     ...          ...                ...
XX.A0095.      1.182061e-05     0.000065          -0.000053
XX.A0096.      1.360653e-06     0.000529          -0.000527
XX.A0097.      9.762485e-07     0.000102          -0.000101
XX.A0098.      4.663780e-06     0.000306          -0.000301
XX.A0099.      1.784453e-06     0.000165          -0.000163

[100 rows x 3 columns]

New p.viz.nb.misfit_comparison() method to display the output from the p.misfits.get_misfit_comparison_table() method in a pretty table.

Copy
p.viz.nb.misfit_comparison(
    reference_data="initial_model",
    other_data=[
        p.inversions.get_simulation_name(
            inverse_problem_configuration="inv", iteration_id=4
        )
    ],
    misfit_configuration="L2-misfit-to-target-model",
    event=p.events.list()[0],
)

New p.viz.seismology.misfit_map() method to display the output from the p.misfits.get_misfit_comparison_table() method on a geographical map.

Copy
p.viz.seismology.misfit_map(
    reference_data="initial_model",
    compare_data=p.inversions.get_simulation_name(
            inverse_problem_configuration="inv", iteration_id=5
        ),
    misfit_configuration="L2-misfit-to-target-model",
    event=p.events.list()[0],
)

Can now delete entities, even if they cannot be instantiated anymore. Useful for certain edge cases.

Fully consistent is_point_in_domain() check for spherical chunk domains.

The map visualization for spherical chunks now also works for domains crossing the international date line.

Copy
sn.domain.dim3.SphericalChunkDomain(
    lat_center=52.0,
    lat_extent=16.0,
    lon_center=170.0,
    lon_extent=66.0,
    radius_in_meter=6371000.0,
).plot()

A few quality-of-life improvements for inversions.

Uppercase variables names during model interpolation from xarray data sets and NetCDF files. Thus lowercase case variables now also work.

Several improvements and small bug fixes in a number of Jupyter notebook inversion widgets.

New inversion action to faciliate model smoothing.

Copy
p.actions.inversion.smooth_model(
    model=gradient,
    smoothing_configuration=sn.ConstantSmoothing(
        smoothing_lengths_in_meters={"VP": 0.01, "RHO": 0.01,},
    ),
    ranks_per_job=4,
    site_name="local",
)

Enable misfit configurations without observed data.

Prevent processing/misfit functions that cannot be deserialized again from being stored in a project in the first place.

New experimental 2-D circular domain.