Bumped to 0.3.4

DESCRIPTION

 Package: RedModRphree
 Title: Utilities Leveraging the R Interface to the PHREEQC Geochemical Solver 
-Version: 0.3.3
+Version: 0.3.4
 Authors@R: c(person(given = "Marco",
                     family = "De Lucia",
                     email = "delucia@gfz-potsdam.de",

R/Rphree_Kinetics.R

 ##  Functions for dealing with simulations with kinetics
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
-### Time-stamp: "Last modified 2021-04-28 17:23:36 delucia"
+### Time-stamp: "Last modified 2021-04-29 20:18:12 delucia"
 
 ##' This function runs the generated input buffer (or a list thereof)
 ##' through \code{phreeqc}, which has been already loaded as
@@ -249,16 +249,14 @@ DistributeMatrix <- function(input, prop, values, minerals, SI=FALSE)
     if (n!=ncol(values)) {
         stopmsg("prop names and values matrix of different dimension.")
     }
-    if (is.character(minerals))
-        {
-            qmin <- gsub("\\(","\\\\(",minerals)
-            qmin <- gsub("\\)","\\\\)",qmin)
-            mins <- match(qmin, qprop)
-        }
-    else
+    if (is.character(minerals)) {
+        qmin <- gsub("\\(","\\\\(",minerals)
+        qmin <- gsub("\\)","\\\\)",qmin)
+        mins <- match(qmin, qprop)
+    } else
         mins <- as.numeric(minerals)
     
-    indAq <- seq(1:n)[-mins]
+    indAq <- seq(1,n)[-mins]
     newinp <- input
 
     for (i in indAq) {
@@ -290,7 +288,7 @@ DistributeMatrix <- function(input, prop, values, minerals, SI=FALSE)
 DistributeKinMatrix <- function(input, prop, values, minerals, kin, dt, ann)
 {
     if (!is.matrix(values)&&!is.data.frame(values))
-        stop("I need a matrix to distribute")
+        stopmsg("I need a matrix to distribute")
     ## correct quoting
     qprop <- gsub("\\(","\\\\(",prop)
     qprop <- gsub("\\)","\\\\)",qprop)
@@ -373,12 +371,13 @@ splitMultiFix <- function(data, procs, base, first, prop, minerals, nmax=200, ve
         ntot <- nrow(datamat)
     }
     
-    if ( ntot  ==1) ## try and call normal Distribute (i.e., ##!is.matrix(data) ||
+    if ( ntot==1 ) {
+        ## try and call normal Distribute (i.e., ##!is.matrix(data) ||
         ## just 1 simulation!)
-    {
+    
         if (verbose)
             msg("just 1 simulation, reverting back to normal Distribute")
-        BigInp <- RepSol(base, 1, first=first)
+        BigInp <- RepSol(base, 1)
         namdis <- colnames(data)
         indminer <- minerals
         nammin <- namdis[indminer]
@@ -401,8 +400,8 @@ splitMultiFix <- function(data, procs, base, first, prop, minerals, nmax=200, ve
         totl <- c(nlong,rep(nsim,procs-1))
         breaks <- c(0,cumsum(totl))
         
-        BigInplong <- RepSol(base, nlong, first=first)
-        BigInp <- RepSol(base, nsim, first=first)
+        BigInplong <- RepSol(base, nlong)
+        BigInp <- RepSol(base, nsim)
         
         newlist <- vector(mode="list", length=procs)
         for (i in seq_along(newlist)) {
@@ -412,14 +411,15 @@ splitMultiFix <- function(data, procs, base, first, prop, minerals, nmax=200, ve
         input <- vector(mode="list", length=procs)
         attr(input,"sims") <- breaks
         
-        input[[1]] <- DistributeMatrix(BigInplong, prop, newlist[[1]], minerals=minerals)
+        input[[1]] <- c(first, DistributeMatrix(BigInplong, prop, newlist[[1]], minerals=minerals))
         input[c(2:procs)] <- lapply(newlist[c(2:procs)],
                                     function(x)
-                                        return(DistributeMatrix(BigInp, prop, x, minerals=minerals))
+                                        return(c(first, DistributeMatrix(BigInp, prop, x, minerals=minerals)))
                                     )
     }  else {
-        BigInp <- RepSol(base, ntot, first=first)
+        BigInp <- RepSol(base, ntot)
         input <- DistributeMatrix(input=BigInp, prop=prop, values=data, minerals=minerals)
+        input <- c(first, BigInp)
     }
     return(input)
 }

R/Rphree_ReactTrans.R

 ##  Functions for dealing with surrogate simulations
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
-### Time-stamp: "Last modified 2021-04-28 18:13:39 delucia"
+### Time-stamp: "Last modified 2021-04-29 18:33:18 delucia"
 
 
 ##' This function takes the current state of a chemical system in form
@@ -619,7 +619,6 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
             parallel::clusterExport(cl=ThisRunCluster, varlist=c("model","surrogate.FUN"), envir = environment())
             msg("All workers are setup with the surrogate model.")
         }
-        
     }
     if (missing(init)) {
         msg("missing initial chemical state")

R/Rphree_Utils.R

 ### Utility functions for RedModRphree
 
 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
-### Time-stamp: "Last modified 2021-04-29 17:36:29 delucia"
+### Time-stamp: "Last modified 2021-04-29 20:17:47 delucia"
 
 ##' Replicates an input buffer containing only one SOLUTION, taking
 ##' care of SOLUTION/KINETICS/PURE identifiers. Eventually insert a
@@ -21,42 +21,35 @@ RepSol <- function(sol, n, first=NULL)
     if (length(dum <- grep("^TITLE|^DATABASE",sol)) > 0 )
         sol <- sol[-dum]
     
-    if (is.character(first))
-        {
-            firstend <- grep("^END",sol)[1]
-            newsol <- c(sol[-firstend],first,"END",rep(sol,n-1))
-        } else {
-            newsol <- rep(sol,n)
-        }
+    newsol <- rep(sol,n)
     
-    if (length(linesol <- grep("^SOLUTION",newsol)) != n)
-        {
-            stopmsg("too many or no SOLUTION defined")
-        } else {
-            newsol[linesol] <- paste("SOLUTION",1:n)
-        }
+    if (length(linesol <- grep("^SOLUTION",newsol)) != n) {
+        stopmsg("too many or no SOLUTION defined")
+    } else {
+        newsol[linesol] <- paste("SOLUTION",1:n)
+    }
     
     linepure <- grep("^PURE|^EQUIL",newsol)
-    if (length(linepure) > 0 )
-        {
-            if (length(linepure) != n)
-                {
-                    stopmsg("too many or no PURE defined")
-                } else {
-                    newsol[linepure] <- paste("PURE",1:n)
-                }
+    if (length(linepure) > 0 ) {
+        if (length(linepure) != n) {
+            stopmsg("too many or no PURE defined")
+        } else {
+            newsol[linepure] <- paste("PURE",1:n)
         }
+    }
     
     linekin <- grep("^KINET",newsol)
-    if (length(linekin) > 0)
-        {
-            if (length(linekin) != n)
-                {
-                    stopmsg("too many or no KINET defined")
-                } else {
-                    newsol[linekin] <- paste("KINETICS",1:n)
-                }
+    if (length(linekin) > 0) {
+        if (length(linekin) != n) {
+            stopmsg("too many or no KINET defined")
+        } else {
+            newsol[linekin] <- paste("KINETICS",1:n)
         }
+    }
+
+    if (is.character(first)) 
+        newsol <- c(first,newsol)
+
     return(newsol)  
 }
 
@@ -174,31 +167,30 @@ AddProp <- function(input, name, values, cat, kinpar=NULL, first=NULL)
 ##' @export
 Distribute <- function(input, prop, values, newname=NULL, first=NULL, wholeline=TRUE)
 {
-
     ## correct quoting of parenthesis
     reg <- gsub("\\(","\\\\(",prop)
     reg <- gsub("\\)","\\\\)",reg)
 
     ## coerce values to character
-    n <- length(val <- as.character(values))
-    nsim <- length(proplines <- grep(paste0("^", reg),input))
+    nval <- length(val <- as.character(values))
+    ## check how many lines contain the property we want to distribute
+    nsim <- length(proplines <- grep(paste0("^", reg, " "),input))
     
     if (n != nsim && nsim != 1) {
         stopmsg("Lengths of simulation and values do not agree")
     }
+
     linesave <- NULL
-    
     if (nsim == 1) {
-        ## deal with comments after the properties (i.e. "as HCO3")
         if (!wholeline) {
-            dum <- unlist(strsplit(gsub(' +',' ',grep(paste0("^", reg), input, value=TRUE))," "))
+            ## deal with comments after the properties (i.e. "as HCO3")
+            dum <- unlist(strsplit(gsub(' +',' ',grep(paste0("^", reg, " "), input, value=TRUE))," "))
             if (length(dum)>2)
                 linesave <- paste(dum[3:length(dum)],collapse=" ")
         }
         ## if the simulation is just 1, then repeat it n times (adding "first")
-        ## grep(paste("^", reg, " ", sep=""), input, value=TRUE)
-        newinp <- RepSol(input, n, first)
-        nsim <- length(proplines <- grep(paste0("^", reg), newinp))
+        newinp <- RepSol(input, n)
+        nsim <- length(proplines <- grep(paste0("^", reg, " "), newinp))
     } else {
         newinp <- input
     }
@@ -206,6 +198,9 @@ Distribute <- function(input, prop, values, newname=NULL, first=NULL, wholeline=
         newname <- prop
 
     newinp[proplines] <- paste(newname, val, linesave)
+
+    if (is.character(first))
+        newinp <- c(first, newinp)
     return(newinp)  
 }
 

demo/demo-eq-surr-RF.R

 ### Licence: LGPL version 2.1
-## Time-stamp: "Last modified 2021-04-28 18:23:19 delucia"
-
+## Time-stamp: "Last modified 2021-04-29 21:21:38 delucia"
+library(RedModRphree)
 require(e1071)
 require(import)
 

demo/demo-equilibrium.R

 ### Licence: LGPL version 2.1
-### Time-stamp: "Last modified 2018-05-06 19:26:31 delucia"
+### Time-stamp: "Last modified 2021-04-29 18:54:08 delucia"
 
 db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
                              package="RedModRphree"), is.db=TRUE)
@@ -39,3 +39,361 @@ sime <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
 MatplotSingle(sim=sime, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]",
               main="Equilibrium simulation, grid 50")
 
+
+
+ebreak=FALSE
+call_pqc=TRUE
+step="time"
+maxtime=21000
+procs=4
+maxsim=50
+writeout=FALSE, 
+reduce=FALSE
+surrogate=FALSE
+
+    step <- match.arg(step)
+    ## prototype of the PHREEQC simulation
+    base     <- setup$base
+    ## SELECTED_OUTPUT and USER_PUNCH sections of the PHREEQC input data
+    first    <- setup$first
+    ## boundary conditions
+    bound    <- setup$bound
+    ## prop are the names of the properties for PHREEQC
+    prop     <- setup$prop
+    immobile <- setup$immobile
+    initsim  <- setup$initsim
+
+    if (procs > 1) {
+        if (Sys.info()[["sysname"]]=="Windows") {
+            ThisRunCluster <- parallel::makePSOCKcluster(procs)
+        } else {
+            ThisRunCluster <- parallel::makeForkCluster(procs)
+        }
+        doParallel::registerDoParallel(ThisRunCluster)
+        ## msg("Registered default doParallel cluster with ", procs, "nodes")
+    }
+
+    n  <- setup$n
+    L  <- setup$len  # m
+    U  <- setup$U    # m3/s
+    dx <- L/n        # m
+    dt <- dx/U * setup$Cf
+    ## msg("Pure advection 1D simulations")
+
+    ## Find out the number of iterations we're about to calculate
+    if (step=="fix_dt"){
+        dt <- setup$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("Courant number=", U*dt/dx, "; dt=",dt)
+        ## msg(maxiter,"iterations required") 
+    }
+    
+    if (step=="iter") {
+        maxiter <- maxtime
+        msg("Courant number=", U*dt/dx, "; dt=",dt,";")
+        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)))
+
+    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... ")
+    phreeqc::phrLoadDatabaseString(db)
+    msg("database loaded.")
+
+    if (procs > 1) {
+        parallel::clusterCall(cl=ThisRunCluster, phreeqc::phrLoadDatabaseString, db)
+        msg("Database loaded on each worker.")
+        if (surrogate) {
+            parallel::clusterExport(cl=ThisRunCluster, varlist=c("model","surrogate.FUN"), envir = environment())
+            msg("All workers are setup with the surrogate model.")
+        }
+    }
+
+    if (missing(init)) {
+        msg("missing initial chemical state")
+        if (!is.character(initsim))
+            runinitsim <- c(first,base)
+        else
+            runinitsim <- initsim
+        
+        msg("running initsim PHREEQC script assuming EQUILIBRIUM and correct PUNCH specifications...")
+        tmpfirstrun <- RunPQC(runinitsim)
+        msg("assuming homogeneous medium")
+        msg("Names gathered after first run:", paste(colnames(tmpfirstrun), collapse=", "))
+        state_C <- matrix(rep(tmpfirstrun,n), byrow=TRUE, ncol=length(prop))
+        colnames(state_C) <- colnames(tmpfirstrun)
+    } else {
+        ## msg("given initial state; assuming correct colnames")
+        state_C <- init
+    }
+    
+    attr(state_C,"immobile") <- immobile
+
+    index_saved <- 1
+    saved_complete[[index_saved]] <- list(C=Act2pH(state_C))
+    index_saved <- index_saved + 1
+    
+    ## transform the pH/pe back into activities for Advection
+    msg("first step, Advection()...")
+    state_T <- AdvectionPQC(pH2Act(state_C), inflow=bound, dx=dx, dt=dt, U=U)
+    msg("Advection done...")
+    Tr <- Ts <- 0
+    
+    ## treat special H+/pH, e-/pe cases
+    state_T <- Act2pH(state_T)
+    
+    if (reduce)
+        Tr <- system.time( reduced <- ReduceState(state_T) )[3]
+    else
+        reduced <- state_T
+    
+    msg("first step, chemistry ...")
+    Ts <- system.time(inplist <- splitMultiFix(data=reduced, procs=procs, base=base, first=first,
+                                               prop=prop, minerals=immobile, nmax=maxsim))[3]
+
+
+data=reduced 
+minerals=immobile
+nmax=maxsim
+
+    if (is.matrix(data))
+        ntot <- dim(data)[1] ## the total number of pqc simulations
+    else {
+        datamat <- matrix(data,ncol=length(prop),byrow=TRUE)
+        colnames(datamat) <- prop
+        data <- datamat
+        ntot <- nrow(datamat)
+    }
+    
+    if ( ntot  ==1) {
+        ## try and call normal Distribute (i.e., ##!is.matrix(data) ||
+        ## just 1 simulation!)
+    
+        if (verbose)
+            msg("just 1 simulation, reverting back to normal Distribute")
+        BigInp <- RepSol(base, 1, first=first)
+        namdis <- colnames(data)
+        indminer <- minerals
+        nammin <- namdis[indminer]
+        namtot <- namdis[-indminer]
+        for (dis in seq_along(nammin)) {
+            BigInp <- Distribute(BigInp, nammin[dis], paste("0.0 ", data[1,nammin[dis] ]))
+        }
+        for (dis in seq_along(namtot)) {
+            BigInp <- Distribute(BigInp, namtot[dis], data[1,namtot[dis] ])
+        }
+        return(BigInp)
+    }
+ 
+BigInp <- RepSol(base, ntot, first=first)
+
+input=BigInp
+values=data
+SI=FALSE
+
+    qprop <- gsub("\\(","\\\\(",prop)
+    qprop <- gsub("\\)","\\\\)",qprop)
+    n <- length(qprop)
+    if (n!=ncol(values)) {
+        stopmsg("prop names and values matrix of different dimension.")
+    }
+    if (is.character(minerals))  {
+            qmin <- gsub("\\(","\\\\(",minerals)
+            qmin <- gsub("\\)","\\\\)",qmin)
+            mins <- match(qmin, qprop)
+        } else
+            mins <- as.numeric(minerals)
+    
+    indAq <- seq(1,n)[-mins]
+    newinp <- input
+
+    for (i in indAq) {
+            newinp <- Distribute(newinp, prop[i], values[,i])
+    }
+
+    for (i in mins) {
+        if (SI)
+            newinp <- Distribute(newinp, prop[i], values[,i])
+        else
+            newinp <- Distribute(newinp, prop[i], paste("0.0 ",values[,i]))
+    }
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+    Tm <- system.time(res_C <- RunPQC(inplist, procs=procs))[3]
+    attr(res_C,"immobile") <- immobile
+    
+    msg("iteration 0, CPU-time = ",round(Tm,3),"[s]")
+
+    ## recompose after the reduction
+    if (reduce)
+        Tr <- Tr + system.time(state_C <- RecomposeState(Act2pH(res_C), reduced))[3]
+    else {
+        state_C <- res_C
+        immobile_col_names <- setup$prop[setup$immobile]
+        immobile_inds <- match( immobile_col_names, colnames(state_C) )
+        attr(state_C,"immobile") <- immobile_inds
+    }
+    
+    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
+    
+    if (ebreak) {
+        msg("Early break invoked, bye.")
+        return(saved_complete[[1]])
+    }
+    time <- dt
+    iter <- 1
+
+    timing[iter,] <- c(dim(reduced)[1],procs,Tm, Tr, Ts)
+
+    msg("going through iterations now!")
+    
+    on.exit({
+        msg("MAIN LOOP interrupted during iteration",iter,"!!")
+        msg("returning last calculated chemistry.")
+        print(traceback())
+        msg(" Bye.")
+        attr(state_C,"timing") <- timing
+        saved_complete$currentT <- Act2pH(state_T)
+        saved_complete$current <- Act2pH(state_C)
+        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
+        time <- time + dt
+        iter <- iter + 1
+
+        ## transform the pH/pe back into activities
+        state_T <- AdvectionPQC(pH2Act(state_C), inflow=bound, dx=dx, dt=dt, U=U)
+
+        ## treat special H+/pH, e-/pe cases
+        state_T <- Act2pH(state_T)
+
+        ## reduction of the problem
+        if(reduce)
+            Tr <- system.time(reduced <- ReduceState(state_T))[3]
+        else
+            reduced <- state_T
+
+        ## this is the place where we switch the simulator to the surrogate model
+        ## for final test for the geochemical surrogate model!
+        
+        if (surrogate) { ## use a surrogate model created by the Surrogate playground script
+            Ts <- 0.0        ## no split and distribute for this guy
+
+            Tm <- system.time( tmpsur <- Selected.Surrogate(reduced, model=model) )[3]
+            bal <- CheckBalance(baleq,reduced,tmpsur)
+            out_bal[[iter]] <- bal
+            totbal <- apply(bal, 1, mae)
+
+            nonok <- which(totbal>tol) ### max(tol,tol2))
+            if (length(nonok)>0){
+                msg("found", length(nonok), "/", nrow(reduced), " offending simulations",
+                    ifelse(call_pqc, ", calling PHREEQC"," "))
+                if (call_pqc) {
+                    reducednonok <- reduced[nonok, ]
+                    Ts <- Ts + system.time(inplist <- splitMultiFix(data=reducednonok, procs=procs, base=base, first=first,
+                                                                    prop=prop, minerals=immobile, nmax=maxsim))[3]
+                    
+                    out_inp[[iter]] <- inplist
+                    if (length(nonok)==1) {
+                        Tm <- Tm + system.time(pqcres_C <- RunPQC(inplist, procs=1))[3]
+                    } else {
+                        Tm <- Tm + system.time(pqcres_C <- RunPQC(input=inplist, procs=procs))[3]
+                    }
+                    ##mat <- apply(state, 2, as.numeric)
+
+                    tmpsur[nonok, ] <- pqcres_C
+                }
+            }
+            res_C <- tmpsur
+            attr(res_C,"immobile") <- immobile
+        } else { ## no surrogate, the original implementation
+            ## 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
+            Tm <- system.time(res_C <- RunPQC(inplist, procs=procs))[3]
+            attr(res_C,"immobile") <- immobile
+        }
+        
+        ## recompose after the reduction
+        if (reduce)
+            Tr <- Tr + system.time(state_C <- RecomposeState(res_C, reduced))[3]
+        else {
+            state_C <- res_C
+            immobile_col_names <- setup$prop[setup$immobile]
+            immobile_inds <- match( immobile_col_names, colnames(state_C) )
+            attr(state_C,"immobile") <- immobile_inds
+        }
+        
+        ## ## transform the pH/pe back into activities
+        ## state_C <- pH2Act(state_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]")
+
+        saved_complete[[index_saved]] <- list(T=Act2pH(state_T), C=Act2pH(state_C))
+        index_saved <- index_saved + 1
+        
+    }