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

Added demos, fixed things, version 0.0.1

parent 57e6d53e
Package: RedModRphree Package: RedModRphree
Title: R infrastructure for training and applying geochemical surrogate models to reactive transport Title: R infrastructure for training and applying geochemical surrogate models to reactive transport
Version: 0.0.0.1 Version: 0.0.1
Authors@R: c(person("Marco", "De Lucia", email = "delucia@gfz-potsdam.de", role = c("aut", "cre")), Authors@R: c(person("Marco", "De Lucia", email = "delucia@gfz-potsdam.de", role = c("aut", "cre")),
person("Janis", "Jatnieks", email = "deltaxzz@gmail.com", role = c("ctb"))) person("Janis", "Jatnieks", email = "deltaxzz@gmail.com", role = c("ctb")))
Description: In the framework of RedMod project, we develop an R infrastructure for training and applying geochemical surrogate models to reactive transport. Description: In the framework of RedMod project, we develop an R infrastructure for training and applying geochemical surrogate models to reactive transport.
Depends: R (>= 3.2.0), parallel, phreeqc, caret, mgcv, randomForest, graphics, methods Depends: R (>= 3.2.0), parallel, plyr, data.table, phreeqc, caret, mgcv, randomForest, graphics, methods
License: LGPL-2.1 License: LGPL-2.1
Encoding: UTF-8 Encoding: UTF-8
LLazyData: true LLazyData: true
......
...@@ -12,6 +12,7 @@ export(DistributeKinMatrix) ...@@ -12,6 +12,7 @@ export(DistributeKinMatrix)
export(DistributeMatrix) export(DistributeMatrix)
export(ElementalBalanceMin) export(ElementalBalanceMin)
export(ExtractPphases) export(ExtractPphases)
export(ExtractSamples)
export(ExtractSpecies) export(ExtractSpecies)
export(ExtractTotals) export(ExtractTotals)
export(FindAllMinNames) export(FindAllMinNames)
...@@ -21,6 +22,8 @@ export(FindLogK) ...@@ -21,6 +22,8 @@ export(FindLogK)
export(FindPhase) export(FindPhase)
export(FlattenList) export(FlattenList)
export(FormulaFromBal) export(FormulaFromBal)
export(Matplot)
export(MatplotSingle)
export(ParseFormula) export(ParseFormula)
export(PlotModsInSample) export(PlotModsInSample)
export(RPhreeFile) export(RPhreeFile)
...@@ -44,7 +47,9 @@ export(splitMultiFix) ...@@ -44,7 +47,9 @@ export(splitMultiFix)
export(splitMultiKin) export(splitMultiKin)
export(stopmsg) export(stopmsg)
import(caret) import(caret)
import(data.table)
import(graphics) import(graphics)
import(mgcv) import(mgcv)
import(phreeqc) import(phreeqc)
import(plyr)
import(randomForest) import(randomForest)
## Functions for simple visualization ### RedModRphree, functions to perform basic visualization
### Licence: LGPL version 2.1
### Time-stamp: "Last modified 2018-05-06 19:30:10 delucia"
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018 ##' @title Display the profiles of some variables along a 1D reactive
### Time-stamp: "Last modified 2018-05-03 17:36:21 delucia" ##' transport simulation
##' @param sim The simulation to display
##' @param sample optional, time frame position to be desplayed. If
##' omitted, the last is showed
##' @param vars optional, character vector containing the names of
##' variables to display (up to 7)
##' @param ... further parameters passed to \code{matplot}, such as
##' main, xlab, ylab...
##' @return NULL
##' @author MDL
##' @export
MatplotSingle <- function(sim, sample=1L, vars=c("Ca","Mg","Calcite","Dolomite"), ...) {
if(missing(sample))
sample <- length(sim)
if (length(vars)>7)
stopmsg("Please specify up to seven variables")
cols <- c("red", "black","blue","olivedrab3", "orange", "light blue", "grey")[1:length(vars)]
frame <- sim[[sample]]$C
matplot(frame[,vars], type="l", lwd=2, lty="solid", col=cols, ...)
legend("topright", vars, text.col=cols, bty="n")
}
##' @title Display the profiles of some variables for 2 reactive
##' transport simulation (one using surrogates)
##' @param sur The simulation with surrogate
##' @param sim The simulation using only phreeqc
##' @param sample optional, time frame position to be desplayed. If
##' omitted, the last is showed
##' @param vars optional, character vector containing the names of
##' variables to display (up to 7)
##' @param ... further parameters passed to \code{matplot}, such as
##' main, xlab, ylab...
##' @return NULL
##' @author MDL
##' @export
Matplot <- function(sur, sim, sample=1L, vars=c("Ca","Mg","Calcite","Dolomite"), ...) {
if(missing(sample))
sample <- length(sim)
if (length(vars)>7)
stopmsg("Please specify up to seven variables")
cols <- c("red", "black","blue","olivedrab3", "orange", "light blue", "grey")[1:length(vars)]
rs <- sur[[sample]]$C
rn <- sim[[sample]]$C
matplot(rn[,vars], type="l", lwd=2, lty="solid", col=cols, ...)
matplot(rs[,vars], type="l", lwd=2, lty="dashed",col=cols, add=TRUE)
legend("topright", vars, text.col=cols, bty="n")
}
##' @title Scatter plots of phreeqc simulations vs surrogate response ##' @title Scatter plots of phreeqc simulations vs surrogate response
##' @param mod list of regression models, one per variable ##' @param mod list of regression models, one per variable
......
## Time-stamp: "Last modified 2018-05-06 18:55:57 delucia"
##' @title Extracts input/output tables from a list of Reactive
##' Transport Simulations
##' @param simlist a list containing the simulations
##' @param var_design optional, a character vector containing the
##' names of the variables to extract for design. If omitted all
##' variables are selected
##' @param var_result optional, a character vector containing the
##' names of the variables to extract for results. If omitted all
##' variables are selected
##' @return a list with 2 elements: design and result
##' @author MDL
##' @export
ExtractSamples <- function(simlist, var_design, var_result)
{
## transform the given simulations in explicit input/output tables for training
totsim <- vector(mode = "list", length = sum(sapply(simlist,length)) - length(simlist))
if (missing(var_design))
var_design <- colnames(simlist[[1]][[2]]$T)
if (missing(var_result))
var_result <- colnames(simlist[[1]][[2]]$C)
k <- 1
for (i in seq_along(simlist)) {
for (j in seq(2, length(simlist[[i]]))) {
totsim[[k]] <- simlist[[i]][[j]]
k = k + 1
}
}
tmpT <- plyr::ldply(totsim, function(x) x$T[, var_design], .id=NULL)
tmpC <- plyr::ldply(totsim, function(x) x$C[, var_result], .id=NULL)
unT <- mgcv::uniquecombs(tmpT)
## unique design combinations
iT <- attr(unT,"index")
index <- match(seq(1,max(iT)), iT)
unC <- tmpC[index,]
saved_index_comp <- as.integer(rownames(unC))
design <- tmpT[index,]
result <- tmpC[index,]
if(!all.equal(unT,tmpT[index,], check.attributes=FALSE))
stopmsg("Something wrong with extracting unique combinations of design and results")
return(list(design=design, result=result))
}
...@@ -3,8 +3,8 @@ ...@@ -3,8 +3,8 @@
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018 ### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2018
### Time-stamp: "Last modified 2018-05-03 17:37:43 delucia" ### Time-stamp: "Last modified 2018-05-03 17:37:43 delucia"
#' @import phreeqc mgcv ##' @import phreeqc mgcv plyr data.table
#' @import graphics caret randomForest ##' @import graphics caret randomForest
NULL NULL
##' Display a message prepending the name of the calling function ##' Display a message prepending the name of the calling function
......
demo-equilibrium 1D reactive transport with equilibrium chemistry demo-equilibrium 1D reactive transport with equilibrium chemistry
demo-kinetics 1D reactive transport with kinetic chemistry demo-kinetics 1D reactive transport with kinetic chemistry
demo-eq-generate-data Generate the training data set from 3 equilibrium simulations
demo-eq-surr-RF Generate the training data set, train randomForest surrogates and perform Reactive Transport simulations using them
### Licence: LGPL version 2.1
## Time-stamp: "Last modified 2018-05-06 19:41:08 delucia"
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat", package="RedModRphree"), is.db=TRUE)
phreeqc::phrLoadDatabaseString(db)
selec <- c(
"SELECTED_OUTPUT",
"-high_precision true",
"-reset false",
"USER_PUNCH",
"-head C Ca Cl Mg pH pe Calcite Dolomite",
"10 PUNCH TOT(\"C(4)\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")")
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pH 9.91 charge",
"pe 4.0 O2(g) -0.68",
"C(4) 1.2279E-04",
"Ca 1.2279E-04",
"Mg 1.0E-12",
"Cl 1.0E-12",
"PURE 1",
"Calcite 0.0 2.07E-04",
"Dolomite 0.0 0.00E+00",
"END")
maxsim <- 50 ## minimum number of single chemical simulations in each time step to spawn parallelization
grid <- 100 ## grid discretization (number of elements)
L <- 0.5 ## total system length [m]
U <- 9.375e-6 ## imposed (constant) Darcy velocity [m3/s]
Cf <- 0.8 ## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound <- c(H=1E-7, e=1E-4, Mg=0.001, Cl=0.002)
prop <- c("C(4)","Ca","Cl","Mg","pH","pe","Calcite","Dolomite")
## simulations with grid 50, 100 and 200
setup <- list(n=50, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim50 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
setup <- list(n=100, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim100 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
setup <- list(n=200, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim200 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
sam <- ExtractSamples(list(sim50, sim100, sim200),
var_design=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"),
var_result=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"))
## data.table::fwrite(sam$design, file = "RedModRphree/inst/extdata/demo_equilibrium_design.csv",
## append = FALSE, quote = "auto", sep = ",", eol = "\n", na = "NA", dec = ".",
## row.names = FALSE, col.names = TRUE)
## data.table::fwrite(sam$result, file = "RedModRphree/inst/extdata/demo_equilibrium_result.csv",
## append = FALSE, quote = "auto", sep = ",", eol = "\n", na = "NA", dec = ".",
## row.names = FALSE, col.names = TRUE)
### Licence: LGPL version 2.1
## Time-stamp: "Last modified 2018-05-06 20:05:39 delucia"
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat", package="RedModRphree"), is.db=TRUE)
phreeqc::phrLoadDatabaseString(db)
selec <- c(
"SELECTED_OUTPUT",
"-high_precision true",
"-reset false",
"USER_PUNCH",
"-head C Ca Cl Mg pH pe Calcite Dolomite",
"10 PUNCH TOT(\"C(4)\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")")
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pH 9.91 charge",
"pe 4.0 O2(g) -0.68",
"C(4) 1.2279E-04",
"Ca 1.2279E-04",
"Mg 1.0E-12",
"Cl 1.0E-12",
"PURE 1",
"Calcite 0.0 2.07E-04",
"Dolomite 0.0 0.00E+00",
"END")
maxsim <- 50 ## minimum number of single chemical simulations in each time step to spawn parallelization
grid <- 100 ## grid discretization (number of elements)
L <- 0.5 ## total system length [m]
U <- 9.375e-6 ## imposed (constant) Darcy velocity [m3/s]
Cf <- 0.8 ## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound <- c(H=1E-7, e=1E-4, Mg=0.001, Cl=0.002)
prop <- c("C(4)","Ca","Cl","Mg","pH","pe","Calcite","Dolomite")
## simulations with grid 50, 100 and 200
setup <- list(n=50, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim50 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
setup <- list(n=100, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim100 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
setup <- list(n=200, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sim200 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
sam <- ExtractSamples(list(sim50, sim100, sim200),
var_design=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"),
var_result=c("C","Ca","Cl","Mg", "pH", "pe", "Calcite", "Dolomite"))
## This code was used to generate the included csv data
## data.table::fwrite(sam$design, file = "RedModRphree/inst/extdata/demo_equilibrium_design.csv",
## append = FALSE, quote = "auto", sep = ",", eol = "\n", na = "NA", dec = ".",
## row.names = FALSE, col.names = TRUE)
## data.table::fwrite(sam$result, file = "RedModRphree/inst/extdata/demo_equilibrium_result.csv",
## append = FALSE, quote = "auto", sep = ",", eol = "\n", na = "NA", dec = ".",
## row.names = FALSE, col.names = TRUE)
## this is how you can import those data without running the whole thing
## sam <- list(
## design = data.matrix(data.table::fread(file = system.file("extdata", "demo_equilibrium_design.csv", package="RedModRphree")))
## result = data.matrix(data.table::fread(file = system.file("extdata", "demo_equilibrium_result.csv", package="RedModRphree")))
## )
baleq <- list(Calcite=c(C=1, Ca=1), Dolomite=c(C=2, Ca=1, Mg=1))
balance <- CheckBalance(baleq, sam$design, sam$result)
set.seed(95014)
fitControl <- trainControl(method = "cv", number=5,
allowParallel = TRUE,
returnData = FALSE,
returnResamp="none",
savePredictions="none",
predictionBounds=rep(TRUE,2))
## this takes around 8 minutes using 4 cores
traintime <- system.time({
mods <- parallel::mclapply(colnames(sam$result), function (x)
return(train(x=sam$design, y=sam$result[,x],trControl = fitControl,
method="parRF")), mc.cores=4)
names(mods) <- colnames(sam$result)
})
PlotModsInSample(mods, sam$result, sam$design)
saveRDS(mods, "models_equilibrium_RF.rds")
## fitControl2 <- trainControl(method = "none",
## allowParallel = TRUE,
## returnData = FALSE,
## returnResamp="none",
## savePredictions="none",
## verboseIter=FALSE,
## predictionBounds=rep(TRUE,2))
## mods2 <- parallel::mclapply(colnames(sam$result), function (x)
## return(train(x=sam$design, y=sam$result[,x],trControl = fitControl2,
## method="parRF")), mc.cores=8)
## names(mods2) <- colnames(sam$result)
## PlotModsInSample(mods2, sam$result, sam$design)
## to use the surrogate in the reactive transport simulations we need a function
MySurr <- function(state, model){
order_for_surrogate <- c("C","Ca","Cl","Mg", "pH","pe", "Calcite","Dolomite")
remember_names <- colnames(state)
rem_attr <- attr(state,"immobile")
ss <- state[, order_for_surrogate]
pred <- lapply(names(model), function(x) as.numeric(predict(model[[x]], ss)))
names(pred) <- names(model)
out <- data.matrix(as.data.frame(pred))
## order back to how it is expected in "state"
out <- out[,remember_names]
attr(out,"immobile") <- rem_attr
return(out)
}
## now run the simulations with surrogates
setup <- list(n=50, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
immobile=c(7,8), initsim=c(selec, base))
sur50_4 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
model=mods, baleq=baleq, tol=1E-4, call_pqc=TRUE)
sur50_5 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
model=mods, baleq=baleq, tol=1E-5, call_pqc=TRUE)
## some visualization
par(mfrow=c(2,1))
Matplot(sur=sur50_4, sim=sim50, xlab="Domain element",
ylab="Ca, Mg [molal], Calcite, Dolomite [mol]",
main="RandomForest surrogate, tol=1E-4")
Matplot(sur=sur50_5, sim=sim50, xlab="Domain element",
ylab="Ca, Mg [molal], Calcite, Dolomite [mol]",
main="RandomForest surrogate, tol=1E-5")
## setup <- list(n=100, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
## immobile=c(7,8), initsim=c(selec, base))
## sur100_4 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
## procs=4, maxsim=maxsim, writeout=FALSE,
## reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
## model=mods, baleq=baleq, tol=1E-4, call_pqc=TRUE)
## sur100_5 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
## procs=4, maxsim=maxsim, writeout=FALSE,
## reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
## model=mods, baleq=baleq, tol=1E-5, call_pqc=TRUE)
## setup <- list(n=200, len=L, U=U, base=base, first=selec, Cf=Cf, bound=bound, prop=prop,
## immobile=c(7,8), initsim=c(selec, base))
## sur200_4 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
## procs=4, maxsim=maxsim, writeout=FALSE,
## reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
## model=mods, baleq=baleq, tol=1E-4, call_pqc=TRUE)
## sur200_5 <- ReactTranspBalanceEq(setup=setup, db=db, step="time", maxtime=21000,
## procs=4, maxsim=maxsim, writeout=FALSE,
## reduce=FALSE, surrogate=TRUE, surrogate.FUN=MySurr,
## model=mods, baleq=baleq, tol=1E-5, call_pqc=TRUE)
## par(mfrow=c(2,1))
## Matplot(50, sur=sur100_4, sim=sim100, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]" ,main="RandomForest surrogate, tol=1E-4")
## Matplot(50, sur=sur100_5, sim=sim100, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]" ,main="RandomForest surrogate, tol=1E-5")
## par(mfrow=c(2,1))
## Matplot(99, sur=sur200_4, sim=sim200, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]" ,main="RandomForest surrogate, tol=1E-4")
## Matplot(99, sur=sur200_5, sim=sim200, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]" ,main="RandomForest surrogate, tol=1E-5")
### Licence: LGPL version 2.1 ### Licence: LGPL version 2.1
### Time-stamp: "Last modified 2016-04-15 17:17:41 delucia" ### Time-stamp: "Last modified 2018-05-06 19:26:31 delucia"
library(RedModRphree)
Matplot <- function(sample, sur, sim, head=FALSE, ...) { db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat",
rs <- sur[[sample]]$C package="RedModRphree"), is.db=TRUE)
rn <- sim[[sample]]$C
if(head){
cat("Surrogate:\n")
print(head(rs))
cat("Simulation:\n")
print(head(rn))
}
matplot(rn[,c("Ca","Mg","Calcite","Dolomite")], type="l", lwd=2, lty="solid", col=c("red", "black","blue","olivedrab3"), ...)
matplot(rs[,c("Ca","Mg","Calcite","Dolomite")], type="l", lwd=2, lty="dashed",col=c("red", "black","blue","olivedrab3"), add=TRUE)
legend("topright", c("Ca","Mg","Calcite","Dolomite"), text.col=c("red", "black","blue","olivedrab3"), bty="n")
}
phreeqc::phrLoadDatabaseString(db)
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",
"-bad_step_max 1000",
"Calcite",
"-m 0.000207",
"-parms 3.20",
"Dolomite",
"-m 0.0",
"-parms 0.32",
"END")
selec <- c("SELECTED_OUTPUT", "-high_precision true", "-reset false",
"USER_PUNCH", "-head C Ca Cl Mg pH pe Calcite Dolomite",
"10 PUNCH TOT(\"C(4)\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), -LA(\"H+\"), -LA(\"e-\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")")
initsim <- c( base <- c("SOLUTION 1", "units mol/kgw", "temp 25.0",
"SELECTED_OUTPUT", "pH 9.91 charge", "pe 4.0 O2(g) -0.68",
"-high_precision true", "C(4) 1.2279E-04", "Ca 1.2279E-04", "Mg 1.0E-12",
"-reset false", "Cl 1.0E-12", "PURE 1", "Calcite 0.0 2.07E-04",
"USER_PUNCH", "Dolomite 0.0 0.00E+00", "END")
"-head C Ca Cl Mg H e O2g Calcite Dolomite",
"10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), ACT(\"H+\"), ACT(\"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")
db <- RPhreeFile( "db/pqc_strip_kin.dat", is.db=TRUE )
phrLoadDatabaseString(db)
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")
prop <- c("C","Ca","Cl","Mg","pH","pe","O2g", "Calcite","Dolomite")
maxsim <- 20 ## minimum number of single chemical simulations in each time step to spawn parallelization
grid <- 50 ## grid discretization (number of elements)
L <- 0.5 ## total system length [m]
U <- 9.375e-6 ## imposed (constant) Darcy velocity [m3/s]
Cf <- 0.9 ## Courant number to add some numerical dispersion
## boundary conditions (concentration of the injected fluid entering into the 1D column)
bound <- c(H=1E-7, e=1E-4, Mg=0.001, Cl=0.002)
setup <- list(n=50, len=L, U=U, base=base, first=selout, Cf=Cf, bound=bound, dt=210,
prop=prop, immobile=c(7,8,9), kin=c(8,9), ann=list(O2g=-0.1675), initsim=initsim)
oo <- ReactTranspBalanceKin(setup=setup, db=db, step="fix_dt", maxtime=2100,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=TRUE, surrogate=FALSE)
of <- ReactTranspBalanceKin(setup=setup, db=db, step="fix_dt", maxtime=2100,
procs=4, maxsim=maxsim, writeout=FALSE,
reduce=FALSE, surrogate=FALSE)
Matplot(101, sur=oo, sim=of, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]" ,main="RandomForest surrogate")
load(file="Rdata/20180411_CalcKin_Surr_EGU.Rdata")
MySurr <- function(state, model){
require(caret)
excl <- which(colnames(state)=="O2g")
remember_state_names <- colnames(state)
des <- state[,-excl]
if (NA %in% des) {
print(des)
stop("Mysurr: NA in data")