Modifying a model to use with EnSRF

Modifying a model to use with EnSRF#

In the following tutorial, one will learn how to modify a model plugin to make it work in batch mode with the EnSRF running mode.

Set-up a working case with independent sampling#

The EnSRF methods rely on generating a random ensemble to approximately solve the inversion equations. There are two possible methods to:

  • batch_sampling: run all samples altogether into a single job; to do so, the control vector is perturbed using virtual species that are propagating all the way to the observation equivalents in the observation operator;

    Basically, for a case when CH4 emissions are optimized, the batch_sampling method will generate an ensemble (flux, CH4_1), …, (flux, CH4_N) of perturbed inputs. Then, perturbations will be propagated through all pyCIF operations, including the CTM. The outputs will be (concs, CH4_1), …, (concs, CH4_N). This will be re-composed into separated monitor.nc files for each samples.

  • parallel_sampling: the ensemble is generated by the master job; then slave jobs are run independently for each perturb control vector. And corresponding monitor.nc files are generated and processed by the master job.

The parallel_sampling mode do not need any modification to work. One can test it directly with a working inversion test case.

Examples of inversion configuration files can be found here. Therein, one can look for EnSRF cases for the dummy model. General inversion configurations can be found in all sub-directories as adjtltest examples.

Implementing the batch_sampling for a specific model#

Creating a restart_inicond component#

At the beginning of a simulation, CIF typically reads a prescribed file with initial conditions for the transported species (inicond). With the EnSRF mode, the first cycle starts from these prescribed initial conditions. However, the subsequent cycles must restart from the posterior concentrations of the previous cycle. To this end, when starting the next cycles, CIF adds a component restart_inicond to the Yaml configuration file controlling the cycle. In addition, it adds a parameter called ensrf_restart_file providing the path to the restart file to use and it also removes any occurence of the inicond component.

This restart_inicond component must be accepted by the model and processed the same way the restart files are processed for forward simulations. Some model routines must be slightly modified to ensure this:

  1. In the function ini_mapper, the restart_inicond component must be added.

    Here is an example for ICON-ART:

    if getattr(model, "ensrf_restart_file", False):
    for spec in acspecies:
        mapper["inputs"][('restart_inicond', spec)] = {
            "input_dates": {model.datei: np.array([[model.datei, model.datei]])},
            "domain": domain,
            "force_dump": True,
            "sparse_data": False,
            "break_fwd_onlyinit_pipe": False,
            "break_adj_onlyinit_pipe": False
        }

else:
    for spec in acspecies:
        mapper["inputs"][('inicond', spec)] = {
            "input_dates": {model.datei: np.array([[model.datei, model.datei]])},
            "domain": domain,
            "force_dump": True,
            "sparse_data": False,
            "loadin_perturb_full_vertical": False,
            "surface_level": -1,
            "break_fwd_onlyinit_pipe": False,
            "break_adj_onlyinit_pipe": False
        }

    Warning

    loadin_perturb_full_vertical and surface_level are two parameters that are important to avoid memory overload during the pre-processing of the ensemble members (called samples in CIF). When a large number of samples are transported, all the 4-D fields (e.g. initial conditions) corresponding to the multiple samples are perturbed and loaded into memory. Setting loading_perturb_full_vertical to False prevents all vertical levels in the data from being perturbed. Instead, only the first sample will be fully loaded (not perturbed because it is the first sample). For the remaining samples, only the surface level is loaded and perturbed. The index of the surface level must be provided using the surface_level setting. This workaround only works if the perturbation is not vertically-resolved and it requires the developers to implement a way of retrieving the full perturbation in the native2inputs routines.

2. In the function native2inputs, the restart_inicond component type must be properly handled. Note that the function to write is dependent on the model. Below, a simple example that moves (as a symbolic link) the restart file generated by the previous simulation to the workdir of the current simulation being prepared:

# Initial conditions from a restart file
elif input_type == 'restart_inicond':
    trids = list(datastore.keys())
    dirorig = datastore[trids[0]]["dirorig"]
    fileorig = datastore[trids[0]]["fileorig"]

    source_file = os.path.join(dirorig, fileorig)
    target_file = "{}/restart_atm_DOM01.nc".format(runsubdir)
    path.link(source_file, target_file)

It is also necessary to add these lines of code in the function outputs2native so the path to restart file is properly passed to the next simulation:

# Fetch end concentrations if input_type is "endconcs"
elif input_type == "endconcs":
    dataout = {}
    fileorig = ddi.strftime("{}/../chain/restart_%Y%m%d%H.nc".format(runsubdir))
    for trid in data2dump:
        dataout[trid] = {"fileorig": fileorig}

Defining a correspondence between perturbed outputs and inputs#

For a given configuration to be computed in batch mode, pyCIF perturbs backwards transforms along the transform pipe (see details here) from the observation vector backwards to the control vector and more generally to inputs.

