diff --git a/R/CreateSupervisor.R b/R/CreateSupervisor.R
index 660f4ecba09c9e98a2ebf4511d81e52822a90ab0..07ab00bb7e0aeb07c75a5853887f0badf29dd793 100644
--- a/R/CreateSupervisor.R
+++ b/R/CreateSupervisor.R
@@ -24,6 +24,7 @@
#' InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
#' sv <- CreateSupervisor(InputsModel)
CreateSupervisor <- function(InputsModel) {
+ if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmInputsModel' (See ?CreateInputsModel.GRiwrm)")
# Create Supervisor environment in the parent of GlobalEnv
e <- new.env(parent = parent.env(globalenv()))
class(e) <- c("Supervisor", class(e))
diff --git a/R/RunModel.GR.R b/R/RunModel.GR.R
index e4fa5562d1a4cfc884e85b2f5f0728e5d3d73470..da2bc072a62b7dc6f720c4df406734898db192fa 100644
--- a/R/RunModel.GR.R
+++ b/R/RunModel.GR.R
@@ -7,7 +7,6 @@
#' @export
#'
RunModel.GR <- function(x, RunOptions, Param, ...) {
- message("RunModel.GR")
if (inherits(x, "SD")) {
# Lag model take one parameter at the beginning of the vector
diff --git a/R/RunModel.GRiwrmInputsModel.R b/R/RunModel.GRiwrmInputsModel.R
index a3d894d28db021f1155a5c07225fbc2cd215b06a..ff6047f1733a38f1a2c35afcf3a5444e9b87d29d 100644
--- a/R/RunModel.GRiwrmInputsModel.R
+++ b/R/RunModel.GRiwrmInputsModel.R
@@ -8,32 +8,25 @@
#' @return \emph{GRiwrmOutputsModel} object which is a list of \emph{OutputsModel} objects (See \code{\link[airGR]{RunModel}}) for each node of the semi-distributed model.
#' @export
RunModel.GRiwrmInputsModel <- function(x, RunOptions, Param, ...) {
- message("RunModel.GRiwrmInputsModel")
- # Run runoff model for each sub-basin
- OutputsModel <- lapply(X = x, FUN = function(IM) {
- RunModel.GR(IM,
- RunOptions = RunOptions[[IM$id]],
- Param = Param[[IM$id]])
- })
- class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+ checkRunModelParameters(x, RunOptions, Param)
- # Loop over sub-basin using SD model
- for(id in getSD_Ids(x)) {
- IM <- x[[id]]
- message("RunModel.GRiwrmInputsModel: Treating sub-basin ", id, "...")
+ OutputsModel <- list()
+ class(OutputsModel) <- c("GRiwrmOutputsModel", class(OutputsModel))
+
+ for(IM in x) {
+ message("RunModel.GRiwrmInputsModel: Treating sub-basin ", IM$id, "...")
# Update x$Qupstream with simulated upstream flows
if(any(IM$UpstreamIsRunoff)) {
- IM <- UpdateQsimUpstream(IM, RunOptions[[id]]$IndPeriod_Run, OutputsModel)
+ IM <- UpdateQsimUpstream(IM, RunOptions[[IM$id]]$IndPeriod_Run, OutputsModel)
}
- # Run the SD model for the sub-basin
- OutputsModel[[id]] <- RunModel.SD(
+ # Run the model for the sub-basin
+ OutputsModel[[IM$id]] <- RunModel.InputsModel(
IM,
- RunOptions = RunOptions[[id]],
- Param = Param[[id]],
- OutputsModel[[id]]
+ RunOptions = RunOptions[[IM$id]],
+ Param = Param[[IM$id]]
)
}
diff --git a/R/RunModel.R b/R/RunModel.R
index a3ea4600cbb87eb55c6b3d3e05dcd367797e6a0b..0ac7f23a9c50d35ca70a0877a1f8f1d3693f3727 100644
--- a/R/RunModel.R
+++ b/R/RunModel.R
@@ -6,6 +6,5 @@
#' @return Either a [list] of OutputsModel object (for GRiwrmInputsModel) or an OutputsModel object (for InputsModel)
#' @export
RunModel <- function(x, ...) {
- message("RunModel")
UseMethod("RunModel", x)
}
diff --git a/R/RunModel.SD.R b/R/RunModel.SD.R
index 95f661f1b32aa1882b299d9d4f08ddbe55c14292..dfff38feefe01349a170c8c2ca367058ae2bb382 100644
--- a/R/RunModel.SD.R
+++ b/R/RunModel.SD.R
@@ -1,15 +1,14 @@
#' Run SD Model from run-off model outputs
#'
#' @inheritParams airGR::RunModel_Lag
-#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel]
-#' @param OutputsModel `OutputsModel` object returned by a GR model by [airGR::RunModel]
+#' @param x `InputsModel` used as `InputsModel` parameter for [airGR::RunModel_Lag]
+#' @param QSimDown a [numeric] corresponding to the runoff of the sub-basin (Typically the `Qsim` outputs of the GR model)
#' @param ... further arguments passed to or from other methods
#'
#' @return `OutputsModel` object. See [airGR::RunModel_Lag]
#' @export
#'
-RunModel.SD <- function(x, RunOptions, Param, OutputsModel, ...) {
- message("RunModel.SD")
- x$OutputsModel <- OutputsModel
+RunModel.SD <- function(x, RunOptions, Param, QsimDown, ...) {
+ x$OutputsModel <- list(Qsim = QsimDown)
RunModel_Lag(x, RunOptions = RunOptions, Param = Param[1])
}
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index f847ecca30a3cf0019d11973b6d98f00851e1603..9f35149d53bc13077d5731405d97d0a61fa2e195 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -14,32 +14,53 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
})
# Run runoff model for each sub-basin
- OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
+ x$OutputsModel <- lapply(X = x$InputsModel, FUN = function(IM) {
RunModel.GR(IM,
RunOptions = RunOptions[[IM$id]],
Param = Param[[IM$id]])
})
- class(OutputsModel) <- append(class(OutputsModel), "GRiwrmOutputsModel")
+ class(x$OutputsModel) <- append(class(x$OutputsModel), "GRiwrmOutputsModel")
+ # Save Qsim for step by step simulation
+ Qsim <- lapply(x$OutputsModel, function(OM) {
+ OM$Qsim
+ })
- # Loop over time steps
- # Loop over sub-basin using SD model
+ # Adapt RunOptions to step by step simulation
for(id in getSD_Ids(x$InputsModel)) {
- IM <- x$InputsModel[[id]]
- message("RunModel.GRiwrmInputsModel: Treating sub-basin ", id, "...")
+ RunOptions[[id]]$IndPeriod_WarmUp <- 0L
+ RunOptions[[id]]$Outputs_Sim <- "StateEnd"
+ }
- # Update InputsModel$Qupstream with simulated upstream flows
- if(any(IM$UpstreamIsRunoff)) {
- IM <- UpdateQsimUpstream(IM, RunOptions[[id]]$IndPeriod_Run, OutputsModel)
- }
+ # Loop over time steps
+ for(iTS in RunOptions[[1]]$IndPeriod_Run) {
+ # Run regulation on the whole basin for the current time step
+ x$ts.index <- iTS
+ x$ts.date <- x$InputsModel[[1]]$DatesR[iTS]
+ doSupervision(x)
- # Run the SD model for the sub-basin
- OutputsModel[[id]] <- RunModel.SD(
- IM,
- RunOptions = RunOptions[[id]],
- Param = Param[[id]],
- OutputsModel[[id]]
- )
+ # Loop over sub-basin using SD model
+ for(id in getSD_Ids(x$InputsModel)) {
+ # Update InputsModel$Qupstream with simulated upstream flows
+ for(i in which(x$InputsModel[[id]]$UpstreamIsRunoff)) {
+ x$InputsModel[[id]]$Qupstream[iTS, i] <-
+ x$OutputsModel[[x$InputsModel[[id]]$UpstreamNodes[i]]]$Qsim[iTS - x$ts.index0]
+ }
+
+ # Run the SD model for the sub-basin and one time step
+ RunOptions[[id]]$IndPeriod_Run <- iTS
+ RunOptions[[id]]$IniStates <- unlist(x$OutputsModel[[id]]$StateEnd)
+ x$OutputsModel[[id]] <- RunModel.SD(
+ x$InputsModel[[id]],
+ RunOptions = RunOptions[[id]],
+ Param = Param[[id]],
+ QsimDown = Qsim[[id]][iTS - x$ts.index0]
+ )
+ Qsim[[id]][iTS - x$ts.index0] <- x$OutputsModel[[id]]$Qsim
+ }
+ }
+ for(id in getSD_Ids(x$InputsModel)) {
+ x$OutputsModel[[id]]$Qsim <- Qsim[[id]]
}
- return(OutputsModel)
+ return(x$OutputsModel)
}
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
new file mode 100644
index 0000000000000000000000000000000000000000..e3bb075b1fa903c07c3df5002c6d35fce018ab52
--- /dev/null
+++ b/R/RunModel_Lag.R
@@ -0,0 +1,108 @@
+RunModel_Lag <- function(InputsModel, RunOptions, Param) {
+ NParam <- 1
+
+ ##Arguments_check
+ if (!inherits(InputsModel, "InputsModel")) {
+ stop("'InputsModel' must be of class 'InputsModel'")
+ }
+ if (!inherits(InputsModel, "SD")) {
+ stop("'InputsModel' must be of class 'SD'")
+ }
+ if (!inherits(RunOptions, "RunOptions")) {
+ stop("'RunOptions' must be of class 'RunOptions'")
+ }
+ if (!is.vector(Param) | !is.numeric(Param)) {
+ stop("'Param' must be a numeric vector")
+ }
+ if (sum(!is.na(Param)) != NParam) {
+ stop(paste("'Param' must be a vector of length", NParam, "and contain no NA"))
+ }
+ if (is.null(InputsModel$OutputsModel)) {
+ stop(
+ "'InputsModel' should contain an 'OutputsModel' key containing the output of the runoff of the downstream subcatchment"
+ )
+ }
+ if (is.null(InputsModel$OutputsModel$Qsim)) {
+ stop(
+ "'InputsModel$OutputsModel' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment"
+ )
+ }
+ if (sum(!is.na(InputsModel$OutputsModel$Qsim)) != length(RunOptions$IndPeriod_Run)) {
+ stop(
+ "'InputsModel$OutputsModel$Qim' should have the same lenght as 'RunOptions$IndPeriod_Run' and contain no NA"
+ )
+ }
+
+ OutputsModel <- InputsModel$OutputsModel
+ OutputsModel$QsimDown <- OutputsModel$Qsim
+
+ if (inherits(InputsModel, "hourly")) {
+ TimeStep <- 60 * 60
+ } else if (inherits(InputsModel, "daily")) {
+ TimeStep <- 60 * 60 * 24
+ } else {
+ stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
+ }
+
+ # propagation time from upstream meshes to outlet
+ PT <- InputsModel$LengthHydro / Param[1L] / TimeStep
+ HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
+
+ NbUpBasins <- length(InputsModel$LengthHydro)
+ LengthTs <- length(OutputsModel$QsimDown)
+ OutputsModel$Qsim <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
+
+ IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
+ if (length(IniSD) > 0) {
+ if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
+ stop(
+ sprintf(
+ "SD initial states has a length of %i and a length of %i is required",
+ length(IniSD),
+ sum(floor(PT)) + NbUpBasins
+ )
+ )
+ }
+ IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+ iStart <- 1
+ if (x > 1) {
+ iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
+ }
+ IniSD[iStart:(iStart + PT[x])]
+ })
+ } else {
+ IniStates <- lapply(seq_len(NbUpBasins), function(x) {
+ rep(0, floor(PT[x] + 1))
+ })
+ }
+
+ for (upstream_basin in seq_len(NbUpBasins)) {
+ Qupstream <- InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin]
+ if (!is.na(InputsModel$BasinAreas[upstream_basin])) {
+ # Upstream flow with area needs to be converted to m3 by time step
+ Qupstream <- Qupstream * InputsModel$BasinAreas[upstream_basin] * 1e3
+ }
+ Qupstream1 <- c(IniStates[[upstream_basin]][-length(IniStates[[upstream_basin]])], Qupstream[1:(LengthTs - floor(PT[upstream_basin]))])
+ Qupstream2 <- IniStates[[upstream_basin]]
+ if(LengthTs - floor(PT[upstream_basin]) - 1 > 0) Qupstream2 <- c(Qupstream2, Qupstream[1:(LengthTs - floor(PT[upstream_basin]) - 1)])
+ OutputsModel$Qsim <- OutputsModel$Qsim +
+ Qupstream1 * HUTRANS[1, upstream_basin] +
+ Qupstream2 * HUTRANS[2, upstream_basin]
+ }
+ # Warning for negative flows
+ if (any(OutputsModel$Qsim < 0)) {
+ warning(length(which(OutputsModel$Qsim < 0)), " time steps with negative flow, set to zero.")
+ OutputsModel$Qsim[OutputsModel$Qsim < 0] <- 0
+ }
+ # Convert back Qsim to mm
+ OutputsModel$Qsim <- OutputsModel$Qsim / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
+
+ if ("StateEnd" %in% RunOptions$Outputs_Sim) {
+ OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
+ LengthTs <- tail(RunOptions$IndPeriod_Run,1)
+ InputsModel$Qupstream[(LengthTs - floor(PT[x])):LengthTs, x]
+ })
+ }
+
+ return(OutputsModel)
+}
diff --git a/R/utils.R b/R/utils.R
index f2fac0a8ad679df4c64acaa2f24475da21786d35..4cf16933935889bd01ffdb97c01e305beb789d47 100644
--- a/R/utils.R
+++ b/R/utils.R
@@ -53,7 +53,6 @@ setDataToLocation <- function(control, supervisor) {
#' @return [NULL]
doSupervision <- function(supervisor) {
for (id in names(supervisor$controllers)) {
- ctrlr <- supervisor$controllers[[id]]
# Read Y from locations in the model
supervisor$controllers[[id]]$Y$v <-
sapply(supervisor$controllers[[id]]$Y$loc, getDataFromLocation, supervisor = supervisor)
@@ -65,3 +64,19 @@ doSupervision <- function(supervisor) {
}
return()
}
+
+#' Check the parameters of RunModel methods
+#'
+#' Stop the execution if an error is detected.
+#'
+#' @param RunOptions a `GRiwrmRunOptions` object (See [CreateRunOptions.GRiwrmInputsModel])
+#' @param Param a [list] of [numeric] containing model parameters of each node of the network
+#'
+#' @return [NULL]
+#'
+checkRunModelParameters <- function(InputsModel, RunOptions, Param) {
+ if(!inherits(InputsModel, "GRiwrmInputsModel")) stop("`InputsModel` parameter must of class 'GRiwrmRunoptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+ if(!inherits(RunOptions, "GRiwrmRunOptions")) stop("Argument `RunOptions` parameter must of class 'GRiwrmRunOptions' (See ?CreateRunOptions.GRiwrmInputsModel)")
+ if(!is.list(Param) || !all(names(InputsModel) %in% names(Param))) stop("Argument `Param` must be a list with names equal to nodes IDs")
+ return()
+}
diff --git a/tests/testthat/test-RunModel.R b/tests/testthat/test-RunModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..b12867f026b01acf99a3508f4b77751c613348a1
--- /dev/null
+++ b/tests/testthat/test-RunModel.R
@@ -0,0 +1,43 @@
+context("RunModel.Supervisor")
+
+test_that("RunModelSupervisor with no regulation should returns same results as RunModel.GRiwrmInputsModel", {
+ data(Severn)
+ # Network configuration
+ nodes <- Severn$BasinsInfo[c(1,2,5), c("gauge_id", "downstream_id", "distance_downstream", "area")]
+ nodes$distance_downstream <- nodes$distance_downstream * 1000 # Conversion km -> m
+ nodes$model <- NA
+ nodes$model[1] <- "RunModel_GR4J"
+ griwrm <- GRiwrm(nodes, list(id = "gauge_id", down = "downstream_id", length = "distance_downstream"))
+ # InputsModel
+ DatesR <- Severn$BasinsObs[[1]]$DatesR
+ PrecipTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$precipitation}))
+ PotEvapTot <- cbind(sapply(Severn$BasinsObs, function(x) {x$peti}))
+ Precip <- ConvertMeteoSD(griwrm, PrecipTot)
+ PotEvap <- ConvertMeteoSD(griwrm, PotEvapTot)
+ Qobs <- cbind(sapply(Severn$BasinsObs, function(x) {x$discharge_spec}))
+ InputsModel <- CreateInputsModel(griwrm, DatesR, Precip, PotEvap, Qobs)
+ # RunOptions
+ IndPeriod_Run <- seq(
+ length(InputsModel[[1]]$DatesR) - 365,
+ length(InputsModel[[1]]$DatesR)
+ )
+ IndPeriod_WarmUp = seq(IndPeriod_Run[1]-366,IndPeriod_Run[1]-1)
+ RunOptions <- CreateRunOptions(
+ InputsModel = InputsModel,
+ IndPeriod_WarmUp = IndPeriod_WarmUp,
+ IndPeriod_Run = IndPeriod_Run
+ )
+ Param <- list("54057" = c(0.727, 175.493, -0.082, 0.029, 4.654))
+ OM_GriwrmInputs <- RunModel(
+ InputsModel,
+ RunOptions = RunOptions,
+ Param = Param
+ )
+ supervisor <- CreateSupervisor(InputsModel)
+ OM_Supervisor <- RunModel(
+ supervisor,
+ RunOptions = RunOptions,
+ Param = Param
+ )
+ expect_equal(OM_Supervisor[["54057"]]$Qsim, OM_GriwrmInputs[["54057"]]$Qsim)
+})