Merge branch '173-run-a-model-from-previous-run-and-changing-inputs' into 'dev'

Resolve "Run a model from previous run and changing inputs"

Closes #173

See merge request !105

NAMESPACE

@@ -21,12 +21,17 @@ S3method(CreateRunOptions,InputsModel)
 S3method(CreateRunOptions,character)
 S3method(RunModel,GR)
 S3method(RunModel,GRiwrmInputsModel)
+S3method(RunModel,GRiwrmOutputsModel)
 S3method(RunModel,InputsModel)
 S3method(RunModel,Supervisor)
+S3method(extractParam,GRiwrmOutputsCalib)
+S3method(extractParam,GRiwrmOutputsModel)
 S3method(isNodeDownstream,GRiwrm)
 S3method(isNodeDownstream,GRiwrmInputsModel)
 S3method(isNodeUpstream,GRiwrm)
 S3method(isNodeUpstream,GRiwrmInputsModel)
+S3method(merge,GRiwrmOutputsModel)
+S3method(merge,OutputsModel)
 S3method(plot,GRiwrm)
 S3method(plot,GRiwrmOutputsModel)
 S3method(plot,OutputsModelReservoir)
@@ -47,6 +52,7 @@ export(RunModel_Reservoir)
 export(as.Qm3s)
 export(extractParam)
 export(getAllNodesProperties)
+export(getNextTimeSteps)
 export(getNoSD_Ids)
 export(getNodeProperties)
 export(getNodeRanking)

R/CreateInputsModel.GRiwrm.R

@@ -116,57 +116,21 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
     warning("The usage of 'Qobs' is deprecated, use 'Qinf' instead")
     Qinf <- Qobs
   }
-  varNames <- c("Precip", "PotEvap", "TempMean", "Qinf", "Qmin",
-                "TempMin", "TempMax", "ZInputs", "HypsoData", "NLayers")
-  names(varNames) <- varNames
-  lapply(varNames, function(varName) {
-    v <- get(varName)
-    if (!is.null(v)) {
-      if (is.matrix(v) || is.data.frame(v)) {
-        if (is.null(colnames(v))) {
-          stop(sprintf(
-            "'%s' must have column names",
-            varName
-          ))
-        } else if (!all(colnames(v) %in% x$id)) {
-          stop(sprintf(
-            "'%s' column names must be included in 'id's of the GRiwrm object",
-            varName
-          ), "\n",
-          sprintf("These columns are not known: %s",
-                  paste(colnames(v)[!colnames(v) %in% x$id], collapse = ", ")))
-        } else if (any(duplicated(colnames(v)))) {
-          stop(sprintf(
-            "'%s' has duplicated column names: '%s'",
-            varName,
-            paste(colnames(v)[duplicated(colnames(v))], collapse = "', '")
-          ))
-        }
-        if (!varName %in% c("ZInputs", "NLayers", "HypsoData") && nrow(v) != length(DatesR)) {
-          stop(sprintf(
-            "'%s' number of rows and the length of 'DatesR' must be equal",
-             varName
-          ))
-        }
-        if (varName %in% c("Precip", "PotEvap", "Qmin")) {
-          if (any(is.na(v))) {
-            stop(sprintf(
-              "`NA` values detected in '%s'. Missing values are not allowed in InputsModel",
-              varName
-            ))
-          }
-          if (any(v < 0)) {
-            stop(sprintf(
-              "'%s' values must be positive or nul. Missing values are not allowed in InputsModel",
-              varName
-            ))
-          }
-        }
-      } else if (!varName %in% c("ZInputs", "NLayers")) {
-        stop(sprintf("'%s' must be a matrix or a data.frame", varName))
-      }
-    }
-  })
+
+  checkInputsModelArguments(
+    x,
+    DatesR,
+    Precip = Precip,
+    PotEvap = PotEvap,
+    TempMean = TempMean,
+    Qinf = Qinf,
+    Qmin = Qmin,
+    TempMin = TempMin,
+    TempMax = TempMax,
+    ZInputs = ZInputs,
+    HypsoData = HypsoData,
+    NLayers = NLayers
+  )
 
   if (is.null(Qinf)) Qinf <- matrix(0, ncol = 0, nrow = length(DatesR))
   if (is.null(Qrelease)) Qrelease <- matrix(0, ncol = 0, nrow = length(DatesR))

R/CreateSupervisor.R

@@ -4,7 +4,7 @@
 #' See [RunModel.Supervisor] and vignettes for examples of use.
 #'
 #' @param InputsModel \[object of type `GRiwrmInputsModel`\] inputs of the model
-#' @param TimeStep [numeric] number of time steps between each supervision
+#' @param TimeStep [integer] number of time steps between each supervision
 #'
 #' @return A `Supervisor` object which is an [environment] containing all the necessary variables to run a supervised simulation, such as:
 #' - `DatesR` [POSIXct]: vector of date from `InputsModel`

R/RunModel.GRiwrmOutputsModel.R

+#' RunModel for hot restart after a previous simulation period
+#'
+#' This function allows to restart a simulation at the end of a previous
+#' simulation period. Parameters `Qinf`, `Qrelease`, and `Qmin` can be
+#' redefined for this new simulation period.
+#'
+#' @details
+#' `IndPeriod_Run` or `DatesR` must must be continuous periods starting the
+#' time step after the last simulation time step of the `GRiwrmOutputsModel`
+#' object provided through the argument `x`.
+#'
+#' `Qinf`, `Qmin`, and `Qrelease` are used for overwriting the corresponding
+#' arguments provided to [CreateInputsModel.GRiwrm] for the period to be simulated.
+#' Therefore, the number of rows of these arguments must correspond to
+#' `IndPeriod_Run` or `DatesR` lengths.
+#'
+#' @inheritParams getNextTimeSteps
+#' @inheritParams RunModel.GRiwrmInputsModel
+#' @inheritParams airGR::CreateRunOptions
+#' @param DatesR (optional) [POSIXt] vector of dates of period to be used for
+#' the model run. See details
+#' @param Qinf (optional) [matrix] or [data.frame] of [numeric] containing
+#'        observed flows. It must be provided only for nodes of type "Direct
+#'        injection" and "Diversion" \[m3 per time step\].
+#'        Column names correspond to node IDs. Negative flows are abstracted from
+#'        the model and positive flows are injected to the model. See details
+#' @param Qmin (optional) [matrix] or [data.frame] of [numeric] containing
+#'        minimum flows to let downstream of a node with a Diversion \[m3 per
+#'        time step\]. Default is zero. Column names correspond to node IDs.
+#'        See details
+#' @param Qrelease (optional) [matrix] or [data.frame] of [numeric] containing
+#'        release flows by nodes using the model `RunModel_Reservoir` \[m3 per
+#'        time step\]. See details
+#' @param merge_outputs [logical] Merge simulation outputs with the one provided
+#' in argument `x`
+#' @param ... Further arguments for compatibility with S3 methods
+#'
+#' @inherit RunModel.GRiwrmInputsModel return
+#' @export
+#'
+RunModel.GRiwrmOutputsModel <- function(x,
+                                        InputsModel,
+                                        RunOptions,
+                                        IndPeriod_Run = which(InputsModel[[1]]$DatesR %in% DatesR),
+                                        DatesR = getNextTimeSteps(x),
+                                        Qinf = NULL,
+                                        Qrelease = NULL,
+                                        Qmin = NULL,
+                                        merge_outputs = TRUE,
+                                        ...) {
+  stopifnot(inherits(x, "GRiwrmOutputsModel"),
+            inherits(InputsModel, "GRiwrmInputsModel"),
+            inherits(RunOptions, "GRiwrmRunOptions"))
+  # Check Run Period
+  next_time_step <- getNextTimeSteps(x)
+  next_index <- which(InputsModel[[1]]$DatesR == next_time_step)
+  if (IndPeriod_Run[1] != next_index) {
+    stop("`IndPeriod_Run` should have its first element equal to ", next_index)
+  }
+
+  # State Initiation
+  for (id in names(RunOptions)) {
+    # Run model for the sub-basin and one time step
+    RunOptions[[id]]$IniResLevels <- NULL
+    RunOptions[[id]]$IniStates <- serializeIniStates(x[[id]]$StateEnd)
+    RunOptions[[id]]$IndPeriod_WarmUp <- 0L
+    RunOptions[[id]]$IndPeriod_Run <- IndPeriod_Run
+  }
+
+  # Inputs change
+  checkInputsModelArguments(
+    attr(InputsModel, "GRiwrm"),
+    InputsModel[[1]]$DatesR[IndPeriod_Run],
+    Qinf = Qinf,
+    Qrelease = Qrelease,
+    Qmin = Qmin
+  )
+  inputs <- list(Qinf = Qinf, Qrelease = Qrelease, Qmin = Qmin)
+  inputs[sapply(inputs, is.null)] <- NULL
+  for (inputArg in names(inputs)) {
+    input <- inputs[[inputArg]]
+    if (length(IndPeriod_Run) != nrow(input)) {
+      stop("The Argument ", inputArg,
+           " must have a number of rows identical to the lenght of `IndPeriod_Run`")
+    }
+    for (id in colnames(input)) {
+      v <- input[, id, drop = TRUE]
+      if (inputArg %in% c("Qrelease", "Qmin")) {
+        if (inputArg == "Qrelease" && !InputsModel[[id]]$isReservoir) {
+            stop("The column ", id, " of the argument `Qrelease` does not refer to a Reservoir node")
+        }
+        if (inputArg == "Qmin" && !InputsModel[[id]]$hasDiversion) {
+          stop("The column ", id, " of the argument `Qmin` does not refer to a Diversion node")
+        }
+        if (is.null(InputsModel[[id]][[inputArg]])) {
+            stop("InputsModel[['", id, "']] should contain a `", inputArg, "` item")
+        }
+        InputsModel[[id]][[inputArg]][IndPeriod_Run] <- v
+      }
+      if (inputArg == "Qinf") {
+        g <- attr(InputsModel, "GRiwrm")
+        if (is.null(InputsModel[[id]])) {
+          # Direct Injection
+          id_down <- g$down[g$id == id]
+          InputsModel[[id_down]]$Qupstream[IndPeriod_Run, id] <- v
+        } else {
+          if (!InputsModel[[id]]$hasDiversion) {
+            stop("The column ", id, " of the argument `Qinf` does not refer to a DirectInjection or a Diversion node")
+          }
+          # Update withdrawal due to Diversion
+          InputsModel[[id]]$Qdiv[IndPeriod_Run] <- -v
+        }
+      }
+    }
+  }
+
+  # Run the model
+  OM <- suppressMessages(
+    RunModel(
+      InputsModel,
+      RunOptions = RunOptions,
+      Param = extractParam(x)
+    )
+  )
+
+  if (merge_outputs) {
+    OM <- merge(x, OM)
+  }
+
+  return(OM)
+}

