From 184b3f3671cf1416186add19f682db4574aed4cf Mon Sep 17 00:00:00 2001 From: Dongdong Kong Date: Tue, 7 Mar 2023 11:16:40 +0800 Subject: [PATCH] update parallel example --- DESCRIPTION | 1 + NAMESPACE | 1 + R/tools.R | 14 ++++ R/vic.R | 172 +++++++++++++++++++++++++++--------------------- R/vic_forcing.R | 5 ++ man/vic.Rd | 8 ++- 6 files changed, 124 insertions(+), 77 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index b7bb437..c0f956e 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -42,6 +42,7 @@ Imports: stats, utils, lubridate, + magrittr, Rcpp, foreach Suggests: diff --git a/NAMESPACE b/NAMESPACE index 12af1dd..6ff0d0d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -11,6 +11,7 @@ export(vic) export(vic_param) export(vic_params) export(vic_version) +import(magrittr) importFrom(Rcpp,sourceCpp) importFrom(lubridate,day) importFrom(lubridate,month) diff --git a/R/tools.R b/R/tools.R index 16da14f..d3e176a 100644 --- a/R/tools.R +++ b/R/tools.R @@ -12,3 +12,17 @@ listk <- function(...) { x <- setNames(list(...), vars) return(x) } + + +`%||%` <- function(x, y) { + if (is.null(x)) { + y + } else { + x + } +} + +set_dimnames <- function(x, value) { + dimnames(x) <- value + x +} diff --git a/R/vic.R b/R/vic.R index d65e910..2e82e61 100644 --- a/R/vic.R +++ b/R/vic.R @@ -248,107 +248,103 @@ vic <- function(forcing, soil, veg, output_info = NULL, veglib = NULL, snowband = NULL, lake = NULL, forcing_veg = NULL, - veg_force_types = c('albedo', 'LAI', 'fcanopy'), + veg_force_types = c("albedo", "LAI", "fcanopy"), parall = FALSE) { - tp1 <- proc.time() - if(!is.list(output_info) || length(output_info) == 0) { - output_info <- list( - fluxes = list(timescale = 'day', aggpar = 1, - outvars = c('OUT_PREC', 'OUT_EVAP', 'OUT_RUNOFF', - 'OUT_BASEFLOW', 'OUT_WDEW', 'OUT_SOIL_MOIST') - ), - snow = list(timescale = 'day', aggpar = 1, - outvars = c('OUT_SWE', 'OUT_SNOW_DEPTH', 'OUT_SNOWF', - 'OUT_SNOW_MELT', 'OUT_SNOW_SURF_TEMP', - 'OUT_SNOW_PACK_TEMP') - ) - ) + if (!is.list(output_info) || length(output_info) == 0) { + output_info <- default_output_info } - if(!is.list(output_info[[1]])) - output_info <- list(output = output_info) + if (!is.list(output_info[[1]])) output_info %<>% list(output = .) + output_info <- deal_output_info(output_info) n_outputs <- length(output_info) out_names <- names(output_info) - if(is.null(veglib)) { + if (is.null(veglib)) { veglib <- veglib_LDAS } else { - if(is.data.frame(veglib)) { - veglib <- veglib[, !sapply(veglib, is.character)] - veglib <- as.matrix(veglib) + if (is.data.frame(veglib)) { + veglib <- veglib[, !sapply(veglib, is.character)] %>% as.matrix() } } - if(is.vector(soil)) + if (is.vector(soil)) { soil <- t(soil) + } soil <- as.matrix(soil) + ncell <- nrow(soil) cellid <- soil[, 2] # Check length of forcing data. minfl <- get_forclen() - for(ft in names(forcing)) { + for (ft in names(forcing)) { fl <- nrow(forcing[[ft]]) - if(is.null(nrow(forcing[[ft]]))) fl <- length(forcing[[ft]]) - if(fl < minfl) { - stop(sprintf('Length of forcing data "%s" (%d) is too short for model require (%d).', - ft, fl, minfl)) + if (is.null(nrow(forcing[[ft]]))) fl <- length(forcing[[ft]]) + if (fl < minfl) { + stop(sprintf( + 'Length of forcing data "%s" (%d) is too short for model require (%d).', + ft, fl, minfl + )) } fc <- ncol(forcing[[ft]]) - if(is.null(fc) && ncell > 1 || fc < ncell) + if (is.null(fc) && ncell > 1 || fc < ncell) { stop(sprintf('Columns of forcing data "%s" (%d) are too few for model require. It should no less than number of gridcells (%d).', ft, fc, ncell)) + } } - if(!is.list(veg)) veg <- list(veg) + if (!is.list(veg)) veg <- list(veg) - if(is.vector(lake)) { + if (is.vector(lake)) { lake <- t(lake) - } else if(!is.null(lake)) { - lake <- as.matrix(lake) + } else if (!is.null(lake)) { + lake %<>% as.matrix() } - if(is.vector(snowband)) { - snowband <- t(snowband) - } else if(!is.null(snowband)) { - snowband <- as.matrix(snowband) + if (is.vector(snowband)) { + snowband %<>% t() + } else if (!is.null(snowband)) { + snowband %<>% as.matrix() } forc_types <- c("PREC", "TEMP", "SW", "LW", "PRESS", "VP", "WIND") # The order of forc_types must not be change. forc_lack <- forc_types[!(forc_types %in% names(forcing))] - celev <- getOption('VIC_global_params')[['nlayers']]*4+10 - if("LW" %in% forc_lack) J <- get_J() + celev <- getOption("VIC_global_params")[["nlayers"]] * 4 + 10 + if ("LW" %in% forc_lack) J <- get_J() - if(ncell == 1) forcing <- data.frame(forcing[forc_types]) - - globalopt <- getOption('VIC_global_params') + if (ncell == 1) forcing <- data.frame(forcing[forc_types]) + + globalopt <- getOption("VIC_global_params") `%dof%` <- ifelse(parall, foreach::`%dopar%`, foreach::`%do%`) tp2 <- proc.time() - tmp_out <- foreach::foreach(i=1:ncell) %dof% { + # 需要传递全局变量 + tmp_out <- foreach::foreach(i = 1:ncell) %dof% { + # tmp_out <- list() + # for (i in 1:ncell) { # running(i, 100) forc_veg <- matrix() # default - if(!is.null(forcing_veg)) { + if (!is.null(forcing_veg)) { forc_veg <- forcing_veg[[i]] attr(forc_veg, "types") <- veg_force_types } - iband <- if (is.null(snowband)) -1 else snowband[i, ] - ilake <- if(is.null(lake)) -1 else lake[i, ] + iband <- snowband[i, ] %||% -1 + ilake <- lake[i, ] %||% -1 - if(is.null(nrow(veg[[i]])) || ncol(veg[[i]]) <= 1){ + if (is.null(nrow(veg[[i]])) || ncol(veg[[i]]) <= 1) { iveg <- matrix(nrow = 0, ncol = 8) } else { iveg <- veg[[i]] } - forc <- sapply(forc_types, function(ft){ - if(ft %in% forc_lack) { + forc <- sapply(forc_types, function(ft) { + if (ft %in% forc_lack) { rep(0, minfl) } else { forcing[[ft]][1:minfl, i] @@ -356,41 +352,48 @@ vic <- function(forcing, soil, veg, }) # Estimate forcing data that not supplied. - if("PRESS" %in% forc_lack) { + if ("PRESS" %in% forc_lack) { elev <- soil[i, celev] - forc[,5] <- 101.3*exp(-elev*9.81/(287*(273.15+forc[,2]-0.5*elev*0.0065))) + forc[, 5] <- cal_Pa(forc[, 2], elev) } - if("LW" %in% forc_lack) { + if ("LW" %in% forc_lack) { lat <- soil[i, 3] - forc[, 4] <- cal_lw(forc[,2], forc[,6], forc[,3], lat, J) + # cal_lw(temp, vp, rdsn. lat, J) + forc[, 4] <- cal_lw(forc[, 2], forc[, 6], forc[, 3], lat, J) # temp, vp, rsds, lat, J } - - p <- listk(globalopt, - forc, soil[i, ], iband, iveg, - ilake, forc_veg, veglib, output_info) + + # p <- listk( + # globalopt, + # forc, soil[i, ], iband, iveg, + # ilake, forc_veg, veglib, output_info + # ) # print(str(p)) - iout <- vic_run_cell(globalopt, - forc, soil[i, ], iband, iveg, - ilake, forc_veg, veglib, output_info) + + iout <- vic_run_cell( + globalopt, + forc, soil[i, ], iband, iveg, + ilake, forc_veg, veglib, output_info + ) + # tmp_out[[i]] = iout iout } tp3 <- proc.time() - # Run time in C - #cpp_time <- rbind( + # cpp_time <- rbind( # c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 1]])), 0), # c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 2]])), 0), # c(rowSums(sapply(tmp_out, function(x)x[[n_outputs + 3]])), 0)) out <- list() - for(i in 1:n_outputs) { - out_header <- attr(tmp_out[[1]][[i]], 'header') - out_time <- attr(tmp_out[[1]][[i]], 'time') + + for (i in 1:n_outputs) { + out_header <- attr(tmp_out[[1]][[i]], "header") + out_time <- attr(tmp_out[[1]][[i]], "time") + iout <- lapply(1:length(out_header), function(h) { - itype <- sapply(1:ncell, function(g) tmp_out[[g]][[i]][, h]) - rownames(itype) <- out_time - colnames(itype) <- cellid + itype <- sapply(1:ncell, function(g) tmp_out[[g]][[i]][, h]) %>% + set_dimnames(list(out_time, cellid)) itype }) names(iout) <- out_header @@ -399,14 +402,35 @@ vic <- function(forcing, soil, veg, tp4 <- proc.time() - time_table <- rbind(tp2-tp1, tp3-tp2, tp4-tp3, tp4-tp1)[,1:3] - #time_table <- rbind(time_table, cpp_time) - #rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time', + time_table <- rbind(tp2 - tp1, tp3 - tp2, tp4 - tp3, tp4 - tp1)[, 1:3] %>% + set_rownames(c("init_time", "run_time", "final_time", "total_time")) %>% + round(4) + + # time_table <- rbind(time_table, cpp_time) + # rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time', # 'init_time(cpp)', 'run_time(cpp)', 'final_time(cpp)') - rownames(time_table) <- c('init_time', 'run_time', 'final_time', 'total_time') - out[['timetable']] <- round(time_table, 4) - out[['global_options']] <- globalopt - out[['output_info']] <- output_info - class(out) <- 'vic_output' + out[["timetable"]] <- time_table + out[["global_options"]] <- globalopt + out[["output_info"]] <- output_info + class(out) <- "vic_output" out } + +#' @import magrittr +default_output_info <- list( + fluxes = list( + timescale = "day", aggpar = 1, + outvars = c( + "OUT_PREC", "OUT_EVAP", "OUT_RUNOFF", + "OUT_BASEFLOW", "OUT_WDEW", "OUT_SOIL_MOIST" + ) + ), + snow = list( + timescale = "day", aggpar = 1, + outvars = c( + "OUT_SWE", "OUT_SNOW_DEPTH", "OUT_SNOWF", + "OUT_SNOW_MELT", "OUT_SNOW_SURF_TEMP", + "OUT_SNOW_PACK_TEMP" + ) + ) +) diff --git a/R/vic_forcing.R b/R/vic_forcing.R index 3df1d7c..8516366 100644 --- a/R/vic_forcing.R +++ b/R/vic_forcing.R @@ -20,6 +20,11 @@ get_J <- function() { J } +# air pressure in kPa +cal_Pa <- function(Tair, elev = 0) { + 101.3 * exp(-elev * 9.81 / (287 * (273.15 + Tair - 0.5 * elev * 0.0065))) +} + cal_lw <- function(temp, vp, rsds, lat, J) { # Calc s lat <- lat * pi/360 diff --git a/man/vic.Rd b/man/vic.Rd index 23a0ffe..3c785e6 100644 --- a/man/vic.Rd +++ b/man/vic.Rd @@ -302,10 +302,12 @@ print(outputs) # Example of parallelization \dontrun{ -library(doParallel) -registerDoParallel(cores=4) +# library(doParallel) +# registerDoParallel(cores=4) +Ipaper::InitCluster(4, FORK = FALSE) + outputs <- vic(forcing, soil, veg, snowband = band, parall = TRUE) -stopImplicitCluster() +# stopImplicitCluster() print(outputs) } }