From 10f48a097acc184d3781483e27058527fcf05a29 Mon Sep 17 00:00:00 2001
From: David <david.dorchies@inrae.fr>
Date: Thu, 5 Dec 2024 03:18:48 +0100
Subject: [PATCH] fix(RunModel.GRiwrmOutputsModel): test and debug final state
Refs #174
---
R/RunModel.Supervisor.R | 12 +-
.../test-RunModel.GRiwrmOutputsModel.R | 135 ++++++++++++------
2 files changed, 99 insertions(+), 48 deletions(-)
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index 5dd28fa..57a59f2 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -159,20 +159,24 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
message(" 100%")
for (id in SD_Ids) {
- StateEnd <- x$OutputsModel[[id]]$StateEnd
if (inherits(x$InputsModel[[id]], "GR")) {
+ # Copy OutputsModel from original run of GR contribution for the whole time series
+ # Except for StateEnd which has been modified
+ StateEnd <- x$OutputsModel[[id]]$StateEnd
x$OutputsModel[[id]] <- OutputsModelGR[[id]]
+ class_StateEnd <- class(x$OutputsModel[[id]]$StateEnd)
+ x$OutputsModel[[id]]$StateEnd <- c(x$OutputsModel[[id]]$StateEnd, StateEnd)
+ class(x$OutputsModel[[id]]$StateEnd) <- class_StateEnd
} else {
+ # Add missing DatesR for non GR models
x$OutputsModel[[id]]$DatesR <- x$InputsModel[[1]]$DatesR[IndPeriod_Run]
}
+ # Copy stored outputs for the whole time series
for (outputVar in outputVars[[id]]) {
x$OutputsModel[[id]][[outputVar]] <- x$storedOutputs[[outputVar]][, id]
}
x$OutputsModel[[id]]$Qsim <-
x$storedOutputs$Qsim_m3[, id] / sum(x$InputsModel[[id]]$BasinAreas, na.rm = TRUE) / 1e3
- class_StateEnd <- class(x$OutputsModel[[id]]$StateEnd)
- x$OutputsModel[[id]]$StateEnd <- c(x$OutputsModel[[id]]$StateEnd, StateEnd)
- class(x$OutputsModel[[id]]$StateEnd) <- class_StateEnd
x$OutputsModel[[id]]$RunOptions$WarmUpQsim_m3 <- OM_WarmUp[[id]]$Qsim_m3
x$OutputsModel[[id]]$RunOptions$WarmUpQsim <- OM_WarmUp[[id]]$Qsim_m3 /
sum(x$InputsModel[[id]]$BasinAreas, na.rm = TRUE) / 1e3
diff --git a/tests/testthat/test-RunModel.GRiwrmOutputsModel.R b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
index ee05413..d0dc612 100644
--- a/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
+++ b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
@@ -1,49 +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(
- # 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
+
+# Setup model
+griwrm <- CreateGRiwrm(rbind(
+ n_derived_rsrvr,
+ data.frame(
+ id = "WD",
+ down = "Dam",
+ length = 0,
+ area = NA,
+ model = NA
)
- for (x in ls(e)) assign(x, get(x, e))
+))
+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)
+# 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")
- )
+# Loop over periods months periods
+dfTS <- data.frame(
+ DatesR = DatesR,
+ yearmonth = format(DatesR, "%Y-%m")
+)
+dfTS <- dfTS[1:(which(dfTS$yearmonth == "1987-01")[1]), ]
+
+test_that("RunModel.GRiwrmOutputsModel works with InputsModel", {
for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
@@ -64,6 +67,50 @@ test_that("RunModel.GRiwrmOutputsModel works", {
Qinf = ym_Qinf)
}
- expect_equal(nrow(attr(OM, "Qm3s")), length(DatesR) - 364)
- expect_equal(length(OM[[1]]$DatesR), length(DatesR) - 364)
+ expect_equal(nrow(attr(OM, "Qm3s")), nrow(dfTS) - 364)
+ expect_equal(length(OM[[1]]$DatesR), nrow(dfTS) - 364)
+})
+
+test_that("RunModel.GRiwrmOutputsModel works with Supervisor", {
+
+ sv <- CreateSupervisor(InputsModel)
+
+ curve <- approx(x = c(31*11 - 365, 30 * 6, 31 * 11, 366 + 30 * 6),
+ y = c(20E6, 90E6, 20E6, 90E6),
+ xout = 1:366)$y
+
+ fn_guide_curve_factory <- function(sv, curve) {
+ function(y) {
+ # How much to release for reaching the filling curve ?
+ deltaV <- sv$OutputsModel$Dam$Vsim - curve[lubridate::yday(sv$ts.date)]
+ # Minimum 500L/s and max 1E6 m3
+ return(max(86400 / 2, min(1E6, deltaV)))
+ }
+ }
+ fn_guide_curve <- fn_guide_curve_factory(sv, curve)
+
+ CreateController(sv, "dam_filling_curve",
+ Y = NULL,
+ U = "Dam",
+ fn_guide_curve)
+
+ for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
+ message("Processing period ", ym)
+ # 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 = sv,
+ RunOptions = RunOptions,
+ IndPeriod_Run = ym_IndPeriod_Run,
+ Qinf = ym_Qinf)
+ }
+ expect_gte(sum(attr(OM, "Qm3s")$WD), sum(Qinf$WD[1:(nrow(dfTS) - 364)]))
})
--
GitLab