R/RunModel.InputsModel.R

@@ -43,6 +43,9 @@ RunModel.InputsModel <- function(x = NULL,
     }
   }
 
+  # Avoiding Error in `FUN_MOD(x, RunOptions, Param)`: NA/NaN/Inf in foreign function call (arg 7)
+  RunOptions$IniStates[is.na(RunOptions$IniStates)] <- 0
+
   FUN_MOD <- match.fun(FUN_MOD)
   if (identical(FUN_MOD, RunModel_Lag)) {
     OutputsModel <- RunModel_Routing(x, RunOptions, Param)

R/RunModel_Reservoir.R

@@ -146,6 +146,7 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
       OutputsModel$RunOptions$WarmUpQdiv_m3 <- Qdiv_m3[iWarmUp]
     }
   }
+  OutputsModel$RunOptions$Param <- Param
   iRun <- length(IndPerWarmUp) + seq(length(RunOptions$IndPeriod_Run))
   OutputsModel$Qsim_m3 <- Qsim_m3[iRun]
   OutputsModel$Vsim <- Vsim[iRun]

R/utils.Calibration.R

@@ -278,6 +278,18 @@ transferGRparams <- function(InputsModel, Param, donor, receiver, default_param
 #' @export
 #'
 extractParam <- function(x) {
-  stopifnot(inherits(x, "GRiwrmOutputsCalib"))
+  UseMethod("extractParam")
+}
+
+#' @export
+#' @rdname extractParam
+extractParam.GRiwrmOutputsCalib <- function(x) {
   lapply(x, "[[", "ParamFinalR")
 }
+
+#' @export
+#' @rdname extractParam
+extractParam.GRiwrmOutputsModel <- function(x) {
+  lapply(x, function(o) o$RunOptions$Param)
+}
+

R/utils.CreateInputsModel.R

@@ -57,3 +57,65 @@ checkQinfQrelease <- function(g, varname, Q) {
   }
   return(Q)
 }
