diff --git a/NAMESPACE b/NAMESPACE index 9b382b193d37a1bcb71f80d6f759e2d74370be2e..0b485d3f26d1a7df051afc63f8ebce54981d7d8c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -21,8 +21,11 @@ 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) diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R index 8d10cac60d97a60a75c139c17941fb926e6750ab..60f11d7f6e57d729bb51808b59150bac6f66a503 100644 --- a/R/CreateInputsModel.GRiwrm.R +++ b/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)) diff --git a/R/RunModel.GRiwrmOutputsModel.R b/R/RunModel.GRiwrmOutputsModel.R new file mode 100644 index 0000000000000000000000000000000000000000..7b1f82b32506648e3ccc6d36cb42e788cb82fea7 --- /dev/null +++ b/R/RunModel.GRiwrmOutputsModel.R @@ -0,0 +1,118 @@ +#' 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 `OutputsModel`. +#' +#' `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 Calibration +#' @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 +#' +#' @inherit RunModel.GRiwrmInputsModel return +#' @export +#' +RunModel.GRiwrmOutputsModel <- function(OutputsModel, + InputsModel, + RunOptions, + IndPeriod_Run = which(InputsModel[[1]]$DatesR %in% DatesR), + DatesR = getNextTimeSteps(OutputsModel), + Qinf = NULL, + Qrelease = NULL, + Qmin = NULL) { + stopifnot(inherits(OutputsModel, "GRiwrmOutputsModel"), + inherits(InputsModel, "GRiwrmInputsModel"), + inherits(RunOptions, "GRiwrmRunOptions")) + # Check Run Period + next_time_step <- getNextTimeSteps(OutputsModel) + 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(OutputsModel[[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]]$isDiversion) { + 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") { + if (is.null(InputsModel[[id]])) { + # Direct Injection + g <- attr(InputsModel, "GRiwrm") + id_down <- g$down[g$id == id] + InputsModel[[id_down]]$Qupstream[IndPeriod_Run, id] <- v + } else { + if (!InputsModel[[id]]$isDiversion) { + stop("The column ", id, " of the argument `Qinf` does not refer to a DirectInjection or a Diversion node") + } + InputsModel[[id]]$Qdiv[IndPeriod_Run] <- v + } + } + } + } + + # Run the model + return(suppressMessages( + RunModel( + InputsModel, + RunOptions = RunOptions, + Param = extractParam(OutputsModel) + ) + )) +} diff --git a/R/RunModel_Reservoir.R b/R/RunModel_Reservoir.R index e0db27aff4e1ef05484b51d9919c409c612c467d..c290116fdf4bed521bfb7dcef846cd942011ed72 100644 --- a/R/RunModel_Reservoir.R +++ b/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] diff --git a/R/utils.Calibration.R b/R/utils.Calibration.R index 2d7157230e9af87d68edcd5fa8a65658cbad8561..96b05f7dba7b049be054a6ca5266c6a1a68a218a 100644 --- a/R/utils.Calibration.R +++ b/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) +} + diff --git a/R/utils.CreateInputsModel.R b/R/utils.CreateInputsModel.R index c0bd112d45971bf608e94b159bd265befe565bf4..6683fec0501dcb48f3682a25fabe8cbd2ee990e7 100644 --- a/R/utils.CreateInputsModel.R +++ b/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)) + } + } + }) +} diff --git a/R/utils.RunModel.R b/R/utils.RunModel.R index f34a136a1e28662475eb01f2107278edf7ae532c..be5ad9375bd447005f1286e745ce36264028438a 100644 --- a/R/utils.RunModel.R +++ b/R/utils.RunModel.R @@ -61,7 +61,9 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) { #' @noRd #' serializeIniStates <- function(IniStates) { - unlist(IniStates) + IniStates <- unlist(IniStates) + IniStates[is.na(IniStates)] <- 0 + return(IniStates) } diff --git a/man/RunModel.GRiwrmOutputsModel.Rd b/man/RunModel.GRiwrmOutputsModel.Rd new file mode 100644 index 0000000000000000000000000000000000000000..3ff2dc4718679ead54b447afa56341f3e99a461f --- /dev/null +++ b/man/RunModel.GRiwrmOutputsModel.Rd @@ -0,0 +1,72 @@ +% 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}( + OutputsModel, + InputsModel, + RunOptions, + IndPeriod_Run = which(InputsModel[[1]]$DatesR \%in\% DatesR), + DatesR = getNextTimeSteps(OutputsModel), + Qinf = NULL, + Qrelease = NULL, + Qmin = NULL +) +} +\arguments{ +\item{OutputsModel}{Object returned by \link{RunModel.GRiwrmInputsModel}, +\link{RunModel.Supervisor}, or \link{RunModel.GRiwrmOutputsModel}} + +\item{InputsModel}{[object of class \emph{InputsModel} or \emph{GRiwrmInputsModel}] see \link{CreateInputsModel}} + +\item{RunOptions}{[object of class \emph{RunOptions} or \emph{GRiwrmRunOptions}] see \link{CreateRunOptions}} + +\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} +} +\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 \code{OutputsModel}. + +\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. +} diff --git a/man/extractParam.Rd b/man/extractParam.Rd index 702831cc13de6f4fc189db5911f15700f62dda3f..4f410a4f592b1bfd43302ca70466f3263304c509 100644 --- a/man/extractParam.Rd +++ b/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}} diff --git a/tests/testthat/helper_1_RunModel.R b/tests/testthat/helper_1_RunModel.R index c8d7fedbc0bf33119c0820ffdc591a7fa8973cb3..57110eb77e2d6273cdcd813668b05a6d64bb3681 100644 --- a/tests/testthat/helper_1_RunModel.R +++ b/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 diff --git a/tests/testthat/test-RunModel.GRiwrmOutputsModel.R b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R new file mode 100644 index 0000000000000000000000000000000000000000..12b64618a44b433a505afff299c2f7e15728f758 --- /dev/null +++ b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R @@ -0,0 +1,52 @@ +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( + `54095` = rep(-1500000, length(DatesR)), + WD = rep(-250000, length(DatesR)) + ) + names(Qinf)[1] <- "54095" + Qrelease <- data.frame(Dam = rep(100000, length(DatesR))) + Qmin <- data.frame("54095" = rep(1000000, 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") + ) + Qm3s <- attr(OM, "Qm3s") + for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) { + OM <- RunModel(OM, + InputsModel = InputsModel, + RunOptions = RunOptions, + IndPeriod_Run = which(dfTS$yearmonth == ym)) + Qm3s <- rbind(Qm3s, attr(OM, "Qm3s")) + } + expect_equal(nrow(Qm3s), length(DatesR) - 364) +})