Skip to content
GitLab
Commit 1ffb5dcc authored by Marco De Lucia's avatar Marco De Lucia
Browse files

Re-added InputFromList. Reworked demo-validate.

parent 5683da32
Show whitespace changes
Inline Side-by-side
NAMESPACE
View file @ 1ffb5dcc
......@@ -24,6 +24,7 @@ export(FlattenList)
export(FormSelectedOutput)
export(FormulaFromBal)
export(FormulaToExpression)
export(InputFromList)
export(Matplot)
export(MatplotSingle)
export(ParseFormula)
......
R/Rphree_Utils.R
View file @ 1ffb5dcc
### Utility functions for RedModRphree
### Marco De Lucia, delucia@gfz-potsdam.de, 2009-2021
### Time-stamp: "Last modified 2021-06-26 17:40:57 delucia"
### Time-stamp: "Last modified 2021-06-30 16:33:18 delucia"
##' Replicates an input buffer containing only one SOLUTION, taking
##' care of SOLUTION/KINETICS/PURE identifiers. Eventually inserts a
......@@ -351,6 +351,55 @@ RGetPhases <- function(lin, which=NULL, delta=TRUE)
}
##' Transforms a parsed solution list in a new input buffer. Given the
##' large number of options offered by PHREEQC, this function must be
##' considered just a stub and the user should adapt it to his needs.
##' In particular, KINETIC blocks are programmatically disregarded.
##'
##' Since the function operates on lists formed from PHREEQC output
##' buffers, it automatically uses molality as unit and adds the line
##' "units mol/kgw" to the script.
##' @title Transform a parsed output list in a new input buffer
##' @param lin a list containing a parsed PHREEQC output, typically
##' resulting from \link{ReadOut}.
##' @param rem Optional, character vector containing the names of the
##' desc properties to include in the new input script. If
##' missing, pH, pe, temp and water are included by default.
##' @param title optional, character. TITLE to appear in the new input
##' buffer.
##' @param signif integer specifying how many significant digits will
##' be written in the buffer. Defaults to 10.
##' @return An input buffer ready for further manipulation and/or a
##' new Rphree call.
##' @author MDL
##' @export
InputFromList <- function(lin, rem, title=NULL, signif=10L)
{
## check that totals, phases, pH and pe are in the list
if (missing(rem))
rem <- c("pH","pe","temp","water")
if (!is.null(title))
title <- paste("TITLE",title)
## Desc
ind_desk <- match(rem,rownames(lin$desc))
desc <- c(paste(rownames(lin$desc)[ind_desk], signif(lin$desc[ind_desk,1],signif), sep=" " ), "units mol/kgw")
## totals
tot <- paste(rownames(lin$tot), signif(lin$tot[,"molal"],signif))
## Equilibrium phases
if ("pphases" %in% names(lin))
phases <- paste(rownames(lin$pphases), 0.0, signif(lin$pphases[,1]))
else
phases <- NULL
## All SOLUTION
buff <- c(title, "SOLUTION 1", desc, tot,"PURE 1", phases, "END")
## print(buff)
return(buff)
}
## ##' Workhorse function for extraction of results from a solution,
## ##' specifically equilibrium phases.
## ##'
......
demo/demo-validate.R
View file @ 1ffb5dcc
### VALIDATE demo
library(RedModRphree)
library(colorspace)
clrs2 <- sequential_hcl(8, "Viridis", alpha = 0.5)
## NB: Du gibst einen Vector mit n Elementen (Farben), in diesem Fall
## 6 und 8, und wenn die Farbe weniger sind als die einzelnen Zonen im
## Plot, sie werden wiederholt.
par(mfrow=c(1,3))
## here rechnen wir ohne Mineralphasen, d.h. alle wässrige Spezien die Fe(2) enthalten
a1 <- Pourbaix(element="Fe(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51), Patm=1,
base=base, aqonly=TRUE, db=llnl.dat, procs = 4)
## here rechnen wir mit Mineralphasen, alle Spezien die Fe(2) enthalten
a2 <- Pourbaix(element="Fe(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51), Patm=1,
base=base, aqonly=FALSE, db=llnl.dat, procs = 4, colors = clrs2)
## here rechnen wir mit Mineralphasen wie vorher, aber schliessen Hematite aus da wir wissen dass
## es nur Fe(3) enthält
a3 <- Pourbaix(element="Fe(2)", pe=seq(-10,12, length=51), pH=seq(2,10, length=51), Patm=1,
base=base, aqonly=FALSE, suppress=c("Hematite"), db=llnl.dat, procs = 4)
str(a3)
## Example from Huang2016
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pressure 1",
"-water 1",
"pH 7",
"pe 4.0",
## "C .1",
## "Cl 0.0001 charge",
"Cu 0.001",
"Fe 0.001",
"S 0.01",
"PURE 1",
"CO2(g) -1 10",
"END")
Hua <- Pourbaix(element = "Cu", pe=seq(-10,12, length=21), pH=seq(2,10, length=21), Patm=1,
base=base, aqonly=FALSE, suppress=c("CO2(g)"), db=llnl.dat, procs = 4)
library(devtools)
setwd("RedModRphree")
load_all()
AddLabel <- function(lab, cex=3, offset=c(0,0)) {
par(xpd=TRUE)
text(grconvertX(0+offset[1], from = "npc", to = "user"),
grconvertY(0+offset[2], from = "npc", to = "user"),
lab, cex=cex, adj=c(1,1))
par(xpd=FALSE)
}
library(RedModRphree)
## Define initial base script
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pressure 1",
"-water 1",
"pH 7",
"pe 4.0",
"Fe 0.001",
"END")
## Color palette
viridis <- colorspace::sequential_hcl(8, "Viridis", alpha = 0.4)
## Calculation grid
pe <- seq(-10, 12, length=101)
pH <- seq(1,12, length=101)
cairo_pdf("FePourbaix.pdf", width = 13, height = 6)
par(mfrow=c(1,2), family="serif",mar=c(3.5,4,0.5,4))
pqcdat <- Pourbaix(element = "Fe", pe = pe, pH = pH, Patm = 1,
base = base, aqonly = TRUE, db = phreeqc.dat,
procs = 4, colors = viridis,
# ann.title = "", ##"Pourbaix diagram for Fe, phreeqc.dat",
ann.sub = "")
text(3,-7,labels = "phreeqc.dat", cex = 2)
AddLabel("a", offset=c(0,-0.05), cex=2.5)
llnldat <- Pourbaix(element = "Fe", pe = pe, pH = pH, Patm = 1, base = base,
aqonly = TRUE, db = llnl.dat, procs = 4,
colors = viridis,
# ann.title = "Pourbaix diagram for Fe, llnl.dat",
ann.sub = "")
text(3,-7,labels = "llnl.dat", cex = 2)
AddLabel("b", offset=c(0,-0.05), cex=2.5)
dev.off()
llnldat <- Pourbaix(element = "Fe(2)", pe = pe, pH = pH, Patm = 1, base = base,
aqonly = TRUE, db = llnl.dat, procs = 4,
colors = viridis)
llnldat <- Pourbaix(element = "Fe(3)", pe = pe, pH = pH, Patm = 1, base = base,
aqonly = TRUE, db = llnl.dat, procs = 4,
colors = viridis)
#dev.off()
par(mfrow=c(1,3), family="serif")
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pressure 1",
"-water 1",
"pH 7",
"pe 4.0",
"Fe 0.001",
"Cl 0.002",
"END")
pqcdat <- Pourbaix(element = "Fe", pe = pe, pH = pH, Patm = 1,
base = base, aqonly = TRUE, db = phreeqc.dat,
procs = 4, colors = viridis,
ann.title = "Pourbaix diagram for Fe, phreeqc.dat",
ann.sub = "")
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat", package="RedModRphree"), is.db=TRUE)
llnldat <- Pourbaix(element = "Fe", pe = pe, pH = pH, Patm = 1, base = base,
aqonly = TRUE, db = llnl.dat, procs = 4,
colors = viridis,
ann.title = "Pourbaix diagram for Fe, llnl.dat",
ann.sub = "")
base <- c("SOLUTION 1",
"units mol/kgw",
"temp 25",
"pressure 1",
"pH 7",
"pe 4.0",
"Ca 0.25",
"Cu 0.001",
"Cl 1 charge",
"Na 0.5",
"Fe 3.0E-05",
"END")
res <- Pourbaix(element="Cu",
pe=seq(-10,12, len = 101),
pH=seq( 2,10, len = 101),
Patm=100,
base=base,
db=llnl.dat,
procs = 2)
##############
## Look at the included files: databases and one simple example PHREEQC script
library(RedModRphree)
library(magrittr) ## enable the "%>%" pipe operator
?Distribute
## Define initial base script
base <- c(
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pressure 1",
"-water 1",
"pH 7",
"pe 4.0",
"END")
## Add some properties
sim <- base %>%
AddProp(name = "Na", val=.1, cat="tot") %>%
AddProp(name = "Cl", val=.1, cat="tot") %>%
AddProp(name = "Calcite", val="0.0 0.1", cat="pphases")
## run the script
phreeqc::phrLoadDatabaseString(phreeqc.dat)
phreeqc::phrSetOutputStringsOn(TRUE)
phreeqc::phrRunString(sim)
out <- phreeqc::phrGetOutputStrings()
rmout <- ReadOut(out)[[1]]
SelOut <- FormSelectedOutput(rmout)
dput(SelOut)
SelOut <- c(SelOut, " -equilibrium_phases Calcite")
sim2 <- Distribute(input=sim, prop="temp", values=seq(25, 100, by=5), first = SelOut)
phreeqc::phrLoadDatabaseString(phreeqc.dat)
phreeqc::phrRunString(sim2)
res <- phreeqc::phrGetSelectedOutput()
out <- phreeqc::phrGetOutputStrings()
rmout <- ReadOut(out)
str(rmout)
SolubCc <- RPinfo(rmout, "tot","Ca")
SolubCc
library(devtools)
setwd("RedModRphree")
load_all()
base <- c("SOLUTION 1",
"units mol/kgw",
"temp 25",
"pressure 1",
"pH 7",
"water 1",
"pH 9.91 charge",
"pe 4.0",
"Ca 0.25",
"Cu 0.001",
"Cl 1 charge",
"Na 0.5",
"Fe 3.0E-05",
"C 1.2279E-04",
"Ca 1.2279E-04",
"Mg 0.001",
"Cl 0.002",
"KINETICS 1",
"-steps 100",
"-bad_step_max 1000",
"Calcite",
"-m 0.000207",
"-parms 3.20",
"Dolomite",
"-m 0.0",
"-parms 0.32",
"END")
## par(mfrow = c(1,2))
cairo_pdf("CuPourbaix.pdf", width = 7, 7, family="serif")
res <- Pourbaix(element="Cu",
pe=seq(-10,12, len = 101),
pH=seq( 2,10, len = 101),
Patm=100,
base=base,
db=llnl.dat,
procs = 4)
dev.off()
library(RedModRphree)
%>%
AddProp(name = "Na", val=0.1, cat="tot") %>%
AddProp(name = "Calcite", val="0.0 10", cat="pphases") %>%
AddProp(name = "Dolomite", val="0.0 10", cat="pphases") %>%
demo("demo-Pourbaix", package="RedModRphree")
## Your first Rphree simulation:
equilln <- Rphree(ex1, db=llndb, sel=mysel, write=TRUE, out="ex1")
equipqc <- Rphree(ex1, db=pqcdb, sel=mysel, write=FALSE)
## understanding the output
str(equilln)
equilln$species
equilln$pphases
## From a solution, create a new script (for restart):
eq2 <- RInputFromList(equilln)
#####################
## Distribute: assign specific values to a property in input
## Only one value:
ext <- Distribute(input=ex1, prop="Cl", values=3)
ext <- Distribute(input=ext, prop="Na", values=3)
equi3pqc <- Rphree(ext, db=pqcdb, sel=mysel)
equi3lln <- Rphree(ext, db=llndb, sel=mysel)
## RPinfo: extract informations from calculated solutions
RPinfo(equi3lln,"tot","C")
RPinfo(equi3pqc,"tot","C")
RPinfo(equi3lln,"pphases","CO2(g)")
RPinfo(equi3pqc,"pphases","CO2(g)")
## Graphically compare the 2 solutions
VizSol(list(equi3lln,equi3pqc),names=c("llnl.dat","phreeqc.dat"),log="y",col=c("light blue","light green"))
#####################
## Distribute more than 1 value
## The argument "values" in function "Distribute" must be of length 1
## OR the same length as the solutions already in the input!
ext_range <- Distribute(ext,"Cl",seq(1,4,0.1))
ext_range <- Distribute(ext_range,"Na",seq(1,4,0.1))
## ext_range now is an input containing 31 solutions. This can be run at once:
equi_extrange_llnl <- Rphree(ext_range, db=llndb,sel=mysel)
equi_extrange_phrq <- Rphree(ext_range, db=pqcdb,sel=mysel)
## the result is a list of lists!
str(equi_extrange_llnl)
## extracting results: RPinfo automagically gets all properties
RPinfo(equi_extrange_llnl,"desc","pH") ## 31 values!
## standard R plot
plot(type="l",seq(1,4,0.1), RPinfo(equi_extrange_llnl,"desc","pH"),
main="pH vs salinity under constant CO2(g) pressure",
xlab="Na, Cl concentration [molal]", ylab="pH", ylim=c(5.44, 5.61))
lines(seq(1,4,0.1), RPinfo(equi_extrange_phrq,"desc","pH"), col="red")
legend("bottomleft",legend=c("phreeqc.dat","llnl.dat"),lwd=2, col=c("red","black"))
## AddProp: add a property to an input!
## remember for "pphases" to "paste" SI and moles!
extmin <- AddProp(ext, name='Calcite', cat="pphases", values="0.0 10")
extmin <- AddProp(extmin, name='Anhydrite', cat="pphases", values="0.0 10")
resextmin_pqc <- Rphree(extmin, db=pqcdb,sel=mysel)
resextmin_lln <- Rphree(extmin, db=llndb,sel=mysel)
VizMinDelta(list(resextmin_pqc,resextmin_lln),phases=c("Anhydrite","Calcite"), names=c("phreeqc.dat","llnl.dat"))
## By default nothing gets written to disc, but you can let Rphree
## write the standard outfile by specifying write=TRUE:
equi_extrange_phrq <- Rphree(ext_range, db=pqcdb, sel=mysel, write=TRUE, out="equi_extrange")
## Now in your current working dir you will have the file "equi_extrange.Rout"
dir(pattern="Rout")
# [1] "equi_extrange.Rout" "ex1.Rout"
## There are helper functions to (partially) parse PHREEQC outfiles:
equi_extrange_fromoutfile <- RReadOut(out="equi_extrange.Rout")
## Currently tot, pphases and partly "desc" are read from outfile.
equi_extrange_phrq[[1]]$pphases
equi_extrange_fromoutfile[[1]]$pphases
library(RedModRphree) ## load RedModRphree
library(magrittr) ## enable the "%>%" pipe operator
## Define initial base script
base <- c(
initsim <- c("SELECTED_OUTPUT",
"-high_precision true",
"-reset false",
"USER_PUNCH",
"-head C Ca Cl Mg H e Calcite Dolomite",
"10 PUNCH TOT(\"C\"), TOT(\"Ca\"), TOT(\"Cl\"), TOT(\"Mg\"), ACT(\"H+\"), ACT(\"e-\"), EQUI(\"Calcite\"), EQUI(\"Dolomite\")",
"SOLUTION 1",
"units mol/kgw",
"temp 25.0",
"pressure 1",
"-water 1",
"pH 7",
"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",
"Calcite 0.0 2.07E-3",
"Dolomite 0.0 0.0",
"END")
## Add some properties: Na, Cl and Calcite
## and put the resulting new script into ``sim''
sim <- base %>%
AddProp(name = "Na", val=.1, cat="tot") %>%
AddProp(name = "Cl", val=.1, cat="tot") %>%
AddProp(name = "Calcite", val='0.0 0.1',
cat="pphases")
## Temperature vector from 25 to 100
tvec <- seq(25, 100, by=5)
## Distribute temperatures to the base script
simT <- Distribute(sim, "temp", values=tvec)
## Load the database
phreeqc::phrLoadDatabaseString(llnl.dat)
## Activate the output as text
phreeqc::phrSetOutputStringsOn(TRUE)
## Run the script
phreeqc::phrRunString(simT)
## Retrieve the results
out <- phreeqc::phrGetOutputStrings()
res <- ReadOut(out)
SelOut <- FormSelectedOutput(res[[1]])
## Add ``-equilibrium_phases'' line to selected output
SelOut <- c(SelOut, " -equilibrium_phases Calcite")
simT <- Distribute(sim, "temp", values=tvec, first=SelOut)
phreeqc::phrLoadDatabaseString(llnl.dat)
res2 <- RunPQC(simT)
NaCl <- seq(0.1, 2, by=0.1)
temp <- seq(25, 100, by=5)
df <- expand.grid(temp=temp, Na=NaCl)
df <- cbind(df, Cl=df$Na)
inplist <- vector(mode="list", length=4)
inplist <- splitMultiFix(data=df, procs=4, base=sim, first=SelOut, prop=colnames(df), minerals=NA, verbose = FALSE)
tmp <- RunPQC(inplist, procs = procs)
procs <- 4
spl <- rep(seq(1, procs), each = nrow(df)/4)
spldf <- split(df, spl)
## Create 320 values to Distribute
NaCl <- seq(0.1, 2, by=0.1)
temp <- seq(25, 100, by=5)
df <- expand.grid(temp=temp, Na=NaCl)
df <- cbind(df, Cl=df$Na)
## Intended
procs <- 4
spl <- rep(seq(1, procs), each = nrow(df)/4)
spldf <- split(df, spl)
## Workhorse function for lapply
.distr <- function (data, base, first) {
Distribute(base, prop="Na", values = data$Na, first=SelOut) %>%
Distribute(prop="Cl", values = data$Cl) %>%
Distribute(prop="temp", values = data$temp)
}
inplist <- lapply(spldf, .distr, base=sim, first=SelOut)
registerDoParallel()
phreeqc::phrLoadDatabaseString(llnl.dat)
res3 <- RunPQC(inplist, procs=4)
df
economics
FindLogK("Calcite", db=llnl.dat, type = c("analytical"))
options(width = 80)
colnames(res3)
resMat <- matrix(-res3$d_Calcite, ncol=length(temp), byrow=TRUE)
plotdf <- data.frame(x=df$temp, y=df$Cl, z=-res3$d_Calcite)
require(plotly)
plot_ly(x = df$temp, y = df$Na, z = -res3[, "d_Calcite"]) %>% add_markers() %>%
layout(scene = list(
xaxis=list(title='T / Celsius'),
yaxis=list(title='NaCl / molal'),
zaxis=list(title='Cc solubility / mol/kgw')))
rgl::plot3d(x = df$temp, y = df$Na, z = -res3[, "d_Calcite"],
xlab='Temperature / °C', ylab='NaCl / molal', zlab='Calcite solubility / mol/kgw', type = "p")
# (4) A globe
lat <- matrix(seq(90, -90, len = 50)*pi/180, 50, 50, byrow = TRUE)
long <- matrix(seq(-180, 180, len = 50)*pi/180, 50, 50)
r <- 6378.1 # radius of Earth in km
x <- r*cos(lat)*cos(long)
y <- r*cos(lat)*sin(long)
z <- r*sin(lat)
library(rgl)
open3d()
persp3d(x, y, z, col = "white",
texture = system.file("textures/worldsmall.png", package = "rgl"),
specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab = "", zlab = "",
normal_x = x, normal_y = y, normal_z = z)
if (!rgl.useNULL())
play3d(spin3d(axis = c(0, 0, 1), rpm = 16), duration = 2.5)
## Not run:
# This looks much better, but is slow because the texture is very big
persp3d(x, y, z, col = "white",
texture = system.file("textures/world.png", package = "rgl"),
specular = "black", axes = FALSE, box = FALSE, xlab = "", ylab = "", zlab = "",
normal_x = x, normal_y = y, normal_z = z)
plot3d(mtcars$mpg,mtcars$wt,mtcars$disp, col="red", size=3)
## visualize
demo("demo-kinetics", package = "RedModRphree")
MatplotSingle(sim=simk, xlab="Domain element", ylab="Ca, Mg [molal], Calcite, Dolomite [mol]", main="1D simulation with grid 50")
### was
library(RedModRphree)
db <- RPhreeFile(system.file("extdata", "phreeqc_kin.dat", package="RedModRphree"), is.db=TRUE)
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",