Skip to content
Snippets Groups Projects

Integrate RcppBTCS with diffusion

Merged Integrate RcppBTCS with diffusion
1-integrate-rcppbtcs-with-diffusion into master
Merged Marco De Lucia requested to merge 1-integrate-rcppbtcs-with-diffusion into master
Compare
and
13 files
+ 795
− 37
Compare changes
  • Side-by-side
  • Inline
Files
13
R/Rphree_Diffusion.R 0 → 100644
+ 454
− 0
## Functions for RT simulations using through RcppBTCS's diffusion
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2022
### Time-stamp: "Last modified 2022-07-13 16:17:43 delucia"
##' This function is somehow equivalent to
##' \code{phreeqc::phrGetSelectedOutput} but it won't call
##' \code{make.names} to preserve special characters in the column
##' names and it also gives the opportunity of selecting one single
##' element from the selected output list it is called upon.
##' @title Modified wrapper to phreeqc method "getSelOutLst"
##' @param n integer, defaults to 1. List index of the selected output
##' that we want
##' @return a selected output array of dim 2.
##' @author MDL
##' @export
GetSelectedOutput <- function(n=1L) {
sel_outs <- .Call("getSelOutLst", PACKAGE = "phreeqc")
if (length(sel_outs)>1)
msg("selout is a list with more than 1 element, returning number ", n)
return(sel_outs[[n]])
}
##' Cfr \code{demo}
##' @title 1D diffusive reactive transport simulations with surface
##' exchange and sorption
##' @param setup a list with several input parameters
##' @param init optional, initial state of the simulation
##' @param maxtime depending on \code{step}, the simulation time
##' @param step one of c("time","iter","fix_dt"), defines the time
##' stepping
##' @param db the chemical database
##' @param maxsim number of simulations above which parallelization is
##' triggered
##' @param ebreak logical, defaults to FALSE. If TRUE, an early break
##' is invoked
##' @param scheme character, one of "FTCS" and "BTCS"
##' @param procs number of processors to be used in parallel execution
##' @param writeout logical, defaults to FALSE. Should we write the
##' phreeqc output on disk?
##' @param surrogate logical, defaults to FALSE. Should we use
##' surrogates instead of phreeqc?
##' @param surrogate.FUN The function which calls the surrogate
##' @param model a list containing all surrogate models, one element
##' per variable
##' @param baleq data structure containing the balance equations
##' @param tol absolute tolerance on the mass balance. If this is
##' trespassed, phreeqc is called instead
##' @param call_pqc logical, defaults to TRUE: if TRUE, calls phreeqc
##' on the surrogate simulations which trespass the tolerance
##' @return a list containing lots of stuff. Each element of the list
##' represent an iteration and has two elements: $T
##' (concentrations after the transport) and $C, concentrations
##' after the chemistry
##' @author MDL
##' @export
DiffRT <- function(setup, init, maxtime,
step=c("time","iter","fix_dt"), db="phreeqc.dat",
maxsim=1, ebreak=FALSE,
scheme="FTCS") {
require(RcppBTCS)
attach(setup)
on.exit({
msg("breaking up before main loop. Traceback:")
print(traceback())
msg("detaching setup, returning state_C")
detach(setup)
msg(" Bye.")
return(state_C)
})
dx <- len/n # m
dt_cfl <- dx^2/alpha/2
msg("Pure diffusion 1D simulations")
## Find out the number of iterations we're about to calculate
if (step=="fix_dt") {
maxiter <- floor(maxtime/dt)+1
msg("Fixed time step of ", dt, ". ", maxiter, " iterations required to reach enditme of", maxtime)
}
if (step=="time") {
maxiter <- floor(maxtime/dt)+1
msg(maxiter,"iterations required")
}
if (step=="iter") {
maxiter <- maxtime
msg("Simulation will end at approx ",floor(maxiter*dt),"secs")
}
## MDL: add "complete" list. This will be a list of lists. Each
## iteration is a 3-components sublist storing state_T (the input
## to chemistry, after transport step), state_C (the input to
## transport, after chemistry), and the whole PHREEQC output.
saved_complete <- vector(mode="list",length=maxiter+1)
pad <- floor(log10(maxiter+1))+1 ## to properly format the step names
spr_string <- paste0("%0", pad,"d")
names(saved_complete) <- paste0("step_",sprintf(spr_string,seq(0,maxiter)))
msg("Loading phreeqc db... ")
phreeqc::phrLoadDatabaseString(db)
msg("database loaded.")
if (missing(init)) {
msg("missing initial chemical state")
if (!is.character(initsim)) {
msg("using base sim to compute initial state")
runinitsim <- RepSol(base, n, first=first)
runinitsim <- runinitsim[-grep("^END",runinitsim)]
} else {
msg("using provided initial script to compute initial state")
runinitsim <- RepSol(initsim, n, first=first)
runinitsim <- runinitsim[-grep("^END",runinitsim)]
}
## msg("running this PHREEQC script:")
## print(runinitsim)
## append the RUN_CELL and make sure there is no other END statement
phreeqc::phrSetOutputStringsOn(TRUE)
phreeqc::phrSetDumpStringsOn(TRUE)
tmpinp <- c(runinitsim, "RUN_CELLS", paste0("-cells 1-",n), "END",
paste0("SAVE solution 1-", n),
paste0("SAVE exchange 1-", n),
paste0("SAVE surface 1-", n), "END",
"DUMP","-all","END")
phreeqc::phrRunString(tmpinp)
tmpfirstrun <- GetSelectedOutput()
## tmpout <- phreeqc::phrGetOutputStrings()
## dump <- phreeqc::phrGetDumpStrings()
if (nrow(tmpfirstrun)>n) {
state_C <- data.matrix(tail(subset(tmpfirstrun, state=="react"), n=n))
} else {
state_C <- data.matrix(tmpfirstrun)
}
rownames(state_C) <- state_C[,"soln"]
} else {
msg("given initial state; assuming correct colnames")
phreeqc::phrRunString(first)
state_C <- data.matrix(init)
rownames(state_C) <- state_C[,"soln"]
## tmpinp <- ""
## tmpout <- ""
}
index_saved <- 1
saved_complete[[index_saved]] <- list(C=state_C) ## , input=tmpinp ,out=tmpout, dump=dump)
index_saved <- index_saved + 1
msg("first step, Diffusion by ", scheme, "...")
state_T <- DiffusionScheme(state_C, inflow=bound, dx=dx, dt=dt, alpha = alpha,
scheme=scheme,
transported=transported)
msg("Diffusion done...")
## remember that total_h, total_o and cb are automatically assumed
modBlock <- c("SOLUTION_MODIFY 1",
"-temp 13",
"-total_h 111.01243359981",
"-total_o 55.506216800086",
"-cb -3.657928589893e-10",
"-totals",
paste(transported, 0.1, sep=" "))
msg("iteration 0...")
## print(state_T)
newinput <- SolutionModify(modBlock, state_T, transported)
## print(newinput)
tmpinp <- c(first, newinput, "RUN_CELLS", paste0("-cells 1-",n),
paste0("SAVE solution 1-", n),
paste0("SAVE exchange 1-", n),
paste0("SAVE surface 1-", n),
"END")
phreeqc::phrRunString(tmpinp)
state_C <- GetSelectedOutput()
msg("dim(state_C)=", paste(dim(state_C), collapse=", "))
## print(c(newinput, "RUN_CELLS", paste("-cells 1-",n),"END"))
saved_complete[[index_saved]] <- list(T = state_T,
C = state_C)
index_saved <- index_saved + 1
if (ebreak) {
msg("Early break invoked, bye.")
return(saved_complete[[1]])
}
simtime <- 0
iter <- 1
msg("going through iterations now!")
on.exit({
msg("MAIN LOOP interrupted during iteration",iter,", simulation time=",round(simtime, 3),"!!")
msg("returning last calculated chemistry.")
print(traceback())
msg(" Bye.")
saved_complete$currentT <- state_T
saved_complete$currentC <- state_C
saved_complete$curinput <- c(newinput, "RUN_CELLS", paste0("-cells 1-",n),"END")
detach(setup)
return(saved_complete)
})
## start of the iterations
while (iter<maxiter) {
## iterations
## TODO: It seems to me that iteration numberins is off by
## one, so that the first run on Rphree is actually one and
## the next ones end
simtime <- simtime + dt
iter <- iter + 1
msg("starting iteration", iter, "/", maxiter,
" - to reach simulation time ", SimTime(simtime))
## transform the pH/pe back into activities
state_T <- DiffusionScheme(state_C, inflow=bound, dx=dx, dt=dt, alpha=alpha,
scheme=scheme, transported=transported)
## new input
newinput <- SolutionModify(modBlock, state_T, transported)
tmpinp <- c(newinput, "RUN_CELLS", paste0("-cells 1-",n),
paste0("SAVE solution 1-", n), paste0("SAVE exchange 1-", n),
paste0("SAVE surface 1-", n),"END")
phreeqc::phrRunString(tmpinp)
## tmpout <- phreeqc::phrGetOutputStrings()
state_C <- GetSelectedOutput()
msg("done iteration", iter, "/", maxiter)
saved_complete[[index_saved]] <- list(T=state_T, C=state_C)
index_saved <- index_saved + 1
}
on.exit()
msg(" detaching setup...")
detach(setup)
msg(" Bye.")
return(saved_complete)
}
##' @title Add a logical block in a PHREEQC input script
##' @param input the input script to be amended
##' @param block a block template. First word on the first line is
##' interpreted as keyword
##' @param before character, defaults to "END". Keyword before which
##' the block will be inserted
##' @return modified input with the added block
##' @author MDL
##' @export
AddBlock <- function(input, block, before="END") {
indsol <- grep("^SOLUTION ", input)
indend <- grep("^END", input)
indbefore <- grep(before,input)
if (length(indsol)!=length(indbefore)) {
stop("Len before != len solution")
}
firstkey <- unlist(strsplit(block[1], " "))[1]
msg("Adding ", length(indsol), "blocks with keyword ", firstkey)
## dimension of new input
newinp <- character(length(input) + length(indsol)*length(block))
## vectorized index magic
offset <- c(0,cumsum(rep(length(block), length(indbefore)-1)))
copy <- mapply(seq, indbefore+offset, indbefore+offset+length(block)-1)
dim(copy) <- NULL
newinp[-copy] <- input
newinp[copy] <- rep(block, length(indsol))
## correct numbering of new block
indkey <- grep(firstkey, newinp)
newinp[indkey] <- paste(firstkey, seq(1, length(indsol)), sep = " ")
return(newinp)
}
##' @title Form a "SOLUTION_MODIFY" input
##' @param block the template "SOLUTION_MODIFY" block
##' @param state named matrix with current concentrations
##' @param prop character vector. Names of properties to Distribute()
##' within the SOLUTION_MODIFY block
##' @return modified input
##' @author MDL
##' @export
SolutionModify <- function(block, state, prop) {
n <- nrow(state)
## msg("N: ",n)
tmp <- rep(block, n)
new <- Distribute(tmp, "-total_o", state[,"total_o"])
new <- Distribute(new, "-total_h", state[,"total_h"])
new <- Distribute(new, "-cb", state[,"cb"])
for (i in prop)
new <- Distribute(new, i, state[,i])
## Fix numbering of SOLUTION_MODIFY
indkey <- grep("SOLUTION_MODIFY", new)
new[indkey] <- paste("SOLUTION_MODIFY", rownames(state), sep = " ")
return(new)
}
##' This function takes the current state of a chemical system in form
##' of a matrix, a vector of concentrations representing boundary
##' conditions (injection), the required time step (dt/s), grid
##' spacing (dx/m) and alpha (diffusion coefficient in m^2/s).
##'
##' @title Diffusion interface to BTCS or FTCS scheme
##' @param conc either a vector or a matrix containing the chemical
##' state to be updated
##' @param inflow vector of fixed concentrations at the inlet
##' @param dx the grid spacing in m
##' @param dt the required time step in s
##' @param alpha diffusion coefficient in m^2/s (constant,
##' homogeneous)
##' @param scheme
##' @param transported
##' @return a matrix containing the state after the advection
##' @author MDL
##' @export
DiffusionScheme <- function(conc, inflow, dx, dt, alpha, scheme, transported) {
totrans <- c(transported, "total_o", "total_h", "cb")
if (is.matrix(conc) | is.data.frame(conc)) {
## cnew is the state, without prepended boundary conditions
cnew <- data.matrix(conc)
if (scheme=="FTCS") {
for (sp in totrans) {
tmp <- as.numeric(c(inflow[sp], cnew[, sp]))
cnew[,sp] <- FTCS(tmp, dx, dt, alpha)
}
} else {
for (sp in totrans) {
## msg("Diffusing names(inflow[sp])", names(inflow[sp]))
tmp <- as.numeric(c(inflow[sp], cnew[, sp]))
cnew[,sp] <- BTCS(tmp, dx, dt, alpha)
}
}
if (typeof(cnew)!="double") {
##mode(res) <- "double"
msg("cnew is not 'double'")
}
return(cnew)
} else {
cj <- c(inflow,conc)
if (scheme=="FTCS") {
cnew <- FTCS(cj, dx, dt, alpha)
} else {
cnew <- BTCS(cj, dx, dt, alpha)
}
if (typeof(cnew)!="double") {
msg("cnew is not double")
}
return(cnew)
}
}
##' @title Call \code{RcppBTCS::BTCS1D()}
##' @param field NumericVector of length n+1 with Dirichlet boundary
##' conditions prepended
##' @param dx the grid spacing in m
##' @param dt the required time step in s
##' @param alpha diffusion coefficient in m^2/s (constant,
##' homogeneous)
##' @return updated concentration vector
##' @author MDL
##' @export
BTCS <- function(field, dx, dt, alpha) {
require(RcppBTCS)
n <- length(field)-1
RcppBTCS::BTCS1D(count_x=n,
grid_size_x=dx*n,
field=field[-1],
alpha=alpha,
dt = dt,
dt_divide=1,
bc_left = field[1],
bc_right= -1)
}
##' @title Call \code{RcppBTCS::RcppFTCS()}
##' @param field NumericVector of length n+1 with Dirichlet boundary
##' conditions prepended
##' @param dx the grid spacing in m
##' @param dt the required time step in s
##' @param alpha diffusion coefficient in m^2/s (constant,
##' homogeneous)
##' @return updated concentration vector
##' @author MDL
##' @export
FTCS <- function(field, dx, dt, alpha) {
require(RcppBTCS)
RcppBTCS::RcppFTCS(n=length(field)-1, length=1,
field=field[-1], alpha=alpha,
bc_left = field[1], timestep = dt)
}
##' @title Run a PHREEQC input and retrieve the DUMP
##' @param input the input script
##' @return the output of \code{phreeqc::phrGetDumpStrings()}
##' @author MDL
##' @export
PQCRunAndDump <- function(input) {
if (!is.list(input)) {
if (is.character(input)) {
out <- phreeqc::phrSetDumpStringsOn(TRUE)
phreeqc::phrRunString(input)
out <- phreeqc::phrGetDumpStrings()
} else {
stopmsg("something wrong with the input, dying!")
}
}
return(out)
}
##' @title Nice formatting of time in s
##' @param ts current time in s
##' @return a string with time converted in days or year as
##' appropriated
##' @author MDL
##' @export
SimTime <- function(ts) {
if (ts <= 3600*24)
return(paste(ts, " [s]"))
else if (ts > 3600*24 && ts <= 3600*24*365.25)
return(paste(round(ts/3600/24, 2), " [d]"))
else
return(paste(round(ts/3600/24/365.25,2), " [y]"))
}