[CHNOSZ-commits] r575 - in pkg/CHNOSZ: . R man tests/testthat vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Jul 24 11:38:48 CEST 2020
Author: jedick
Date: 2020-07-24 11:38:47 +0200 (Fri, 24 Jul 2020)
New Revision: 575
Added:
pkg/CHNOSZ/tests/testthat/test-mix.R
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/mix.R
pkg/CHNOSZ/man/diagram.Rd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Add test-mix.R
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2020-07-24 00:40:29 UTC (rev 574)
+++ pkg/CHNOSZ/DESCRIPTION 2020-07-24 09:38:47 UTC (rev 575)
@@ -1,6 +1,6 @@
Date: 2020-07-24
Package: CHNOSZ
-Version: 1.3.6-48
+Version: 1.3.6-49
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/R/diagram.R
===================================================================
--- pkg/CHNOSZ/R/diagram.R 2020-07-24 00:40:29 UTC (rev 574)
+++ pkg/CHNOSZ/R/diagram.R 2020-07-24 09:38:47 UTC (rev 575)
@@ -719,9 +719,23 @@
out2D <- list(namesx=pn$namesx, namesy=pn$namesy)
} # end if(nd==2)
} # end if(plot.it)
+
# even if plot=FALSE, return the diagram clipped to the water stability region (for testing) 20200719
if(isTRUE(limit.water) & !is.null(H2O.predominant)) predominant <- H2O.predominant
- outstuff <- list(plotvar=plotvar, plotvals=plotvals, names=names, predominant=predominant)
+ # make a matrix with the affinities of predominant species 20200724
+ # (for calculating affinities of metastable species - multi-metal.Rmd example)
+ predominant.values <- NA
+ if(!identical(predominant, NA)) {
+ predominant.values <- eout$values[[1]]
+ predominant.values[] <- NA
+ for(ip in na.omit(unique(as.vector(predominant)))) {
+ ipp <- predominant == ip
+ ipp[is.na(ipp)] <- FALSE
+ predominant.values[ipp] <- eout$values[[ip]][ipp]
+ }
+ }
+
+ outstuff <- list(plotvar = plotvar, plotvals = plotvals, names = names, predominant = predominant, predominant.values = predominant.values)
# include the balance name and coefficients if we diagrammed affinities 20200714
if(eout.is.aout) outstuff <- c(list(balance = balance, n.balance = n.balance), outstuff)
out <- c(eout, outstuff)
Modified: pkg/CHNOSZ/R/mix.R
===================================================================
--- pkg/CHNOSZ/R/mix.R 2020-07-24 00:40:29 UTC (rev 574)
+++ pkg/CHNOSZ/R/mix.R 2020-07-24 09:38:47 UTC (rev 575)
@@ -82,14 +82,26 @@
}
}
# Solve for the mole fractions of each species that give the required mixture
- frac1 <- frac2 <- numeric()
+ isingular <- frac1 <- frac2 <- numeric()
for(i in 1:nrow(combs)) {
# The stoichiometric matrix
A <- matrix(as.numeric(c(s1[i, ibal], s2[i, ibal])), nrow = 2, byrow = TRUE)
- x <- solve(t(A), parts)
- frac1 <- c(frac1, x[1])
- frac2 <- c(frac2, x[2])
+ x <- tryCatch(solve(t(A), parts), error = function(e) e)
+ # Check if we have a singular combination (e.g. FeV and FeVO4)
+ if(inherits(x, "error")) {
+ isingular <- c(isingular, i)
+ } else {
+ frac1 <- c(frac1, x[1])
+ frac2 <- c(frac2, x[2])
+ }
}
+ if(length(isingular) > 0) {
+ if(length(isingular) > 1) stxt <- "s" else stxt <- ""
+ message(paste0("mix: removing ", length(isingular), " combination", stxt, " with a singular stoichiometric matrix"))
+ combs <- combs[-isingular, ]
+ s1 <- s1[-isingular, ]
+ s2 <- s2[-isingular, ]
+ }
# Note that some of frac1 or frac2 might be < 0 ... we remove these combinations below
# Make a new species data frame starting with sum of scaled formation reactions
Modified: pkg/CHNOSZ/man/diagram.Rd
===================================================================
--- pkg/CHNOSZ/man/diagram.Rd 2020-07-24 00:40:29 UTC (rev 574)
+++ pkg/CHNOSZ/man/diagram.Rd 2020-07-24 09:38:47 UTC (rev 575)
@@ -111,7 +111,7 @@
A new plot is started unless \code{add} is TRUE.
If \code{plot.it} is FALSE, no plot will be generated but all the intermediate computations will be performed and the results returned.
-Line or field labels use the names of the species as provided in \code{eout}; formatting is applied to chemical formulas unless \code{format.names} is FALSE.
+Line or field labels use the names of the species as provided in \code{eout}; formatting is applied to chemical formulas (see \code{\link{expr.species}}) unless \code{format.names} is FALSE.
Set \code{names} to TRUE or NULL to plot the names, or FALSE, NA, or \code{""} to prevent plotting the names, or a character argument to replace the default species names.
Alternatively, supply a numeric value to \code{names} to indicate a subset of default names that should or shouldn't be plotted (positive and negative indices, respectively).
Use \code{col.names} and \code{cex.names} to change the colors and size of the labels.
@@ -173,12 +173,13 @@
}
\section{Value}{
- \code{diagram} returns an \code{\link{invisible}} list containing, first, the contents of \code{eout}, i.e. the provided output of \code{\link{affinity}} or \code{\link{equilibrate}}.
- To this are added the name of the plotted variable in \code{plotvar}, the plotted values in \code{plotvals}, and the names used for labeling the plot in \code{names}.
+ \code{diagram} returns an \code{\link{invisible}} list containing, first, the contents of \code{eout}, i.e. the output of \code{\link{affinity}} or \code{\link{equilibrate}} supplied in the function call.
+ To this are added the names of the plotted variable in \code{plotvar}, the labels used for species (which may be \code{\link{plotmath}} expressions if \code{format.names} is TRUE) in \code{names}, and the values used for plotting in a list named \code{plotvals}.
For 1-D diagrams, \code{plotvals} usually corresponds to the chemical activities of the species (i.e. \code{eout$loga.equil}), or, if \code{alpha} is \code{TRUE}, their mole fractions (degrees of formation).
- For 2-D diagrams, the output also contains \code{predominant}, giving the numbers (from the \code{\link{species}} definition) of the predominant (aka maximum-affinity) species at each grid point.
- The rows and columns of \code{predominant} correspond to the x- and y-variables, respectively.
- Finally, the output for 2-D diagrams contains a \code{lines} component, giving the x- and y-coordinates of the field boundaries computed using \code{\link{contourLines}}; the values are padded to equal length with NAs to faciliate exporting the results using \code{\link{write.csv}}.
+ For 2-D diagrams, \code{plotvals} corresponds to the values of affinity (from \code{eout$values}) divided by the respective balancing coefficients for each species.
+ For 2-D diagrams, the output also contains the matrices \code{predominant}, which identifies the predominant species in \code{eout$species} at each grid point, and \code{predominant.values}, which has the affinities of the predominant species.
+ The rows and columns of these matrices correspond to the \emph{x} and \emph{y} variables, respectively.
+ Finally, the output for 2-D diagrams contains a \code{lines} component, giving the \emph{x} and \emph{y} coordinates of the field boundaries computed using \code{\link{contourLines}}; the values are padded to equal length with NAs to faciliate exporting the results using \code{\link{write.csv}}.
}
\seealso{
Added: pkg/CHNOSZ/tests/testthat/test-mix.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-mix.R (rev 0)
+++ pkg/CHNOSZ/tests/testthat/test-mix.R 2020-07-24 09:38:47 UTC (rev 575)
@@ -0,0 +1,66 @@
+context("mix")
+
+test_that("mix() and mosaic() yield same affinities for a bimetallic mineral", {
+
+ # Test that we get the same values for affinity of a bimetallic mineral formed from different single-metal basis species using:
+ # 1. affinity() of the bimetallic mineral minus mix() for two single-metal systems in correct proportions
+ # 2. mosaic() with two sets of changing basis species corresponding to the single-metal systems
+ # 20200724
+
+ plot.it <- FALSE
+ res <- 50
+ ## Uncomment to show plots
+ #plot.it <- TRUE
+ #par(mfrow = c(2, 2))
+
+ # Set up system
+ pH <- c(0, 7, res)
+ Eh <- c(-1, 1, res)
+ basis(c("copper", "iron", "H2S", "H2O", "e-", "H+"))
+ iFe.cr <- info(c("iron", "ferrous-oxide", "magnetite", "hematite"))
+ iFe.aq <- info(c("Fe+2", "FeOH+"))
+ iCu.cr <- info(c("copper", "cuprite", "tenorite"))
+ iCu.aq <- info(c("Cu+", "CuOH"))
+
+ # METHOD 1 (mix)
+
+ # Fe-O-H diagram
+ # Use logact = 0 for everything so we can compare results with mosaic()
+ species(c(iFe.aq, iFe.cr), 0)
+ # Increase the temperature so we get a ferrous oxide field
+ aFe <- affinity(pH = pH, Eh = Eh, T = 300)
+ dFe <- diagram(aFe, fill = fill(aFe), plot.it = plot.it)
+
+ # Cu-O-H diagram
+ species(c(iCu.aq, iCu.cr), 0)
+ aCu <- affinity(aFe) # argument recall
+ dCu <- diagram(aCu, fill = fill(aCu), plot.it = plot.it)
+
+ # Combine the diagrams for a 1:5 mixture of Fe:Cu
+ aFeCu15 <- mix(dFe, dCu, parts = c(1, 5))
+ dFeCu15 <- diagram(aFeCu15, fill = fill(aFeCu15), min.area = 0.01, plot.it = FALSE)
+ # Calculate affinity of bornite (Cu5FeS4)
+ species("bornite")
+ abornite <- affinity(aFe) # argument recall
+ # Because of the 1) equal stoichiometry and 2) same basis species
+ # used for bornite and the mix() diagram, subtracting these makes sense
+ abornite_vs_predominant.values <- abornite$values[[1]] - dFeCu15$predominant.values
+ if(plot.it) image(abornite_vs_predominant.values)
+
+ # METHOD 2 (mosaic)
+
+ # Make a mosaic diagram to calculate affinity of bornite with changing basis species
+ # (which correspond to the Fe-O-H and Cu-O-H diagrams)
+ iFe <- c(iFe.cr, iFe.aq)
+ iCu <- c(iCu.cr, iCu.aq)
+ # TODO: allow numeric values for bases in mosaic()
+ #mbornite <- mosaic(list(iFe, iCu), pH = pH, Eh = Eh, T = 300, predominant = list(dFe$predominant, dCu$predominant))
+ Fe <- info(iFe)$name
+ Cu <- info(iCu)$name
+ mbornite <- mosaic(list(Fe, Cu), pH = pH, Eh = Eh, T = 300, blend = FALSE)
+ if(plot.it) image(mbornite$A.species$values[[1]])
+
+ # The test: values calculated both ways should be equal
+ expect_equal(abornite_vs_predominant.values, mbornite$A.species$values[[1]], tol = 1e-5, scale = 1)
+
+})
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-24 00:40:29 UTC (rev 574)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2020-07-24 09:38:47 UTC (rev 575)
@@ -155,7 +155,6 @@
"FeO2-" = -368.2, # SSWS97
"FeO4-2" = -322.2 # Mis73
)
-Fe.aq <- 1000 * c("Fe+2" = -78.90, "FeOH+" = -277.4)
V.aq <- 1000 * c("VO+2" = -446.4, "VO2+" = -587.0, "VO3-" = -783.6, "VO4-3" = -899.0,
"V2O7-4" = -1719, "HVO4" = -745.1, "HVO4-2" = -974.8,
@@ -162,7 +161,6 @@
"VOH2O2+3" = -523.4, "VO2H2O2+" = -746.3, "V2HO7-3" = 1792.2, "V2H3O7-" = -1863.8,
"HV10O28-5" = -7702, "H2V10O28-4" = -7723
)
-V.aq <- 1000 * c("VO+2" = -446.4, "HVO4-2" = -974.8)
# Gibbs energies of formation (J / mol) for solids from Wagman et al., 1982
Fe3O4 <- 1000 * -1015.4 # magnetite
More information about the CHNOSZ-commits
mailing list