Commit be82c72c authored by Marco De Lucia's avatar Marco De Lucia
Browse files

MDL: some fixes for CRAN submission

parent 7673975e
## Functions for dealing with surrogate simulations
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
### Time-stamp: "Last modified 2021-04-28 15:30:47 delucia"
### Time-stamp: "Last modified 2021-04-28 15:54:48 delucia"
##' Todo
##' 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 spcing
##' (dx/m) and U (Darcy velocity/m^3/s). Note that it is
##' responsibility of user to check that the required dt does not
##' trespass the CFL condition.
##'
##' @title Simple explicit 1D forward-Euler advection, multispecies
##' @param conc matrix containing the chemical state
##' @param conc either a vector or a matrix containing the chemical
##' state to be updated
##' @param inflow vector of concentrations entering the inlet
##' @param dx the grid spacing in m
##' @param dt the required time step in s
##' @param U Darcy velocity in m^3/s
##' @param U Darcy velocity in m^3/s
##' @return a matrix containing the state after the advection
##' @author MDL
##' @export
AdvectionPQC <- function(conc, inflow=rep(0,ncol(conc)), dx, dt, U)
{
if (is.matrix(conc)) {
## msg("conc is matrix")
## CurrConc <<- conc
## actually the first iteration does not have the attributes attached
if ( "immobile" %in% names(attributes(conc)) ) {
immobile <- attr(conc,"immobile")
......@@ -146,8 +150,6 @@ Act2pH <- function(state){
}
##' TODO
##'
##' @title Recompose a chemical state after reduction
##' @param state the state to be recomposed
##' @param reduced the initial reduced state containing the "index" attribute
......@@ -168,13 +170,13 @@ RecomposeState <- function(state, reduced)
return(res)
}
##' This functions requires the package \code{mgcv}
##'
##' TODO
##' Applies \code{mgcv::uniquecombs} to the matrix representing the
##' current chemical state.
##' @title Finds unique problems in a geochemical state
##' @param data the matrix
##' @return a matrix plus an attribute containing the indeces of the
##' compressed matrix
##' @return a matrix with attached attribute containing the indeces of
##' the compressed matrix for later restoring and the "immobile"
##' indices
##' @author MDL
##' @export
ReduceState <- function(data)
......@@ -187,7 +189,7 @@ ReduceState <- function(data)
##' TODO
##' Cfr \code{demo}
##' @title 1D Reactive transport simulations using phreeqc or a
##' surrogate, equilibrium simulations
##' @param setup a list with several input parameters
......@@ -291,11 +293,12 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
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)))
out_inp <<- vector(mode="list",length=length(maxiter))
out_res <<- vector(mode="list",length=length(maxiter))
out_bal <<- vector(mode="list",length=length(maxiter))
out_of_tol <<- vector(mode="list",length=length(maxiter))
utils::globalVariables(c("out_inp","out_res","out_bal"))
out_inp <- vector(mode="list",length=length(maxiter))
out_res <- vector(mode="list",length=length(maxiter))
out_bal <- vector(mode="list",length=length(maxiter))
timing <- matrix(NA,ncol=5,nrow=maxiter)
......@@ -332,7 +335,6 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
attr(state_C,"immobile") <- immobile
## ssC <<- state_C
index_saved <- 1
saved_complete[[index_saved]] <- list(C=Act2pH(state_C))
index_saved <- index_saved + 1
......@@ -372,8 +374,8 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
saved_complete[[index_saved]] <- list(T=Act2pH(state_T), C=Act2pH(state_C))
index_saved <- index_saved + 1
out_inp[[1]] <<- inplist
out_res[[1]] <<- res_C
out_inp[[1]] <- inplist
out_res[[1]] <- res_C
if (ebreak) {
msg("Early break invoked, bye.")
......@@ -426,7 +428,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
Tm <- system.time( tmpsur <- Selected.Surrogate(reduced, model=model) )[3]
bal <- CheckBalance(baleq,reduced,tmpsur)
out_bal[[iter]] <<- bal
out_bal[[iter]] <- bal
totbal <- apply(bal, 1, mae)
nonok <- which(totbal>tol) ### max(tol,tol2))
......@@ -455,7 +457,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
## new input
Ts <- system.time(inplist <- splitMultiFix(data=reduced, procs=procs, base=base, first=first,
prop=prop, minerals=immobile, nmax=maxsim))[3]
out_inp[[iter]] <<- inplist
out_inp[[iter]] <- inplist
Tm <- system.time(res_C <- RunPQC(inplist, procs=procs))[3]
attr(res_C,"immobile") <- immobile
}
......@@ -473,7 +475,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
## ## transform the pH/pe back into activities
## state_C <- pH2Act(state_C)
out_res[[iter]] <<- res_C
out_res[[iter]] <- res_C
timing[iter,] <- c( dim(reduced)[1], procs, Tm, Tr, Ts)
msg("done iteration", iter, "/", maxiter," CPU-time ",round(Tm,3),"[s]")
......@@ -601,10 +603,12 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
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)))
utils::globalVariables(c("out_inp","out_res","out_bal"))
out_inp <<- vector(mode="list",length=length(maxiter))
out_res <<- vector(mode="list",length=length(maxiter))
out_bal <<- vector(mode="list",length=length(maxiter))
out_inp <- vector(mode="list",length=length(maxiter))
out_res <- vector(mode="list",length=length(maxiter))
out_bal <- vector(mode="list",length=length(maxiter))
timing <- matrix(NA,ncol=5,nrow=maxiter)
msg("Loading phreeqc db... ")
......@@ -675,8 +679,8 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
saved_complete[[index_saved]] <- list(T=Act2pH(state_T), C=Act2pH(state_C))
index_saved <- index_saved + 1
out_inp[[1]] <<- inplist
out_res[[1]] <<- res_C
out_inp[[1]] <- inplist
out_res[[1]] <- res_C
if (ebreak) {
msg("Early break invoked, bye.")
......@@ -730,7 +734,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
Tm <- system.time( tmpsur <- Selected.Surrogate(reduced, model=model) )[3]
## print(head(tmpsur))
bal <- CheckBalance(baleq,reduced,tmpsur)
out_bal[[iter]] <<- bal
out_bal[[iter]] <- bal
totbal <- apply(bal, 1, mae)
## cat(":: Surrogate Balance Summary:\n")
......@@ -747,7 +751,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
Ts <- Ts + system.time(inplist <- splitMultiKin(data=reducednonok, procs=procs, base=base, first=first, ann=ann,
prop=prop, minerals=immobile, kin=kin, dt=dt, nmax=maxsim))[3]
out_inp[[iter]] <<- inplist
out_inp[[iter]] <- inplist
if (length(nonok)==1) {
Tm <- Tm + system.time(tmp2 <- RunPQC(inplist, procs=1))[3]
pqcres_C <- tmp2[,prop]
......@@ -765,7 +769,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
## new input
Ts <- system.time(inplist <- splitMultiKin(data=reduced, procs=procs, base=base, first=first, ann=ann,
prop=prop, minerals=immobile, kin=kin, dt=dt, nmax=maxsim))[3]
out_inp[[iter]] <<- inplist
out_inp[[iter]] <- inplist
Tm <- system.time(tmp3 <- RunPQC(inplist, procs=procs))[3]
res_C <- tmp3[,prop]
attr(res_C,"immobile") <- immobile
......@@ -784,7 +788,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
## transform the pH/pe back into activities
state_C <- pH2Act(state_C)
out_res[[iter]] <<- res_C
out_res[[iter]] <- res_C
timing[iter,] <- c( dim(reduced)[1], procs, Tm, Tr, Ts)
msg("done iteration", iter, "/", maxiter," CPU-time ",round(Tm,3),"[s]")
......
## Functions for dealing with surrogate simulations
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
### Time-stamp: "Last modified 2021-04-28 15:32:05 delucia"
### Time-stamp: "Last modified 2021-04-28 15:44:10 delucia"
##' Computes the average of absolute values of a vector
##' @title Average of absolute values
......@@ -94,9 +94,8 @@ ExtractSpecies <- function(flatlist, species) {
##' @param flatlist the list as input to search in
##' @param pphases character vector containing the names of the Pure
##' Phases whose values are to be extracted
##' @param species vector containing the names of the species to
##' search for
##' @return a matrix containing as many columns as specified pphases
##' @export
##' @author MDL
ExtractPphases <- function(flatlist, pphases) {
tmp <- matrix(NA, len <- length(flatlist), speclen <- length(pphases))
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment