Skip to content

Commit

Permalink
Merge branch '179-runmodel-outputsmodel-results-are-not-correct' into…
Browse files Browse the repository at this point in the history 
… 'dev'

Resolve "RunModel.OutputsModel results are not correct"

Closes #179

See merge request in-wop/airGRiwrm!107
  • Loading branch information
Dorchies David committed Jan 1, 2025
2 parents 0e1f5f9 + 8817f48 commit fd1262c
Show file tree
Hide file tree
Showing 5 changed files with 150 additions and 66 deletions.
2 changes: 1 addition & 1 deletion R/RunModel.GRiwrmOutputsModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
Expand Down
2 changes: 1 addition & 1 deletion R/RunModel.Supervisor.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) {
Expand Down
26 changes: 25 additions & 1 deletion R/utils.RunModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
}

Expand All @@ -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)
}

Expand Down Expand Up @@ -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)
}

Expand All @@ -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)
}

Expand Down
90 changes: 46 additions & 44 deletions tests/testthat/helper_1_RunModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,8 @@ setupRunModel <-
Qinf = NULL,
Qrelease = NULL,
Qmin = NULL,
IsHyst = FALSE) {
IsHyst = FALSE,
ParamMichel = getDefaultParamMichel()) {

data(Severn)

Expand Down Expand Up @@ -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())
Expand Down Expand Up @@ -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
)
)
}
96 changes: 77 additions & 19 deletions tests/testthat/test-RunModel.GRiwrmOutputsModel.R
View file
Original file line number Diff line number Diff line change
@@ -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,
Expand All @@ -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)),
Expand All @@ -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,
Expand All @@ -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),
Expand All @@ -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)
Expand Down

0 comments on commit fd1262c

Please sign in to comment.