diff --git a/R/RunModel.GRiwrmOutputsModel.R b/R/RunModel.GRiwrmOutputsModel.R
index cf8d00daae411a0ad29deccf41e662f5fabfd53e..878e614c2dced2902480983869edbabbba82a701 100644
--- a/R/RunModel.GRiwrmOutputsModel.R
+++ b/R/RunModel.GRiwrmOutputsModel.R
@@ -73,7 +73,7 @@ RunModel.GRiwrmOutputsModel <- function(x,
for (id in names(RunOptions)) {
# Run model for the sub-basin and one time step
RunOptions[[id]]$IniResLevels <- NULL
- RunOptions[[id]]$IniStates <- serializeIniStates(x[[id]]$StateEnd)
+ RunOptions[[id]]$IniStates <- serializeIniStates(x[[id]]$StateEnd, InputsModel[[id]])
RunOptions[[id]]$IndPeriod_WarmUp <- 0L
RunOptions[[id]]$IndPeriod_Run <- IndPeriod_Run
}
diff --git a/R/RunModel.Supervisor.R b/R/RunModel.Supervisor.R
index 57a59f2cf0112f237578a40a04b3cae268c67c47..29591227d920cbad99303c0b707a06a457555497 100644
--- a/R/RunModel.Supervisor.R
+++ b/R/RunModel.Supervisor.R
@@ -120,7 +120,7 @@ RunModel.Supervisor <- function(x, RunOptions, Param, ...) {
# Loop over sub-basin using SD model
for (id in SD_Ids) {
# Run model for the sub-basin and one time step
- RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd)
+ RunOptions[[id]]$IniStates <- serializeIniStates(x$OutputsModel[[id]]$StateEnd, x$InputsModel[[id]])
RunOptions[[id]]$IndPeriod_Run <- iTS
# Route upstream flows for SD nodes
if (x$InputsModel[[id]]$isReservoir) {
diff --git a/R/utils.RunModel.R b/R/utils.RunModel.R
index dbf187bfb894844598d9b09ad1382b72f0005228..71b0512fbd7f586e2f0911881dd1793677e1d388 100644
--- a/R/utils.RunModel.R
+++ b/R/utils.RunModel.R
@@ -49,6 +49,7 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
dfQsim <- cbind(data.frame(DatesR = InputsModel[[1]]$DatesR[IndPeriod_Run]),
do.call(cbind,lQsim) / attr(InputsModel, "TimeStep"))
dfQsim <- as.Qm3s(dfQsim)
+ rownames(dfQsim) <- NULL
return(dfQsim)
}
@@ -60,8 +61,27 @@ OutputsModelQsim <- function(InputsModel, OutputsModel, IndPeriod_Run) {
#' @return A vector as in `RunOptions$IniStates`
#' @noRd
#'
-serializeIniStates <- function(IniStates) {
+serializeIniStates <- function(IniStates, InputsModel) {
+ if (!is.list(IniStates)) return(IniStates)
+ ObjectClass <- class(InputsModel)
+ if (!"CemaNeige" %in% ObjectClass && any(is.na(IniStates$CemaNeigeLayers$G))) {
+ IniStates$CemaNeigeLayers$G <- NULL
+ }
+ if (!"CemaNeige" %in% ObjectClass && any(is.na(IniStates$CemaNeigeLayers$eTG))) {
+ IniStates$CemaNeigeLayers$eTG <- NULL
+ }
+ if (!"CemaNeige" %in% ObjectClass && any(is.na(IniStates$CemaNeigeLayers$Gthr))) {
+ IniStates$CemaNeigeLayers$Gthr <- NULL
+ }
+ if (!"CemaNeige" %in% ObjectClass && any(is.na(IniStates$CemaNeigeLayers$Glocmax))) {
+ IniStates$CemaNeigeLayers$Glocmax <- NULL
+ }
+ IniStates$Store$Rest <- rep(NA, 3)
IniStates <- unlist(IniStates)
+ IniStates[is.na(IniStates) & !grepl("SD", names(IniStates))] <- 0
+ if ("monthly" %in% ObjectClass) {
+ IniStates <- IniStates[seq_len(NState)]
+ }
return(IniStates)
}
@@ -128,6 +148,9 @@ merge.OutputsModel <- function(x, y, ...) {
for (item in items) {
y[[item]] <- c(x[[item]], y[[item]])
}
+ # We keep original warm-up data
+ if (!is.null(x$RunOptions$WarmUpQsim)) y$RunOptions$WarmUpQsim <- x$RunOptions$WarmUpQsim
+ if (!is.null(x$RunOptions$WarmUpQsim_m3)) y$RunOptions$WarmUpQsim_m3 <- x$RunOptions$WarmUpQsim_m3
return(y)
}
@@ -140,6 +163,7 @@ merge.GRiwrmOutputsModel <- function(x, y, ...) {
})
attributes(y) <- y_attributes
attr(y, "Qm3s") <- rbind(attr(x, "Qm3s"), attr(y, "Qm3s"))
+ rownames(attr(y, "Qm3s")) <- NULL
return(y)
}
diff --git a/tests/testthat/helper_1_RunModel.R b/tests/testthat/helper_1_RunModel.R
index 57110eb77e2d6273cdcd813668b05a6d64bb3681..81423de9c271f38cbe6f014694ac94922716c4ff 100644
--- a/tests/testthat/helper_1_RunModel.R
+++ b/tests/testthat/helper_1_RunModel.R
@@ -18,7 +18,8 @@ setupRunModel <-
Qinf = NULL,
Qrelease = NULL,
Qmin = NULL,
- IsHyst = FALSE) {
+ IsHyst = FALSE,
+ ParamMichel = getDefaultParamMichel()) {
data(Severn)
@@ -52,49 +53,6 @@ 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())
@@ -132,3 +90,47 @@ setupRunOptions <- function(InputsModel) {
IndPeriod_Run = IndPeriod_Run)
return(environment())
}
+
+getDefaultParamMichel <- function() {
+ 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
+ )
+ )
+}
diff --git a/tests/testthat/test-RunModel.GRiwrmOutputsModel.R b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
index d0dc61299b1dc0ab8726a872cba870a1e2745375..134a5a02f1a6595bd264b4ee09decad85400cfdd 100644
--- a/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
+++ b/tests/testthat/test-RunModel.GRiwrmOutputsModel.R
@@ -1,5 +1,62 @@
skip_on_cran()
+data(Severn)
+
+test_that("Single node returns same result as RunModel.GRiwrmInputsModel", {
+ nodes <- loadSevernNodes()
+ nodes <- nodes[nodes$id == "54029", ]
+ nodes$down <- NA_character_
+ nodes$length <- NA_real_
+ e <- setupRunModel(
+ runRunOptions = FALSE,
+ griwrm = CreateGRiwrm(nodes)
+ )
+ for (x in ls(e)) assign(x, get(x, e))
+ ROref <- CreateRunOptions(
+ InputsModel,
+ IndPeriod_WarmUp = 1:364,
+ IndPeriod_Run = 365:366
+ )
+ OMref <- RunModel(
+ InputsModel,
+ RunOptions = ROref,
+ Param = ParamMichel
+ )
+ ROwarmUp <- CreateRunOptions(
+ InputsModel,
+ IndPeriod_WarmUp = 1:364,
+ IndPeriod_Run = 365L
+ )
+ OMwarmUp <- RunModel(
+ InputsModel,
+ RunOptions = ROwarmUp,
+ Param = ParamMichel
+ )
+ ROhotStart <- CreateRunOptions(
+ InputsModel,
+ IniStates = lapply(OMwarmUp, "[[", "StateEnd"),
+ IndPeriod_WarmUp = 0L,
+ IndPeriod_Run = 366L
+ )
+ ROhotStart$`54029`$IniResLevels <- NULL
+ # State Initiation
+ ROtest <- ROwarmUp
+ for (id in names(ROtest)) {
+ # Run model for the sub-basin and one time step
+ ROtest[[id]]$IniResLevels <- NULL
+ ROtest[[id]]$IniStates <- serializeIniStates(OMwarmUp[[id]]$StateEnd, InputsModel[[id]])
+ ROtest[[id]]$IndPeriod_WarmUp <- 0L
+ ROtest[[id]]$IndPeriod_Run <- 366L
+ }
+ expect_equal(ROtest$`54029`, ROhotStart$`54029`)
+ OMtest <- RunModel(OMwarmUp,
+ InputsModel = InputsModel,
+ RunOptions = ROwarmUp,
+ IndPeriod_Run = 366L
+ )
+ expect_equal(OMtest$`54029`, OMref$`54029`)
+})
+
# Setup model
griwrm <- CreateGRiwrm(rbind(
n_derived_rsrvr,
@@ -11,8 +68,7 @@ griwrm <- CreateGRiwrm(rbind(
model = NA
)
))
-data(Severn)
-DatesR <- Severn$BasinsObs[[1]]$DatesR
+DatesR <- Severn$BasinsObs[[1]]$DatesR
Qinf <- data.frame(
# Diversion to the dam
`54095` = rep(-1E6, length(DatesR)),
@@ -26,40 +82,42 @@ Qrelease <- data.frame(Dam = rep(100E3, length(DatesR)))
Qmin <- data.frame("54095" = rep(3E6, length(DatesR)))
names(Qmin) <- "54095"
e <- setupRunModel(
+ runRunModel = FALSE,
griwrm = griwrm,
Qinf = Qinf,
Qrelease = Qrelease,
- Qmin = Qmin,
- runRunOptions = FALSE
+ Qmin = Qmin
)
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
+# Simulation periods up to 31/12/1986
dfTS <- data.frame(
DatesR = DatesR,
yearmonth = format(DatesR, "%Y-%m")
)
-dfTS <- dfTS[1:(which(dfTS$yearmonth == "1987-01")[1]), ]
+dfTS <- dfTS[1:(which(dfTS$yearmonth == "1987-01")[1] - 1), ]
-test_that("RunModel.GRiwrmOutputsModel works with InputsModel", {
+# Run simulation in "normal" mode
+Param <- c(ParamMichel[names(ParamMichel) %in% griwrm$id], list(Dam = c(100E6, 1)))
+ROref <- CreateRunOptions(InputsModel, IndPeriod_WarmUp = 1:364, IndPeriod_Run = 365:nrow(dfTS))
+OMref <- RunModel(InputsModel, ROref, Param)
- for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
+# Set up initial conditions
+ROO <- CreateRunOptions(InputsModel, IndPeriod_WarmUp = 1:364, IndPeriod_Run = 365L)
+OM <- RunModel(InputsModel, ROO, Param)
+test_that("RunModel.GRiwrmOutputsModel works with InputsModel", {
+ for (ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
# 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
- }
+ # 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 = InputsModel,
RunOptions = RunOptions,
@@ -69,10 +127,10 @@ test_that("RunModel.GRiwrmOutputsModel works with InputsModel", {
expect_equal(nrow(attr(OM, "Qm3s")), nrow(dfTS) - 364)
expect_equal(length(OM[[1]]$DatesR), nrow(dfTS) - 364)
+ expect_equal(attr(OM, "Qm3s"), attr(OMref, "Qm3s"))
})
test_that("RunModel.GRiwrmOutputsModel works with Supervisor", {
-
sv <- CreateSupervisor(InputsModel)
curve <- approx(x = c(31*11 - 365, 30 * 6, 31 * 11, 366 + 30 * 6),
@@ -94,7 +152,7 @@ test_that("RunModel.GRiwrmOutputsModel works with Supervisor", {
U = "Dam",
fn_guide_curve)
- for(ym in unique(dfTS$yearmonth[dfTS$DatesR > OM[[1]]$DatesR])) {
+ 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)