+
+#' Check the parameters provided to CreateInputsModel.GRiwrm
+#'
+#' @param x GRiwrm
+#' @param DatesR DatesR
+#' @param ... Parameters to check
+#'
+#' @return Nothing
+#' @noRd
+#'
+checkInputsModelArguments <- function(x, DatesR, ...) {
+  dots <- list(...)
+
+  lapply(names(dots), function(varName) {
+    v <- dots[[varName]]
+    if (!is.null(v)) {
+      if (is.matrix(v) || is.data.frame(v)) {
+        if (is.null(colnames(v))) {
+          stop(sprintf(
+            "'%s' must have column names",
+            varName
+          ))
+        } else if (!all(colnames(v) %in% x$id)) {
+          stop(sprintf(
+            "'%s' column names must be included in 'id's of the GRiwrm object",
+            varName
+          ), "\n",
+          sprintf("These columns are not known: %s",
+                  paste(colnames(v)[!colnames(v) %in% x$id], collapse = ", ")))
+        } else if (any(duplicated(colnames(v)))) {
+          stop(sprintf(
+            "'%s' has duplicated column names: '%s'",
+            varName,
+            paste(colnames(v)[duplicated(colnames(v))], collapse = "', '")
+          ))
+        }
+        if (!varName %in% c("ZInputs", "NLayers", "HypsoData") && nrow(v) != length(DatesR)) {
+          stop(sprintf(
+            "'%s' number of rows and the length of 'DatesR' must be equal",
+            varName
+          ))
+        }
+        if (varName %in% c("Precip", "PotEvap", "Qmin")) {
+          if (any(is.na(v))) {
+            stop(sprintf(
+              "`NA` values detected in '%s'. Missing values are not allowed in InputsModel",
+              varName
+            ))
+          }
+          if (any(v < 0)) {
+            stop(sprintf(
+              "'%s' values must be positive or nul. Missing values are not allowed in InputsModel",
+              varName
+            ))
+          }
+        }
+      } else if (!varName %in% c("ZInputs", "NLayers")) {
+        stop(sprintf("'%s' must be a matrix or a data.frame", varName))
+      }
+    }
+  })
+}

R/utils.RunModel.R

@@ -61,7 +61,8 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
 #' @noRd
 #'
 serializeIniStates <- function(IniStates) {
-  unlist(IniStates)
+  IniStates <- unlist(IniStates)
+  return(IniStates)
 }
 
 
@@ -88,3 +89,56 @@ calcOverAbstraction <- function(O, WarmUp) {
   }
   return(O)
 }
