From cb8b0be58bbb39079ea3355876117318ab2144f8 Mon Sep 17 00:00:00 2001
From: Dorchies David <david.dorchies@inrae.fr>
Date: Thu, 17 Jun 2021 14:39:56 +0200
Subject: [PATCH] feat(RunModel_Lag): handling of warm-up simulated flows
- Implementation and tests
Refs #123
---
R/RunModel_Lag.R | 18 +++++++++++++++---
tests/testthat/test-RunModel_Lag.R | 25 ++++++++++++++++++++-----
2 files changed, 35 insertions(+), 8 deletions(-)
diff --git a/R/RunModel_Lag.R b/R/RunModel_Lag.R
index 473e7758..29712446 100644
--- a/R/RunModel_Lag.R
+++ b/R/RunModel_Lag.R
@@ -69,9 +69,21 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
IniSD[iStart:(iStart + PT[x])]
})
} else {
- IniStates <- lapply(seq_len(NbUpBasins), function(x) {
- rep(0, floor(PT[x] + 1))
- })
+ IniStates <- lapply(
+ seq_len(NbUpBasins),
+ function(iUpBasins) {
+ iWarmUp <- seq.int(
+ max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
+ max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
+ )
+ ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
+ if(length(ini) != floor(PT[iUpBasins] + 1)) {
+ # If warm-up period is not enough long complete beginning with first value
+ ini <- c(rep(ini[1], floor(PT[iUpBasins] + 1) - length(ini)), ini)
+ }
+ return(ini)
+ }
+ )
}
# message("Initstates: ", paste(IniStates, collapse = ", "))
diff --git a/tests/testthat/test-RunModel_Lag.R b/tests/testthat/test-RunModel_Lag.R
index c7cc53ec..d71189c1 100644
--- a/tests/testthat/test-RunModel_Lag.R
+++ b/tests/testthat/test-RunModel_Lag.R
@@ -165,7 +165,7 @@ Qm3GR4Only <- OutputsGR4JOnly$Qsim * InputsModel$BasinAreas[2L] * 1e3
test_that("1 input with lag of 1 time step delay out gives an output delayed of one time step", {
OutputsSD <- RunModel(InputsModel, RunOptions, Param = ParamSD, FUN_MOD = RunModel_GR4J)
Qm3SdSim <- OutputsSD$Qsim_m3
- Qm3UpstLagObs <- c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * InputsModel$BasinAreas[1] * 1e3
+ Qm3UpstLagObs <- Qupstream[Ind_Run - 1] * InputsModel$BasinAreas[1] * 1e3
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
@@ -174,7 +174,7 @@ test_that("1 input with lag of 0.5 time step delay out gives an output delayed o
Param = c(InputsModel$LengthHydro * 1e3 / (12 * 3600), Param),
FUN_MOD = RunModel_GR4J)
Qm3SdSim <- OutputsSD$Qsim_m3
- Qm3UpstLagObs <- (Qupstream[Ind_Run] + c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])) / 2 * InputsModel$BasinAreas[1] * 1e3
+ Qm3UpstLagObs <- (Qupstream[Ind_Run] + Qupstream[Ind_Run - 1]) / 2 * InputsModel$BasinAreas[1] * 1e3
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
@@ -233,7 +233,7 @@ InputsCrit <- CreateInputsCrit(
InputsModel = InputsModel,
RunOptions = RunOptions,
VarObs = "Q",
- Obs = (c(0, Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]]) * BasinAreas[1L] +
+ Obs = (Qupstream[Ind_Run - 1] * BasinAreas[1L] +
BasinObs$Qmm[Ind_Run] * BasinAreas[2L]) / sum(BasinAreas)
)
@@ -283,7 +283,7 @@ test_that("1 no area input with lag of 1 time step delay out gives an output del
expect_false(any(is.na(OutputsSD$Qsim)))
Qm3SdSim <- OutputsSD$Qsim_m3
- Qm3UpstLagObs <- c(0, InputsModel$Qupstream[Ind_Run[1:(length(Ind_Run) - 1)]])
+ Qm3UpstLagObs <- InputsModel$Qupstream[Ind_Run - 1]
expect_equal(Qm3SdSim - Qm3GR4Only, Qm3UpstLagObs)
})
@@ -338,6 +338,21 @@ test_that("Error on Wrong length of iniState$SD", {
InputsModel = IM, IndPeriod_Run = Ind_Run2,
IndPeriod_WarmUp = 0L,
IniStates = OutputsModel1$StateEnd)
- expect_error(RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
+ expect_error(
+ RunModel(InputsModel = IM, RunOptions = RunOptions2, Param = PSDini, FUN_MOD = RunModel_GR4J)
)
})
+
+test_that("First Qupstream time steps must be repeated if warm-up period is too short", {
+ IM2 <- IM[2558:3557]
+ IM2$BasinAreas[3] <- 0
+ IM2$Qupstream <- matrix(rep(1:1000, 2), ncol = 2)
+ RunOptions2 <- CreateRunOptions(FUN_MOD = RunModel_GR4J,
+ InputsModel = IM2, IndPeriod_Run = seq_len(1000),
+ IndPeriod_WarmUp = 0L)
+ OM2 <- RunModel(InputsModel = IM2,
+ RunOptions = RunOptions2,
+ Param = PSDini,
+ FUN_MOD = RunModel_GR4J)
+ expect_equal(OM2$Qsim_m3[1:3], rep(2,3))
+})
--
GitLab