max/poet
1
0
Fork 0
mirror of https://git.gfz-potsdam.de/naaice/poet.git synced 2025-12-16 12:54:50 +01:00
Code Issues Packages Projects Releases Wiki Activity

testing staged changes

Browse Source
This commit is contained in:
Marco De Lucia 2023-08-17 14:29:09 +02:00
parent a5c8982b98
commit b961c26a85
2 changed files with 203 additions and 195 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

390
R_lib/kin_r_library.R
View File

@ -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 ### 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 ### 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 ### this program; if not, write to the Free Software Foundation, Inc., 51
### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. ### Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
## Simple function to check file extension. It is needed to check if ## Simple function to check file extension. It is needed to check if
## the GridFile is SUM (MUFITS format) or rds/RData ## the GridFile is SUM (MUFITS format) or rds/RData
FileExt <- function(x) { FileExt <- function(x) {
pos <- regexpr("\\.([[:alnum:]]+)$", x) pos <- regexpr("\\.([[:alnum:]]+)$", x)
ifelse(pos > -1L, substring(x, pos + 1L), "") ifelse(pos > -1L, substring(x, pos + 1L), "")
} }
master_init <- function(setup) { 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 ## Setup the directory where we will store the results
verb <- FALSE verb <- FALSE
if (local_rank == 0) { if (local_rank == 0) {
verb <- TRUE ## verbosity loading MUFITS results verb <- TRUE ## verbosity loading MUFITS results
if (!dir.exists(fileout)) { if (!dir.exists(fileout)) {
dir.create(fileout) dir.create(fileout)
msgm("created directory ", 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 { } else {
msgm("dir ", fileout, " already exists, I will overwrite!")
} }
if (!exists("store_result")) {
msgm("store_result doesn't exist!") setup$iter <- 1
} else { setup$maxiter <- setup$iterations
msgm("store_result is ", store_result) setup$timesteps <- setup$timesteps
setup$simulation_time <- 0
if (is.null(setup[["store_result"]])) {
setup$store_result <- TRUE
} }
} else {
if (setup$store_result) {
} if (is.null(setup[["out_save"]])) {
setup$out_save <- seq(1, setup$iterations)
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)
} }
}
return(setup)
return(setup)
} }
## This function, called only by master, stores on disk the last ## 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 ## function for the workers to compute chemistry through PHREEQC
slave_chemistry <- function(setup, data) { slave_chemistry <- function(setup, data) {
base <- setup$base base <- setup$base
first <- setup$first first <- setup$first
prop <- setup$prop prop <- setup$prop
immobile <- setup$immobile immobile <- setup$immobile
kin <- setup$kin kin <- setup$kin
ann <- setup$ann ann <- setup$ann
iter <- setup$iter iter <- setup$iter
timesteps <- setup$timesteps timesteps <- setup$timesteps
dt <- timesteps[iter] dt <- timesteps[iter]
state_T <- data ## not the global field, but the work-package state_T <- data ## not the global field, but the work-package
## treat special H+/pH, e-/pe cases ## treat special H+/pH, e-/pe cases
state_T <- RedModRphree::Act2pH(state_T) state_T <- RedModRphree::Act2pH(state_T)
## reduction of the problem ## reduction of the problem
if (setup$reduce) { if (setup$reduce) {
reduced <- ReduceStateOmit(state_T, omit = setup$ann) reduced <- ReduceStateOmit(state_T, omit = setup$ann)
} else { } else {
reduced <- state_T reduced <- state_T
} }
## form the PHREEQC input script for the current work package ## form the PHREEQC input script for the current work package
inplist <- SplitMultiKin( inplist <- SplitMultiKin(
data = reduced, procs = 1, base = base, first = first, data = reduced, procs = 1, base = base, first = first,
ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt ann = ann, prop = prop, minerals = immobile, kin = kin, dt = dt
) )
## if (local_rank==1 & iter==1) ## if (local_rank==1 & iter==1)
## RPhreeWriteInp("FirstInp", inplist) ## RPhreeWriteInp("FirstInp", inplist)
tmpC <- RunPQC(inplist, procs = 1, second = TRUE) tmpC <- RunPQC(inplist, procs = 1, second = TRUE)
## recompose after the reduction ## recompose after the reduction
if (setup$reduce) { if (setup$reduce) {
state_C <- RecomposeState(tmpC, reduced) state_C <- RecomposeState(tmpC, reduced)
} else { } else {
state_C <- tmpC state_C <- tmpC
} }
## the next line is needed since we don't need all columns of ## the next line is needed since we don't need all columns of
## PHREEQC output ## PHREEQC output
return(state_C[, prop]) return(state_C[, prop])
} }
## This function, called by master ## This function, called by master
master_chemistry <- function(setup, data) { master_chemistry <- function(setup, data) {
state_T <- setup$state_T state_T <- setup$state_T
msgm(" chemistry iteration", setup$iter) msgm(" chemistry iteration", setup$iter)
## treat special H+/pH, e-/pe cases ## treat special H+/pH, e-/pe cases
state_T <- RedModRphree::Act2pH(state_T) state_T <- RedModRphree::Act2pH(state_T)
## reduction of the problem ## reduction of the problem
if (setup$reduce) { if (setup$reduce) {
reduced <- ReduceStateOmit(state_T, omit = setup$ann) reduced <- ReduceStateOmit(state_T, omit = setup$ann)
} else { } else {
reduced <- state_T reduced <- state_T
} }
### inject data from workers ## inject data from workers
res_C <- data res_C <- data
rownames(res_C) <- NULL rownames(res_C) <- NULL
## print(res_C) ## print(res_C)
if (nrow(res_C) > nrow(reduced)) { if (nrow(res_C) > nrow(reduced)) {
res_C <- res_C[seq(2, nrow(res_C), by = 2), ] res_C <- res_C[seq(2, nrow(res_C), by = 2), ]
} }
## recompose after the reduction ## recompose after the reduction
if (setup$reduce) { if (setup$reduce) {
state_C <- RecomposeState(res_C, reduced) state_C <- RecomposeState(res_C, reduced)
} else { } else {
state_C <- res_C state_C <- res_C
} }
setup$state_C <- state_C setup$state_C <- state_C
setup$reduced <- reduced setup$reduced <- reduced
return(setup) return(setup)
} }
## Adapted version for "reduction" ## Adapted version for "reduction"
ReduceStateOmit <- function(data, omit = NULL, sign = 6) { 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) colnames(red)[now + 1] <- names(omit)
if (is.list(omit)) {
indomi <- match(names(omit), colnames(data))
datao <- data[, -indomi]
} else {
datao <- data
}
datao <- signif(datao, sign) ret <- red[, colnames(data)]
red <- mgcv::uniquecombs(datao) } else {
inds <- attr(red, "index") ret <- red
now <- ncol(red) }
rownames(ret) <- NULL
attr(ret, "index") <- inds
## reattach the omitted column(s) return(ret)
## 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)
} }
@ -221,67 +224,70 @@ ReduceStateOmit <- function(data, omit = NULL, sign = 6) {
## Attach the name of the calling function to the message displayed on ## Attach the name of the calling function to the message displayed on
## R's stdout ## R's stdout
msgm <- function(...) { msgm <- function(...) {
if (local_rank == 0) { if (local_rank == 0) {
fname <- as.list(sys.call(-1))[[1]] fname <- as.list(sys.call(-1))[[1]]
prefix <- paste0("R: ", fname, " ::") prefix <- paste0("R: ", fname, " ::")
cat(paste(prefix, ..., "\n")) cat(paste(prefix, ..., "\n"))
} }
invisible() invisible()
} }
## Function called by master R process to store on disk all relevant ## Function called by master R process to store on disk all relevant
## parameters for the simulation ## parameters for the simulation
StoreSetup <- function(setup) { 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 to_store <- vector(mode = "list", length = 4)
tmpbuff <- file(filesim, "r") ## names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
setupfile <- readLines(tmpbuff) names(to_store) <- c("Sim", "Transport", "DHT", "Cmdline")
close.connection(tmpbuff)
## 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( ## to_store$Chemistry <- list(
# snapshots = setup$snapshots, ## nprocs = n_procs,
# gridfile = setup$gridfile, ## wp_size = work_package_size,
# phase = setup$phase, ## base = setup$base,
# density = setup$density, ## first = setup$first,
# dt_differ = setup$dt_differ, ## init = setup$initsim,
# prolong = setup$prolong, ## db = db,
# maxiter = setup$maxiter, ## kin = setup$kin,
# saved_iter = setup$iter_output, ## ann = setup$ann)
# out_save = setup$out_save )
to_store$Transport <- setup$diffusion if (dht_enabled) {
# to_store$Chemistry <- list( to_store$DHT <- list(
# nprocs = n_procs, enabled = dht_enabled,
# wp_size = work_package_size, log = dht_log
# base = setup$base, ## signif = dht_final_signif,
# first = setup$first, ## proptype = dht_final_proptype
# init = setup$initsim, )
# db = db, } else {
# kin = setup$kin, to_store$DHT <- FALSE
# ann = setup$ann) }
if (dht_enabled) { to_store$Cmdline <- commandArgs(trailingOnly=FALSE)
to_store$DHT <- list( saveRDS(to_store, file = paste0(fileout, "/setup.rds"))
enabled = dht_enabled, msgm("initialization stored in ", paste0(fileout, "/setup.rds"))
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"))
} }
GetWorkPackageSizesVector <- function(n_packages, package_size, len) { GetWorkPackageSizesVector <- function(n_packages, package_size, len) {
ids <- rep(1:n_packages, times=package_size, each = 1)[1:len] ids <- rep(1:n_packages, times=package_size, each = 1)[1:len]
return(as.integer(table(ids))) return(as.integer(table(ids)))
} }

8
src/SimParams.cpp
View File

@ -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 ** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam) ** 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 ** 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 ** 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"}]; // simparams.print_progressbar = cmdl[{"P", "progress"}];
/*Parse DHT arguments*/ /* Parse DHT arguments */
chem_params.use_dht = cmdl["dht"]; chem_params.use_dht = cmdl["dht"];
chem_params.use_interp = cmdl["interp"]; chem_params.use_interp = cmdl["interp"];
// cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n'; // cout << "CPP: DHT is " << ( dht_enabled ? "ON" : "OFF" ) << '\n';