+
+
+#' Get the next time steps date/time of a simulation
+#'
+#' @param x Object returned by [RunModel.GRiwrmInputsModel],
+#' [RunModel.Supervisor], or [RunModel.GRiwrmOutputsModel]
+#' @param TimeStep [integer] number of time steps to get after the end of the
+#' simulation
+#'
+#' @return A [POSIXct] containing the date/time of the time steps following
+#' the end of the simulation.
+#' @export
+#'
+getNextTimeSteps <- function(x, TimeStep = 1L) {
+  stopifnot(inherits(x, "GRiwrmOutputsModel"),
+            is.integer(TimeStep))
+  last_date <- dplyr::last(x[[1]]$DatesR)
+  first_date <- last_date + attr(x, "TimeStep")
+  return(seq(first_date, length.out = TimeStep, by = attr(x, "TimeStep")))
+}
+
+
+#' Merge Two outputs of airGR simulations
+#'
+#' @param x,y **OutputsModel** objects from [airGR::RunModel],
+#' [RunModel.GRiwrmInputsModel], [RunModel.GRiwrmOutputsModel], [RunModel.Supervisor]
+#' @param ... Not used
+#'
+#' @return An object **OutputsModel** with merged times series of simulation
+#' results.
+#' @rdname merge.OutputsModel
+#' @export
+#'
+merge.OutputsModel <- function(x, y, ...) {
+  items <- names(x)
+  items <- items[!grepl("RunOptions|StateEnd", items)]
+  for (item in items) {
+    y[[item]] <- c(x[[item]], y[[item]])
+  }
+  return(y)
+}
+
+#' @rdname merge.OutputsModel
+#' @export
+merge.GRiwrmOutputsModel <- function(x, y, ...) {
+  y_attributes <- attributes(y)
+  y <- lapply(setNames(nm = names(y)), function(id) {
+    merge(x[[id]], y[[id]])
+  })
+  attributes(y) <- y_attributes
+  attr(y, "Qm3s") <- rbind(attr(x, "Qm3s"), attr(y, "Qm3s"))
+  return(y)
+}

man/CreateSupervisor.Rd

@@ -9,7 +9,7 @@ CreateSupervisor(InputsModel, TimeStep = 1L)
 \arguments{
 \item{InputsModel}{[object of type \code{GRiwrmInputsModel}] inputs of the model}
 
-\item{TimeStep}{\link{numeric} number of time steps between each supervision}
+\item{TimeStep}{\link{integer} number of time steps between each supervision}
 }
 \value{
 A \code{Supervisor} object which is an \link{environment} containing all the necessary variables to run a supervised simulation, such as:

man/RunModel.GRiwrmOutputsModel.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/RunModel.GRiwrmOutputsModel.R
+\name{RunModel.GRiwrmOutputsModel}
+\alias{RunModel.GRiwrmOutputsModel}
+\title{RunModel for hot restart after a previous simulation period}
+\usage{
+\method{RunModel}{GRiwrmOutputsModel}(
+  x,
+  InputsModel,
+  RunOptions,
+  IndPeriod_Run = which(InputsModel[[1]]$DatesR \%in\% DatesR),
+  DatesR = getNextTimeSteps(x),
+  Qinf = NULL,
+  Qrelease = NULL,
+  Qmin = NULL,
+  merge_outputs = TRUE,
+  ...
+)
+}
+\arguments{
+\item{x}{Object returned by \link{RunModel.GRiwrmInputsModel},
+\link{RunModel.Supervisor}, or \link{RunModel.GRiwrmOutputsModel}}
+
+\item{InputsModel}{[object of class \emph{InputsModel}] see \code{\link[airGR]{CreateInputsModel}} for details}
+
+\item{RunOptions}{[object of class \emph{GRiwrmRunOptions}] see \link{CreateRunOptions.GRiwrmInputsModel} for details}
+
+\item{IndPeriod_Run}{[numeric] index of period to be used for the model run [-]. See details}
+
+\item{DatesR}{(optional) \link{POSIXt} vector of dates of period to be used for
+the model run. See details}
+
+\item{Qinf}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+observed flows. It must be provided only for nodes of type "Direct
+injection" and "Diversion" [m3 per time step].
+Column names correspond to node IDs. Negative flows are abstracted from
+the model and positive flows are injected to the model. See details}
+
+\item{Qrelease}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+release flows by nodes using the model \code{RunModel_Reservoir} [m3 per
+time step]. See details}
+
+\item{Qmin}{(optional) \link{matrix} or \link{data.frame} of \link{numeric} containing
+minimum flows to let downstream of a node with a Diversion [m3 per
+time step]. Default is zero. Column names correspond to node IDs.
+See details}
+
+\item{merge_outputs}{\link{logical} Merge simulation outputs with the one provided
+in argument \code{x}}
+
+\item{...}{Further arguments for compatibility with S3 methods}
+}
+\value{
+An object of class \emph{GRiwrmOutputsModel}.
+This object is a \link{list} of \emph{OutputsModel} objects produced by \link{RunModel.InputsModel}
+for each node of the semi-distributed model.
+
+It also contains the following attributes (see \link{attr}):
+\itemize{
+\item "Qm3s": a \link{data.frame} containing the dates of simulation and one column by node
+with the simulated flows in cubic meters per seconds (See \link{plot.Qm3s})
+\item "GRiwrm":  a copy of the \emph{GRiwrm} object produced by \link{CreateGRiwrm} and used for the simulation
+\item "TimeStep":  time step of the simulation in seconds
+}
+}
+\description{
+This function allows to restart a simulation at the end of a previous
+simulation period. Parameters \code{Qinf}, \code{Qrelease}, and \code{Qmin} can be
+redefined for this new simulation period.
+}
+\details{
+\code{IndPeriod_Run} or \code{DatesR} must must be continuous periods starting the
+time step after the last simulation time step of the \code{GRiwrmOutputsModel}
+object provided through the argument \code{x}.
+
+\code{Qinf}, \code{Qmin}, and \code{Qrelease} are used for overwriting the corresponding
+arguments provided to \link{CreateInputsModel.GRiwrm} for the period to be simulated.
+Therefore, the number of rows of these arguments must correspond to
+\code{IndPeriod_Run} or \code{DatesR} lengths.
+}

