Skip to content

GitLab

Commit 4b353575 authored by Marco De Lucia's avatar Marco De Lucia
Browse files

Renamed splitMultiFix in SplitMulti; demo-equilibrium not working

parent 2a58aebd
Hide whitespace changes
Inline Side-by-side
NAMESPACE
View file @ 4b353575
......@@ -43,6 +43,8 @@ export(RecomposeState)
export(ReduceState)
export(RepSol)
export(RunPQC)
export(SplitMulti)
export(SplitMultiKin)
export(StoichiometricMatrix)
export(SuppressSim)
export(TrailSpaces)
......@@ -50,8 +52,6 @@ export(TranspAsPert)
export(mae)
export(msg)
export(pH2Act)
export(splitMultiFix)
export(splitMultiKin)
export(stopmsg)
export(strrev)
import(doParallel)
......
R/Rphree_Kinetics.R
View file @ 4b353575
......@@ -139,7 +139,7 @@ DistributeKin <- function(input, prop=NULL, values, ident="-m0")
##' in parallel
##' @author MDL
##' @export
splitMultiKin <- function(data, procs, base, first, prop, minerals, kin, dt, ann, nmax=200, verbose=FALSE)
SplitMultiKin <- function(data, procs, base, first, prop, minerals, kin, dt, ann, nmax=200, verbose=FALSE)
{
if (is.matrix(data))
ntot <- dim(data)[1] ## the total number of pqc simulations
......@@ -360,7 +360,7 @@ DistributeKinMatrix <- function(input, prop, values, minerals, kin, dt, ann)
##' @return a list with n=procs inputs
##' @author MDL
##' @export
splitMultiFix <- function(data, procs, base, first, prop, minerals, nmax=200, verbose=FALSE)
SplitMulti <- function(data, procs, base, first, prop, minerals, nmax=200, verbose=FALSE)
{
if (is.matrix(data))
ntot <- dim(data)[1] ## the total number of pqc simulations
......
R/Rphree_ReactTrans.R
View file @ 4b353575
......@@ -352,7 +352,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
reduced <- state_T
msg("first step, chemistry ...")
Ts <- system.time(inplist <- splitMultiFix(data=reduced, procs=procs, base=base, first=first,
Ts <- system.time(inplist <- SplitMulti(data=reduced, procs=procs, base=base, first=first,
prop=prop, minerals=immobile, nmax=maxsim))[3]
Tm <- system.time(res_C <- RunPQC(inplist, procs=procs))[3]
......@@ -434,7 +434,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
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,
Ts <- Ts + system.time(inplist <- SplitMulti(data=reducednonok, procs=procs, base=base, first=first,
prop=prop, minerals=immobile, nmax=maxsim))[3]
out_inp[[iter]] <- inplist
......@@ -452,7 +452,7 @@ ReactTranspBalanceEq <- function(setup, init, maxtime, step=c("time","iter","fix
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,
Ts <- system.time(inplist <- SplitMulti(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]
......@@ -656,7 +656,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
reduced <- state_T
msg("first step, chemistry ...")
Ts <- system.time(inplist <- splitMultiKin(data=reduced, procs=procs, base=base, first=first, ann=ann,
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]
Tm <- system.time(tmp <- RunPQC(inplist, procs=procs))[3]
......@@ -743,7 +743,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
ifelse(call_pqc, ", calling PHREEQC"," "))
if (call_pqc) {
reducednonok <- reduced[nonok, ]
Ts <- Ts + system.time(inplist <- splitMultiKin(data=reducednonok, procs=procs, base=base, first=first, ann=ann,
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
......@@ -762,7 +762,7 @@ ReactTranspBalanceKin <- function(setup, init, maxtime, step=c("time","iter","fi
attr(res_C,"immobile") <- immobile
} else { ## no surrogate, the original implementation
## new input
Ts <- system.time(inplist <- splitMultiKin(data=reduced, procs=procs, base=base, first=first, ann=ann,
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
Tm <- system.time(tmp3 <- RunPQC(inplist, procs=procs))[3]
......
demo/demo-equilibrium.R
View file @ 4b353575
### Licence: LGPL version 2.1
### Time-stamp: "Last modified 2021-04-29 18:54:08 delucia"
### Time-stamp: "Last modified 2021-06-28 00:03:14 delucia"
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
package="RedModRphree"), is.db=TRUE)
......@@ -38,362 +38,3 @@ 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
}
on.exit()
if (procs > 1) {
parallel::stopCluster(ThisRunCluster)
msg("Parallelization cluster stopped")
}
msg("assign out_balance, out_input and out_results to global environment for possible later use")
assign("out_balance", out_bal, envir = .GlobalEnv)
assign("out_results", out_res, envir = .GlobalEnv)
assign("out_input", out_inp, envir = .GlobalEnv)
msg("total time of chemistry: ",round(sum(timing[,3]),3),"seconds")
msg("total number of simulations: ",sum(timing[,1]))
msg(" Bye.")
attr(saved_complete,"timing") <- timing
return(saved_complete)
}
man/splitMultiFix.Rd → man/SplitMulti.Rd
View file @ 4b353575
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/Rphree_Kinetics.R
\name{splitMultiFix}
\alias{splitMultiFix}
\name{SplitMulti}
\alias{SplitMulti}
\title{Splits a data (matrix) in n simulations to run in parallel}
\usage{
splitMultiFix(
SplitMulti(
data,
procs,
base,
......
man/splitMultiKin.Rd → man/SplitMultiKin.Rd
View file @ 4b353575
% Generated by roxygen2: do not edit by hand
% Please edit documentation in R/Rphree_Kinetics.R
\name{splitMultiKin}
\alias{splitMultiKin}
\name{SplitMultiKin}
\alias{SplitMultiKin}
\title{Splits a named matrix of data (a column is a variable) in n
simulations to run in parallel}
\usage{
splitMultiKin(
SplitMultiKin(
data,
procs,
base,
......
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