diff --git a/R/CreateInputsModel.GRiwrm.R b/R/CreateInputsModel.GRiwrm.R
index 55d73e60fc5752c839b01a746652ce4b8579849b..d3e8b3244fe0977a498c0d2ff3b70388af542e2a 100644
--- a/R/CreateInputsModel.GRiwrm.R
+++ b/R/CreateInputsModel.GRiwrm.R
@@ -16,6 +16,9 @@
#' @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
+#' @param Qrelease (optional) [matrix] or [data.frame] of [numeric] containing
+#' release flows by nodes using the model `RunModel_Reservoir` \[m3 per
+#' time step\]
#' @param PrecipScale (optional) named [vector] of [logical] indicating if the
#' mean of the precipitation interpolated on the elevation layers must be
#' kept or not, required to create CemaNeige module inputs, default `TRUE`
@@ -69,6 +72,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
PotEvap = NULL,
Qobs = NULL,
Qmin = NULL,
+ Qrelease = NULL,
PrecipScale = TRUE,
TempMean = NULL, TempMin = NULL,
TempMax = NULL, ZInputs = NULL,
@@ -120,31 +124,14 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
}
})
- directFlowIds <- x$id[is.na(x$model) | x$model == "Diversion" | x$model == "RunModel_Reservoir"]
- if (length(directFlowIds) > 0) {
- err <- FALSE
- if (is.null(Qobs)) {
- err <- TRUE
- } else {
- Qobs <- as.matrix(Qobs)
- if (is.null(colnames(Qobs))) {
- err <- TRUE
- } else if (!all(directFlowIds %in% colnames(Qobs))) {
- err <- TRUE
- }
- }
- if (err) stop(sprintf("'Qobs' column names must at least contain %s", paste(directFlowIds, collapse = ", ")))
- }
- if (!all(colnames(Qobs) %in% directFlowIds)) {
- warning(
- "The following columns in 'Qobs' are ignored since they don't match with ",
- "Direction Injection (model=`NA`), ",
- "Reservoir (model=\"RunModelReservoir\"), ",
- "or Diversion nodes (model=\"Diversion\"): ",
- paste(setdiff(colnames(Qobs), directFlowIds), collapse = ", ")
- )
- Qobs <- Qobs[, directFlowIds]
- }
+ if (is.null(Qobs)) Qobs <- matrix(0, ncol = 0, nrow = length(DatesR))
+ if (is.null(Qrelease)) Qrelease <- matrix(0, ncol = 0, nrow = length(DatesR))
+ l <- updateQObsQrelease(g = x, Qobs = Qobs, Qrelease = Qrelease)
+ Qobs <- l$Qobs
+ Qrelease <- l$Qrelease
+ checkQobsQrelease(x, "Qobs", Qobs)
+ checkQobsQrelease(x, "Qrelease", Qrelease)
+
diversionRows <- getDiversionRows(x)
if (length(diversionRows) > 0) {
warn <- FALSE
@@ -171,7 +158,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
warning(
sprintf(
"'Qmin' would include the following columns %s.\n Zero values are applied by default.",
- paste(directFlowIds, collapse = ", ")
+ paste(x$id[diversionRows], collapse = ", ")
)
)
}
@@ -188,16 +175,6 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
InputsModel <- CreateEmptyGRiwrmInputsModel(x)
- # Qobs completion for at least filling Qupstream of all nodes by zeros
- Qobs0 <- matrix(0, nrow = length(DatesR), ncol = nrow(x))
- colnames(Qobs0) <- x$id
- if (is.null(Qobs)) {
- Qobs <- Qobs0
- } else {
- Qobs0[, colnames(Qobs)] <- Qobs
- Qobs <- Qobs0
- }
-
for(id in getNodeRanking(x)) {
message("CreateInputsModel.GRiwrm: Processing sub-basin ", id, "...")
@@ -216,6 +193,7 @@ CreateInputsModel.GRiwrm <- function(x, DatesR,
NLayers = getInputBV(NLayers, id, 5),
Qobs = Qobs,
Qmin = getInputBV(Qmin, id),
+ Qrelease = Qrelease,
IsHyst = IsHyst
)
}
@@ -251,7 +229,7 @@ CreateEmptyGRiwrmInputsModel <- function(griwrm) {
#'
#' @return \emph{InputsModel} object for one.
#' @noRd
-CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
+CreateOneGRiwrmInputsModel <- function(id, griwrm, DatesR, ..., Qobs, Qmin, Qrelease, IsHyst) {
np <- getNodeProperties(id, griwrm)
if (np$Diversion) {
@@ -272,7 +250,20 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
if(length(UpstreamNodeRows) > 0) {
# Sub-basin with hydraulic routing
- Qupstream <- as.matrix(Qobs[ , griwrm$id[UpstreamNodeRows], drop=FALSE])
+ Qupstream <- NULL
+ Qupstream <- as.matrix(cbind(
+ Qobs[ , colnames(Qobs)[colnames(Qobs) %in% griwrm$id[UpstreamNodeRows]], drop = FALSE],
+ Qrelease[ , colnames(Qrelease)[colnames(Qrelease) %in% griwrm$id[UpstreamNodeRows]], drop = FALSE]
+ ))
+ # Qupstream completion with zeros for all upstream nodes
+ Qupstream0 <- matrix(0, nrow = length(DatesR), ncol = length(UpstreamNodeRows))
+ colnames(Qupstream0) <- griwrm$id[UpstreamNodeRows]
+ if (is.null(Qupstream) || ncol(Qupstream) == 0) {
+ Qupstream <- Qupstream0
+ } else {
+ Qupstream0[, colnames(Qupstream)] <- Qupstream
+ Qupstream <- Qupstream0
+ }
upstreamDiversion <- which(
sapply(griwrm$id[UpstreamNodeRows],
function(id) {
@@ -308,6 +299,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
# Set model inputs with the **airGR** function
InputsModel <- CreateInputsModel(
FUN_MOD,
+ DatesR = DatesR,
...,
Qupstream = Qupstream,
LengthHydro = LengthHydro,
@@ -352,7 +344,8 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
InputsModel$diversionOutlet <- diversionOutlet
InputsModel$Qdiv <- -Qobs[, id]
InputsModel$Qmin <- Qmin
- } else if(np$Reservoir) {
+ }
+ if (np$Reservoir) {
# If an upstream node is ungauged and the donor is downstream then we are in an ungauged reduced network
iUpstreamUngaugedNodes <- which(griwrm$id %in% griwrm$id[UpstreamNodeRows] &
griwrm$model == "Ungauged")
@@ -361,7 +354,7 @@ CreateOneGRiwrmInputsModel <- function(id, griwrm, ..., Qobs, Qmin, IsHyst) {
InputsModel$isUngauged <- any(griwrm$donor[iUpstreamUngaugedNodes] == InputsModel$gaugedId)
}
# Fill reservoir release with Qobs
- InputsModel$Qrelease <- Qobs[, id]
+ InputsModel$Qrelease <- Qrelease[, id]
}
# Add class for S3 process (Prequel of HYCAR-Hydro/airgr#60)
diff --git a/R/RunModel.InputsModel.R b/R/RunModel.InputsModel.R
index f36e4b0e4a08f524678c1c1107cdba6b2aa888fb..8fb7a70e4c1271bb9dd5077bc7d57cd071af732a 100644
--- a/R/RunModel.InputsModel.R
+++ b/R/RunModel.InputsModel.R
@@ -70,7 +70,7 @@ RunModel.InputsModel <- function(x = NULL,
OutputsModel$RunOptions$WarmUpQsim_m3 <-
OutputsModel$RunOptions$WarmUpQsim * sum(x$BasinAreas, na.rm = TRUE) * 1e3
}
- if (x$hasDiversion) {
+ if (x$hasDiversion && !x$isReservoir) {
OutputsModel <- RunModel_Diversion(x, RunOptions, OutputsModel)
}
return(OutputsModel)
diff --git a/R/RunModel_Reservoir.R b/R/RunModel_Reservoir.R
index 0fff58a478736696bf4147ef5187d906cea377d1..3e20530576207bc9c1744900d5c410282e5530c6 100644
--- a/R/RunModel_Reservoir.R
+++ b/R/RunModel_Reservoir.R
@@ -68,10 +68,17 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
iPerTot <- seq(length(IndPerTot))
Vsim <- rep(0, length(IndPerTot))
Qsim_m3 <- Vsim
+ if (InputsModel$hasDiversion) {
+ Qdiv_m3 <- Vsim
+ }
# Time series volume and release calculation
for(i in iPerTot) {
Vsim[i] <- V0 + Qinflows_m3[i]
+ if (InputsModel$hasDiversion) {
+ Qdiv_m3[i] <- min(Vsim[i] + InputsModel$Qmin[IndPerTot[i]], InputsModel$Qdiv[IndPerTot[i]])
+ Vsim[i] <- Vsim[i] - Qdiv_m3[i]
+ }
Qsim_m3[i] <- min(Vsim[i], InputsModel$Qrelease[IndPerTot[i]])
Vsim[i] <- Vsim[i] - Qsim_m3[i]
if (Vsim[i] > Vmax) {
@@ -86,10 +93,16 @@ RunModel_Reservoir <- function(InputsModel, RunOptions, Param) {
iWarmUp <- seq(length(RunOptions$IndPeriod_WarmUp))
OutputsModel$RunOptions$WarmUpQsim_m3 <- Qsim_m3[iWarmUp]
OutputsModel$RunOptions$WarmUpVsim <- Vsim[iWarmUp]
+ if (InputsModel$hasDiversion) {
+ OutputsModel$RunOptions$WarmUpQdiv_m3 <- Qdiv_m3[iWarmUp]
+ }
}
iRun <- length(IndPerWarmUp) + seq(length(RunOptions$IndPeriod_Run))
OutputsModel$Qsim_m3 <- Qsim_m3[iRun]
OutputsModel$Vsim <- Vsim[iRun]
+ if (InputsModel$hasDiversion) {
+ OutputsModel$Qdiv_m3 <- Qdiv_m3[iRun]
+ }
OutputsModel$StateEnd$Reservoir <- list(V = Vsim[length(Vsim)])
class(OutputsModel) <- c("OutputsModelReservoir", class(OutputsModel))
return(OutputsModel)
diff --git a/R/utils.CreateInputsModel.R b/R/utils.CreateInputsModel.R
new file mode 100644
index 0000000000000000000000000000000000000000..c5a4ab9481327ec8df388d96361540e45ea50f34
--- /dev/null
+++ b/R/utils.CreateInputsModel.R
@@ -0,0 +1,59 @@
+updateQObsQrelease <- function(g, Qobs, Qrelease) {
+ reservoirIds <- g$id[!is.na(g$model) & g$model == "RunModel_Reservoir"]
+ # Fill Qrelease with Qobs
+ warn_ids <- NULL
+ for(id in reservoirIds) {
+ if (!id %in% colnames(Qrelease)) {
+ if (id %in% colnames(Qobs)) {
+ if (!any(g$id == id & (!is.na(g$model) & g$model == "Diversion"))) {
+ if (is.null(Qrelease)) {
+ Qrelease = Qobs[, id, drop = FALSE]
+ } else {
+ Qrelease = cbind(Qrelease, Qobs[, id, drop = FALSE])
+ }
+ Qobs <- Qobs[, colnames(Qobs) != id, drop = FALSE]
+ warn_ids = c(warn_ids, id)
+ }
+ }
+ }
+ }
+ if (!is.null(warn_ids)) {
+ warning("Use of the `Qobs` parameter for reservoir releases is deprecated, please use `Qrelease` instead.\n",
+ "Processing `Qrelease <- cbind(Qrelease, Qobs[, c(\"", paste(warn_ids, collapse = "\", `"), "\"])`...")
+ }
+ return(list(Qobs = Qobs, Qrelease = Qrelease))
+}
+
+checkQobsQrelease <- function(g, varname, Q) {
+ if (varname == "Qobs") {
+ directFlowIds <- g$id[is.na(g$model) | g$model == "Diversion"]
+ } else {
+ directFlowIds <- g$id[!is.na(g$model) & g$model == "RunModel_Reservoir"]
+ }
+ if (length(directFlowIds) > 0) {
+ err <- FALSE
+ if (is.null(Q)) {
+ err <- TRUE
+ } else {
+ Q <- as.matrix(Q)
+ if (is.null(colnames(Q))) {
+ err <- TRUE
+ } else if (!all(directFlowIds %in% colnames(Q))) {
+ err <- TRUE
+ }
+ }
+ if (err) stop(sprintf("'%s' column names must at least contain %s", varname, paste(directFlowIds, collapse = ", ")))
+ }
+ if (!all(colnames(Q) %in% directFlowIds)) {
+ warning(
+ sprintf("The following columns in '%s' are ignored since they don't match with ", varname),
+ ifelse(varname == "Qobs",
+ c("Direction Injection (model=`NA`), ",
+ "or Diversion nodes (model=\"Diversion\"): "),
+ "Reservoir nodes (model=\"RunModelReservoir\"): "),
+ paste(setdiff(colnames(Q), directFlowIds), collapse = ", ")
+ )
+ Q <- Q[, directFlowIds]
+ }
+ return(Q)
+}
diff --git a/man/CreateInputsModel.GRiwrm.Rd b/man/CreateInputsModel.GRiwrm.Rd
index 26b9a66bc318a276b31f98ce25216ad896e2e9db..08a947ba73830f918da711b6fbd399a949f0b7df 100644
--- a/man/CreateInputsModel.GRiwrm.Rd
+++ b/man/CreateInputsModel.GRiwrm.Rd
@@ -11,6 +11,7 @@
PotEvap = NULL,
Qobs = NULL,
Qmin = NULL,
+ Qrelease = NULL,
PrecipScale = TRUE,
TempMean = NULL,
TempMin = NULL,
@@ -45,6 +46,10 @@ the model and positive flows are injected to the model}
minimum flows to let downstream of a node with a Diversion [m3 per
time step]. Default is zero. Column names correspond to node IDs}
+\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]}
+
\item{PrecipScale}{(optional) named \link{vector} of \link{logical} indicating if the
mean of the precipitation interpolated on the elevation layers must be
kept or not, required to create CemaNeige module inputs, default \code{TRUE}
diff --git a/tests/testthat/helper_1_RunModel.R b/tests/testthat/helper_1_RunModel.R
index 6db2a92ac21ebb1bded8fdebe629a116cfa74f1a..ef1d0b7cf85fa9ebf32c16f8739a224e2616a20d 100644
--- a/tests/testthat/helper_1_RunModel.R
+++ b/tests/testthat/helper_1_RunModel.R
@@ -16,6 +16,7 @@ setupRunModel <-
runRunModel = TRUE,
griwrm = NULL,
Qobs2 = NULL,
+ Qrelease = NULL,
IsHyst = FALSE) {
data(Severn)
@@ -50,6 +51,49 @@ setupRunModel <-
TempMean <- NULL
}
+ # Calibration parameters
+ ParamMichel <- list(
+ `54057` = c(
+ 0.781180650559296,
+ 74.4247133147015,
+ -1.26590474908235,
+ 0.96012365697022,
+ 2.51101785373787
+ ),
+ `54032` = c(
+ 0.992743594649893,
+ 1327.19309205366,
+ -0.505902143697436,
+ 7.91553615210291,
+ 2.03604525989572
+ ),
+ `54001` = c(
+ 1.03,
+ 24.7790862245877,
+ -1.90430150145153,
+ 21.7584023961971,
+ 1.37837837837838
+ ),
+ `54095` = c(
+ 256.844150254651,
+ 0.0650458497009288,
+ 57.523675209819,
+ 2.71809513102128
+ ),
+ `54002` = c(
+ 419.437754485522,
+ 0.12473266292168,
+ 13.0379482833606,
+ 2.12230907892238
+ ),
+ `54029` = c(
+ 219.203385553954,
+ 0.389211590110934,
+ 48.4242150713452,
+ 2.00300300300301
+ )
+ )
+
# set up inputs
if (!runInputsModel)
return(environment())
@@ -57,6 +101,7 @@ setupRunModel <-
suppressWarnings(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
TempMean = TempMean,
Qobs = Qobs2,
+ Qrelease = Qrelease,
IsHyst = IsHyst))
# RunOptions
@@ -83,49 +128,5 @@ setupRunOptions <- function(InputsModel) {
RunOptions <- CreateRunOptions(InputsModel,
IndPeriod_WarmUp = IndPeriod_WarmUp,
IndPeriod_Run = IndPeriod_Run)
-
- # Calibration parameters
- ParamMichel <- list(
- `54057` = c(
- 0.781180650559296,
- 74.4247133147015,
- -1.26590474908235,
- 0.96012365697022,
- 2.51101785373787
- ),
- `54032` = c(
- 0.992743594649893,
- 1327.19309205366,
- -0.505902143697436,
- 7.91553615210291,
- 2.03604525989572
- ),
- `54001` = c(
- 1.03,
- 24.7790862245877,
- -1.90430150145153,
- 21.7584023961971,
- 1.37837837837838
- ),
- `54095` = c(
- 256.844150254651,
- 0.0650458497009288,
- 57.523675209819,
- 2.71809513102128
- ),
- `54002` = c(
- 419.437754485522,
- 0.12473266292168,
- 13.0379482833606,
- 2.12230907892238
- ),
- `54029` = c(
- 219.203385553954,
- 0.389211590110934,
- 48.4242150713452,
- 2.00300300300301
- )
- )
-
return(environment())
}
diff --git a/tests/testthat/test-CreateInputsModel.R b/tests/testthat/test-CreateInputsModel.R
index f8695a76068709e601d1413ed788cacfc7b087c6..fa6acf9d2a2575f7ebcadebd14c8e00929b125f4 100644
--- a/tests/testthat/test-CreateInputsModel.R
+++ b/tests/testthat/test-CreateInputsModel.R
@@ -272,3 +272,12 @@ test_that("Node with upstream nodes having area = NA should return correct Basin
expect_equal(sum(InputsModel$`54001`$BasinAreas),
g$area[g$id == "54001"])
})
+
+test_that("Use of Qobs for Qrelease should raise a warning", {
+ g <- CreateGRiwrm(n_rsrvr)
+ e <- setupRunModel(griwrm = g, runInputsModel = FALSE)
+ for(x in ls(e)) assign(x, get(x, e))
+ expect_warning(CreateInputsModel(griwrm, DatesR, Precip, PotEvap,
+ TempMean = TempMean,
+ Qobs = Qobs_rsrvr))
+})
diff --git a/tests/testthat/test-CreateSupervisor.R b/tests/testthat/test-CreateSupervisor.R
index 3ba3e848c467626651b4147bec5d774bcda3d480..a484f708895b9a646ce0fdf56eeaac600e360bbb 100644
--- a/tests/testthat/test-CreateSupervisor.R
+++ b/tests/testthat/test-CreateSupervisor.R
@@ -54,7 +54,7 @@ test_that("CreateSupervisor using reservoir and diversion", {
))
g <- CreateGRiwrm(nodes)
# Add Qobs for the 2 new nodes and create InputsModel
- Qobs <- matrix(data = rep(0, 2*length(DatesR)), ncol = 2)
+ Qobs <- matrix(data = 0, ncol = 2, nrow = length(DatesR))
colnames(Qobs) <- c("54029", "Reservoir")
InputsModel <- suppressWarnings(
CreateInputsModel(g, DatesR, Precip, PotEvap, Qobs)
diff --git a/tests/testthat/test-RunModel_Reservoir.R b/tests/testthat/test-RunModel_Reservoir.R
index bd227a0947ed96808fe4934c159e356ca97355c9..7417935bf54e330d38bc31df006271aeac34de6a 100644
--- a/tests/testthat/test-RunModel_Reservoir.R
+++ b/tests/testthat/test-RunModel_Reservoir.R
@@ -98,3 +98,40 @@ test_that("Calibration with ungauged node and reservoir filled by a diversion wo
colnames(Qobs2) <- c("Dam", "54095")
expect_dam(n_derived_rsrvr, Qobs2)
})
+
+test_that("Diversion on a reservoir works #146", {
+ e <- setupRunModel(runInputsModel = FALSE)
+ for (x in ls(e)) assign(x, get(x, e))
+ Qrelease <- data.frame(Dam = rep(3508465, length(DatesR)))
+ Param <- c(ParamMichel[names(ParamMichel) %in% griwrm$id], list(Dam = c(10E6, 1)))
+ e <- setupRunModel(runRunModel = FALSE,
+ griwrm = CreateGRiwrm(n_rsrvr),
+ Qrelease = Qrelease)
+ for (x in ls(e)) assign(x, get(x, e))
+ OM_resOnly <- RunModel(InputsModel,
+ RunOptions = RunOptions,
+ Param = Param)
+ nodes <- rbind(
+ n_rsrvr,
+ data.frame(
+ id = "Dam",
+ down = NA,
+ length = NA,
+ area = NA,
+ model = "Diversion"
+ )
+ )
+ Qobs2 <- Qrelease * 0.1
+ g <- CreateGRiwrm(nodes)
+ e <- setupRunModel(griwrm = g,
+ runRunModel = FALSE,
+ Qobs2 = Qobs2,
+ Qrelease = Qrelease)
+ for (x in ls(e)) assign(x, get(x, e))
+
+ OM <- RunModel(InputsModel,
+ RunOptions = RunOptions,
+ Param = Param)
+ expect_true(max(OM$Dam$Vsim) - min(OM$Dam$Vsim) > 0)
+ expect_false(all(OM$Dam$Vsim == OM_resOnly$Dam$Vsim))
+})