Running#

Restarting a crashed simulation#

If the reload_results argument is set to true, a previously interrupted response-functions simulation will skip already computed response functions and resume from where it stopped.

The observation operator plugin autorestart input argument can also be used to resume the simulation that was interrupted.

Exploit results#

The following code examples can be used to read and use the response function results:

\(\mathbf{H}\) matrix#

Reading the \(\mathbf{H}\) matrix (with control vector coordinates) from the observation vector decomposition.

Note

The response function simulation must be run with the dump_obsvect_decompostion argument set to true for the CIF to produce the required output files.

import os

import numpy as np
import pandas as pd
import xarray as xr
from pycif.utils.sparse_array import to_dense_dataset

def read_decomp(path, region_names=None, use_sparse_dtype=False):
    with xr.open_dataset(path) as ds:
        if 'sparse_sim_tl' in ds:
            ds = to_dense_dataset(ds)

        ds = ds.isel(vresol=0).drop_vars('vresol')  # Dropping vertical levels

        # Replacing hresol with region names
        if region_names is not None and 'region_id' in ds:
            ds['hresol'] = (['hresol'], region_names[ds.region_id.values])

        decomp = ds.sim_tl.to_series().unstack(['date', 'hresol'])

    if use_sparse_dtype:
        decomp = decomp.astype(pd.SparseDtype("float", 0.0))

    return decomp

def read_h_matrix(workdir, region_names=None):
    h_dict = {}
    dirname = os.path.join(workdir, "base_functions/decomposition")

    for component in os.listdir(dirname):
        comp_dirname = os.path.join(dirname, component)
          for parameter in os.listdir(comp_dirname):
              print(f"reading H matrix for {component}/{parameter}")
              path = os.path.join(comp_dirname, parameter, "obsvect_decomp.nc")
              h_dict[(component, parameter)] = read_decomp(path, region_names, True)

    df = pd.concat(
        h_dict.values(),
        keys=h_dict.keys(),
        names=['component', 'parameter'],
        axis='columns',
    )

    return df


with xr.open_dataset("region_mask.nc") as ds_mask:
    region_names = ds_mask.regions_name.values

h_matrix = read_h_matrix("path/to/workdir", region_names)

Variational inv#

  • TODO: Output

  • TODO: Rerun with another B, R