Commit 51258480 authored by Dorchies David's avatar Dorchies David
Browse files

feat(RunModel_Lag): add warm-up period simulation handling

Refs #132
parent bd42c74a
RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
NParam <- 1
##Arguments_check
## argument check
if (!inherits(InputsModel, "InputsModel")) {
stop("'InputsModel' must be of class 'InputsModel'")
}
......@@ -21,35 +21,49 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
if (is.null(QcontribDown$Qsim)) {
stop("'QcontribDown' should contain a key 'Qsim' containing the output of the runoff of the downstream subcatchment")
}
OutputsModel <- QcontribDown
OutputsModel$QsimDown <- OutputsModel$Qsim
if (length(QcontribDown$Qsim) != length(RunOptions$IndPeriod_Run)) {
stop("Time series Qsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
}
if (!is.null(QcontribDown$WarmUpQsim) &&
length(QcontribDown$WarmUpQsim) != length(RunOptions$IndPeriod_WarmUp) &&
RunOptions$IndPeriod_WarmUp != 0L) {
stop("Time series WarmUpQsim in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_WarmUp'")
}
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
OutputsModel <- list()
class(OutputsModel) <- c("OutputsModel", class(OutputsModel))
OutputsModel$QsimDown <- QcontribDown
if (length(QcontribDown) != length(RunOptions$IndPeriod_Run)) {
stop("'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
}
} else {
stop("'QcontribDown' must be a numeric vector or a 'OutputsModel' object")
}
if (length(OutputsModel$QsimDown) != length(RunOptions$IndPeriod_Run)) {
stop("Time series in 'QcontribDown' should have the same lenght as 'RunOptions$IndPeriod_Run'")
# data set up
NbUpBasins <- length(InputsModel$LengthHydro)
if (identical(RunOptions$IndPeriod_WarmUp, 0L)) {
RunOptions$IndPeriod_WarmUp <- NULL
}
IndPeriod1 <- c(RunOptions$IndPeriod_WarmUp, RunOptions$IndPeriod_Run)
LInputSeries <- as.integer(length(IndPeriod1))
IndPeriod2 <- (length(RunOptions$IndPeriod_WarmUp)+1):LInputSeries
if (inherits(InputsModel, "hourly")) {
TimeStep <- 60 * 60
} else if (inherits(InputsModel, "daily")) {
TimeStep <- 60 * 60 * 24
} else {
stop("'InputsModel' should be of class \"daily\" or \"hourly\"")
if (inherits(QcontribDown, "OutputsModel")) {
OutputsModel <- QcontribDown
if (is.null(OutputsModel$WarmUpQsim)) {
OutputsModel$WarmUpQsim <- rep(NA, length(RunOptions$IndPeriod_WarmUp))
}
QsimDown <- c(OutputsModel$WarmUpQsim, OutputsModel$Qsim)
} else if (is.vector(QcontribDown) && is.numeric(QcontribDown)) {
OutputsModel <- list()
class(OutputsModel) <- c("OutputsModel", class(RunOptions)[-1])
QsimDown <- c(rep(NA, length(RunOptions$IndPeriod_WarmUp)),
QcontribDown)
}
# propagation time from upstream meshes to outlet
PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / TimeStep
## propagation time from upstream meshes to outlet
PT <- InputsModel$LengthHydro * 1e3 / Param[1L] / RunOptions$FeatFUN_MOD$TimeStep
HUTRANS <- rbind(1 - (PT - floor(PT)), PT - floor(PT))
NbUpBasins <- length(InputsModel$LengthHydro)
LengthTs <- length(OutputsModel$QsimDown)
OutputsModel$Qsim_m3 <- OutputsModel$QsimDown * InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
## set up initial states
IniSD <- RunOptions$IniStates[grep("SD", names(RunOptions$IniStates))]
if (length(IniSD) > 0) {
if (sum(floor(PT)) + NbUpBasins != length(IniSD)) {
......@@ -66,15 +80,15 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
if (x > 1) {
iStart <- iStart + sum(floor(PT[1:x - 1]) + 1)
}
IniSD[iStart:(iStart + PT[x])]
as.vector(IniSD[iStart:(iStart + PT[x])])
})
} else {
IniStates <- lapply(
seq_len(NbUpBasins),
function(iUpBasins) {
iWarmUp <- seq.int(
from = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)] - floor(PT[iUpBasins])),
to = max(1, RunOptions$IndPeriod_WarmUp[length(RunOptions$IndPeriod_WarmUp)])
from = max(1, IndPeriod1[1] - floor(PT[iUpBasins]) - 1),
to = max(1, IndPeriod1[1] - 1)
)
ini <- InputsModel$Qupstream[iWarmUp, iUpBasins]
if (length(ini) != floor(PT[iUpBasins] + 1)) {
......@@ -85,19 +99,30 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
}
)
}
# message("Initstates: ", paste(IniStates, collapse = ", "))
# message("IniStates: ", paste(IniStates, collapse = ", "))
## Lag model computation
Qsim_m3 <- QsimDown *
InputsModel$BasinAreas[length(InputsModel$BasinAreas)] * 1e3
for (upstream_basin in seq_len(NbUpBasins)) {
Qupstream <- c(IniStates[[upstream_basin]],
InputsModel$Qupstream[RunOptions$IndPeriod_Run, upstream_basin])
InputsModel$Qupstream[IndPeriod1, upstream_basin])
# message("Qupstream[", upstream_basin, "]: ", paste(Qupstream, collapse = ", "))
OutputsModel$Qsim_m3 <- OutputsModel$Qsim_m3 +
Qupstream[2:(1 + LengthTs)] * HUTRANS[1, upstream_basin] +
Qupstream[1:LengthTs] * HUTRANS[2, upstream_basin]
Qsim_m3 <- Qsim_m3 +
Qupstream[2:(1 + LInputSeries)] * HUTRANS[1, upstream_basin] +
Qupstream[1:LInputSeries] * HUTRANS[2, upstream_basin]
}
## OutputsModel
OutputsModel$Qsim_m3 <- Qsim_m3[IndPeriod2]
if ("Qsim" %in% RunOptions$Outputs_Sim) {
# Convert back Qsim to mm
OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
# message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
}
# Convert back Qsim to mm
OutputsModel$Qsim <- OutputsModel$Qsim_m3 / sum(InputsModel$BasinAreas, na.rm = TRUE) / 1e3
# message("Qsim: ", paste(OutputsModel$Qsim, collapse = ", "))
# Warning for negative flows or NAs only in extended outputs
if (length(RunOptions$Outputs_Sim) > 2) {
......@@ -112,12 +137,20 @@ RunModel_Lag <- function(InputsModel, RunOptions, Param, QcontribDown) {
}
if ("StateEnd" %in% RunOptions$Outputs_Sim) {
OutputsModel$StateEnd$SD <- lapply(seq(NbUpBasins), function(x) {
SD <- lapply(seq(NbUpBasins), function(x) {
lastTS <- RunOptions$IndPeriod_Run[length(RunOptions$IndPeriod_Run)]
InputsModel$Qupstream[(lastTS - floor(PT[x])):lastTS, x]
})
if(is.null(OutputsModel$StateEnd)) {
OutputsModel$StateEnd <- CreateIniStates(RunModel_Lag, InputsModel, SD = SD)
} else {
OutputsModel$StateEnd$SD <- SD
}
# message("StateEnd: ", paste(OutputsModel$StateEnd$SD, collapse = ", "))
}
if ("WarmUpQsim" %in% RunOptions$Outputs_Sim) {
OutputsModel$WarmUpQsim <- Qsim_m3[seq_len(length(RunOptions$IndPeriod_WarmUp))]
}
if ("Param" %in% RunOptions$Outputs_Sim) {
OutputsModel$Param <- c(Param, OutputsModel$Param)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment