diff --git a/R_lib/kin_r_library.R b/R_lib/kin_r_library.R index c759bf9e1..273fe6c4f 100644 --- a/R_lib/kin_r_library.R +++ b/R_lib/kin_r_library.R @@ -1,4 +1,6 @@ -### Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) +## Time-stamp: "Last modified 2023-08-15 11:58:23 delucia" + +### Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) ### ### POET is free software; you can redistribute it and/or modify it under the ### terms of the GNU General Public License as published by the Free Software @@ -13,51 +15,52 @@ ### this program; if not, write to the Free Software Foundation, Inc., 51 ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + ## Simple function to check file extension. It is needed to check if ## the GridFile is SUM (MUFITS format) or rds/RData FileExt <- function(x) { - pos <- regexpr("\\.([[:alnum:]]+)$", x) - ifelse(pos > -1L, substring(x, pos + 1L), "") + pos <- regexpr("\\.([[:alnum:]]+)$", x) + ifelse(pos > -1L, substring(x, pos + 1L), "") } master_init <- function(setup) { - msgm("Process with rank 0 reading GRID properties") + msgm("Process with rank 0 reading GRID properties") - ## Setup the directory where we will store the results - verb <- FALSE - if (local_rank == 0) { - verb <- TRUE ## verbosity loading MUFITS results - if (!dir.exists(fileout)) { - dir.create(fileout) - msgm("created directory ", fileout) + ## Setup the directory where we will store the results + verb <- FALSE + if (local_rank == 0) { + verb <- TRUE ## verbosity loading MUFITS results + if (!dir.exists(fileout)) { + dir.create(fileout) + msgm("created directory ", fileout) + } else { + msgm("dir ", fileout, " already exists, I will overwrite!") + } + if (!exists("store_result")) { + msgm("store_result doesn't exist!") + } else { + msgm("store_result is ", store_result) + } } else { - msgm("dir ", fileout, " already exists, I will overwrite!") + } - if (!exists("store_result")) { - msgm("store_result doesn't exist!") - } else { - msgm("store_result is ", store_result) + + setup$iter <- 1 + setup$maxiter <- setup$iterations + setup$timesteps <- setup$timesteps + setup$simulation_time <- 0 + + if (is.null(setup[["store_result"]])) { + setup$store_result <- TRUE } - } else { - - } - - setup$iter <- 1 - setup$maxiter <- setup$iterations - setup$timesteps <- setup$timesteps - setup$simulation_time <- 0 - - if (is.null(setup[["store_result"]])) { - setup$store_result <- TRUE - } - - if (setup$store_result) { - if (is.null(setup[["out_save"]])) { - setup$out_save <- seq(1, setup$iterations) + + if (setup$store_result) { + if (is.null(setup[["out_save"]])) { + setup$out_save <- seq(1, setup$iterations) + } } - } - - return(setup) + + return(setup) } ## This function, called only by master, stores on disk the last @@ -95,125 +98,125 @@ master_iteration_end <- function(setup) { ## function for the workers to compute chemistry through PHREEQC slave_chemistry <- function(setup, data) { - base <- setup$base - first <- setup$first - prop <- setup$prop - immobile <- setup$immobile - kin <- setup$kin - ann <- setup$ann - - iter <- setup$iter - timesteps <- setup$timesteps - dt <- timesteps[iter] - - state_T <- data ## not the global field, but the work-package - - ## treat special H+/pH, e-/pe cases - state_T <- RedModRphree::Act2pH(state_T) - - ## reduction of the problem - if (setup$reduce) { - reduced <- ReduceStateOmit(state_T, omit = setup$ann) - } else { - reduced <- state_T - } - - ## form the PHREEQC input script for the current work package - inplist <- SplitMultiKin( - data = reduced, procs = 1, base = base, first = first, - ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt - ) - - ## if (local_rank==1 & iter==1) - ## RPhreeWriteInp("FirstInp", inplist) - - tmpC <- RunPQC(inplist, procs = 1, second = TRUE) - - ## recompose after the reduction - if (setup$reduce) { - state_C <- RecomposeState(tmpC, reduced) - } else { - state_C <- tmpC - } - - ## the next line is needed since we don't need all columns of - ## PHREEQC output - return(state_C[, prop]) + base <- setup$base + first <- setup$first + prop <- setup$prop + immobile <- setup$immobile + kin <- setup$kin + ann <- setup$ann + + iter <- setup$iter + timesteps <- setup$timesteps + dt <- timesteps[iter] + + state_T <- data ## not the global field, but the work-package + + ## treat special H+/pH, e-/pe cases + state_T <- RedModRphree::Act2pH(state_T) + + ## reduction of the problem + if (setup$reduce) { + reduced <- ReduceStateOmit(state_T, omit = setup$ann) + } else { + reduced <- state_T + } + + ## form the PHREEQC input script for the current work package + inplist <- SplitMultiKin( + data = reduced, procs = 1, base = base, first = first, + ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt + ) + + ## if (local_rank==1 & iter==1) + ## RPhreeWriteInp("FirstInp", inplist) + + tmpC <- RunPQC(inplist, procs = 1, second = TRUE) + + ## recompose after the reduction + if (setup$reduce) { + state_C <- RecomposeState(tmpC, reduced) + } else { + state_C <- tmpC + } + + ## the next line is needed since we don't need all columns of + ## PHREEQC output + return(state_C[, prop]) } ## This function, called by master master_chemistry <- function(setup, data) { - state_T <- setup$state_T - - msgm(" chemistry iteration", setup$iter) - - ## treat special H+/pH, e-/pe cases - state_T <- RedModRphree::Act2pH(state_T) - - ## reduction of the problem - if (setup$reduce) { - reduced <- ReduceStateOmit(state_T, omit = setup$ann) - } else { - reduced <- state_T - } - - ### inject data from workers - res_C <- data - - rownames(res_C) <- NULL - - ## print(res_C) - - if (nrow(res_C) > nrow(reduced)) { - res_C <- res_C[seq(2, nrow(res_C), by = 2), ] - } - - ## recompose after the reduction - if (setup$reduce) { - state_C <- RecomposeState(res_C, reduced) - } else { - state_C <- res_C - } - - setup$state_C <- state_C - setup$reduced <- reduced - - return(setup) + state_T <- setup$state_T + + msgm(" chemistry iteration", setup$iter) + + ## treat special H+/pH, e-/pe cases + state_T <- RedModRphree::Act2pH(state_T) + + ## reduction of the problem + if (setup$reduce) { + reduced <- ReduceStateOmit(state_T, omit = setup$ann) + } else { + reduced <- state_T + } + + ## inject data from workers + res_C <- data + + rownames(res_C) <- NULL + + ## print(res_C) + + if (nrow(res_C) > nrow(reduced)) { + res_C <- res_C[seq(2, nrow(res_C), by = 2), ] + } + + ## recompose after the reduction + if (setup$reduce) { + state_C <- RecomposeState(res_C, reduced) + } else { + state_C <- res_C + } + + setup$state_C <- state_C + setup$reduced <- reduced + + return(setup) } ## Adapted version for "reduction" ReduceStateOmit <- function(data, omit = NULL, sign = 6) { - require(mgcv) + require(mgcv) + + rem <- colnames(data) + if (is.list(omit)) { + indomi <- match(names(omit), colnames(data)) + datao <- data[, -indomi] + } else { + datao <- data + } + + datao <- signif(datao, sign) + red <- mgcv::uniquecombs(datao) + inds <- attr(red, "index") + now <- ncol(red) + + + ## reattach the omitted column(s) + ## FIXME: control if more than one ann is present + if (is.list(omit)) { + red <- cbind(red, rep(data[1, indomi], nrow(red))) - rem <- colnames(data) - if (is.list(omit)) { - indomi <- match(names(omit), colnames(data)) - datao <- data[, -indomi] - } else { - datao <- data - } + colnames(red)[now + 1] <- names(omit) - datao <- signif(datao, sign) - red <- mgcv::uniquecombs(datao) - inds <- attr(red, "index") - now <- ncol(red) - - - ## reattach the omitted column(s) - ## FIXME: control if more than one ann is present - if (is.list(omit)) { - red <- cbind(red, rep(data[1, indomi], nrow(red))) - - colnames(red)[now + 1] <- names(omit) - - ret <- red[, colnames(data)] - } else { - ret <- red - } - rownames(ret) <- NULL - attr(ret, "index") <- inds - return(ret) + ret <- red[, colnames(data)] + } else { + ret <- red + } + rownames(ret) <- NULL + attr(ret, "index") <- inds + return(ret) } @@ -221,67 +224,70 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) { ## Attach the name of the calling function to the message displayed on ## R's stdout msgm <- function(...) { - if (local_rank == 0) { - fname <- as.list(sys.call(-1))[[1]] - prefix <- paste0("R: ", fname, " ::") - cat(paste(prefix, ..., "\n")) - } - invisible() + if (local_rank == 0) { + fname <- as.list(sys.call(-1))[[1]] + prefix <- paste0("R: ", fname, " ::") + cat(paste(prefix, ..., "\n")) + } + invisible() } ## Function called by master R process to store on disk all relevant ## parameters for the simulation StoreSetup <- function(setup) { - to_store <- vector(mode = "list", length = 3) - # names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") - names(to_store) <- c("Sim", "Transport", "DHT") - ## read the setup R file, which is sourced in kin.cpp - tmpbuff <- file(filesim, "r") - setupfile <- readLines(tmpbuff) - close.connection(tmpbuff) + to_store <- vector(mode = "list", length = 4) + ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT") + names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline") + + ## read the setup R file, which is sourced in kin.cpp + tmpbuff <- file(filesim, "r") + setupfile <- readLines(tmpbuff) + close.connection(tmpbuff) + + to_store$Sim <- setupfile + + ## to_store$Flow <- list( + ## snapshots = setup$snapshots, + ## gridfile = setup$gridfile, + ## phase = setup$phase, + ## density = setup$density, + ## dt_differ = setup$dt_differ, + ## prolong = setup$prolong, + ## maxiter = setup$maxiter, + ## saved_iter = setup$iter_output, + ## out_save = setup$out_save ) - to_store$Sim <- setupfile + to_store$Transport <- setup$diffusion - # to_store$Flow <- list( - # snapshots = setup$snapshots, - # gridfile = setup$gridfile, - # phase = setup$phase, - # density = setup$density, - # dt_differ = setup$dt_differ, - # prolong = setup$prolong, - # maxiter = setup$maxiter, - # saved_iter = setup$iter_output, - # out_save = setup$out_save ) + ## to_store$Chemistry <- list( + ## nprocs = n_procs, + ## wp_size = work_package_size, + ## base = setup$base, + ## first = setup$first, + ## init = setup$initsim, + ## db = db, + ## kin = setup$kin, + ## ann = setup$ann) - to_store$Transport <- setup$diffusion - # to_store$Chemistry <- list( - # nprocs = n_procs, - # wp_size = work_package_size, - # base = setup$base, - # first = setup$first, - # init = setup$initsim, - # db = db, - # kin = setup$kin, - # ann = setup$ann) + if (dht_enabled) { + to_store$DHT <- list( + enabled = dht_enabled, + log = dht_log + ## signif = dht_final_signif, + ## proptype = dht_final_proptype + ) + } else { + to_store$DHT <- FALSE + } - if (dht_enabled) { - to_store$DHT <- list( - enabled = dht_enabled, - log = dht_log - #signif = dht_final_signif, - #proptype = dht_final_proptype - ) - } else { - to_store$DHT <- FALSE - } - - saveRDS(to_store, file = paste0(fileout, "/setup.rds")) - msgm("initialization stored in ", paste0(fileout, "/setup.rds")) + to_store$Cmdline <- commandArgs(trailingOnly=FALSE) + saveRDS(to_store, file = paste0(fileout, "/setup.rds")) + msgm("initialization stored in ", paste0(fileout, "/setup.rds")) } GetWorkPackageSizesVector <- function(n_packages, package_size, len) { - ids <- rep(1:n_packages, times=package_size, each = 1)[1:len] - return(as.integer(table(ids))) + ids <- rep(1:n_packages, times=package_size, each = 1)[1:len] + return(as.integer(table(ids))) } diff --git a/src/SimParams.cpp b/src/SimParams.cpp index fd77dde20..0a81c5b73 100644 --- a/src/SimParams.cpp +++ b/src/SimParams.cpp @@ -1,8 +1,10 @@ +// Time-stamp: "Last modified 2023-08-15 12:12:36 delucia" + /* ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of ** Potsdam) ** -** Copyright (C) 2018-2022 Marco De Lucia, Max Luebke (GFZ Potsdam) +** Copyright (C) 2018-2023 Marco De Lucia, Max Luebke (GFZ Potsdam) ** ** POET is free software; you can redistribute it and/or modify it under the ** terms of the GNU General Public License as published by the Free Software @@ -144,9 +146,9 @@ int SimParams::parseFromCmdl(char *argv[], RInsidePOET &R) { simparams.print_progressbar = cmdl[{"P", "progress"}]; - simparams.print_progressbar = cmdl[{"P", "progress"}]; + // simparams.print_progressbar = cmdl[{"P", "progress"}]; - /*Parse DHT arguments*/ + /* Parse DHT arguments */ chem_params.use_dht = cmdl["dht"]; chem_params.use_interp = cmdl["interp"]; // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';