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

BREAKING CHANGE: integrate 'tug' as diffusion module

Browse Source 
It is now possible to run a simulation for one iteration in sequential
mode without the use of Rmufits.

According scripts are provided.

refactor: TransportSim renamed to DiffusionModule

refactor: parsing of R input script is now done outside of simulation
modules (except ChemSim)
This commit is contained in:
Max Lübke 2022-10-25 09:06:48 +02:00
parent 6483b9c249
commit 28b59ff6c3
22 changed files with 1467 additions and 486 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
.gitmodules vendored Normal file
View File

@ -0,0 +1,3 @@
[submodule "ext/tug"]
path = ext/tug
url = ../../sec34/tug.git

4
CMakeLists.txt
View File

@ -4,7 +4,7 @@ cmake_minimum_required(VERSION 3.9)
project(POET CXX C)
# specify the C++ standard
set(CMAKE_CXX_STANDARD 14)
set(CMAKE_CXX_STANDARD 17)
set(CMAKE_CXX_STANDARD_REQUIRED True)
set(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE)
@ -23,6 +23,8 @@ add_subdirectory(R_lib)
add_subdirectory(data)
add_subdirectory(app)
add_subdirectory(ext/tug EXCLUDE_FROM_ALL)
option(BUILD_DOC "Build documentation with doxygen" OFF)
if(BUILD_DOC)

8
README.md
View File

@ -3,6 +3,14 @@
Time-stamp: "Last modified 2021-02-08 13:46:00 mluebke"
-->
**Po**tsdamer R**e**active **D**iffusion
# Forked Project
POED is a forked project of POET, which aims to replace the advection transport
by diffusion. The following README is also applicable for this project.
# POET
POET is a coupled reactive transport simulator implementing a parallel

736
R_lib/kin_r_library.R
View File

