[CHNOSZ-commits] r326 - in pkg/CHNOSZ: . R inst man tests/testthat
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Sep 23 05:55:31 CEST 2018
Author: jedick
Date: 2018-09-23 05:55:29 +0200 (Sun, 23 Sep 2018)
New Revision: 326
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/diagram.R
pkg/CHNOSZ/R/subcrt.R
pkg/CHNOSZ/R/water.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/subcrt.Rd
pkg/CHNOSZ/tests/testthat/test-subcrt.R
pkg/CHNOSZ/tests/testthat/test-swap.basis.R
Log:
subcrt(): add water density threshold for properties for HKF species
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/DESCRIPTION 2018-09-23 03:55:29 UTC (rev 326)
@@ -1,6 +1,6 @@
-Date: 2018-09-22
+Date: 2018-09-23
Package: CHNOSZ
-Version: 1.1.3-33
+Version: 1.1.3-34
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
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 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/R/diagram.R 2018-09-23 03:55:29 UTC (rev 326)
@@ -166,10 +166,6 @@
# some additional steps for affinity values, but not for equilibrated activities
if(eout.is.aout) {
for(i in 1:length(pv)) {
- # change any NAs in the plotvals to -Inf, so that
- # they don't get on the plot, but permit others to
- # (useful for making mineral stability diagrams beyond transition temperatures of one or more minerals)
- pv[[i]][is.na(pv[[i]])] <- -Inf
# TODO: see vignette for an explanation for how this is normalizing
# the formulas in a predominance calculation
if(normalize) pv[[i]] <- (pv[[i]] + eout$species$logact[i] / n.balance[i]) - log10(n.balance[i])
Modified: pkg/CHNOSZ/R/subcrt.R
===================================================================
--- pkg/CHNOSZ/R/subcrt.R 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/R/subcrt.R 2018-09-23 03:55:29 UTC (rev 326)
@@ -67,7 +67,8 @@
}
# warn for too high temperatures for Psat 20171110
- if(identical(P, "Psat") & any(T > 647.067)) warning("attempting calculation at P = 'Psat' for some T > Tcritical; set P = 1 (or higher)")
+ warnings <- character()
+ if(identical(P, "Psat") & any(T > 647.067)) warnings <- c(warnings, "P = 'Psat' undefined for T > Tcritical")
# gridding?
do.grid <- FALSE
@@ -248,10 +249,10 @@
return(subcrt(species=newspecies, coeff=newcoeff, state=newstate,
property=property, T=outvert(T, "K"), P=P, grid=grid, convert=convert, logact=logact, exceed.Ttr=FALSE))
} else if(identical(action.unbalanced,'warn'))
- warning(paste('reaction was unbalanced, missing', as.chemical.formula(miss)),call.=FALSE)
+ warnings <- c(warnings, paste('reaction was unbalanced, missing', as.chemical.formula(miss)))
} else {
if(identical(action.unbalanced,'warn'))
- warning(paste('reaction was unbalanced, missing', as.chemical.formula(miss)),call.=FALSE)
+ warnings <- c(warnings, paste('reaction was unbalanced, missing', as.chemical.formula(miss)))
}
}
}
@@ -281,6 +282,14 @@
hkfstuff <- hkf(eosprop, parameters = param, T = T, P = P, H2O.props=H2O.props)
p.aq <- hkfstuff$aq
H2O.PT <- hkfstuff$H2O
+ # set properties to NA for density below 0.35 g/cm3 (near-critical isochore; threshold used in SUPCRT92) 20180922
+ ilowrho <- H2O.PT$rho < 350
+ ilowrho[is.na(ilowrho)] <- FALSE
+ if(any(ilowrho)) {
+ for(i in 1:length(p.aq)) p.aq[[i]][ilowrho, ] <- NA
+ if(sum(ilowrho)==1) ctext <- "condition" else ctext <- "conditions"
+ warnings <- c(warnings, paste0("below density threshold for applicability of revised HKF equations (", sum(ilowrho), " T,P ", ctext, ")"))
+ }
# calculate activity coefficients if ionic strength is not zero
if(any(IS != 0)) {
if(grepl("Helgeson", thermo$opt$nonideal)) p.aq <- nonideal(iphases[isaq], p.aq, newIS, T, P, H2O.PT$A_DH, H2O.PT$B_DH)
@@ -556,6 +565,11 @@
# add names to the output
names(out$out) <- as.character(reaction$name)
}
+ # add warnings to output 20180922
+ if(length(warnings) > 0) {
+ out <- c(out, list(warnings=warnings))
+ for(warn in warnings) warning(warn)
+ }
return(out)
}
Modified: pkg/CHNOSZ/R/water.R
===================================================================
--- pkg/CHNOSZ/R/water.R 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/R/water.R 2018-09-23 03:55:29 UTC (rev 326)
@@ -117,9 +117,11 @@
# assemble additional properties: V, rho, Psat, E, kT
if(any(iprop > 23)) {
mwH2O <- 18.0152 # SUP92.f
- V=mwH2O/rho.out
- rho=rho.out*1000
- Psat=P.out
+ V <- mwH2O/rho.out
+ rho <- rho.out*1000
+ # rho==0 should be NA 20180923
+ rho[rho==0] <- NA
+ Psat <- P.out
E <- V*w.out$alpha
kT <- V*w.out$beta
# A and B parameters in Debye-Huckel equation:
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/inst/NEWS 2018-09-23 03:55:29 UTC (rev 326)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-33 (2018-09-22)
+CHANGES IN CHNOSZ 1.1.3-34 (2018-09-23)
---------------------------------------
THERMODYNAMIC DATA
@@ -45,37 +45,39 @@
- Change abbreviation of grossular to Grs.
-OTHER CHANGES
+DIAGRAMS
- Lines in 1-D diagram()s can optionally be drawn as splines using the
method for splinefun() given in the 'spline.method' argument (the
default of NULL means no splines).
+- Add 'srt' argument to diagram() (rotation of line labels).
+
+- Export thermo.axis(), as it is useful for adding major and minor tick
+ marks after (above) other plot elements such as legends.
+
+- diagram(): rename 'what' argument to 'type'.
+
+- digram(): add new type of diagram, 'saturation', which is used to
+ plot saturation lines for minerals (where their affinity equals
+ zero).
+
+OTHER CHANGES
+
- Add dumpdata() for returning/writing all packaged thermodynamic data
(including default database and optional data files).
-- Add 'srt' argument to diagram() (rotation of line labels).
-
- TODO: fix overly long message for info("SiO2").
- In equilibrate(), accept a length > 1 'normalize' argument to
normalize the chemical formulas of only the selected species.
-- Export thermo.axis(), as it is useful for adding major and minor tick
- marks after (above) other plot elements such as legends.
-
- Add C implementation of counting occurrences of all letters in a
string (src/count_letters.c) to speed up operation of count.aa().
- read.fasta(): add support for file connections created using
archive::archive_read (https://github.com/jimhester/archive).
-- diagram(): rename 'what' argument to 'type'.
-
-- digram(): add new type of diagram, 'saturation', which is used to
- plot saturation lines for minerals (where their affinity equals
- zero).
-
- Add demo/bison.R (average oxidation state of carbon of metagenome-
derived proteins in different microbial phyla at Bison Pool)
@@ -99,6 +101,10 @@
Thanks to Feng Lai for reporting the incorrect behavior caused by
this bug.
+- For species in the revised HKF model, subcrt() now sets properties to
+ NA where the density of H2O is less than 0.35 g/cm3, avoiding the
+ output of bogus values in this region. Thanks to Evgeniy Bastrakov.
+
CHANGES IN CHNOSZ 1.1.3 (2017-11-13)
------------------------------------
Modified: pkg/CHNOSZ/man/subcrt.Rd
===================================================================
--- pkg/CHNOSZ/man/subcrt.Rd 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/man/subcrt.Rd 2018-09-23 03:55:29 UTC (rev 326)
@@ -83,9 +83,6 @@
\code{subcrt} is modeled after the functionality of the \acronym{SUPCRT92} package (Johnson et al., 1992).
Certain features of \acronym{SUPCRT92} are not available here, for example, calculations as a function of density of \H2O instead of pressure, or calculations of temperatures of univariant curves (i.e. for which \code{logK} is zero).
-The informative messages produced by \code{SUPCRT92} when temperature or pressure limits of the equations of state are exceeded generally are not reproduced here.
-However, \code{NA}s may be produced in the output of \code{subcrt} if the requisite thermodynamic or electrostatic properties of water can not be calculated at given conditions.
-Specifically, \code{NA}s are produced for calculations at \samp{Psat} when the temperature exceeds the critical temperature of \H2O.
For calculations below 273.16 K, the pressure should be set to 1, as \Psat is not defined in these conditions.
@@ -94,9 +91,15 @@
}
\section{Warning}{
-Although \acronym{SUPCRT92} prohibits calculations above 350 \degC at \Psat (\dQuote{beyond range of applicability of aqueous species equations}), there is no corresponding limit in place in \code{subcrt} (or \code{\link{hkf}}).
-Therefore, CHNOSZ can perform calculations up to the critical temperature (373.917 \degC) at \Psat, but these settings represent untested extrapolations.
-Unexpected results are evident in the discontinuity in the value of \logK at \Psat shown in \code{\link{demos}("NaCl")}.
+Although \acronym{SUPCRT92} prohibits calculations above 350 \degC at \Psat (\dQuote{beyond range of applicability of aqueous species equations}), CHNOSZ does not impose this limitation, and allows calculations up to the critical temperature (373.917 \degC) at \Psat.
+Interpret calculations between 350 \degC and the critical temperature at \Psat at your own risk.
+The discontinuity in the value of \logK at \Psat that is apparent in \code{\link{demos}("NaCl")} demonstrates one unexpected result.
+
+\code{NA}s are produced for calculations at \samp{Psat} when the temperature exceeds the critical temperature of \H2O.
+In addition, properties of species using the revised HKF equations are set to \code{NA} wherever the density of \H2O < 0.35 g/cm\S{3} (threshold just above the critical isochore; Johnson et al., 1992).
+Both of these situations produce warnings, which are stored in the \samp{warnings} element of the return value.
+
+\code{NA}s are also output if the T, P conditions are otherwise beyond the capabilities of the water equations of state derived from SUPCRT92 (H2O92D.f), but the messages about this are produced by \code{\link{water.SUPCRT92}} rather than \code{subcrt}.
}
\value{
Modified: pkg/CHNOSZ/tests/testthat/test-subcrt.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/tests/testthat/test-subcrt.R 2018-09-23 03:55:29 UTC (rev 326)
@@ -184,6 +184,13 @@
expect_equal(unique(s4$out$logK), 0)
})
+test_that("properties of HKF species below 0.35 g/cm3 are NA and give a warning", {
+ wtext <- "below density threshold for applicability of revised HKF equations \\(2 T,P conditions\\)"
+ expect_warning(s1 <- subcrt(c("Na+", "quartz"), T=450, P=c(400, 450, 500)), wtext)
+ expect_equal(sum(is.na(s1$out$`Na+`$logK)), 2)
+ expect_equal(sum(is.na(s1$out$quartz$logK)), 0)
+})
+
# references
# Amend, J. P. and Shock, E. L. (2001)
Modified: pkg/CHNOSZ/tests/testthat/test-swap.basis.R
===================================================================
--- pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2018-09-22 05:38:00 UTC (rev 325)
+++ pkg/CHNOSZ/tests/testthat/test-swap.basis.R 2018-09-23 03:55:29 UTC (rev 326)
@@ -26,7 +26,7 @@
expect_error(basis.logact(ep), "element\\(s\\) O not found in basis")
})
-test_that("equil.potentials - basis.logact - element.mu makes a roundtrip at 25 and 100 degrees C", {
+test_that("equil.potentials - basis.logact - element.mu makes a roundtrip at 25 and 99 degrees C", {
basis(c("graphite", "H2", "O2"), c("cr", "gas", "gas"))
ispecies <- info(c("ethane", "propane", "acetic acid", "propanoic acid"))
# at 25 degrees C
@@ -37,10 +37,10 @@
basis(names(bl25), bl25)
# element.mu() calculates the chemical potentials of the elements from the current setting of basis species
expect_equal(element.mu(), ep25)
- # at 100 degrees C
- w100 <- run.wjd(ispecies, as.chemical.formula(colMeans(i2A(ispecies))), T=100)
- ep100 <- equil.potentials(w100)
- bl100 <- basis.logact(ep100, T=100)
- basis(names(bl100), bl100)
- expect_equal(element.mu(T=100), ep100)
+ # at 99 degrees C
+ w99 <- run.wjd(ispecies, as.chemical.formula(colMeans(i2A(ispecies))), T=99)
+ ep99 <- equil.potentials(w99)
+ bl99 <- basis.logact(ep99, T=99)
+ basis(names(bl99), bl99)
+ expect_equal(element.mu(T=99), ep99)
})
More information about the CHNOSZ-commits
mailing list