[CHNOSZ-commits] r218 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Tue Sep 26 09:06:47 CEST 2017
Author: jedick
Date: 2017-09-26 09:06:46 +0200 (Tue, 26 Sep 2017)
New Revision: 218
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/DEW.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/DEW.Rd
pkg/CHNOSZ/man/eos.Rd
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/man/util.water.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-DEW.R
pkg/CHNOSZ/tests/testthat/test-EOSregress.R
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
enable g-function for DEW and IAPWS-95
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-26 07:06:46 UTC (rev 218)
@@ -1,6 +1,6 @@
-Date: 2017-09-25
+Date: 2017-09-26
Package: CHNOSZ
-Version: 1.1.0-16
+Version: 1.1.0-17
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/DEW.R
===================================================================
--- pkg/CHNOSZ/R/DEW.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/R/DEW.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -199,3 +199,49 @@
A * exp(B) * density ^ (A - 1)
}
+
+### testing functions ###
+# These unexported functions are included for testing purposes only.
+# In CHNOSZ, the g function and omega(P,T) are calculated via hkf().
+
+# 'Returns the value of omega at the input P and T.
+# The value returned is 'in units of cal/mol and NOT multiplied by 10^-5.
+#'pressure - Pressure to calculate at, in bars
+#'temperature- Temperature to calculate at, in Celsius
+#'density - Density of water to calculate omega at, in g/cm^3.
+#'wref - The value of omega at standard pressure and temperature, in units of cal/mol.
+#'Z - The charge of the species
+calculateOmega <- function(pressure, temperature, density, wref, Z) {
+ # 'These equations are given by Shock et al. (1992)
+ eta <- 166027 # 'Value in units of Angstroms cal mol^-1
+ # 'Defines the electrostatic radius at reference pressure and temperature
+ reref <- Z * Z / (wref / eta + Z / 3.082)
+ # 'This represents the pressure and temperature dependent solvent function
+ g <- calculateG(pressure, temperature, density)
+ # 'Defines the electrostatic radius at the input P and T
+ re <- reref + abs(Z) * g
+ omega <- eta * (Z * Z / re - Z / (3.082 + g))
+ # 'If species is hydrogen, the species is neutral, or the pressure is above 6 kb,
+ # 'this equation is not necessary because omega is very close to wref.
+ if(Z==0) omega[] <- wref
+ omega[pressure > 6000] <- wref
+}
+
+# 'Returns the value of the g function. If the density is greater than 1 g/cm^3, then zero is returned.
+# 'pressure - The pressure to calculate at, in bars
+# 'temperature- The temperature to calculate at, in celsius
+# 'density - The density of water at which to calculate g at, in g/cm^3
+calculateG <- function(pressure, temperature, density) {
+ T <- temperature
+ P <- pressure
+ a_g <- -2.037662 + 0.005747 * T - 6.557892E-06 * T * T
+ b_g <- 6.107361 - 0.01074377 * T + 1.268348E-05 * T * T
+ # 'Calculates the difference function in the case where we need to calculate at Psat conditions
+ f <- (((T - 155) / 300)^4.8 + 36.66666 * ((T - 155) / 300)^16) *
+ (-1.504956E-10 * (1000 - P)^3 + 5.017997E-14 * (1000 - P)^4)
+ f[P > 1000 | T < 155 | T > 355] <- 0
+ g <- a_g * (1 - density)^b_g - f
+ # use g = 0 for density >= 1
+ g[density >= 1] <- 0
+ g
+}
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -30,7 +30,7 @@
notcontrib <- ! contrib %in% c('n', 's', 'o')
if(TRUE %in% notcontrib) stop(paste("contrib must be in c('n', 's', 'o'); got", c2s(contrib[notcontrib])))
# get water properties
- # rho - for subcrt() output and derivatives of omega (if needed)
+ # rho - for subcrt() output and g function
# Born functions and epsilon - for HKF calculations
H2O.props <- c(H2O.props, "QBorn", "XBorn", "YBorn", "epsilon")
if(grepl("SUPCRT", thermo$opt$water)) {
@@ -74,7 +74,7 @@
dwdP <- dwdT <- d2wdT2 <- numeric(length(T))
Z <- PAR$Z
omega.PT <- rep(PAR$omega, length(T))
- if(!identical(Z, 0) & !identical(PAR$name, "H+") & grepl("SUPCRT", thermo$opt$water)) {
+ if(!identical(Z, 0) & !identical(PAR$name, "H+")) {
# compute derivatives of omega: g and f functions (Shock et al., 1992; Johnson et al., 1992)
rhohat <- H2O.PT$rho/1000 # just converting kg/m3 to g/cm3
g <- gfun(rhohat, convert(T, "C"), P, H2O.PT$alpha, H2O.PT$daldT, H2O.PT$beta)
@@ -234,46 +234,54 @@
ifg <- ifg & !is.na(ifg)
# Eq. 32
g[ifg] <- g[ifg] - f[ifg]
+ # at P > 6000 bar (in DEW calculations), g is zero 20170926
+ g[P > 6000] <- 0
## now we have g at P, T
+ # put the results in their right place (where rhohat < 1)
+ out$g[idoit] <- g
## the rest is to get its partial derivatives with pressure and temperature
## after Johnson et al., 1992
# alpha - coefficient of isobaric expansivity (K^-1)
# daldT - temperature derivative of coefficient of isobaric expansivity (K^-2)
# beta - coefficient of isothermal compressibility (bar^-1)
- # Eqn. 76
- d2fdT2 <- (0.0608/300*((Tc-155)/300)^2.8 + af1/375*((Tc-155)/300)^14) * (af2*(1000-P)^3 + af3*(1000-P)^4)
- # Eqn. 75
- dfdT <- (0.016*((Tc-155)/300)^3.8 + 16*af1/300*((Tc-155)/300)^15) * (af2*(1000-P)^3 + af3*(1000-P)^4)
- # Eqn. 74
- dfdP <- -(((Tc-155)/300)^4.8 + af1*((Tc-155)/300)^16) * (3*af2*(1000-P)^2 + 4*af3*(1000-P)^3)
- d2bdT2 <- 2 * bg3 # Eqn. 73
- d2adT2 <- 2 * ag3 # Eqn. 72
- dbdT <- bg2 + 2*bg3*Tc # Eqn. 71
- dadT <- ag2 + 2*ag3*Tc # Eqn. 70
- # Eqn. 69
- dgadT <- bg*rhohat*alpha*(1-rhohat)^(bg-1) + log(1-rhohat)*g/ag*dbdT
- D <- rhohat
- # transcribed from SUPCRT92/reac92.f
- dDdT <- -D * alpha
- dDdP <- D * beta
- dDdTT <- -D * (daldT - alpha^2)
- Db <- (1-D)^bg
- dDbdT <- -bg*(1-D)^(bg-1)*dDdT + log(1-D)*Db*dbdT
- dDbdTT <- -(bg*(1-D)^(bg-1)*dDdTT + (1-D)^(bg-1)*dDdT*dbdT +
- bg*dDdT*(-(bg-1)*(1-D)^(bg-2)*dDdT + log(1-D)*(1-D)^(bg-1)*dbdT)) +
- log(1-D)*(1-D)^bg*d2bdT2 - (1-D)^bg*dbdT*dDdT/(1-D) + log(1-D)*dbdT*dDbdT
- d2gdT2 <- ag*dDbdTT + 2*dDbdT*dadT + Db*d2adT2
- d2gdT2[ifg] <- d2gdT2[ifg] - d2fdT2[ifg]
- dgdT <- g/ag*dadT + ag*dgadT # Eqn. 67
- dgdT[ifg] <- dgdT[ifg] - dfdT[ifg]
- dgdP <- -bg*rhohat*beta*g*(1-rhohat)^-1 # Eqn. 66
- dgdP[ifg] <- dgdP[ifg] - dfdP[ifg]
- # phew! done with those derivatives
- # put the results in their right place (where rhohat < 1)
- out$g[idoit] <- g
- out$dgdT[idoit] <- dgdT
- out$d2gdT2[idoit] <- d2gdT2
- out$dgdP[idoit] <- dgdP
+ # if these are NULL or NA (for IAPWS-95 and DEW), we skip the calculation
+ if(is.null(alpha)) alpha <- NA
+ if(is.null(daldT)) daldT <- NA
+ if(is.null(beta)) beta <- NA
+ if(!all(is.na(alpha)) & !all(is.na(daldT)) & !all(is.na(beta))) {
+ # Eqn. 76
+ d2fdT2 <- (0.0608/300*((Tc-155)/300)^2.8 + af1/375*((Tc-155)/300)^14) * (af2*(1000-P)^3 + af3*(1000-P)^4)
+ # Eqn. 75
+ dfdT <- (0.016*((Tc-155)/300)^3.8 + 16*af1/300*((Tc-155)/300)^15) * (af2*(1000-P)^3 + af3*(1000-P)^4)
+ # Eqn. 74
+ dfdP <- -(((Tc-155)/300)^4.8 + af1*((Tc-155)/300)^16) * (3*af2*(1000-P)^2 + 4*af3*(1000-P)^3)
+ d2bdT2 <- 2 * bg3 # Eqn. 73
+ d2adT2 <- 2 * ag3 # Eqn. 72
+ dbdT <- bg2 + 2*bg3*Tc # Eqn. 71
+ dadT <- ag2 + 2*ag3*Tc # Eqn. 70
+ # Eqn. 69
+ dgadT <- bg*rhohat*alpha*(1-rhohat)^(bg-1) + log(1-rhohat)*g/ag*dbdT
+ D <- rhohat
+ # transcribed from SUPCRT92/reac92.f
+ dDdT <- -D * alpha
+ dDdP <- D * beta
+ dDdTT <- -D * (daldT - alpha^2)
+ Db <- (1-D)^bg
+ dDbdT <- -bg*(1-D)^(bg-1)*dDdT + log(1-D)*Db*dbdT
+ dDbdTT <- -(bg*(1-D)^(bg-1)*dDdTT + (1-D)^(bg-1)*dDdT*dbdT +
+ bg*dDdT*(-(bg-1)*(1-D)^(bg-2)*dDdT + log(1-D)*(1-D)^(bg-1)*dbdT)) +
+ log(1-D)*(1-D)^bg*d2bdT2 - (1-D)^bg*dbdT*dDdT/(1-D) + log(1-D)*dbdT*dDbdT
+ d2gdT2 <- ag*dDbdTT + 2*dDbdT*dadT + Db*d2adT2
+ d2gdT2[ifg] <- d2gdT2[ifg] - d2fdT2[ifg]
+ dgdT <- g/ag*dadT + ag*dgadT # Eqn. 67
+ dgdT[ifg] <- dgdT[ifg] - dfdT[ifg]
+ dgdP <- -bg*rhohat*beta*g*(1-rhohat)^-1 # Eqn. 66
+ dgdP[ifg] <- dgdP[ifg] - dfdP[ifg]
+ # phew! done with those derivatives
+ out$dgdT[idoit] <- dgdT
+ out$d2gdT2[idoit] <- d2gdT2
+ out$dgdP[idoit] <- dgdP
+ }
return(out)
}
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/R/water.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -9,10 +9,11 @@
# set water option
if(length(property)==1 & any(property %in% c("SUPCRT", "SUPCRT92", "IAPWS", "IAPWS95", "DEW"))) {
thermo <- get("thermo")
+ owater <- thermo$opt$water
thermo$opt$water <- property
assign("thermo", thermo, "CHNOSZ")
- message(paste("water: thermo$opt$water was set to", property))
- return(invisible(property))
+ message(paste("water: setting thermo$opt$water to", property))
+ return(invisible(owater))
}
# make T and P equal length
if(!identical(P, "Psat")) {
@@ -23,15 +24,12 @@
T[T == 273.15] <- 273.16
wopt <- get("thermo")$opt$water
if(grepl("SUPCRT", wopt)) {
- # get the values of the properties using SUPCRT92
+ # get properties using SUPCRT92
w.out <- water.SUPCRT92(property, T, P)
}
if(grepl("IAPWS", wopt)) {
- # here we get properties using IAPWS-95
+ # get properties using IAPWS-95
w.out <- water.IAPWS95(property, T, P)
- # normalize the names to use upper case (expected by subcrt())
- iprop <- match(tolower(property), tolower(water.props("IAPWS95")))
- colnames(w.out) <- water.props("IAPWS95")[iprop]
}
if(grepl("DEW", wopt)) {
# use the Deep Earth Water (DEW) model
@@ -291,7 +289,7 @@
return(p)
}
### main loop; init dataframe output and density holders
- w.out <- NULL
+ w.out <- data.frame(matrix(nrow=length(T), ncol=length(property)))
my.rho <- NULL
# get densities unless only Psat is requested
if(!identical(tolower(property), "psat")) {
@@ -302,18 +300,23 @@
rho <- function() my.rho
}
for(i in 1:length(property)) {
- if(tolower(property[i]) %in% c("e", "kt")) {
+ if(tolower(property[i]) %in% c("e", "kt", "alpha", "daldt", "beta")) {
# E and kT aren't here yet... set them to NA
- warning("water.IAPWS95: values of ", property[i], " are NA\n", call.=FALSE)
- inew <- rep(NA, length(T))
+ # also set alpha, daldT, and beta (for derivatives of g function) to NA 20170926
+ warning("water.IAPWS95: values of ", property[i], " are NA", call.=FALSE)
+ wnew <- rep(NA, length(T))
} else {
message(paste(" ", property[i], sep=""), appendLF=FALSE)
- inew <- get(tolower(property[i]))()
+ wnew <- get(tolower(property[i]))()
}
- wnew <- data.frame(inew)
- if(i > 1) w.out <- cbind(w.out, wnew) else w.out <- wnew
+ w.out[, i] <- wnew
}
message("")
+ # use uppercase property names (including properties available in SUPCRT that might be NA here)
+ wprop <- unique(c(water.props("SUPCRT"), water.props("IAPWS")))
+ iprop <- match(tolower(property), tolower(wprop))
+ property[!is.na(iprop)] <- wprop[na.omit(iprop)]
+ colnames(w.out) <- property
return(w.out)
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-26 07:06:46 UTC (rev 218)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-16 (2017-09-25)
+CHANGES IN CHNOSZ 1.1.0-17 (2017-09-26)
---------------------------------------
MAJOR CHANGES:
@@ -49,6 +49,9 @@
- In water.* functions, rename the "diel" variable to "epsilon".
+- calculations of the g function (but not its derivatives) are now
+ enabled for DEW and IAPWS-95.
+
CHANGES IN CHNOSZ 1.1.0 (2017-05-04)
------------------------------------
Modified: pkg/CHNOSZ/man/DEW.Rd
===================================================================
--- pkg/CHNOSZ/man/DEW.Rd 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/man/DEW.Rd 2017-09-26 07:06:46 UTC (rev 218)
@@ -38,7 +38,7 @@
}
\seealso{
-\code{\link{water.DEW}}
+\code{\link{water.DEW}}; use \code{water("DEW")} to activate these equations for the main functions in CHNOSZ.
}
\examples{
Modified: pkg/CHNOSZ/man/eos.Rd
===================================================================
--- pkg/CHNOSZ/man/eos.Rd 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/man/eos.Rd 2017-09-26 07:06:46 UTC (rev 218)
@@ -34,7 +34,7 @@
\code{H2O.props} is an optional argument that lists the properties of water that should be returned; it is used by \code{\link{subcrt}} so that the time-consuming water calculations are only performed once.
The temperature and pressure derivatives of the \code{omega} parameter for charged species (where \code{Z != 0}, but not for the aqueous proton, H+) are calculated using the \eqn{g}{g}- and \eqn{f}{f}-functions described by Shock et al., 1992 and Johnson et al., 1992.
-This option is turned off when the IAPWS-95 equations are activated (see \code{\link{water}}).
+If the IAPWS-95 or DEW equations are activated (see \code{\link{water}}), only the \eqn{g}{g}-function (applicable to \samp{G}), but not its derivatives (needed for \samp{H}, \samp{S}, \samp{Cp}, and \samp{V}), is calculated.
The parameters in the \code{cgl} equations of state for crystalline, gas and liquid species (except liquid water) include \code{V}, \code{a}, \code{b}, \code{c}, \code{d}, \code{e}, \code{f} and \code{lambda}.
The terms denoted by \code{a}, \code{b} and \code{c} correspond to the Maier-Kelley equation for heat capacity (Maier and Kelley, 1932); the additional terms are useful for representing heat capacities of minerals (Robie and Hemingway, 1995) and gaseous or liquid organic species (Helgeson et al., 1998).
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/man/subcrt.Rd 2017-09-26 07:06:46 UTC (rev 218)
@@ -262,7 +262,7 @@
## Subzero (degrees C) calculations
# uncomment the following to try IAPWS95 instead of SUPCRT92
-water("IAPWS95")
+#water("IAPWS95")
# the limit for H2O92D.f (from SUPCRT92) is currently -20 deg C
# but we go to -30 knowing properties will become NA
sb <- subcrt(c("H2O", "Na+"), T=seq(-30, 10), P=1)$out
Modified: pkg/CHNOSZ/man/util.water.Rd
===================================================================
--- pkg/CHNOSZ/man/util.water.Rd 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/man/util.water.Rd 2017-09-26 07:06:46 UTC (rev 218)
@@ -63,12 +63,13 @@
water.IAPWS95("V", T, "Psat")
# the calculated Psat, while not exact, should be close enough for most
# geochemical calculations to select liquid or vapor
-water("IAPWS95")
+ow <- water("IAPWS95")
thermo$opt$IAPWS.sat <<- "vapor"
V.vapor <- subcrt("water", T=convert(445:455, "C"))$out[[1]]$V
thermo$opt$IAPWS.sat <<- "liquid" # the default
V.liquid <- subcrt("water", T=convert(445:455, "C"))$out[[1]]$V
stopifnot(all(V.vapor > V.liquid))
+water(ow)
}
\references{
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/man/water.Rd 2017-09-26 07:06:46 UTC (rev 218)
@@ -40,7 +40,7 @@
}
-Use e.g. \code{water("DEW")} to set the option.
+Use e.g. \code{water("DEW")} to set the option; the previous option (at the time of the function call) is returned invisibly.
Subsequent calculations with \code{water}, or other functions such as \code{subcrt} and \code{affinity}, will use that option.
The allowed \code{property}s for \code{water} are one or more of those given below, depending on the computational option; availability is shown by an asterisk.
@@ -132,15 +132,14 @@
## comparing the formulations
T <- convert(c(25, 100, 200, 300), "K")
# IAPWS-95
-water("IAPWS95")
+ow <- water("IAPWS95")
water(water.props(), T=T)
# SUPCRT92 (the default)
-water("SUPCRT92")
+water(ow)
water(water.props(), T=T)
## calculating Q Born function
# after Table 22 of Johnson and Norton, 1991
-water("SUPCRT92")
T <- rep(c(375, 400, 425, 450, 475), each=5)
P <- rep(c(250, 300, 350, 400, 450), 5)
w <- water("QBorn", T=convert(T, "K"), P=P)
Modified: pkg/CHNOSZ/tests/testthat/test-DEW.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-DEW.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/tests/testthat/test-DEW.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -1,9 +1,13 @@
context("DEW")
+# since this is alphabetically the first test,
+# we need to load the 'thermo' object (for running tests in R CMD check)
+suppressMessages(data(thermo))
+
test_that("density of water is calculated correctly", {
pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
- # density from R functions
+ # density from R translation of DEW macro functions
RDensity <- calculateDensity(pressure, temperature)
# density from DEW spreadsheet
DEWDensity <- c(1.108200, 0.597623, 1.196591, 0.798331, 1.321050, 1.000735, 1.578116, 1.287663)
@@ -13,7 +17,7 @@
test_that("Gibbs energy of water is calculated correctly", {
pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
- # Gibbs energies from R functions
+ # Gibbs energies from R translation of DEW macro functions
RGibbs <- calculateGibbsOfWater(pressure, temperature)
# Gibbs energies from DEW spreadsheet
DEWGibbs <- c(-56019.85419280258, -84262.028821198, -54155.004480575895, -81210.38766217149,
@@ -24,7 +28,7 @@
test_that("dielectric constant of water is calculated correctly", {
pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
- # epsilon from R functions
+ # epsilon from R translation of DEW macro functions
Repsilon <- calculateEpsilon(calculateDensity(pressure, temperature), temperature)
# epsilon from DEW spreadsheet
DEWepsilon <- c(65.63571, 6.10465, 72.40050, 8.97800, 82.16244, 12.13131, 103.12897, 16.97266)
@@ -34,7 +38,7 @@
test_that("Born coefficient Q is calculated correctly", {
pressure <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
temperature <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
- # Q from R functions
+ # Q from R translation of DEW macro functions
RQ <- calculateQ(calculateDensity(pressure, temperature), temperature)
# Q from DEW spreadsheet
DEWQ <- c(0.32319817, 14.50286092, 0.19453478, 3.12650897,
@@ -42,7 +46,30 @@
expect_equal(RQ, DEWQ)
})
+test_that("g function is calculated correctly", {
+ pressure <- c(1000, 1000, 5000, 5000, 10000)
+ temperature <- c(100, 1000, 100, 1000, 100)
+ # properties of water from DEW implementation in CHNOSZ
+ ow <- water("DEW")
+ # note that values returned for alpha, daldT, beta are NA
+ w <- water(c("rho", "alpha", "daldT", "beta", "Psat"), T=convert(temperature, "K"), P=pressure)
+ # g from CHNOSZ functions
+ Rg <- gfun(w$rho/1000, temperature, pressure, w$alpha, w$daldT, w$beta)$g
+ # g from R translation of DEW macro functions (not used in CHNOSZ)
+ DEWg <- calculateG(pressure, temperature, w$rho/1000)
+ expect_equal(Rg, DEWg)
+ water(ow)
+})
+
test_that("Gibbs energies of species are calculated correctly", {
+ P <- c(5000, 5000, 10000, 10000, 20000, 20000, 50000, 50000)
+ T <- c(100, 1000, 100, 1000, 100, 1000, 100, 1000)
+ ow <- water("DEW")
+ RG_HCl <- subcrt("HCl", P=P, T=T)$out[[1]]$G
DEWG_HCl <- c(-28784.99, -58496.85, -26520.94, -55276.92, -21928.89, -50337.19, -8014.34, -36746.87)
+ # TODO: expect_equal(RG_HCl, DEWG_HCl)
+ RG_Cl <- subcrt("Cl-", P=P, T=T)$out[[1]]$G
DEWG_Cl <- c(-30054.59, -22839.35, -27910.68, -28094.07, -23568.45, -27959.67, -10443.07, -18744.93)
+ # TODO: expect_equal(RG_Cl, DEWG_Cl)
+ water(ow)
})
Modified: pkg/CHNOSZ/tests/testthat/test-EOSregress.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-EOSregress.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/tests/testthat/test-EOSregress.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -1,9 +1,5 @@
context("EOSregress")
-# since this is the first test (alphabetically)
-# we need to load the 'thermo' object (for R CMD check)
-suppressMessages(data(thermo))
-
test_that("EOSvar stops with unknown variables", {
expect_error(EOSvar("TX", T=25, P=1), "can't find a variable named TX")
# why can't the test find these?
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-25 13:31:36 UTC (rev 217)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-26 07:06:46 UTC (rev 218)
@@ -28,48 +28,49 @@
})
test_that("standard Gibbs energies of reactions involving aqueous species are consistent with the literature", {
- # choose units of Joules for subcrt()
- E.units("J")
- # from Amend and Shock, 2001 Table 3
+
+ # from Amend and Shock, 2001 [AS01] Table 3
T <- c(2, 18, 25, 37, 45, 55, 70, 85, 100, 115, 150, 200)
# standard Gibbs energies in kJ/mol
- DG0.H2O <- c(78.25, 79.34, 79.89, 80.90, 81.63, 82.59, 84.13, 85.78, 87.55, 89.42, 94.22, 102.21)
+ E.units("J")
+
# H2O(l) = H+ + OH-
+ AS01.H2O <- c(78.25, 79.34, 79.89, 80.90, 81.63, 82.59, 84.13, 85.78, 87.55, 89.42, 94.22, 102.21)
sout.H2O <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T=T)$out
- # from Amend and Shock, 2001 Table 4.3 Reaction A1
- DG0.A1 <- c(-263.94, -263.45, -263.17, -262.62, -262.20, -261.63, -260.67, -259.60, -258.44, -257.18, -253.90, -248.44)
- # H2(aq) + 0.5O2(aq) = H2O(l)
+ # tolerances set to lowest order of magnitute to pass
+ expect_equal(sout.H2O$G/1000, AS01.H2O, tolerance=1e-4)
+
+ # AS01 Table 4.3 Reaction A1: H2(aq) + 0.5O2(aq) = H2O(l)
+ AS01.A1 <- c(-263.94, -263.45, -263.17, -262.62, -262.20, -261.63, -260.67, -259.60, -258.44, -257.18, -253.90, -248.44)
sout.A1 <- subcrt(c("H2", "O2", "H2O"), "aq", c(-1, -0.5, 1), T=T)$out
- # from Amend and Shock, 2001 Table 5.1 NO(g) and NO(aq)
+ expect_equal(sout.A1$G/1000, AS01.A1, tolerance=1e-4)
+
+ # AS01 Table 5.1 NO(g) -> NO(aq)
DG0.NO.g <- c(91.39, 88.04, 86.57, 84.03, 82.33, 80.20, 76.99, 73.75, 70.50, 67.23, 59.53, 48.38)
DG0.NO.aq <- c(104.56, 102.87, 102.06, 100.58, 99.53, 98.15, 95.98, 93.67, 91.25, 88.72, 82.46, 72.77)
- # NO(g) = NO(aq) ... greater tolerance, our values for NO(aq) differ slightly from AS01
sout.NO <- subcrt("NO", c("gas", "aq"), c(-1, 1), T=T)$out
- # from Amend and Shock, 2001 Table 5.1 Reaction B10
- DG0.B10 <- c(-268.85, -268.01, -267.50, -266.46, -265.66, -264.55, -262.66, -260.54, -258.19, -255.63, -248.81, -237.06)
- # NH3(aq) + 1.5O2(aq) = H+ + NO2- + H2O(l)
+ # higher tolerance, our values for NO(aq) differ slightly from AS01
+ expect_equal(sout.NO$G/1000, DG0.NO.aq - DG0.NO.g, tolerance=1e-2)
+
+ # AS01 Table 5.1 Reaction B10: NH3(aq) + 1.5O2(aq) = H+ + NO2- + H2O(l)
+ AS01.B10 <- c(-268.85, -268.01, -267.50, -266.46, -265.66, -264.55, -262.66, -260.54, -258.19, -255.63, -248.81, -237.06)
sout.B10 <- subcrt(c("NH3", "O2", "H+", "NO2-", "H2O"), "aq", c(-1, -1.5, 1, 1, 1), T=T)$out
- # from Amend and Shock, 2001 Table 6.3 Reaction C7
- DG0.C7 <- c(-1695.30, -1686.90, -1682.80, -1675.30, -1670.00, -1663.10, -1652.00, -1640.30, -1628.00, -1615.20, -1583.50, -1533.00)
- # 5S2O3-2 + H2O(l) + 4O2(aq) = 6SO4-2 + 2H+ + 4S(s)
+ expect_equal(sout.B10$G/1000, AS01.B10, tolerance=1e-3)
+
+ # AS01 Table 6.3 Reaction C7: 5S2O3-2 + H2O(l) + 4O2(aq) = 6SO4-2 + 2H+ + 4S(s)
+ AS01.C7 <- c(-1695.30, -1686.90, -1682.80, -1675.30, -1670.00, -1663.10, -1652.00, -1640.30, -1628.00, -1615.20, -1583.50, -1533.00)
s.C7 <- subcrt(c("S2O3-2", "H2O", "O2", "SO4-2", "H+", "S"), c("aq", "liq", "aq", "aq", "aq", "cr"), c(-5, -1, -4, 6, 2, 4), T=T)
sout.C7 <- s.C7$out
- ## from Amend and Shock, 2001 Table 8.3 Reaction E12
+ expect_equal(sout.C7$G/1000, AS01.C7, tolerance=1e-4)
+ # we can also check that sulfur has expected phase transitions
+ expect_equal(s.C7$state$sulfur, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3))
+
+ ## AS01 Table 8.3 Reaction E12: 4(2-)propanol(aq) + 3CO2(aq) + 2H2O(l) = 3CH4(aq) + 4lactic acid(aq)
#DG0.E12 <- c(132.52, 132.26, 132.29, 132.49, 132.74, 133.15, 133.98, 135.04, 136.31, 137.79, 141.97, 149.53)
- ## 4(2-)propanol(aq) + 3CO2(aq) + 2H2O(l) = 3CH4(aq) + 4lactic acid(aq)
#sout.E12 <- subcrt(c("2-propanol", "CO2", "H2O", "CH4", "lactic acid"), c(-4, -3, -2, 3, 4), T=T)$out
- # now the tests, tolerances set to lowest order of magnitute to pass
- expect_equal(sout.H2O$G/1000, DG0.H2O, tolerance=1e-4)
- expect_equal(sout.A1$G/1000, DG0.A1, tolerance=1e-4)
- # greater tolerance, our values for NO(aq) differ slightly from AS01
- expect_equal(sout.NO$G/1000, DG0.NO.aq - DG0.NO.g, tolerance=1e-2)
- expect_equal(sout.B10$G/1000, DG0.B10, tolerance=1e-3)
- # we can check that sulfur has expected phase transitions
- expect_equal(s.C7$state$sulfur, c(1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 3, 3))
- expect_equal(sout.C7$G/1000, DG0.C7, tolerance=1e-4)
- # this one is on hold until the HKF parameters of 2-propanol can be located
+ ## this one is on hold until the HKF parameters of 2-propanol can be located
#expect_equal(sout.E12$G/1000, DG0.E12, 1e-4)
- # also todo: COS(g) in Amend and Helgeson, 2001 Table 7.2?
+
# return to our favourite units
E.units("cal")
})
@@ -86,12 +87,12 @@
})
test_that("calculations using IAPWS-95 are possible", {
- thermo$opt$water <<- "IAPWS95"
+ ow <- water("IAPWS95")
sb <- subcrt(c("H2O", "Na+"), T=c(-30, -20, 0, 10), P=1)$out
# the test is not a demanding numerical comparison, more that we got numbers and no error
expect_that(all(sb$`Na+`$G < sb$water$G), is_true())
# clean up
- suppressMessages(data(thermo))
+ water(ow)
})
# references
More information about the CHNOSZ-commits
mailing list