man/extractParam.Rd

@@ -2,9 +2,15 @@
 % Please edit documentation in R/utils.Calibration.R
 \name{extractParam}
 \alias{extractParam}
+\alias{extractParam.GRiwrmOutputsCalib}
+\alias{extractParam.GRiwrmOutputsModel}
 \title{Extract calibrated parameters}
 \usage{
 extractParam(x)
+
+\method{extractParam}{GRiwrmOutputsCalib}(x)
+
+\method{extractParam}{GRiwrmOutputsModel}(x)
 }
 \arguments{
 \item{x}{A \emph{GRiwrmOutputsModel} object returned by \link{Calibration.GRiwrmInputsModel}}

man/getNextTimeSteps.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.RunModel.R
+\name{getNextTimeSteps}
+\alias{getNextTimeSteps}
+\title{Get the next time steps date/time of a simulation}
+\usage{
+getNextTimeSteps(x, TimeStep = 1L)
+}
+\arguments{
+\item{x}{Object returned by \link{RunModel.GRiwrmInputsModel},
+\link{RunModel.Supervisor}, or \link{RunModel.GRiwrmOutputsModel}}
+
+\item{TimeStep}{\link{integer} number of time steps to get after the end of the
+simulation}
+}
+\value{
+A \link{POSIXct} containing the date/time of the time steps following
+the end of the simulation.
+}
+\description{
+Get the next time steps date/time of a simulation
+}

man/merge.OutputsModel.Rd

+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/utils.RunModel.R
+\name{merge.OutputsModel}
+\alias{merge.OutputsModel}
+\alias{merge.GRiwrmOutputsModel}
+\title{Merge Two outputs of airGR simulations}
+\usage{
+\method{merge}{OutputsModel}(x, y, ...)
+
+\method{merge}{GRiwrmOutputsModel}(x, y, ...)
+}
+\arguments{
+\item{x, y}{\strong{OutputsModel} objects from \link[airGR:RunModel]{airGR::RunModel},
+\link{RunModel.GRiwrmInputsModel}, \link{RunModel.GRiwrmOutputsModel}, \link{RunModel.Supervisor}}
+
+\item{...}{Not used}
+}
+\value{
+An object \strong{OutputsModel} with merged times series of simulation
+results.
+}
+\description{
+Merge Two outputs of airGR simulations
+}

tests/testthat/helper_1_RunModel.R

@@ -17,6 +17,7 @@ setupRunModel <-
            griwrm = NULL,
            Qinf = NULL,
            Qrelease = NULL,
+           Qmin = NULL,
            IsHyst = FALSE) {
 
     data(Severn)
@@ -100,8 +101,9 @@ setupRunModel <-
     InputsModel <-
       suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
                                          TempMean = TempMean,
-                                         Qobs = Qinf,
+                                         Qinf = Qinf,
                                          Qrelease = Qrelease,
+                                         Qmin = Qmin,
                                          IsHyst = IsHyst))
 
     # RunOptions

