Source code for pycif.plugins.models.satwetch4.run

import pandas as pd
import datetime
import os
import shutil
import subprocess
import glob
import numpy as np

from logging import info, debug




def run(self, runsubdir, mode, workdir, ddi, nbproc=1, do_simu=True,
        approx_transf=False, ref_fwd_dir="", overlap=False, datastore=None, **kwargs):
    """Run the model in forward, tangent-linear or adjoint mode.
    This includes:

        - executing the model external executable
        - updating ``adj_refdir``
        - moving files needed for chained simulations to
          "{}/../".format(runsubdir)

    Note:
        For model for which the adjoint is not coded, make sure to return
        a clear error if the ``run`` function is called
        in ``adj`` mode and with ``do_simu = True``

    Args:
        self: the model Plugin
        runsubdir (str): working directory for the current run
        mode (str): forward or backward
        workdir (str): pyCIF working directory
        do_simu (bool): re-run or not existing simulation

    """

    T0 = 273
    ucf = 10**-9  # unit conversion facrot to convert ugCH4 m-2 s-1 to kgCH4 m-2 s-1
    fw = self.meteo_inundated_fraction[ddi]
    Clabile = self.meteo_C_labile[ddi]
    T = self.meteo_soil_temperature[ddi]
    K = self.satwetch4_model_param_k[ddi]  # in ugCH4/m2/s
    Q10value = self.satwetch4_model_param_q10[ddi]

    info("Temperature's time Dimention: {}".format(T.time.values))
    info("Inundation Fraction's time Dimention: {} ".format(fw.time.values))
    info("Clabile's time Dimention: {} ".format(Clabile.time.values))
    info("K's time Dimention: {} ".format(K.time.values))
    info("Q10's time Dimention: {} ".format(Q10value.time.values))

    if mode == 'fwd':
        Q10fun = np.power(Q10value, (T0/T))
        flux = ucf * fw * K * Clabile * np.power(Q10fun, ((T - T0)/10))

        self.flux_out = flux

    elif mode == 'tl':
        xmod_in = datastore["inputs"]
        xmod_out = datastore["outputs"]

        # hard coding only to parameters in CV. Please Suggest me a better way
        # I assume the RHS exists if they're in control vector and mode is tl
        dQ10value = self.satwetch4_model_param_q10_tl[ddi]
        dK = self.satwetch4_model_param_k_tl[ddi]

        Q10fun = np.power(Q10value, (T0 / T))
        flux = ucf * fw * K * Clabile * np.power(Q10fun, ((T - T0)/10))
        dflux = ucf * fw * Clabile * (
            K * (T0 / T * (T - T0) / 10) * np.power(Q10value,
                                                    T0 / T * (T - T0) / 10 - 1) * dQ10value
            + np.power(Q10fun, ((T - T0) / 10)) * dK
        )

        dflux.values = np.where((T == 0) | np.isnan(T), 0, dflux.values)
        flux.values = np.where((T == 0) | np.isnan(T), 0, flux.values)

        self.flux_out = flux
        self.dflux_out = dflux

    else:  # if adj
        xmod_in = datastore["inputs"]
        xmod_out = datastore["outputs"]

        for trid_out in xmod_out:
            # I assume this is adjoint flux i.e d_observable/d_flux or d_costfunction/d_flux
            adjoint_flux = self.flux_adjoint

        Q10fun = np.power(Q10value, (T0 / T))

        adjoint_Q10value = ucf * fw * Clabile * K * \
            (T0 / T) * (T - T0) / 10 * np.power(Q10value,
                                                T0 / T * (T - T0) / 10 - 1) * adjoint_flux
        adjoint_K = ucf * fw * Clabile * \
            np.power(Q10fun, ((T - T0) / 10)) * adjoint_flux

        adjoint_Q10value.values = np.where(
            (T == 0) | np.isnan(T), 0, adjoint_Q10value)
        adjoint_K.values = np.where((T == 0) | np.isnan(T), 0, adjoint_K)

        # Initialize dictionary for later use in native2inputs_adj
        if not hasattr(self, "satwetch4_model_param_k_adjoint"):
            self.satwetch4_model_param_k_adjoint = {}

        if not hasattr(self, "satwetch4_model_param_q10_adjoint"):
            self.satwetch4_model_param_q10_adjoint = {}

        self.satwetch4_model_param_k_adjoint[ddi] = adjoint_K
        self.satwetch4_model_param_q10_adjoint[ddi] = adjoint_Q10value

    pass