@ -1,4 +1,4 @@
### Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
### Copyright (C) 2018-2022 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
@ -21,177 +21,243 @@ FileExt <- function (x)
ifelse(pos > -1L, substring(x, pos + 1L), "")
}
master_init <- function(setup) {
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)
} 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 {
}
# NOTE: Removed the mufits part since we will use diffusion with unconditional
# stable BTCS
## now that we know how many iterations we're gonna have, we can
## setup the optional output
#if (local_rank == 0) {
# if (is.null(setup$iter_output)) {
# ## "iter_output" is not specified: store all iterations
# setup$out_save <- seq(1, maxiter)
# msgm("setup$iter_output unspecified, storing all iterations")
# } else if (setup$iter_output == "all") {
# ## "iter_output" is "all": store all iterations
# setup$out_save <- seq(1, maxiter)
# msgm("storing all iterations")
# } else if (is.numeric(setup$iter_output)) {
# msgm("storing iterations:", paste(setup$iter_output, collapse = ", "))
# setup$out_save <- as.integer(setup$iter_output)
# } else if (setup$iter_output == "last") {
# msgm("storing only the last iteration")
# setup$out_save <- maxiter
# } else { ## fallback to "all"
# setup$out_save <- seq(1, maxiter)
# msgm("invalid setup$iter_output: storing all iterations")
# }
#}
setup$iter <- 1
setup$maxiter <- setup$iterations
#setup$timesteps <- setup$iterations
setup$simulation_time <- 0
# setup$dt_differ <- dt_differ
# if (nrow(setup$bound) == 1) {
# boundmatAct <- t(RedModRphree::pH2Act(setup$bound))
# msg("formed correct matrix from setup$bound:")
# print(boundmatAct)
# } else {
# boundmatAct <- RedModRphree::pH2Act(setup$bound)
# }
# setup$boundmatAct <- boundmatAct
return(setup)
}
### This function is only called by the master process. It sets up
### the data structures for the coupled simulations and performs the
### first timestep or "iteration zero", possibly involving the phreeqc
### calculation of the initial state (if it is not already provided as
### matrix) and the first Advection iteration
master_init <- function(setup)
{
msgm("Process with rank 0 reading GRID and MUFITS_sims")
# master_init <- function(setup)
# {
# msgm("Process with rank 0 reading GRID and MUFITS_sims")
#
# ## MDL TODO: actually grid and all the MUFITS snapshots should be
# ## read and stored only by the master process!!
#
# ## 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 {
#
# }
#
# ## to enhance flexibility, gridfile can be now given in SUM or,
# ## already processed, in rds format. We check for extension first.
# gridext <- FileExt(setup$gridfile)
# if (gridext=="SUM") {
# msgm("Generating grid from MUFITS SUM file...")
# mufits_grid <- ReadGrid(setup$gridfile, verbose=verb)
# } else {
# msgm("Reading grid from rds file...")
# mufits_grid <- readRDS(setup$gridfile)
# }
#
# ## load all the snapshots at once. The setup argument "snapshots"
# ## now can be either a *directory* where the .SUM files are stored
# ## or a rds file.
# if (dir.exists(setup$snapshots)) {
# msgm("<snapshots> points to a directory; reading from SUM files in there...")
# mufits_sims <- LoadMufitsSumRes(dir=setup$snapshots, verbose=verb)
# } else {
# msgm("<snapshots> points to a file. Reading as rds...")
# mufits_sims <- readRDS(setup$snapshots)
# }
#
# ## Since this function is evaluated by the R process called from
# ## the C++ program, we need to make available all these variables
# ## to the R parent frame!
# assign("mufits_grid", mufits_grid, pos=parent.frame())
# assign("mufits_sims", mufits_sims, pos=parent.frame())
# ## cat(local_rank, "assignement complete\n")
#
#
# ## we calculate the *coupling* iterations
# nstep <- length(mufits_sims)
# msgm("using", nstep,"flow snapshots")
# ## cat(local_rank, "nstep:", nstep, "\n")
# ## cat(local_rank, "names:", paste(names(mufits_sims), collate=","), "\n")
#
# ## we have these snapshots; we output the results of the coupling
# ## after each timesteps
# timesteps <- diff(sapply(mufits_sims, function(x) {return(x$info)}))
# ## cat(local_rank, "timesteps:", paste0(timesteps, collate=" "), "\n")
#
# dt_differ <- abs(max(timesteps) - min(timesteps)) != 0
#
# maxiter <- length(timesteps)
#
#
# ## steady state after last flow snapshot? It is controlled by
# ## "setup$prolong", containing an integer which represents the
# ## number of further iterations
# setup$steady <- FALSE
# if ("prolong" %in% names(setup)) {
#
# last_dt <- timesteps[length(timesteps)]
# timesteps <- c(timesteps, rep(last_dt, setup$prolong))
# msgm("simulation prolonged for", setup$prolong, "additional iterations")
# ## we set this flag to TRUE to check afterwards which snapshot we need to use at each iteration
# setup$steady <- TRUE
# setup$last_snapshot <- maxiter
# maxiter <- maxiter + setup$prolong
#
# }
#
# ## now that we know how many iterations we're gonna have, we can
# ## setup the optional output
# if (local_rank==0) {
# if (is.null(setup$iter_output)) {
# ## "iter_output" is not specified: store all iterations
# setup$out_save <- seq(1,maxiter)
# msgm("setup$iter_output unspecified, storing all iterations")
# } else if (setup$iter_output=="all") {
# ## "iter_output" is "all": store all iterations
# setup$out_save <- seq(1,maxiter)
# msgm("storing all iterations")
# } else if (is.numeric(setup$iter_output)) {
# msgm("storing iterations:", paste(setup$iter_output, collapse=", "))
# setup$out_save <- as.integer(setup$iter_output)
# } else if (setup$iter_output=="last") {
# msgm("storing only the last iteration")
# setup$out_save <- maxiter
# } else {## fallback to "all"
# setup$out_save <- seq(1,maxiter)
# msgm("invalid setup$iter_output: storing all iterations")
# }
# }
#
# setup$iter <- 1
# setup$maxiter <- maxiter
# setup$timesteps <- timesteps
# setup$simulation_time <- 0
# setup$dt_differ <- dt_differ
#
# if (nrow(setup$bound)==1) {
# boundmatAct <- t(RedModRphree::pH2Act(setup$bound))
# msg("formed correct matrix from setup$bound:")
# print(boundmatAct)
# } else {
# boundmatAct <- RedModRphree::pH2Act(setup$bound)
# }
#
# setup$boundmatAct <- boundmatAct
# return(setup)
# }
## MDL TODO: actually grid and all the MUFITS snapshots should be
## read and stored only by the master process!!
## 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 {
}
## to enhance flexibility, gridfile can be now given in SUM or,
## already processed, in rds format. We check for extension first.
gridext <- FileExt(setup$gridfile)
if (gridext=="SUM") {
msgm("Generating grid from MUFITS SUM file...")
mufits_grid <- ReadGrid(setup$gridfile, verbose=verb)
} else {
msgm("Reading grid from rds file...")
mufits_grid <- readRDS(setup$gridfile)
}
## load all the snapshots at once. The setup argument "snapshots"
## now can be either a *directory* where the .SUM files are stored
## or a rds file.
if (dir.exists(setup$snapshots)) {
msgm("<snapshots> points to a directory; reading from SUM files in there...")
mufits_sims <- LoadMufitsSumRes(dir=setup$snapshots, verbose=verb)
} else {
msgm("<snapshots> points to a file. Reading as rds...")
mufits_sims <- readRDS(setup$snapshots)
}
## Since this function is evaluated by the R process called from
## the C++ program, we need to make available all these variables
## to the R parent frame!
assign("mufits_grid", mufits_grid, pos=parent.frame())
assign("mufits_sims", mufits_sims, pos=parent.frame())
## cat(local_rank, "assignement complete\n")
## we calculate the *coupling* iterations
nstep <- length(mufits_sims)
msgm("using", nstep,"flow snapshots")
## cat(local_rank, "nstep:", nstep, "\n")
## cat(local_rank, "names:", paste(names(mufits_sims), collate=","), "\n")
## we have these snapshots; we output the results of the coupling
## after each timesteps
timesteps <- diff(sapply(mufits_sims, function(x) {return(x$info)}))
## cat(local_rank, "timesteps:", paste0(timesteps, collate=" "), "\n")
dt_differ <- abs(max(timesteps) - min(timesteps)) != 0
maxiter <- length(timesteps)
## steady state after last flow snapshot? It is controlled by
## "setup$prolong", containing an integer which represents the
## number of further iterations
setup$steady <- FALSE
if ("prolong" %in% names(setup)) {
last_dt <- timesteps[length(timesteps)]
timesteps <- c(timesteps, rep(last_dt, setup$prolong))
msgm("simulation prolonged for", setup$prolong, "additional iterations")
## we set this flag to TRUE to check afterwards which snapshot we need to use at each iteration
setup$steady <- TRUE
setup$last_snapshot <- maxiter
maxiter <- maxiter + setup$prolong
}
## now that we know how many iterations we're gonna have, we can
## setup the optional output
if (local_rank==0) {
if (is.null(setup$iter_output)) {
## "iter_output" is not specified: store all iterations
setup$out_save <- seq(1,maxiter)
msgm("setup$iter_output unspecified, storing all iterations")
} else if (setup$iter_output=="all") {
## "iter_output" is "all": store all iterations
setup$out_save <- seq(1,maxiter)
msgm("storing all iterations")
} else if (is.numeric(setup$iter_output)) {
msgm("storing iterations:", paste(setup$iter_output, collapse=", "))
setup$out_save <- as.integer(setup$iter_output)
} else if (setup$iter_output=="last") {
msgm("storing only the last iteration")
setup$out_save <- maxiter
} else {## fallback to "all"
setup$out_save <- seq(1,maxiter)
msgm("invalid setup$iter_output: storing all iterations")
}
}
setup$iter <- 1
setup$maxiter <- maxiter
setup$timesteps <- timesteps
setup$simulation_time <- 0
setup$dt_differ <- dt_differ
if (nrow(setup$bound)==1) {
boundmatAct <- t(RedModRphree::pH2Act(setup$bound))
msg("formed correct matrix from setup$bound:")
print(boundmatAct)
} else {
boundmatAct <- RedModRphree::pH2Act(setup$bound)
}
setup$boundmatAct <- boundmatAct
return(setup)
}
## This function is called by all processes.
## Since the worker processes need state_C in worker.cpp -> worker_function()
## even the worker need to run this code.
## TODO: worker.cpp -> worker_function() --> state_C is only needed to obtain headers
## and size of dataframe ... if this will be adjusted, this function will be
## unnecessary for workers
init_chemistry <- function(setup)
{
## setup the chemistry
if (!is.matrix(setup$initsim)) {
msgm("initial state defined through PHREEQC simulation, assuming homogeneous medium!")
tmpfirstrun <- RunPQC(setup$initsim, second=TRUE)
## if (local_rank==0){
## print(tmpfirstrun)
## }
remcolnames <- colnames(tmpfirstrun)
if (nrow(tmpfirstrun) > 1) {
## msgm("Iter 0 selecting second row")
firstrun <- matrix(tmpfirstrun[2,], nrow=1, byrow=TRUE)
colnames(firstrun) <- remcolnames
} else {
firstrun <- tmpfirstrun
}
state_C <- matrix(rep(firstrun,setup$n), byrow=TRUE, ncol=length(setup$prop))
colnames(state_C) <- colnames(firstrun)
## if (local_rank==0)
## saveRDS(state_C, "initstate.rds")
} else {
msgm("given initial state")
state_C <- setup$initsim
}
## save state_T and state_C
setup$state_C <- state_C
return(setup)
}
# ## This function is called by all processes.
# ## Since the worker processes need state_C in worker.cpp -> worker_function()
# ## even the worker need to run this code.
# ## TODO: worker.cpp -> worker_function() --> state_C is only needed to obtain headers
# ## and size of dataframe ... if this will be adjusted, this function will be
# ## unnecessary for workers
# init_chemistry <- function(setup)
# {
# ## setup the chemistry
# if (!is.matrix(setup$initsim)) {
# msgm("initial state defined through PHREEQC simulation, assuming homogeneous medium!")
# tmpfirstrun <- RunPQC(setup$initsim, second=TRUE)
#
# ## if (local_rank==0){
# ## print(tmpfirstrun)
# ## }
# remcolnames <- colnames(tmpfirstrun)
# if (nrow(tmpfirstrun) > 1) {
# ## msgm("Iter 0 selecting second row")
# firstrun <- matrix(tmpfirstrun[2,], nrow=1, byrow=TRUE)
# colnames(firstrun) <- remcolnames
# } else {
# firstrun <- tmpfirstrun
# }
# state_C <- matrix(rep(firstrun,setup$n), byrow=TRUE, ncol=length(setup$prop))
# colnames(state_C) <- colnames(firstrun)
#
# ## if (local_rank==0)
# ## saveRDS(state_C, "initstate.rds")
# } else {
# msgm("given initial state")
# state_C <- setup$initsim
# }
# ## save state_T and state_C
# setup$state_C <- state_C
#
# return(setup)
# }
## relic code, may be useful in future
finit <- function(setup)
@ -321,8 +387,8 @@ master_iteration_end <- function(setup) {
iter <- setup$iter
## MDL Write on disk state_T and state_C after every iteration
## comprised in setup$out_save
if (store_result) {
if (iter %in% setup$out_save) {
#if (store_result) {
# if (iter %in% setup$out_save) {
nameout <- paste0(fileout, '/iter_', sprintf('%03d', iter), '.rds')
info <- list(tr_req_dt = as.integer(setup$requested_dt),
tr_allow_dt = setup$allowed_dt,
@ -332,132 +398,132 @@ master_iteration_end <- function(setup) {
simtime=as.integer(setup$simulation_time),
tr_info=info), file=nameout)
msgm("results stored in <", nameout, ">")
}
}
# }
#}
msgm("done iteration", iter, "/", setup$maxiter)
setup$iter <- setup$iter + 1
return(setup)
}
## master: compute advection
master_advection <- function(setup) {
msgm("requested dt=", setup$requested_dt)
concmat <- RedModRphree::pH2Act(setup$state_C)
inner_timesteps <- setup$inner_timesteps
Cf <- setup$Cf
iter <- setup$iter
## MDL: not used at the moment, so commenting it out
## excl <- setup$no_transport_conn
## msgm("excl"); print(excl)
immobile <- setup$immobile
boundmat <- setup$boundmatAct
## setting the current flow snapshot - we have to check if
## we are in "prolongation" or not
if (setup$steady) {
if (iter > setup$last_snapshot)
snap <- mufits_sims[[setup$last_snapshot]]
else
snap <- mufits_sims[[iter]]
} else
snap <- mufits_sims[[iter]]
## conc is a matrix with colnames
if (is.matrix(concmat)) {
if (length(immobile)>0)
val <- concmat[,-immobile]
else
val <- concmat
## sort the columns of the matrix matching the names of boundary
sptr <- colnames(val)
inflow <- boundmat[, colnames(boundmat)!="cbound", drop=FALSE]
spin <- colnames(inflow)
ind <- match(spin,sptr)
if (TRUE %in% is.na(ind)) {
msgm("Missing boundary conditions for some species")
val <- cbind(val,matrix(0,ncol=sum(is.na(ind)),nrow=nrow(val) ))
colnames(val) <- c(sptr,spin[is.na(ind)])
}
sptr <- colnames(val)
ind <- match(sptr,spin)
## msgm("ind"); print(ind)
cnew <- val
msgm("Computing transport of ", ncol(val), " species")
## if (local_rank==0) {
## saveRDS(list(conc=newconc, times=times, activeCells=grid$cell$CELLID[-exclude_cell], fluxv=fluxv),
## file=paste0("TranspUpwind_", local_rank, ".RData"))
## }
## cnew[ boundmat[, 1] ] <- boundmat[, 2]
## MDL 20200227: Here was the bug: "excl" refers to
## CONNECTIONS and not GRID_CELLS!!
## cells where we won't update the concentrations: union of
## boundary and no-flow cells
## if (length(excl) == 0)
## exclude_cell <- sort(c(boundmat[, 1]))
## else
## exclude_cell <- sort(c(boundmat[, 1], excl))
exclude_cell <- sort(c(boundmat[, 1]))
## msgm("mufits_grid$cell$CELLID[-exclude_cell]:")
## print(mufits_grid$cell$CELLID[-exclude_cell])
for (i in seq(1, ncol(val))) {
## form a 2-column matrix with cell id and boundary
## concentration for those elements
newboundmat <- boundmat[,c(1, ind[i]+1)]
## vector with the old concentrations
concv <- val[,i]
## apply boundary conditions to the concentration vector
## (they should stay constant but better safe than sorry)
concv[newboundmat[, 1]] <- newboundmat[, 2]
## call the function
cnew[,i] <- CppTransportUpwindIrr(concv = concv,
times = inner_timesteps,
activeCells = mufits_grid$cell$CELLID[-exclude_cell],
fluxv = setup$fluxv,
listconn = mufits_grid$listconn,
porVol = setup$cellporvol)
}
colnames(cnew) <- colnames(val)
if ( length(immobile) > 0) {
res <- cbind(cnew, concmat[,immobile])
} else {
res <- cnew
}
## check for negative values. This SHOULD NOT OCCUR and may be
## the result of numerical dispersion (or bug in my transport
## code!)
if (any(res <0 )) {
rem_neg <- which(res<0, arr.ind=TRUE)
a <- nrow(rem_neg)
res[res < 0 ] <- 0
msgm("-> ", a, "concentrations were negative")
print(rem_neg)
}
## msgm("state_T after iteration", setup$iter, ":")
## print(head(res))
} else {
msgm("state_C at iteration", setup$iter, " is not a Matrix, doing nothing!")
}
## retransform concentrations H+ and e- into pH and pe
state_T <- RedModRphree::Act2pH(res)
setup$state_T <- state_T
return(setup)
}
#master_advection <- function(setup) {
# msgm("requested dt=", setup$requested_dt)
# concmat <- RedModRphree::pH2Act(setup$state_C)
# inner_timesteps <- setup$inner_timesteps
# Cf <- setup$Cf
# iter <- setup$iter
#
# ## MDL: not used at the moment, so commenting it out
# ## excl <- setup$no_transport_conn
# ## msgm("excl"); print(excl)
# immobile <- setup$immobile
# boundmat <- setup$boundmatAct
#
# ## setting the current flow snapshot - we have to check if
# ## we are in "prolongation" or not
# if (setup$steady) {
# if (iter > setup$last_snapshot)
# snap <- mufits_sims[[setup$last_snapshot]]
# else
# snap <- mufits_sims[[iter]]
# } else
# snap <- mufits_sims[[iter]]
#
# ## conc is a matrix with colnames
# if (is.matrix(concmat)) {
# if (length(immobile)>0)
# val <- concmat[,-immobile]
# else
# val <- concmat
#
# ## sort the columns of the matrix matching the names of boundary
# sptr <- colnames(val)
# inflow <- boundmat[, colnames(boundmat)!="cbound", drop=FALSE]
# spin <- colnames(inflow)
# ind <- match(spin,sptr)
# if (TRUE %in% is.na(ind)) {
# msgm("Missing boundary conditions for some species")
# val <- cbind(val,matrix(0,ncol=sum(is.na(ind)),nrow=nrow(val) ))
# colnames(val) <- c(sptr,spin[is.na(ind)])
# }
#
# sptr <- colnames(val)
# ind <- match(sptr,spin)
# ## msgm("ind"); print(ind)
#
# cnew <- val
#
# msgm("Computing transport of ", ncol(val), " species")
# ## if (local_rank==0) {
# ## saveRDS(list(conc=newconc, times=times, activeCells=grid$cell$CELLID[-exclude_cell], fluxv=fluxv),
# ## file=paste0("TranspUpwind_", local_rank, ".RData"))
# ## }
# ## cnew[ boundmat[, 1] ] <- boundmat[, 2]
#
# ## MDL 20200227: Here was the bug: "excl" refers to
# ## CONNECTIONS and not GRID_CELLS!!
#
# ## cells where we won't update the concentrations: union of
# ## boundary and no-flow cells
# ## if (length(excl) == 0)
# ## exclude_cell <- sort(c(boundmat[, 1]))
# ## else
# ## exclude_cell <- sort(c(boundmat[, 1], excl))
# exclude_cell <- sort(c(boundmat[, 1]))
#
# ## msgm("mufits_grid$cell$CELLID[-exclude_cell]:")
# ## print(mufits_grid$cell$CELLID[-exclude_cell])
#
# for (i in seq(1, ncol(val))) {
# ## form a 2-column matrix with cell id and boundary
# ## concentration for those elements
# newboundmat <- boundmat[,c(1, ind[i]+1)]
# ## vector with the old concentrations
# concv <- val[,i]
# ## apply boundary conditions to the concentration vector
# ## (they should stay constant but better safe than sorry)
# concv[newboundmat[, 1]] <- newboundmat[, 2]
#
# ## call the function
# cnew[,i] <- CppTransportUpwindIrr(concv = concv,
# times = inner_timesteps,
# activeCells = mufits_grid$cell$CELLID[-exclude_cell],
# fluxv = setup$fluxv,
# listconn = mufits_grid$listconn,
# porVol = setup$cellporvol)
# }
#
# colnames(cnew) <- colnames(val)
#
# if ( length(immobile) > 0) {
# res <- cbind(cnew, concmat[,immobile])
# } else {
# res <- cnew
# }
#
# ## check for negative values. This SHOULD NOT OCCUR and may be
# ## the result of numerical dispersion (or bug in my transport
# ## code!)
# if (any(res <0 )) {
# rem_neg <- which(res<0, arr.ind=TRUE)
# a <- nrow(rem_neg)
# res[res < 0 ] <- 0
# msgm("-> ", a, "concentrations were negative")
# print(rem_neg)
# }
# ## msgm("state_T after iteration", setup$iter, ":")
# ## print(head(res))
# } else {
# msgm("state_C at iteration", setup$iter, " is not a Matrix, doing nothing!")
# }
#
# ## retransform concentrations H+ and e- into pH and pe
# state_T <- RedModRphree::Act2pH(res)
#
# setup$state_T <- state_T
#
# return(setup)
#}
## function for the workers to compute chemistry through PHREEQC
@ -480,7 +546,7 @@ slave_chemistry <- function(setup, data)
## 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)
@ -507,10 +573,10 @@ slave_chemistry <- function(setup, data)
## }
## form the PHREEQC input script for the current work package
inplist <- splitMultiKin(data=reduced, procs=1, base=base, first=first,
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)
## if (local_rank==1 & iter==1)
## RPhreeWriteInp("FirstInp", inplist)
tmpC <- RunPQC(inplist, procs=1, second=TRUE)
@ -622,56 +688,50 @@ msgm <- function (...)
## 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")
to_store <- vector(mode="list", length=5)
names(to_store) <- c("Sim", "Flow", "Transport", "Chemistry", "DHT")
## read the setup R file, which is sourced in kin.cpp
tmpbuff <- file(filesim, "r")
setupfile <- readLines(tmpbuff)
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$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$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$Transport <- list(
boundary = setup$bound,
Cf = setup$Cf,
prop = setup$prop,
immobile = setup$immobile,
reduce = setup$reduce )
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
}
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)
saveRDS(to_store, file=paste0(fileout,'/setup.rds'))
msgm("initialization stored in ", paste0(fileout,'/setup.rds'))
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"))
}