tests/testthat/test-RunModel.GRiwrmOutputsModel.R

+skip_on_cran()
+test_that("RunModel.GRiwrmOutputsModel works", {
+  # Setup model
+  griwrm <- CreateGRiwrm(rbind(
+    n_derived_rsrvr,
+    data.frame(
+      id = "WD",
+      down = "Dam",
+      length = 0,
+      area = NA,
+      model = NA
+    )
+  ))
+  data(Severn)
+  DatesR <-  Severn$BasinsObs[[1]]$DatesR
+  Qinf <- data.frame(
+    # Diversion to the dam
+    `54095` = rep(-1E6, length(DatesR)),
+    # Withdrawal in the dam
+    WD = rep(-250000, length(DatesR))
+  )
+  names(Qinf)[1] <- "54095"
+  # Release of the dam back to the river
+  Qrelease <- data.frame(Dam = rep(100E3, length(DatesR)))
+  # Diversion limited by fixed minimum flow to let in the river
+  Qmin <- data.frame("54095" = rep(3E6, length(DatesR)))
+  names(Qmin) <- "54095"
+  e <- setupRunModel(
+    griwrm = griwrm,
+    Qinf = Qinf,
+    Qrelease = Qrelease,
+    Qmin = Qmin,
+    runRunOptions = FALSE
+  )
+  for (x in ls(e)) assign(x, get(x, e))
+
+  # Set up initial conditions
+  RunOptions <- CreateRunOptions(InputsModel, IndPeriod_WarmUp = 1:364, IndPeriod_Run = 365L)
+  Param <- c(ParamMichel[names(ParamMichel) %in% griwrm$id], list(Dam = c(100E6, 1)))
+  OM <- RunModel(InputsModel, RunOptions, Param)
+
+  # Loop over periods months periods
+  dfTS <- data.frame(
+    DatesR = DatesR,
+    yearmonth = format(DatesR, "%Y-%m")
+  )
+
+  for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
+
+    # Preparing extract of Qinf for the current run
+    ym_IndPeriod_Run <- which(dfTS$yearmonth == ym)
+    ym_Qinf <- Qinf[ym_IndPeriod_Run, , drop = FALSE]
+    ym_Qrelease <- Qrelease[ym_IndPeriod_Run, , drop = FALSE]
+
+    # 50% Restriction on reservoir withdrawals if remaining less than 90 days of water
+    nb_remain_days <- OM$Dam$StateEnd$Reservoir$V / (-ym_Qinf$`WD`[1] + ym_Qrelease$Dam[1])
+    if (nb_remain_days < 180) {
+      ym_Qinf$`WD` <- -(max(0, OM$Dam$StateEnd$Reservoir$V - sum(ym_Qrelease$Dam))) / 365
+    }
+    OM <- RunModel(OM,
+                   InputsModel = InputsModel,
+                   RunOptions = RunOptions,
+                   IndPeriod_Run = ym_IndPeriod_Run,
+                   Qinf = ym_Qinf)
+  }
+
+  expect_equal(nrow(attr(OM, "Qm3s")), length(DatesR) - 364)
+  expect_equal(length(OM[[1]]$DatesR), length(DatesR) - 364)
+})

vignettes/V06_Modelling_regulated_diversion.Rmd

@@ -195,7 +195,7 @@ with and without low-flow support at station 54001:
 dfQdiv <- data.frame(DatesR = OM_div[[1]]$DatesR,
                      Diverted_flow = OM_div$`54001`$Qdiv_m3 / 86400)
 
-oldpar <- par(mfrow=c(3,1), mar = c(2.5,4,1,1))
+oldpar <- par(mfrow = c(3,1), mar = c(2.5,4,1,1))
 plot.Qm3s(dfQdiv)
 
 # Plot natural and influenced flow at station "54001" and "54029"