[CHNOSZ-commits] r227 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Sep 29 11:47:24 CEST 2017
Author: jedick
Date: 2017-09-29 11:47:24 +0200 (Fri, 29 Sep 2017)
New Revision: 227
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/cgl.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/util.misc.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
implement SUPCRT92 handling of quartz and coesite volumes
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-29 09:47:24 UTC (rev 227)
@@ -1,6 +1,6 @@
Date: 2017-09-29
Package: CHNOSZ
-Version: 1.1.0-25
+Version: 1.1.0-26
Title: Thermodynamic Calculations for Geobiochemistry
Author: Jeffrey Dick
Maintainer: Jeffrey Dick <j3ffdick at gmail.com>
Modified: pkg/CHNOSZ/R/cgl.R
===================================================================
--- pkg/CHNOSZ/R/cgl.R 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/R/cgl.R 2017-09-29 09:47:24 UTC (rev 227)
@@ -12,12 +12,21 @@
prop <- eargs$prop
EOS.Prop <- eargs$Prop
+ # the number of T, P conditions
+ ncond <- max(c(length(T), length(P)))
+
# initialize output list
out <- list()
+ # loop over each species
for(k in 1:nrow(parameters)) {
- # loop over each species
+ # the parameters for *this* species
PAR <- parameters[k, ]
- w <- NULL
+ # start with NA values
+ values <- data.frame(matrix(NA, ncol = length(prop), nrow=ncond))
+ colnames(values) <- EOS.Prop
+ # additional calculations for quartz and coesite
+ qtz <- quartz_coesite(PAR, T, P)
+ isqtz <- !identical(qtz$V, 0)
for(i in 1:length(prop)) {
PROP <- prop[i]
# a test for availability of the EoS parameters
@@ -33,30 +42,33 @@
else p <- PAR$a + PAR$b * T + PAR$c * T^-2 + PAR$d * T^-0.5 + PAR$e * T^2 + PAR$f * T^PAR$lambda
}
if(PROP == "v") {
- p <- rep(PAR$V, length(T))
+ if(isqtz) p <- qtz$V
+ else p <- rep(PAR$V, ncond)
}
if(PROP %in% c("e", "kt")) {
- p <- rep(NA, length(T))
+ p <- rep(NA, ncond)
warning("cgl: E and/or kT of cr, gas and/or liq species are NA.")
}
if(PROP == "g") {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$G + PAR$Cp * (T - Tr - T * log(T/Tr)) else {
+ if(!hasEOS) p <- PAR$Cp * (T - Tr - T * log(T/Tr)) else {
# Gibbs energy integral: the value at Tref plus heat capacity terms
- p <- PAR$G + PAR$a * (T - Tr - T * log(T/Tr)) -
+ p <- PAR$a * (T - Tr - T * log(T/Tr)) -
PAR$b * (T - Tr)^2 / 2 - PAR$c * (1/T + T/Tr^2 - 2/Tr) / 2 -
PAR$d * (T^0.5 - 0.5 * T * Tr^-0.5 - 0.5 * Tr^0.5) / -0.25 -
PAR$e * (T^3 - 3 * T * Tr^2 + 2 * Tr^3) / 6
}
- # entropy and volume terms
- if(!is.na(PAR$S)) p <- p - PAR$S * (T - Tr)
- if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
# use additional heat capacity term if it's defined
if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
if(PAR$lambda == -1) p <- p + PAR$f * (log(T/Tr) - T * (1/Tr - 1/T))
else p <- p + PAR$f * ( T^(PAR$lambda + 1) - (PAR$lambda + 1) * T * Tr^PAR$lambda +
PAR$lambda * Tr^(PAR$lambda + 1) ) / ( PAR$lambda * (PAR$lambda + 1) )
}
+ # entropy and volume terms
+ if(!is.na(PAR$S)) p <- p - PAR$S * (T - Tr)
+ if(isqtz) p <- p + qtz$G
+ else if(!is.na(PAR$V)) p <- p + convert(PAR$V * (P - Pr), "calories")
+ p <- PAR$G + p
}
if(PROP == "h") {
# use constant Cp if the EoS parameters are not available
@@ -64,32 +76,83 @@
p <- PAR$a * (T - Tr) + PAR$b * (T^2 - Tr^2) / 2 +
PAR$c * (1/T - 1/Tr) / -1 + PAR$d * (T^0.5 - Tr^0.5) / 0.5 +
PAR$e * (T^3 - Tr^3) / 3
- # SUPCRT seems to ignore this term? ... 20070802
- # + convert(PAR$V*(P-Pr),'calories')
}
- p <- PAR$H + p
if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
if(PAR$lambda == -1) p <- p + PAR$f * log(T/Tr)
else p <- p - PAR$f * ( T^(PAR$lambda + 1) - Tr^(PAR$lambda + 1) ) / (PAR$lambda + 1)
}
+ if(isqtz) p <- p + qtz$H
+ ## SUPCRT seems to ignore this term? ... 20070802
+ #else p <- p + convert(PAR$V*(P-Pr),'calories')
+ p <- PAR$H + p
}
if(PROP=="s") {
# use constant Cp if the EoS parameters are not available
- if(!hasEOS) p <- PAR$Cp * log(T/Tr) else
- p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
- PAR$c * (T^-2 - Tr^-2) / -2 + PAR$e * (T^2 - Tr^2) / 2 +
- PAR$d * (T^-0.5 - Tr^-0.5) / -0.5
- p <- PAR$S + p
+ if(!hasEOS) p <- PAR$Cp * log(T/Tr) else {
+ p <- PAR$a * log(T / Tr) + PAR$b * (T - Tr) +
+ PAR$c * (T^-2 - Tr^-2) / -2 + PAR$e * (T^2 - Tr^2) / 2 +
+ PAR$d * (T^-0.5 - Tr^-0.5) / -0.5
+ }
if(!is.na(PAR$f) & !is.na(PAR$lambda)) if(PAR$f != 0) {
p <- p + PAR$f * (T^PAR$lambda - Tr^PAR$lambda) / PAR$lambda
}
+ p <- PAR$S + p + qtz$S
}
- wnew <- data.frame(p)
- if(i > 1) w <- cbind(w, wnew) else w <- wnew
+ values[, i] <- p
}
- colnames(w) <- EOS.Prop
- out[[k]] <- w
+ out[[k]] <- values
}
return(out)
}
+### unexported function ###
+
+# calculate GHS and V corrections for quartz and coesite 20170929
+# (these are the only mineral phases which SUPCRT92 applies a variable volume)
+quartz_coesite <- function(PAR, T, P) {
+ # the corrections are 0 for anything other than quartz and coesite
+ if(!PAR$name %in% c("quartz", "coesite")) return(list(G=0, H=0, S=0, V=0))
+ ncond <- max(c(length(T), length(P)))
+ # get Tr, Pr and TtPr (transition temperature at Pr)
+ Pr <- get("thermo")$opt$Pr # 1 bar
+ Tr <- get("thermo")$opt$Tr # 298.15 K
+ TtPr <- 848 # Kelvin
+ # constants from SUP92D.f
+ aa <- 549.824
+ ba <- 0.65995
+ ca <- -0.4973e-4
+ VPtTta <- 23.348
+ VPrTtb <- 23.72
+ Stran <- 0.342
+ # constants from REAC92D.f
+ VPrTra <- 22.688 # VPrTr(a-quartz)
+ Vdiff <- 2.047 # VPrTr(a-quartz) - VPrTr(coesite)
+ #k <- 38.5 # dPdTtr(a/b-quartz)
+ k <- 38.45834 # calculated in CHNOSZ: dPdTtr(info("quartz"))
+ # code adapted from REAC92D.f
+ qphase <- gsub("cr", "", PAR$state)
+ if(qphase == 2) {
+ Pstar <- P
+ Sstar <- rep(0, ncond)
+ V <- rep(VPrTtb, ncond)
+ } else {
+ Pstar <- Pr + k * (T - TtPr)
+ Sstar <- rep(Stran, ncond)
+ V <- VPrTra + ca*(P-Pr) + (VPtTta - VPrTra - ca*(P-Pr))*(T-Tr) / (TtPr + (P-Pr)/k - Tr)
+ }
+ Pstar[T < TtPr] <- Pr
+ Sstar[T < TtPr] <- 0
+ if(PAR$name == "coesite") {
+ VPrTra <- VPrTra - Vdiff
+ VPrTtb <- VPrTtb - Vdiff
+ V <- V - Vdiff
+ }
+ GVterm <- convert(1, "calories") * (VPrTra * (P - Pstar) + VPrTtb * (Pstar - Pr) -
+ 0.5 * ca * (2 * Pr * (P - Pstar) - (P^2 - Pstar^2)) -
+ ca * k * (T - Tr) * (P - Pstar) +
+ k * (ba + aa * ca * k) * (T - Tr) * log((aa + P/k) / (aa + Pstar/k)))
+ SVterm <- convert(1, "calories") * (-k * (ba + aa * ca * k) *
+ log((aa + P/k) / (aa + Pstar/k)) + ca * k * (P - Pstar)) - Sstar
+ list(G=GVterm, S=SVterm, H=GVterm + T*SVterm, V=V)
+}
+
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/R/subcrt.R 2017-09-29 09:47:24 UTC (rev 227)
@@ -10,7 +10,7 @@
#source("util.data.R")
subcrt <- function(species, coeff=1, state=NULL, property=c('logK','G','H','S','V','Cp'),
- T=seq(273.15,623.15,25), P='Psat', grid=NULL, convert=TRUE, check.Ttr=TRUE, exceed.Ttr=FALSE,
+ T=seq(273.15,623.15,25), P='Psat', grid=NULL, convert=TRUE, exceed.Ttr=FALSE,
logact=NULL, action.unbalanced='warn', IS=0) {
# revise the call if the states have
@@ -311,8 +311,6 @@
for(i in 1:length(iscgl)) {
# not if we're not cgl
if(!iscgl[i]) next
- # not if check.Ttr is FALSE (e.g. subcrt is called by dPdTtr)
- if(!check.Ttr) next
# name and state
myname <- reaction$name[i]
mystate <- reaction$state[i]
Modified: pkg/CHNOSZ/R/util.misc.R
===================================================================
--- pkg/CHNOSZ/R/util.misc.R 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/R/util.misc.R 2017-09-29 09:47:24 UTC (rev 227)
@@ -6,14 +6,17 @@
# calculate dP/dT for a phase transition
# (argument is index of the lower-T phase)
thermo <- get("thermo")
- t1 <- subcrt(ispecies, P=0, T=thermo$obigt$z.T[ispecies], convert=FALSE, check.Ttr=FALSE)
- t2 <- subcrt(ispecies+1, P=0, T=thermo$obigt$z.T[ispecies], convert=FALSE, check.Ttr=FALSE)
+ pars <- info(c(ispecies, ispecies+1), check.it=FALSE)
+ # the special handling for quartz and coesite interferece with this function,
+ # so we convert to uppercase names to prevent cgl() from calling quartz_coesite()
+ pars$name <- toupper(pars$name)
+ props <- cgl(c("G", "S", "V"), pars, P=0, T=thermo$obigt$z.T[ispecies])
# if these aren't the same mineral all we can say is zero
# actually, should be infinity ... the volume change is zero
- if(as.character(t1$species$name) != as.character(t2$species$name)) return(Inf)
+ if(as.character(pars$name[1]) != as.character(pars$name[2])) return(Inf)
# we really hope the G's are the same ...
- #if(abs(t2$out[[1]]$G - t2$out[[1]]$G) > 0.1) warning('dP.dT: inconsistent values of G for different phases of ',ispecies,call.=FALSE)
- dP.dT <- convert( ( t2$out[[1]]$S - t1$out[[1]]$S ) / ( t2$out[[1]]$V - t1$out[[1]]$V ), 'cm3bar' )
+ #if(abs(props[[2]]$G - props[[1]]$G) > 0.1) warning('dP.dT: inconsistent values of G for different phases of ',ispecies,call.=FALSE)
+ dP.dT <- convert( ( props[[2]]$S - props[[1]]$S ) / ( props[[2]]$V - props[[1]]$V ), 'cm3bar' )
return(dP.dT)
}
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/inst/NEWS 2017-09-29 09:47:24 UTC (rev 227)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.0-24 (2017-09-29)
+CHANGES IN CHNOSZ 1.1.0-26 (2017-09-29)
---------------------------------------
MAJOR CHANGES:
@@ -13,6 +13,10 @@
- Add demo DEW.R.
+- Implement SUPCRT92's handling of variable volume for quartz and
+ coesite. Calculations for other minerals still assume constant
+ volume of each phase.
+
- water.lines() now works for diagrams of Eh, pe, logfO2, logaO2,
logfH2, or logaH2 vs pH, T, or P. It is possible to have T or P on
either the x- or y-axis.
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/man/subcrt.Rd 2017-09-29 09:47:24 UTC (rev 227)
@@ -10,8 +10,8 @@
subcrt(species, coeff = 1, state = NULL,
property = c("logK","G","H","S","V","Cp"),
T = seq(273.15,623.15,25), P = "Psat", grid = NULL,
- convert = TRUE, check.Ttr = TRUE, exceed.Ttr = FALSE,
- logact = NULL, action.unbalanced = "warn", IS = 0)
+ convert = TRUE, exceed.Ttr = FALSE, logact = NULL,
+ action.unbalanced = "warn", IS = 0)
}
\arguments{
@@ -22,7 +22,6 @@
\item{T}{numeric, temperature(s) of the calculation}
\item{P}{numeric, pressure(s) of the calculation, or character, \samp{Psat}}
\item{grid}{character, type of \code{P}\eqn{\times}{x}\code{T} grid to produce (NULL, the default, means no gridding)}
- \item{check.Ttr}{logical, check if temperatures are beyond transition temperatures?}
\item{exceed.Ttr}{logical, calculate Gibbs energies of mineral phases and other species beyond their transition temperatures?}
\item{logact}{numeric, logarithms of activities of species in reaction}
\item{convert}{logical, are input and output units of T and P those of the user (\code{TRUE}) (see \code{\link{T.units}}), or are they Kelvin and bar (\code{FALSE})?}
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-29 06:26:36 UTC (rev 226)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-29 09:47:24 UTC (rev 227)
@@ -111,7 +111,7 @@
expect_equal(CHNOSZ_1bar$G, S92_1bar$G, tolerance = 1e-5)
expect_equal(CHNOSZ_1bar$H, S92_1bar$H, tolerance = 1e-5)
expect_equal(CHNOSZ_1bar$S, S92_1bar$S, tolerance = 1e-2)
-# expect_equal(CHNOSZ_1bar$V, S92_1bar$V, tolerance = 1e-2)
+ expect_equal(CHNOSZ_1bar$V, S92_1bar$V, tolerance = 1e-2)
# output from SUPCRT92 for reaction specified as "1 QUARTZ" run at 500 bar
# (SUPCRT shows phase transition at 587.811 deg C)
@@ -125,6 +125,10 @@
590 -214653 -208668 25.4 23.7
")
CHNOSZ_500bar <- subcrt("quartz", T=seq(585, 590), P=500)$out[[1]]
+ expect_equal(CHNOSZ_500bar$G, S92_500bar$G, tolerance = 1e-5)
+ expect_equal(CHNOSZ_500bar$H, S92_500bar$H, tolerance = 1e-4)
+ expect_equal(CHNOSZ_500bar$S, S92_500bar$S, tolerance = 1e-3)
+ expect_equal(CHNOSZ_500bar$V, S92_500bar$V, tolerance = 1e-2)
})
# references
More information about the CHNOSZ-commits
mailing list