For any parameter in the observation vector, e.g., (‘concs’, ‘CH4’), pyCIF will create virtual equivalents to be computed in parallel: (‘concs’, ‘CH4__sample#000’), (‘concs’, ‘CH4__sample#001’), …, (‘concs’, ‘CH4__sample#00N’). The mapper of every transform in the transform pipe is thus perturbed according to the perturbed inputs of the transform successors. To perturb inputs depending on perturbed outputs, there are two possibilities:

  1. inputs and outputs have the same name, thus perturbations are directly propagated

  2. input names differ from output names In this case, pyCIF needs extra information to know which inputs should propagate the perturbations. This is done by defining the variable outputs2inputs in the transform mapper, defined in the function ini_mapper. This is valid for any transform, and any model.

    The variable outputs2inputs is defined as a dictionary for which each entry is an unperturbed output, e.g., (‘concs’, ‘CH4’), and the values for each entry are the corresponding inputs to be perturbed.

    For instance, CHIMERE runs with the output (‘concs’, ‘CH4’) and the inputs (‘flux’, ‘CH4’), (‘inicond’, ‘CH4’), (‘latcond’, ‘CH4’), (‘topcond’, ‘CH4’), (‘meteo’, ‘’).

    To produce perturbed outputs, the model will need perturbed inputs for (‘flux’, ‘CH4’), (‘inicond’, ‘CH4’), (‘latcond’, ‘CH4’), (‘topcond’, ‘CH4’), but not for (‘meteo’, ‘’). Therefore, the variable outputs2inputs will look like:

    # Save propagation of perturbations
inputs = ["inicond", "latcond", "topcond", "restart_inicond"]
mapper["outputs2inputs"] = {
    (outcomp, s):
    [(cmp, s) for cmp in inputs]
    + ([("flux", s)] if s in model.chemistry.emis_species.attributes
       else [])
    for s in model.chemistry.outspecies.attributes
    for outcomp in model.output_components
}

    Warning

    The variable outputs2inputs can have a complex behaviour in the case of complex chemical schemes, or when models are flexible in terms of needed inputs to transport a given species. Defining this variable thus need thorough thought before starting. Do not hesitate to contact the core team to help you.

Perturbing the behaviour of the model#

To allow for the batch_sampling to work, pyCIF should be able to modify the behaviour of the model. A dedicated function perturb_model should be added to the model plugin to do so.

Usually, the perturb_model function should modify the underlying chemistry plugin to update, e.g., acspecies, emis_species, etc.

One can find an example of such a function here:

pycif.plugins.models.chimere.perturb_model(self, nsamples, transf_mapper)[source]

or here:

pycif.plugins.models.iconart.perturb_model(self, nsamples, transf_mapper)[source]

Improving the native2inputs routines to handle multiple tracers efficiently#

Previously, CIF processed each virtual tracer (e.g., CH4__sample#00i) individually in the native2inputs routines, leading to long pre-processing times. To optimize performance, CIF can now supply data for multiple tracers simultaneously, allowing developers to modify routines for improved efficiency. The following scenarios may arise:

  1. No perturbed components (default behavior, no EnSRF mode).

    In this case, no modifications are needed. Each non-perturbed tracer is processed independently, with only its corresponding data provided.

  2. At least one (component, tracer) pair is perturbed

    An ensemble of virtual tracers is generated for the perturbed tracer. Two sub-cases exist:

    1. The called native2inputs routine belongs to the perturbed component

      A unique perturbation is applied to each sample (e.g., CH4__sample#00i). All sample data is available in the datastore. Developers can process them collectively for efficiency.

    2. The called native2inputs routine does not belong to the perturbed component.

      To avoid redundant computations, only the reference tracer’s data (e.g., CH4) is stored. The developer must duplicate this data for each model-transported tracer.

    Some examples are provided in the native2inputs routines of ICON-ART.

Here is an example of how to detect 1) an ensemble of tracers and 2) if the data being processed is perturbed, at the beginning of an native2inputs routine:

def make_fluxes(self, datastore, ddi, ddf, runsubdir, mode):
    """Native to input function for the ICON-ART flux module"""

    # Inputs and outputs
    trids_in = list(datastore.keys())
    trids_out = [("flux", s) for s in self.chemistry.emis_species]

    # Flags to detect the ensemble
    is_ensemble = False
    is_perturbed_comp = False

    # Loop over the species that must be emitted
    for trid_out in trids_out:

        trid_in = trid_out
        emspec_ref = emspec = trid_out[1]

        # If ensemble, check what the datastore contains
        if "__sample#" in emspec:
            is_ensemble = True
            is_perturbed_comp = True

            emspec_ref = emspec.split("__sample#")[0]
            sample_id = emspec.split("__sample#")[1]
            if int(sample_id) > 0:
                continue

            if trid_in not in datastore:
                trid_in = ("flux", emspec_ref)
                is_perturbed_comp = False

        if trid_in not in datastore:
            continue

        tracer = datastore[trid_in]
        tracer_data = tracer["data"][ddi]

        ...
Contents