[CHNOSZ-commits] r216 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Mon Sep 25 13:40:05 CEST 2017
Author: jedick
Date: 2017-09-25 13:40:05 +0200 (Mon, 25 Sep 2017)
New Revision: 216
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/DEW.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-water.R
Log:
add DEW option to water()
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-25 11:40:05 UTC (rev 216)
@@ -1,6 +1,6 @@
Date: 2017-09-25
Package: CHNOSZ
-Version: 1.1.0-14
+Version: 1.1.0-15
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-25 11:40:05 UTC (rev 216)
@@ -2,6 +2,9 @@
# calculate thermodynamic properties using equations of state
# 11/17/03 jmd
+## if this file is interactively sourced, the following is also needed to provide unexported functions:
+#source("util.args.R")
+
hkf <- function(property = NULL, T = 298.15, P = 1, parameters = NULL,
contrib = c('n', 's', 'o'), H2O.props="rho") {
# calculate G, H, S, Cp, V, kT, and/or E using
@@ -18,6 +21,11 @@
property <- eargs$prop
EOS.props <- eargs$props
EOS.Props <- eargs$Prop
+ # make T and P equal length
+ if(!identical(P, "Psat")) {
+ if(length(P) < length(T)) P <- rep(P, length.out = length(T))
+ if(length(T) < length(P)) T <- rep(T, length.out = length(P))
+ }
# nonsolvation, solvation, and origination contribution
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])))
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/R/subcrt.R 2017-09-25 11:40:05 UTC (rev 216)
@@ -2,7 +2,7 @@
# calculate standard molal thermodynamic propertes
# 20060817 jmd
-## if we interactively source this file, the following are also needed to provide unexported functions:
+## if this file is interactively sourced, the following are also needed to provide unexported functions:
#source("util.args.R")
#source("util.character.R")
#source("info.R")
@@ -291,6 +291,7 @@
# we're not using the HKF, but still want water
H2O.PT <- water(c("rho", eosprop), T = T, P = P)
}
+
# crystalline, gas, liquid (except water) species
iscgl <- reaction$state %in% c('liq','cr','gas','cr1','cr2','cr3',
'cr4','cr5','cr6','cr7','cr8','cr9') & reaction$name != 'water'
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/R/water.R 2017-09-25 11:40:05 UTC (rev 216)
@@ -13,36 +13,41 @@
}
# turn 273.15 K to 273.16 K (needed for water.SUPCRT92 at Psat)
T[T == 273.15] <- 273.16
- # this tells us to do the calculations using code taken from SUPCRT
- do.supcrt <- get("thermo")$opt$water != "IAPWS95"
- if(do.supcrt) {
+ wopt <- get("thermo")$opt$water
+ if(grepl("SUPCRT", wopt)) {
# get the values of the properties using SUPCRT92
w.out <- water.SUPCRT92(property, T, P)
- return(w.out)
- } else {
+ }
+ if(grepl("IAPWS", wopt)) {
# here we 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]
- return(w.out)
}
+ if(grepl("DEW", wopt)) {
+ # use the Deep Earth Water (DEW) model
+ w.out <- water.DEW(property, T, P)
+ }
+ w.out
}
water.props <- function(formulation=get("thermo")$opt$water) {
# return the names of properties that are available in SUPCRT92 or IAPWS95
# added 20130212 jmd
- if(formulation=="SUPCRT92")
+ if(grepl("SUPCRT", formulation))
props <- c("A", "G", "S", "U", "H", "Cv", "Cp",
"Speed", "alpha", "beta", "epsilon", "visc",
"tcond", "surten", "tdiff", "Prndtl", "visck", "albe",
"ZBorn", "YBorn", "QBorn", "daldT", "XBorn",
"V", "rho", "Psat", "E", "kT")
- else if(formulation=="IAPWS95")
+ if(grepl("IAPWS", formulation))
props <- c("A", "G", "S", "U", "H", "Cv", "Cp",
"Speed", "epsilon",
"YBorn", "QBorn", "XBorn", "NBorn", "UBorn",
"V", "rho", "Psat", "de.dT", "de.dP", "P")
+ if(grepl("DEW", formulation))
+ props <- c("G", "epsilon", "QBorn", "V", "rho")
return(props)
}
@@ -306,6 +311,12 @@
# get water properties from DEW model for use by subcrt() 20170925
water.DEW <- function(property, T = 373.15, P = 1000) {
+ # we can't use Psat here
+ if(identical(P, "Psat")) stop("Psat isn't supported in this implementation of the DEW model")
+ # use uppercase property names (including H, S, etc., so we get them from the SUPCRT properties)
+ wprop <- water.props("SUPCRT")
+ iprop <- match(tolower(property), tolower(wprop))
+ property[!is.na(iprop)] <- wprop[na.omit(iprop)]
# convert temperature units
pressure <- P
temperature <- convert(T, "C")
@@ -317,11 +328,19 @@
# calculate rho if it's needed for any other properties
if(any(c("rho", "V", "QBorn", "epsilon") %in% property)) rho <- calculateDensity(pressure, temperature)
# fill in columns with values
- if("rho" %in% property) out$rho <- rho
+ if("rho" %in% property) out$rho <- rho*1000 # use kg/m^3 (like SUPCRT)
if("V" %in% property) out$V <- 18.01528/rho
if("G" %in% property) out$G <- calculateGibbsOfWater(pressure, temperature)
if("QBorn" %in% property) out$QBorn <- calculateQ(rho, temperature)
if("epsilon" %in% property) out$epsilon <- calculateEpsilon(rho, temperature)
+ # use SUPCRT-calculated values below 100 degC and/or below 1000 bar
+ ilow <- T < 373.15 | P < 1000
+ if(any(ilow)) {
+ out[ilow, ] <- water.SUPCRT92(property, T=T[ilow], P=P[ilow])
+ iPrTr <- T==get("thermo")$opt$Tr & P==get("thermo")$opt$Pr
+ if(sum(iPrTr)==sum(ilow)) message(paste("water.DEW: using SUPCRT calculations for Pr,Tr"))
+ if(sum(iPrTr)==0) message(paste("water.DEW: using SUPCRT calculations for", sum(ilow), "low-T or low-P condition(s)"))
+ if(sum(ilow) > sum(iPrTr)) message(paste("water.DEW: using SUPCRT calculations for Pr,Tr and", sum(ilow)-1, "other low-T or low-P condition(s)"))
+ }
out
}
-
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-25 11:40:05 UTC (rev 216)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-14 (2017-09-25)
+CHANGES IN CHNOSZ 1.1.0-15 (2017-09-25)
---------------------------------------
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
Modified: pkg/CHNOSZ/man/DEW.Rd
===================================================================
--- pkg/CHNOSZ/man/DEW.Rd 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/man/DEW.Rd 2017-09-25 11:40:05 UTC (rev 216)
@@ -30,6 +30,7 @@
The spreadsheet provides multiple options for some calculations; here the default equations for density of water (Zhang and Duan, 2005), dielectric constant (Sverjensky et al., 2014) and Gibbs energy of water (integral of volume, equation created by Brandon Harrison) are used.
Comments in the original code indicate that \code{calculateGibbsOfWater} is valid for 100 \le \T \le 1000 \degC and \P \ge 1000 bar.
+Likewise, the power function fit of the dielectric constant (epsilon) is valid for 100 \le \T \le 1200 \degC and \P \ge 1000 bar (Sverjensky et al., 2014).
}
\value{
@@ -37,6 +38,7 @@
}
\seealso{
+\code{\link{water.DEW}}
}
\examples{
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/man/water.Rd 2017-09-25 11:40:05 UTC (rev 216)
@@ -36,7 +36,7 @@
\item{\samp{IAPWS95}}{Thermodynamic properties are calculated using an implementation in \R code of the \acronym{IAPWS-95} formulation (Wagner and Pruss, 2002), and electrostatic properties are calculated using the equations of Archer and Wang, 1990. See \code{\link{IAPWS95}} and more information below.}
- \item{\samp{DEW}}{Thermodynamic and electrostatic properties are calculated using the Deep Earth Water (\acronym{DEW}) model (Sverjensky et al., 2014). The defaults for \code{T} and \code{P} reflect the minimum values for applicability of the model. See \code{\link{DEW}}.}
+ \item{\samp{DEW}}{Thermodynamic and electrostatic properties are calculated using the Deep Earth Water (\acronym{DEW}) model (Sverjensky et al., 2014). The defaults for \code{T} and \code{P} reflect the minimum values for applicability of the model; calculations for lower \code{T} and/or \code{P} points fall back to using \samp{SUPCRT92}. See \code{\link{DEW}}.}
}
@@ -73,7 +73,7 @@
\code{NBorn} \tab N Born function \tab bar\eqn{^{-2}}{^-2} \tab * \tab NA \tab NA \cr
\code{UBorn} \tab U Born function \tab bar\eqn{^{-1}}{^-1} K\eqn{^{-1}}{^-1} \tab * \tab NA \tab NA \cr
\code{V} \tab Volume \tab cm\eqn{^3} mol\eqn{^{-1}}{^-1} \tab * \tab * \tab * \cr
- \code{rho} \tab Density \tab kg cm\eqn{^3} \tab * \tab * \tab * \cr
+ \code{rho} \tab Density \tab kg m\eqn{^3} \tab * \tab * \tab * \cr
\code{Psat} \tab Saturation vapor pressure \tab bar \tab * \tab * \tab NA \cr
\code{E} \tab Isobaric expansivity \tab cm\eqn{^3} K\eqn{^{-1}}{^-1} \tab NA \tab * \tab NA \cr
\code{kT} \tab Isothermal compressibility \tab cm\eqn{^3} bar\eqn{^{-1}}{^-1} \tab NA \tab * \tab NA \cr
Modified: pkg/CHNOSZ/tests/testthat/test-water.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-water.R 2017-09-25 09:13:16 UTC (rev 215)
+++ pkg/CHNOSZ/tests/testthat/test-water.R 2017-09-25 11:40:05 UTC (rev 216)
@@ -16,6 +16,11 @@
expect_equal(water.SUPCRT92("alpha", T, P)[, 1] * 1e6, c(268.06, 625.55), tolerance=1e-2)
})
+test_that("water.DEW() gives expected erros and messages", {
+ expect_error(water.DEW("G", P="Psat"), "Psat isn\'t supported")
+ expect_message(water.DEW("G", T=c(298.15, 373.15), P=c(1, 1000)), "using SUPCRT calculations")
+})
+
# reference
# Fine, R. A. and Millero, F. J. (1973)
More information about the CHNOSZ-commits
mailing list