[CHNOSZ-commits] r223 - in pkg/CHNOSZ: . R man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Wed Sep 27 17:35:15 CEST 2017
Author: jedick
Date: 2017-09-27 17:35:15 +0200 (Wed, 27 Sep 2017)
New Revision: 223
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/DEW.R
pkg/CHNOSZ/R/hkf.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/man/util.water.Rd
pkg/CHNOSZ/man/water.Rd
pkg/CHNOSZ/tests/testthat/test-DEW.R
pkg/CHNOSZ/tests/testthat/test-subcrt.R
Log:
enable g-function in volume calculations with DEW
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/DESCRIPTION 2017-09-27 15:35:15 UTC (rev 223)
@@ -1,6 +1,6 @@
Date: 2017-09-27
Package: CHNOSZ
-Version: 1.1.0-21
+Version: 1.1.0-22
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-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/R/DEW.R 2017-09-27 15:35:15 UTC (rev 223)
@@ -152,7 +152,7 @@
ZD05_R * TK * density * delta / m
}
-# 'Calculates the partial derivative of density with respect to pressure, i.e. (d(rho)/dP)_T
+# 'Calculates the partial derivative of density with respect to pressure, i.e. (d(rho)/dP)_T, in units of g^3/cm^3/bar
# 'density - The density of water, in g/cm^3
# 'temperature - The temperature of water, in Celsius
calculate_drhodP <- function(density, temperature) {
Modified: pkg/CHNOSZ/R/hkf.R
===================================================================
--- pkg/CHNOSZ/R/hkf.R 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/R/hkf.R 2017-09-27 15:35:15 UTC (rev 223)
@@ -34,15 +34,17 @@
# Born functions and epsilon - for HKF calculations
H2O.props <- c(H2O.props, "QBorn", "XBorn", "YBorn", "epsilon")
if(grepl("SUPCRT", thermo$opt$water)) {
- # using H2O92D.f from SUPCRT92
- # (alpha, daldT, beta - for partial derivatives of omega (g function))
+ # using H2O92D.f from SUPCRT92: alpha, daldT, beta - for partial derivatives of omega (g function)
H2O.props <- c(H2O.props, "alpha", "daldT", "beta")
}
if(grepl("IAPWS", thermo$opt$water)) {
- # using IAPWS-95
- # (NBorn, UBorn - for compressibility, expansibility)
- H2O.props <- c(H2O.props, 'NBorn', 'UBorn')
+ # using IAPWS-95: NBorn, UBorn - for compressibility, expansibility
+ H2O.props <- c(H2O.props, "NBorn", "UBorn")
}
+ if(grepl("DEW", thermo$opt$water)) {
+ # using DEW model: get beta to calculate dgdP
+ H2O.props <- c(H2O.props, "beta")
+ }
H2O <- water(H2O.props, T = c(thermo$opt$Tr, T), P = c(thermo$opt$Pr, P))
H2O.PrTr <- H2O[1, ]
H2O.PT <- H2O[-1, ]
@@ -248,23 +250,23 @@
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. 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
+ if(!all(is.na(alpha)) & !all(is.na(daldT))) {
# 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
+ #dDdP <- D * beta
dDdTT <- -D * (daldT - alpha^2)
Db <- (1-D)^bg
dDbdT <- -bg*(1-D)^(bg-1)*dDdT + log(1-D)*Db*dbdT
@@ -275,11 +277,13 @@
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
+ }
+ if(!all(is.na(beta))) {
+ dgdP <- -bg*rhohat*beta*g*(1-rhohat)^-1 # Eqn. 66
+ dgdP[ifg] <- dgdP[ifg] - dfdP[ifg]
out$dgdP[idoit] <- dgdP
}
return(out)
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/R/water.R 2017-09-27 15:35:15 UTC (rev 223)
@@ -9,11 +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
+ oldwat <- thermo$opt$water
thermo$opt$water <- property
assign("thermo", thermo, "CHNOSZ")
message(paste("water: setting thermo$opt$water to", property))
- return(invisible(owater))
+ return(invisible(oldwat))
}
# make T and P equal length
if(!identical(P, "Psat")) {
@@ -53,7 +53,7 @@
"YBorn", "QBorn", "XBorn", "NBorn", "UBorn",
"V", "rho", "Psat", "de.dT", "de.dP", "P")
if(grepl("DEW", formulation))
- props <- c("G", "epsilon", "QBorn", "V", "rho")
+ props <- c("G", "epsilon", "QBorn", "V", "rho", "beta")
return(props)
}
@@ -337,13 +337,15 @@
out <- as.data.frame(out)
colnames(out) <- property
# calculate rho if it's needed for any other properties
- if(any(c("rho", "V", "QBorn", "epsilon") %in% property)) rho <- calculateDensity(pressure, temperature)
+ if(any(c("rho", "V", "QBorn", "epsilon", "beta") %in% property)) rho <- calculateDensity(pressure, temperature)
# fill in columns with values
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)
+ # divide drhodP by rho to get units of bar^-1
+ if("beta" %in% property) out$beta <- calculate_drhodP(rho, temperature) / rho
# use SUPCRT-calculated values below 100 degC and/or below 1000 bar
ilow <- T < 373.15 | P < 1000
if(any(ilow)) {
Modified: pkg/CHNOSZ/man/util.water.Rd
===================================================================
--- pkg/CHNOSZ/man/util.water.Rd 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/man/util.water.Rd 2017-09-27 15:35:15 UTC (rev 223)
@@ -63,13 +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
-ow <- water("IAPWS95")
+oldwat <- 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)
+water(oldwat)
}
\references{
Modified: pkg/CHNOSZ/man/water.Rd
===================================================================
--- pkg/CHNOSZ/man/water.Rd 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/man/water.Rd 2017-09-27 15:35:15 UTC (rev 223)
@@ -132,10 +132,10 @@
## comparing the formulations
T <- convert(c(25, 100, 200, 300), "K")
# IAPWS-95
-ow <- water("IAPWS95")
+oldwat <- water("IAPWS95")
water(water.props(), T=T)
# SUPCRT92 (the default)
-water(ow)
+water(oldwat)
water(water.props(), T=T)
## calculating Q Born function
Modified: pkg/CHNOSZ/tests/testthat/test-DEW.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-DEW.R 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/tests/testthat/test-DEW.R 2017-09-27 15:35:15 UTC (rev 223)
@@ -50,7 +50,7 @@
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")
+ oldwat <- 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
@@ -58,13 +58,13 @@
# 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)
+ water(oldwat)
})
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")
+ oldwat <- water("DEW")
add.obigt("DEW_aq")
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)
@@ -72,14 +72,14 @@
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)
expect_equal(RG_Cl, DEWG_Cl, tolerance = 1e-7)
- water(ow)
+ water(oldwat)
})
test_that("Delta G, logK, and Delta V of reactions are calculated correctly", {
- # these are reactions corresponding to Fig. 1b of Sverjensky et al., 2014 (Nat. Geosci.)
- # the properties are calculated using parameters from the DEW spreadsheet,
- # which are not necessarily identical those that were used for the paper
- ow <- water("DEW")
+ # These are reactions corresponding to Fig. 1b of Sverjensky et al., 2014 (Nat. Geosci.).
+ # The properties are calculated using parameters from the DEW spreadsheet,
+ # which are not necessarily identical those that were used for the paper.
+ oldwat <- water("DEW")
add.obigt("DEW_aq", c("CO2", "HCO3-", "CO3-2", "acetic acid", "acetate", "methane"))
T <- 600
P <- c(5000, 50000)
@@ -109,12 +109,18 @@
35.45364304557209, 33.632014510923746)
expect_equal(c(R1$logK, R2$logK, R3$logK, R4$logK, R5$logK), DEW_logK, tolerance=1e-3)
- # TODEW: add domega/dP calculation for DEW in CHNOSZ
- ## Delta V values calculated using the DEW spreadsheet
- #DEW_DV <- c(-45.26925983499276, -14.640599169742725,
- # -47.95180846799733, -9.469432706749927,
- # -27.042087540311922, -6.836267057722694,
- # -18.1550937649195, 5.513800665649967,
- # -37.37077435045512, -45.08570662275889)
- #expect_equal(c(R1$V, R2$V, R3$V, R4$V, R5$V), DEW_DV, tolerance=1e-3)
+ # Delta V values calculated using the DEW spreadsheet
+ DEW_DV <- c(-45.26925983499276, -14.640599169742725,
+ -47.95180846799733, -9.469432706749927,
+ -27.042087540311922, -6.836267057722694,
+ -18.1550937649195, 5.513800665649967,
+ -37.37077435045512, -45.08570662275889)
+ # for the aqueous species we're getting very close results
+ # (at P=5000 bar this depends on calculating drhodP -> beta -> dgdP -> dwdP -> V correctly, which is not tested above)
+ expect_equal(c(R1$V, R2$V, R3$V), DEW_DV[1:6], tolerance=1e-15)
+ # TODO: why does DEW spreadsheet use V (O2,g) == 24.465?
+ #expect_equal(c(R4$V, R5$V), DEW_DV[7:10])
+
+ # reset computational option for water
+ water(oldwat)
})
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-27 13:22:42 UTC (rev 222)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2017-09-27 15:35:15 UTC (rev 223)
@@ -87,12 +87,12 @@
})
test_that("calculations using IAPWS-95 are possible", {
- ow <- water("IAPWS95")
+ oldwat <- 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
- water(ow)
+ water(oldwat)
})
# references
More information about the CHNOSZ-commits
mailing list