67
app/poet.cpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -21,9 +21,9 @@
#include <Rcpp.h>
#include <poet/ChemSim.hpp>
#include <poet/Grid.hpp>
#include <poet/TransportSim.hpp>
#include <poet/RRuntime.hpp>
#include <poet/SimParams.hpp>
#include <poet/DiffusionModule.hpp>
#include <cstring>
#include <iostream>
@ -65,6 +65,7 @@ int main(int argc, char *argv[]) {
RRuntime R(argc, argv);
/*Loading Dependencies*/
// TODO: kann raus
std::string r_load_dependencies = "suppressMessages(library(Rmufits));"
"suppressMessages(library(RedModRphree));"
"source('../R_lib/kin_r_library.R');"
@ -84,32 +85,40 @@ int main(int argc, char *argv[]) {
cout << "CPP: R Init (RInside) on process " << world_rank << endl;
bool dt_differ;
if (world_rank == 0) {
// get timestep vector from grid_init function ...
// HACK: we disable master_init and dt_differ propagation here for testing
// purposes
//
// bool dt_differ;
if (world_rank == 0) { // get timestep vector from
// grid_init function ... //
std::string master_init_code = "mysetup <- master_init(setup=setup)";
R.parseEval(master_init_code);
dt_differ = R.parseEval("mysetup$dt_differ");
// ... and broadcast it to every other rank unequal to 0
MPI_Bcast(&dt_differ, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
}
/* workers will only read the setup DataFrame defined by input file */
else {
R.parseEval("mysetup <- setup");
MPI_Bcast(&dt_differ, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
// dt_differ = R.parseEval("mysetup$dt_differ");
// // ... and broadcast it to every other rank unequal to 0
// MPI_Bcast(&dt_differ, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
// }
// /* workers will only read the setup DataFrame defined by input file */
// else {
// R.parseEval("mysetup <- setup");
// MPI_Bcast(&dt_differ, 1, MPI_C_BOOL, 0, MPI_COMM_WORLD);
}
params.setDtDiffer(dt_differ);
// params.setDtDiffer(dt_differ);
// initialize chemistry on all processes
std::string init_chemistry_code = "mysetup <- init_chemistry(setup=mysetup)";
R.parseEval(init_chemistry_code);
// TODO: einlesen einer initial matrix (DataFrame)
// std::string init_chemistry_code = "mysetup <-
// init_chemistry(setup=mysetup)"; R.parseEval(init_chemistry_code);
Grid grid(R);
grid.init();
// TODO: Grid anpassen
params.initVectorParams(R, grid.getCols());
R.parseEvalQ("mysetup <- setup");
Grid grid(R, poet::GridParams(R));
// grid.init_from_R();
params.initVectorParams(R, grid.getSpeciesCount());
// MDL: store all parameters
if (world_rank == 0) {
@ -126,18 +135,24 @@ int main(int argc, char *argv[]) {
/* THIS IS EXECUTED BY THE MASTER */
if (world_rank == 0) {
ChemMaster master(params, R, grid);
TransportSim trans(R);
DiffusionModule diffusion(poet::DiffusionParams(R), grid);
// diffusion.initialize();
sim_start = MPI_Wtime();
/* Iteration Count is dynamic, retrieving value from R (is only needed by
* master for the following loop) */
uint32_t maxiter = R.parseEval("mysetup$maxiter");
uint32_t maxiter = R.parseEval("mysetup$iterations");
/* SIMULATION LOOP */
for (uint32_t iter = 1; iter < maxiter + 1; iter++) {
cout << "CPP: Evaluating next time step" << endl;
R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
// cout << "CPP: Evaluating next time step" << endl;
// R.parseEvalQ("mysetup <- master_iteration_setup(mysetup)");
double dt = Rcpp::as<double>(
R.parseEval("mysetup$timesteps[" + std::to_string(iter) + "]"));
cout << "CPP: Next time step is " << dt << "[s]" << endl;
/* displaying iteration number, with C++ and R iterator */
cout << "CPP: Going through iteration " << iter << endl;
@ -147,11 +162,13 @@ int main(int argc, char *argv[]) {
cout << "CPP: Calling Advection" << endl;
/* run transport */
trans.run();
// TODO: transport to diffusion
diffusion.simulate(dt);
cout << "CPP: Chemistry" << endl;
/* Fallback for sequential execution */
// TODO: use new grid
if (world_size == 1) {
master.ChemSim::run();
}
@ -186,7 +203,7 @@ int main(int argc, char *argv[]) {
R["simtime"] = sim_end - sim_start;
R.parseEvalQ("profiling$simtime <- simtime");
trans.end();
diffusion.end();
if (world_size == 1) {
master.ChemSim::end();

2
data/CMakeLists.txt
View File

@ -1 +1 @@
install(FILES SimDol2D.R DESTINATION data)
install(FILES SimDol2D.R SimDol2D_diffu.R SimDol1D_diffu.R DESTINATION data)

196
data/SimDol1D_diffu.R Normal file
View File

@ -0,0 +1,196 @@
## chemical database
# NOTE: not needed now
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
package = "RedModRphree"
), is.db = TRUE)
phreeqc::phrLoadDatabaseString(db)
## only the directory
# demodir <- system.file("extdata", "demo_rtwithmufits", package="Rmufits")
n <- 10
m <- 10
prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")
grid <- list(
n_cells = c(n),
s_cells = c(n)
)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"water 1",
"pH 9.91 charge",
"pe 4.0",
"C 1.2279E-04",
"Ca 1.2279E-04",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"SELECTED_OUTPUT", "-high_precision true", "-reset false",
"-time", "-soln", "-temperature true", "-water true",
"-pH", "-pe", "-totals C Ca Cl Mg",
"-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
)
#
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "SOLUTION 1",
# "units mol/kgw",
# "temp 25.0", "water 1",
# "pH 9.91 charge",
# "pe 4.0",
# "C 1.2279E-04",
# "Ca 1.2279E-04",
# "Cl 1E-12",
# "Mg 1E-12",
# "PURE 1",
# "O2g -0.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
init_diffu <- c(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pe" = 4,
"pH" = 7
)
alpha_diffu <- c(
"C" = 1E-6,
"Ca" = 1E-6,
"Cl" = 1E-6,
"Mg" = 1E-6,
"pe" = 1E-8,
"pH" = 1E-8
)
vecinj_diffu <- list(
list(
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001,
"pe" = 4,
"pH" = 7
)
)
inner_index <- c(5)
inner_vecinj_index <- rep(1, 3)
vecinj_inner <- cbind(inner_index, inner_vecinj_index)
# constant_boundary <- list(
# n_cell = seq(1, n),
# i_vecinj = rep(1, n)
# )
#
# closed_boundary <- list(
# n_cell = seq(1, n),
# i_vecinj = rep(0, n)
# )
boundary <- list(
"N" = rep(1, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
vecinj_inner = vecinj_inner,
boundary = boundary,
alpha = alpha_diffu
)
# cbound <- seq(from = 0, to = n - 1)
## setup boundary conditions for transport - we have already read the
## GRID with the following code:
## grid <- Rmufits::ReadGrid(paste0(demodir,"/d2ascii.run.GRID.SUM"))
## cbound <- which(grid$cell$ACTNUM == 2)
## dput(cbound)
# cbound <- c(
# 1L, 50L, 100L, 150L, 200L, 250L, 300L, 350L, 400L, 450L, 500L,
# 550L, 600L, 650L, 700L, 750L, 800L, 850L, 900L, 950L, 1000L,
# 1050L, 1100L, 1150L, 1200L, 1250L, 1300L, 1350L, 1400L, 1450L,
# 1500L, 1550L, 1600L, 1650L, 1700L, 1750L, 1800L, 1850L, 1900L,
# 1950L, 2000L, 2050L, 2100L, 2150L, 2200L, 2250L, 2300L, 2350L,
# 2400L, 2450L, 2451L, 2452L, 2453L, 2454L, 2455L, 2456L, 2457L,
# 2458L, 2459L, 2460L, 2461L, 2462L, 2463L, 2464L, 2465L, 2466L,
# 2467L, 2468L, 2469L, 2470L, 2471L, 2472L, 2473L, 2474L, 2475L,
# 2476L, 2477L, 2478L, 2479L, 2480L, 2481L, 2482L, 2483L, 2484L,
# 2485L, 2486L, 2487L, 2488L, 2489L, 2490L, 2491L, 2492L, 2493L,
# 2494L, 2495L, 2496L, 2497L, 2498L, 2499L, 2500L
# )
# boundinit <- matrix(rep(init[-c(7, 8)], length(cbound)), byrow = TRUE, nrow = length(cbound))
# myboundmat <- cbind(cbound, boundinit)
# myboundmat[cbound == 1, c(2:7)] <- vecinj
# colnames(myboundmat) <- c("cbound", names(vecinj))
# TODO: dt and iterations
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
prop = prop,
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
s_grid = c(n),
n_grid = c(n),
dt = 1,
iterations = 1,
timesteps = rep(1, 10)
)

2
data/SimDol2D.R
View File

@ -100,6 +100,8 @@ myboundmat <- cbind(cbound,boundinit)
myboundmat[cbound==1, c(2:7)] <- vecinj
colnames(myboundmat) <- c("cbound", names(vecinj))
# TODO: dt and iterations
setup <- list(n=2500,
bound=myboundmat,
base=base,

190
data/SimDol2D_diffu.R Normal file
View File

@ -0,0 +1,190 @@
#################################################################
## Section 1 ##
## Grid initialization ##
#################################################################
n <- 5
m <- 5
types <- c("scratch", "phreeqc", "rds")
# initsim <- c("SELECTED_OUTPUT", "-high_precision true",
# "-reset false",
# "USER_PUNCH",
# "-head C Ca Cl Mg pH pe O2g Calcite Dolomite",
# "10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"O2g\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
# "SOLUTION 1",
# "units mol/kgw",
# "temp 25.0", "water 1",
# "pH 9.91 charge",
# "pe 4.0",
# "C 1.2279E-04",
# "Ca 1.2279E-04",
# "Cl 1E-12",
# "Mg 1E-12",
# "PURE 1",
# "O2g -0.6788 10.0",
# "Calcite 0.0 2.07E-3",
# "Dolomite 0.0 0.0",
# "END")
# needed if init type is set to "scratch"
# prop <- c("C", "Ca", "Cl", "Mg", "pH", "pe", "O2g", "Calcite", "Dolomite")
init_cell <- list(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pH" = 9.91,
"pe" = 4,
"O2g" = 10,
"Calcite" = 2.07e-4,
"Dolomite" = 0
)
grid <- list(
n_cells = c(n, m),
s_cells = c(1,1),
type = types[1],
init_cell = as.data.frame(init_cell),
props = names(init_cell),
init_script = NULL
)
##################################################################
## Section 2 ##
## Diffusion parameters and boundary conditions ##
##################################################################
init_diffu <- c(
"C" = 1.2279E-4,
"Ca" = 1.2279E-4,
"Cl" = 0,
"Mg" = 0,
"pe" = 4,
"pH" = 7
)
alpha_diffu <- c(
"C" = 1E-4,
"Ca" = 1E-4,
"Cl" = 1E-4,
"Mg" = 1E-4,
"pe" = 1E-4,
"pH" = 1E-4
)
vecinj_diffu <- list(
list(
"C" = 0,
"Ca" = 0,
"Cl" = 0.002,
"Mg" = 0.001,
"pe" = 4,
"pH" = 9.91
)
)
inner_index <- c(5, 15, 25)
inner_vecinj_index <- rep(1, 3)
vecinj_inner <- cbind(inner_index, inner_vecinj_index)
boundary <- list(
"N" = rep(1, n),
"E" = rep(0, n),
"S" = rep(0, n),
"W" = rep(0, n)
)
diffu_list <- names(alpha_diffu)
diffusion <- list(
init = init_diffu,
vecinj = do.call(rbind.data.frame, vecinj_diffu),
#vecinj_inner = vecinj_inner,
vecinj_index = boundary,
alpha = alpha_diffu
)
##################################################################
## Section 3 ##
## Phreeqc simulation ##
##################################################################
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
package = "RedModRphree"
), is.db = TRUE)
phreeqc::phrLoadDatabaseString(db)
# NOTE: This won't be needed in the future either. Could also be done in a. pqi
# file
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"water 1",
"pH 9.91 charge",
"pe 4.0",
"C 1.2279E-04",
"Ca 1.2279E-04",
"Mg 0.001",
"Cl 0.002",
"PURE 1",
"O2g -0.1675 10",
"KINETICS 1",
"-steps 100",
"-step_divide 100",
"-bad_step_max 2000",
"Calcite", "-m 0.000207",
"-parms 0.0032",
"Dolomite",
"-m 0.0",
"-parms 0.00032",
"END"
)
selout <- c(
"SELECTED_OUTPUT", "-high_precision true", "-reset false",
"-time", "-soln", "-temperature true", "-water true",
"-pH", "-pe", "-totals C Ca Cl Mg",
"-kinetic_reactants Calcite Dolomite", "-equilibrium O2g"
)
# TODO: dt and iterations
setup <- list(
# bound = myboundmat,
base = base,
first = selout,
# initsim = initsim,
# Cf = 1,
grid = grid,
diffusion = diffusion,
prop = names(init_cell),
immobile = c(7, 8, 9),
kin = c(8, 9),
ann = list(O2g = -0.1675),
# phase = "FLUX1",
# density = "DEN1",
reduce = FALSE,
# snapshots = demodir, ## directory where we will read MUFITS SUM files
# gridfile = paste0(demodir, "/d2ascii.run.GRID.SUM")
# init = init,
# vecinj = vecinj,
# cinj = c(0,1),
# boundary = boundary,
# injections = FALSE,
s_grid = c(n, m),
n_grid = c(n, m),
dt = 1,
iterations = 1,
timesteps = rep(1, 10)
)
# not needed yet
# signif_vector <- c(7, 7, 7, 7, 7, 7, 7, 5, 5)
# prop_type <- c("act", "act", "act", "act", "logact", "logact", "ignore", "act", "act")

1
ext/tug Submodule

@ -0,0 +1 @@
Subproject commit 0abd4e04ae40724f9ddb62da9f7c6ebd234330bb

14
include/poet/ChemSim.hpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2021 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
@ -22,14 +22,15 @@
#define CHEMSIM_H
#include "DHT_Wrapper.hpp"
#include "Grid.hpp"
#include "RRuntime.hpp"
#include "SimParams.hpp"
#include "Grid.hpp"
#include <cstdint>
#include <mpi.h>
#include <string>
#include <vector>
/** Number of data elements that are kept free at each work package */
#define BUFFER_OFFSET 5
@ -52,7 +53,10 @@ namespace poet {
* containing basic parameters for simulation.
*
*/
constexpr const char *CHEMISTRY_MODULE_NAME = "state_c";
class ChemSim {
public:
/**
* @brief Construct a new ChemSim object
@ -192,6 +196,10 @@ protected:
*
*/
double chem_t = 0.f;
poet::StateMemory *state;
uint32_t n_cells_per_prop;
std::vector<std::string> prop_names;
};
/**

53
include/poet/TransportSim.hpp → include/poet/DiffusionModule.hpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -18,10 +18,17 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#ifndef TRANSPORT_SIM_H
#define TRANSPORT_SIM_H
#ifndef DIFFUSION_MODULE_H
#define DIFFUSION_MODULE_H
#include "RRuntime.hpp"
#include "poet/SimParams.hpp"
#include <array>
#include <bits/stdint-uintn.h>
#include <poet/Grid.hpp>
#include <string>
#include <tug/BoundaryCondition.hpp>
#include <tug/Diffusion.hpp>
#include <vector>
namespace poet {
/**
@ -30,8 +37,11 @@ namespace poet {
* Offers simple methods to run an iteration and end the simulation.
*
*/
class TransportSim {
public:
constexpr const char *DIFFUSION_MODULE_NAME = "state_t";
class DiffusionModule {
public:
/**
* @brief Construct a new TransportSim object
*
@ -39,7 +49,7 @@ class TransportSim {
*
* @param R RRuntime object
*/
TransportSim(RRuntime &R);
DiffusionModule(poet::DiffusionParams diffu_args, Grid &grid_);
/**
* @brief Run simulation for one iteration
@ -47,7 +57,7 @@ class TransportSim {
* This will simply call the R function 'master_advection'
*
*/
void run();
void simulate(double dt);
/**
* @brief End simulation
@ -64,12 +74,31 @@ class TransportSim {
*/
double getTransportTime();
private:
private:
/**
* @brief Instance of RRuntime
*
*/
RRuntime &R;
// RRuntime &R;
enum { DIM_1D = 1, DIM_2D };
void initialize(poet::DiffusionParams args);
Grid &grid;
uint8_t dim;
uint32_t prop_count;
tug::diffusion::TugInput diff_input;
std::vector<double> alpha;
std::vector<uint32_t> index_constant_cells;
std::vector<std::string> prop_names;
std::vector<tug::bc::BoundaryCondition> bc_vec;
poet::StateMemory *state;
uint32_t n_cells_per_prop;
/**
* @brief time spent for transport
@ -77,6 +106,6 @@ class TransportSim {
*/
double transport_t = 0.f;
};
} // namespace poet
} // namespace poet
#endif // TRANSPORT_SIM_H
#endif // DIFFUSION_MODULE_H

86
include/poet/Grid.hpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -22,9 +22,29 @@
#define GRID_H
#include "RRuntime.hpp"
#include "poet/SimParams.hpp"
#include <Rcpp.h>
#include <array>
#include <bits/stdint-intn.h>
#include <bits/stdint-uintn.h>
#include <cstdint>
#include <map>
#include <string>
#include <vector>
#define MAX_DIM 2
namespace poet {
enum { GRID_X_DIR, GRID_Y_DIR, GRID_Z_DIR };
using StateMemory = struct s_StateMemory {
std::vector<double> mem;
std::vector<std::string> props;
};
constexpr const char *GRID_MODULE_NAME = "grid_init";
/**
* @brief Class describing the grid
*
@ -35,7 +55,10 @@ namespace poet {
*
*/
class Grid {
public:
private:
using prop_vec = std::vector<std::string>;
public:
/**
* @brief Construct a new Grid object
*
@ -44,8 +67,9 @@ class Grid {
*
* @param R
*/
Grid(RRuntime &R) : R(R){};
Grid(RRuntime &R, poet::GridParams grid_args);
~Grid();
/**
* @brief Init the grid
*
@ -54,21 +78,24 @@ class Grid {
* is done by the R runtime. This may change in the future.
*
*/
void init();
void init_from_R();
/**
* @brief Returns the number of elements for each gridcell
*
* @return unsigned int Number of elements
*/
unsigned int getCols();
auto getGridDimension() -> uint8_t;
auto getTotalCellCount() -> uint32_t;
auto getGridCellsCount(uint8_t direction) -> uint32_t;
auto getDomainSize(uint8_t dircetion) -> uint32_t;
/**
* @brief Returns the number of grid cells
*
* @return unsigned int Number of grid cells
*/
unsigned int getRows();
StateMemory *registerState(std::string module_name,
std::vector<std::string> props);
StateMemory *getStatePointer(std::string module_name);
void deregisterState(std::string module_name);
auto getSpeciesCount() -> uint32_t;
auto getPropNames() -> prop_vec;
auto getSpeciesByName(std::string name,
std::string module_name = poet::GRID_MODULE_NAME)
-> std::vector<double>;
/**
* @brief Shuffle the grid and export it to C memory area
@ -112,24 +139,25 @@ class Grid {
*/
void exportWP(double *buffer);
private:
private:
/**
* @brief Instance of RRuntime
*
*/
RRuntime R;
/**
* @brief Number of columns of grid
*
*/
unsigned int ncol;
// uint32_t species_count;
/**
* @brief Number of rows of grid
*
*/
unsigned int nrow;
std::uint8_t dim;
std::array<std::uint32_t, MAX_DIM> d_spatial;
std::array<std::uint32_t, MAX_DIM> n_cells;
std::map<std::string, StateMemory *> state_register;
// std::map<std::string, std::vector<std::string>> prop_register;
prop_vec prop_names;
std::map<std::string, std::vector<double>> grid_init;
/**
* @brief Get a SkeletonDataFrame
@ -142,5 +170,5 @@ class Grid {
*/
Rcpp::DataFrame getSkeletonDataFrame(unsigned int rows);
};
} // namespace poet
#endif // GRID_H
} // namespace poet
#endif // GRID_H

42
include/poet/SimParams.hpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -21,10 +21,14 @@
#ifndef PARSER_H
#define PARSER_H
#include <bits/stdint-uintn.h>
#include <string>
#include <vector>
#include "RRuntime.hpp"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
#include "Rcpp/DataFrame.h"
#include "argh.hpp" // Argument handler https://github.com/adishavit/argh
#include <Rcpp.h>
// BSD-licenced
/** Return value if no error occured */
@ -70,6 +74,32 @@ typedef struct {
bool store_result;
} t_simparams;
using GridParams = struct s_GridParams {
std::vector<uint32_t> n_cells;
std::vector<double> s_cells;
std::string type;
Rcpp::DataFrame init_df;
std::string init_sim;
std::vector<std::string> props;
s_GridParams(poet::RRuntime &R);
};
using DiffusionParams = struct s_DiffusionParams {
std::vector<std::string> prop_names;
Rcpp::NumericVector alpha;
Rcpp::NumericMatrix vecinj_inner;
Rcpp::DataFrame vecinj;
Rcpp::DataFrame vecinj_index;
s_DiffusionParams(poet::RRuntime &R);
};
/**
* @brief Reads information from program arguments and R runtime
*
@ -81,7 +111,7 @@ typedef struct {
*
*/
class SimParams {
public:
public:
/**
* @brief Construct a new SimParams object
*
@ -205,7 +235,7 @@ class SimParams {
*/
std::string getOutDir();
private:
private:
/**
* @brief Validate program parameters and flags
*
@ -270,5 +300,5 @@ class SimParams {
*/
std::string out_dir;
};
} // namespace poet
#endif // PARSER_H
} // namespace poet
#endif // PARSER_H

2
src/CMakeLists.txt
View File

@ -8,5 +8,5 @@ find_library(CRYPTO_LIBRARY crypto)
add_library(poet_lib ${poet_lib_SRC})
target_include_directories(poet_lib PUBLIC ${PROJECT_SOURCE_DIR}/include)
target_link_libraries(poet_lib PUBLIC
MPI::MPI_C ${MATH_LIBRARY} ${CRYPTO_LIBRARY} RRuntime)
MPI::MPI_C ${MATH_LIBRARY} ${CRYPTO_LIBRARY} RRuntime tug)
target_compile_definitions(poet_lib PUBLIC STRICT_R_HEADERS OMPI_SKIP_MPICXX)

20
src/ChemMaster.cpp
View File

@ -23,6 +23,7 @@
#include <iostream>
#include <poet/ChemSim.hpp>
#include <poet/Grid.hpp>
using namespace poet;
using namespace std;
@ -39,15 +40,15 @@ ChemMaster::ChemMaster(SimParams &params, RRuntime &R_, Grid &grid_)
/* allocate memory */
workerlist = (worker_struct *)calloc(world_size - 1, sizeof(worker_struct));
send_buffer = (double *)calloc((wp_size * (grid.getCols())) + BUFFER_OFFSET,
sizeof(double));
mpi_buffer =
(double *)calloc(grid.getRows() * grid.getCols(), sizeof(double));
send_buffer = (double *)calloc(
(wp_size * (grid.getSpeciesCount())) + BUFFER_OFFSET, sizeof(double));
mpi_buffer = (double *)calloc(grid.getGridCellsCount(GRID_Y_DIR) *
grid.getGridCellsCount(GRID_X_DIR),
sizeof(double));
/* calculate distribution of work packages */
R.parseEvalQ(
"wp_ids <- distribute_work_packages(len=nrow(mysetup$state_C), "
"package_size=work_package_size)");
R.parseEvalQ("wp_ids <- distribute_work_packages(len=length(mysetup$prop), "
"package_size=work_package_size)");
// we only sort once the vector
R.parseEvalQ("ordered_ids <- order(wp_ids)");
@ -175,7 +176,7 @@ void ChemMaster::sendPkgs(int &pkg_to_send, int &count_pkgs,
workerlist[p].send_addr = work_pointer;
/* push work pointer to next work package */
end_of_wp = local_work_package_size * grid.getCols();
end_of_wp = local_work_package_size * grid.getSpeciesCount();
work_pointer = &(work_pointer[end_of_wp]);
// fill send buffer starting with work_package ...
@ -364,7 +365,8 @@ void ChemMaster::end() {
/* do some cleanup */
free(timings);
if (dht_enabled) free(dht_perfs);
if (dht_enabled)
free(dht_perfs);
}
double ChemMaster::getSendTime() { return this->send_t; }

57
src/ChemSim.cpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -18,11 +18,17 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "poet/DiffusionModule.hpp"
#include <poet/ChemSim.hpp>
#include <poet/Grid.hpp>
#include <Rcpp.h>
#include <algorithm>
#include <bits/stdint-uintn.h>
#include <iostream>
#include <string>
#include <vector>
using namespace Rcpp;
using namespace poet;
@ -34,6 +40,22 @@ ChemSim::ChemSim(SimParams &params, RRuntime &R_, Grid &grid_)
this->world_size = tmp.world_size;
this->wp_size = tmp.wp_size;
this->out_dir = params.getOutDir();
this->prop_names = this->grid.getPropNames();
this->n_cells_per_prop = this->grid.getTotalCellCount();
this->state =
this->grid.registerState(poet::CHEMISTRY_MODULE_NAME, this->prop_names);
auto &field = this->state->mem;
field.resize(this->n_cells_per_prop * this->prop_names.size());
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
std::vector<double> prop_vec =
this->grid.getSpeciesByName(this->prop_names[i]);
std::copy(prop_vec.begin(), prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
}
}
void ChemSim::run() {
@ -42,10 +64,43 @@ void ChemSim::run() {
/* start time measuring */
chem_a = MPI_Wtime();
std::vector<double> &field = this->state->mem;
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
try {
std::vector<double> t_prop_vec = this->grid.getSpeciesByName(
this->prop_names[i], poet::DIFFUSION_MODULE_NAME);
std::copy(t_prop_vec.begin(), t_prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
} catch (...) {
continue;
}
}
// HACK: transfer the field into R data structure serving as input for phreeqc
R["TMP_T"] = field;
R.parseEvalQ("mysetup$state_T <- setNames(data.frame(matrix(TMP_T, "
"ncol=length(mysetup$grid$props), nrow=" +
std::to_string(this->n_cells_per_prop) +
")), mysetup$grid$props)");
R.parseEvalQ(
"result <- slave_chemistry(setup=mysetup, data=mysetup$state_T)");
R.parseEvalQ("mysetup <- master_chemistry(setup=mysetup, data=result)");
// HACK: copy R data structure back to C++ field
Rcpp::DataFrame result = R.parseEval("mysetup$state_C");
for (uint32_t i = 0; i < this->prop_names.size(); i++) {
std::vector<double> c_prop =
Rcpp::as<std::vector<double>>(result[this->prop_names[i].c_str()]);
std::copy(c_prop.begin(), c_prop.end(),
field.begin() + (i * this->n_cells_per_prop));
}
/* end time measuring */
chem_b = MPI_Wtime();

8
src/ChemWorker.cpp
View File

@ -41,19 +41,19 @@ ChemWorker::ChemWorker(SimParams &params, RRuntime &R_, Grid &grid_,
this->dht_file = params.getDHTFile();
mpi_buffer = (double *)calloc((wp_size * (grid.getCols())) + BUFFER_OFFSET,
mpi_buffer = (double *)calloc((wp_size * (grid.getSpeciesCount())) + BUFFER_OFFSET,
sizeof(double));
mpi_buffer_results =
(double *)calloc(wp_size * (grid.getCols()), sizeof(double));
(double *)calloc(wp_size * (grid.getSpeciesCount()), sizeof(double));
if (world_rank == 1)
cout << "CPP: Worker: DHT usage is " << (dht_enabled ? "ON" : "OFF")
<< endl;
if (dht_enabled) {
int data_size = grid.getCols() * sizeof(double);
int data_size = grid.getSpeciesCount() * sizeof(double);
int key_size =
grid.getCols() * sizeof(double) + (dt_differ * sizeof(double));
grid.getSpeciesCount() * sizeof(double) + (dt_differ * sizeof(double));
int dht_buckets_per_process =
dht_size_per_process / (1 + data_size + key_size);

188
src/DiffusionModule.cpp Normal file
View File

@ -0,0 +1,188 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2022 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
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "poet/SimParams.hpp"
#include "tug/BoundaryCondition.hpp"
#include "tug/Diffusion.hpp"
#include <Rcpp.h>
#include <algorithm>
#include <cstdint>
#include <poet/DiffusionModule.hpp>
#include <poet/Grid.hpp>
#include <array>
#include <bits/stdint-uintn.h>
#include <cassert>
#include <mpi.h>
#include <string>
#include <vector>
using namespace poet;
constexpr std::array<uint8_t, 4> borders = {
tug::bc::BC_SIDE_LEFT, tug::bc::BC_SIDE_RIGHT, tug::bc::BC_SIDE_TOP,
tug::bc::BC_SIDE_BOTTOM};
inline const char *convert_bc_to_config_file(uint8_t in) {
switch (in) {
case tug::bc::BC_SIDE_TOP:
return "N";
case tug::bc::BC_SIDE_RIGHT:
return "E";
case tug::bc::BC_SIDE_BOTTOM:
return "S";
case tug::bc::BC_SIDE_LEFT:
return "W";
}
return "";
}
DiffusionModule::DiffusionModule(poet::DiffusionParams diffu_args, Grid &grid_)
: grid(grid_) {
this->diff_input.setGridCellN(grid_.getGridCellsCount(GRID_X_DIR),
grid_.getGridCellsCount(GRID_Y_DIR));
this->diff_input.setDomainSize(grid_.getDomainSize(GRID_X_DIR),
grid_.getDomainSize(GRID_Y_DIR));
this->dim = grid_.getGridDimension();
this->n_cells_per_prop = grid_.getTotalCellCount();
this->initialize(diffu_args);
}
void DiffusionModule::initialize(poet::DiffusionParams args) {
// const poet::DiffusionParams args(this->R);
// name of props
this->prop_names = args.prop_names;
this->prop_count = this->prop_names.size();
this->state =
this->grid.registerState(DIFFUSION_MODULE_NAME, this->prop_names);
auto &field = this->state->mem;
// get alphas - we cannot assume alphas are provided in same order as initial
// input
std::vector<double> local_alpha(this->prop_count);
for (uint32_t i = 0; i < this->prop_count; i++) {
local_alpha[i] = args.alpha[this->prop_names[i]];
}
this->alpha = local_alpha;
// initialize field
field.reserve(this->n_cells_per_prop * this->prop_count);
for (uint32_t i = 0; i < this->prop_count; i++) {
std::vector<double> prop_vec =
grid.getSpeciesByName(this->prop_names[i]);
std::copy(prop_vec.begin(), prop_vec.end(),
field.begin() + (i * this->n_cells_per_prop));
if (this->dim == this->DIM_2D) {
tug::bc::BoundaryCondition bc(this->grid.getGridCellsCount(GRID_X_DIR),
this->grid.getGridCellsCount(GRID_Y_DIR));
this->bc_vec.push_back(bc);
} else {
tug::bc::BoundaryCondition bc(this->grid.getGridCellsCount(GRID_X_DIR));
this->bc_vec.push_back(bc);
}
}
// apply inner grid constant cells
// NOTE: opening a scope here for distinguish variable names
if (args.vecinj_inner.rows() != 0){
// get indices of constant grid cells
Rcpp::NumericVector indices_const_cells = args.vecinj_inner(Rcpp::_, 0);
this->index_constant_cells =
Rcpp::as<std::vector<uint32_t>>(indices_const_cells);
// get indices to vecinj for constant cells
Rcpp::NumericVector vecinj_indices = args.vecinj_inner(Rcpp::_, 1);
// apply inner constant cells for every concentration
for (int i = 0; i < this->prop_count; i++) {
std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
for (int j = 0; j < indices_const_cells.size(); j++) {
tug::bc::boundary_condition bc = {tug::bc::BC_TYPE_CONSTANT,
bc_vec[vecinj_indices[j] - 1]};
uint32_t x = this->index_constant_cells[j] %
this->grid.getGridCellsCount(GRID_X_DIR);
uint32_t y = (this->dim == this->DIM_1D
? 0
: this->index_constant_cells[j] /
this->grid.getGridCellsCount(GRID_Y_DIR));
curr_bc.setInnerBC(bc, x, y);
}
}
}
// apply boundary conditions to each ghost node
uint8_t bc_size = (this->dim == this->DIM_1D ? 2 : 4);
for (uint8_t i = 0; i < bc_size; i++) {
const auto &side = borders[i];
std::vector<uint32_t> vecinj_i = Rcpp::as<std::vector<uint32_t>>(
args.vecinj_index[convert_bc_to_config_file(side)]);
for (int i = 0; i < this->prop_count; i++) {
std::vector<double> bc_vec = args.vecinj[this->prop_names[i]];
tug::bc::BoundaryCondition &curr_bc = *(this->bc_vec.begin() + i);
for (int j = 0; j < vecinj_i.size(); j++) {
if (vecinj_i[j] == 0) {
continue;
}
curr_bc(side, j) = {tug::bc::BC_TYPE_CONSTANT, bc_vec[vecinj_i[j] - 1]};
}
}
}
}
void DiffusionModule::simulate(double dt) {
double sim_a_transport, sim_b_transport;
sim_b_transport = MPI_Wtime();
auto *field = this->state->mem.data();
this->diff_input.setTimestep(dt);
for (int i = 0; i < this->prop_count; i++) {
auto *in_field = &field[i * this->n_cells_per_prop];
std::vector<double> in_alpha(this->n_cells_per_prop, this->alpha[i]);
this->diff_input.setBoundaryCondition(this->bc_vec[i]);
if (this->dim == this->DIM_1D) {
tug::diffusion::BTCS_1D(this->diff_input, in_field, in_alpha.data());
} else {
tug::diffusion::ADI_2D(this->diff_input, in_field, in_alpha.data());
}
}
sim_a_transport = MPI_Wtime();
transport_t += sim_a_transport - sim_b_transport;
}
void DiffusionModule::end() {
// R["simtime_transport"] = transport_t;
// R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
}
double DiffusionModule::getTransportTime() { return this->transport_t; }

187
src/Grid.cpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -18,20 +18,193 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "poet/SimParams.hpp"
#include <Rcpp.h>
#include <algorithm>
#include <bits/stdint-intn.h>
#include <bits/stdint-uintn.h>
#include <cassert>
#include <new>
#include <numeric>
#include <poet/Grid.hpp>
#include <stdexcept>
#include <string>
#include <vector>
using namespace poet;
using namespace Rcpp;
void Grid::init() {
R.parseEval("GRID_TMP <- mysetup$state_C");
this->ncol = R.parseEval("ncol(GRID_TMP)");
this->nrow = R.parseEval("nrow(GRID_TMP)");
enum { INIT_SCRATCH, INIT_PHREEQC, INIT_RDS };
static inline int8_t get_type_id(std::string type) {
if (type == "scratch") {
return INIT_SCRATCH;
}
if (type == "phreeqc") {
return INIT_PHREEQC;
}
if (type == "rds") {
return INIT_RDS;
}
return -1;
}
unsigned int Grid::getCols() { return this->ncol; }
Grid::Grid(RRuntime &R, poet::GridParams grid_args) : R(R) {
// get dimension of grid
// this->dim = R.parseEval("length(mysetup$grid$n_cells)");
this->dim = grid_args.n_cells.size();
unsigned int Grid::getRows() { return this->nrow; }
// assert that dim is 1 or 2
assert(this->dim <= MAX_DIM && this->dim > 0);
// also assert that vetor spatial discretization is equal in size of
// dimensions
assert(this->dim == grid_args.s_cells.size());
this->n_cells.fill(0);
this->d_spatial.fill(0);
// store grid dimensions
std::copy(grid_args.n_cells.begin(), grid_args.n_cells.end(),
this->n_cells.begin());
std::copy(grid_args.s_cells.begin(), grid_args.s_cells.end(),
this->d_spatial.begin());
auto init_type = get_type_id(grid_args.type);
assert(init_type >= 0);
uint32_t vec_size = this->getTotalCellCount();
if (init_type != INIT_PHREEQC) {
this->prop_names = grid_args.props;
}
for (auto const &prop : this->prop_names) {
std::vector<double> prop_vec(vec_size);
// size of the vector shall be 1
if (init_type == INIT_SCRATCH) {
double prop_val =
((Rcpp::NumericVector)grid_args.init_df[prop.c_str()])[0];
std::fill(prop_vec.begin(), prop_vec.end(), prop_val);
} else if (init_type == INIT_RDS) {
prop_vec = Rcpp::as<std::vector<double>>(grid_args.init_df[prop.c_str()]);
}
// TODO: define behaviour for phreeqc script input
// NOTE: if there are props defined more than once the first written value
// will be used
this->grid_init.insert(std::pair(prop, prop_vec));
}
};
Grid::~Grid() {
for (auto &[key, val] : this->state_register) {
delete val;
}
}
poet::StateMemory *Grid::registerState(std::string module_name,
std::vector<std::string> props) {
poet::StateMemory *mem = new poet::StateMemory;
mem->props = props;
const auto [it, success] = this->state_register.insert({module_name, mem});
if (!success) {
delete mem;
throw std::bad_alloc();
}
return mem;
}
poet::StateMemory *Grid::getStatePointer(std::string module_name) {
const auto it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
throw std::range_error(
std::string("Requested module " + module_name + " was not found"));
}
return it->second;
}
void Grid::deregisterState(std::string module_name) {
const auto it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
throw std::range_error(
std::string("Requested module " + module_name + " was not found"));
}
auto *p = it->second;
delete p;
this->state_register.erase(it);
}
auto Grid::getDomainSize(uint8_t dircetion) -> uint32_t {
return this->d_spatial.at(dircetion);
}
auto Grid::getGridDimension() -> uint8_t { return this->dim; }
auto Grid::getTotalCellCount() -> uint32_t {
uint32_t sum = 1;
for (auto const &dim : this->n_cells) {
sum *= dim;
}
return sum;
}
auto Grid::getGridCellsCount(uint8_t direction) -> uint32_t {
return this->n_cells.at(direction);
}
auto Grid::getSpeciesCount() -> uint32_t { return this->prop_names.size(); }
auto Grid::getPropNames() -> prop_vec { return this->prop_names; }
auto poet::Grid::getSpeciesByName(std::string name, std::string module_name)
-> std::vector<double> {
// if module name references to initial grid
if (module_name == poet::GRID_MODULE_NAME) {
auto const it = this->grid_init.find(name);
if (it == this->grid_init.end()) {
throw std::range_error(
std::string("Species " + name + " was not found for initial grid"));
}
return it->second;
}
// else find module with name
auto const it = this->state_register.find(module_name);
if (it == this->state_register.end()) {
throw std::range_error(
std::string("Requested module " + module_name + " was not found"));
}
poet::StateMemory const *module_memory = it->second;
auto const &prop_vector = module_memory->props;
auto const find_res = std::find(prop_vector.begin(), prop_vector.end(), name);
if (find_res == prop_vector.end()) {
throw std::range_error(std::string(
"Species " + name + " was not found for module " + module_name));
}
uint32_t prop_index = find_res - prop_vector.begin();
uint32_t begin_vec = prop_index * this->getTotalCellCount(),
end_vec = ((prop_index + 1) * this->getTotalCellCount());
return std::vector<double>(module_memory->mem.begin() + begin_vec,
module_memory->mem.begin() + end_vec);
}
void Grid::shuffleAndExport(double *buffer) {
R.parseEval("tmp <- shuffle_field(mysetup$state_T, ordered_ids)");

43
src/SimParams.cpp
View File

@ -2,7 +2,7 @@
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (GFZ Potsdam)
** Copyright (C) 2018-2022 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
@ -18,16 +18,49 @@
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include "poet/RRuntime.hpp"
#include <bits/stdint-uintn.h>
#include <poet/SimParams.hpp>
#include <Rcpp.h>
#include <iostream>
#include <string>
#include <vector>
using namespace poet;
using namespace std;
using namespace Rcpp;
poet::GridParams::s_GridParams(poet::RRuntime &R) {
this->n_cells =
Rcpp::as<std::vector<uint32_t>>(R.parseEval("mysetup$grid$n_cells"));
this->s_cells =
Rcpp::as<std::vector<double>>(R.parseEval("mysetup$grid$s_cells"));
this->type = Rcpp::as<std::string>(R.parseEval("mysetup$grid$type"));
this->init_df =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$grid$init_cell"));
this->props =
Rcpp::as<std::vector<std::string>>(R.parseEval("mysetup$grid$props"));
// this->init_sim =
// Rcpp::as<std::string>(R.parseEval("mysetup$grid$init_sim"));
}
poet::DiffusionParams::s_DiffusionParams(poet::RRuntime &R) {
this->prop_names = Rcpp::as<std::vector<std::string>>(
R.parseEval("names(mysetup$diffusion$init)"));
this->alpha =
Rcpp::as<Rcpp::NumericVector>(R.parseEval("mysetup$diffusion$alpha"));
if (Rcpp::as<bool>(R.parseEval("exists('mysetup$diffusion$vecinj_inner')"))) {
this->vecinj_inner = Rcpp::as<Rcpp::NumericMatrix>(
R.parseEval("mysetup$diffusion$vecinj_inner"));
}
this->vecinj =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj"));
this->vecinj_index =
Rcpp::as<Rcpp::DataFrame>(R.parseEval("mysetup$diffusion$vecinj_index"));
}
SimParams::SimParams(int world_rank_, int world_size_) {
this->simparams.world_rank = world_rank_;
this->simparams.world_size = world_size_;
@ -169,9 +202,8 @@ void SimParams::initVectorParams(RRuntime &R, int col_count) {
// MDL: new output on signif_vector and prop_type
if (signif_vector_exists) {
cout << "CPP: using problem-specific rounding digits: " << endl;
R.parseEval(
"print(data.frame(prop=prop, type=prop_type, "
"digits=signif_vector))");
R.parseEval("print(data.frame(prop=prop, type=prop_type, "
"digits=signif_vector))");
} else {
cout << "CPP: using DHT default rounding digits = "
<< simparams.dht_significant_digits << endl;
@ -206,7 +238,8 @@ std::list<std::string> SimParams::validateOptions(argh::parser cmdl) {
/* loop over all flags and compare to given flaglist*/
for (auto &flag : cmdl.flags()) {
if (!(flaglist.find(flag) != flaglist.end())) retList.push_back(flag);
if (!(flaglist.find(flag) != flaglist.end()))
retList.push_back(flag);
}
/* and loop also over params and compare to given paramlist */

44
src/TransportSim.cpp
View File

@ -1,44 +0,0 @@
/*
** Copyright (C) 2018-2021 Alexander Lindemann, Max Luebke (University of
** Potsdam)
**
** Copyright (C) 2018-2021 Marco De Lucia (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
** Foundation; either version 2 of the License, or (at your option) any later
** version.
**
** POET is distributed in the hope that it will be useful, but WITHOUT ANY
** WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
** A PARTICULAR PURPOSE. See the GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License along with
** this program; if not, write to the Free Software Foundation, Inc., 51
** Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
*/
#include <poet/TransportSim.hpp>
#include <mpi.h>
using namespace poet;
TransportSim::TransportSim(RRuntime &R_) : R(R_) {}
void TransportSim::run() {
double sim_a_transport, sim_b_transport;
sim_b_transport = MPI_Wtime();
R.parseEvalQ("mysetup <- master_advection(setup=mysetup)");
sim_a_transport = MPI_Wtime();
transport_t += sim_a_transport - sim_b_transport;
}
void TransportSim::end() {
R["simtime_transport"] = transport_t;
R.parseEvalQ("profiling$simtime_transport <- simtime_transport");
}
double TransportSim::getTransportTime() { return this->transport_t; }