Add bryan case vapor function

src/func_table.cpp

@@ -8,20 +8,22 @@ namespace kintera {
 /////  func1 registry  /////
 
 user_func1 func1_table_cpu[] = {
-    h2o_ideal,     h2o_ideal_ddT,     nh3_ideal,       nh3_ideal_ddT,
-    nh3_h2s_lewis, nh3_h2s_lewis_ddT, h2s_ideal,       h2s_ideal_ddT,
-    h2s_antoine,   h2s_antoine_ddT,   ch4_ideal,       ch4_ideal_ddT,
-    so2_antoine,   so2_antoine_ddT,   co2_antoine,     co2_antoine_ddT,
-    kcl_lodders,   kcl_lodders_ddT,   na_h2s_visscher, na_h2s_visscher_ddT};
+    h2o_ideal,       h2o_ideal_ddT,      h2o_bryan,     h2o_bryan_ddT,
+    nh3_ideal,       nh3_ideal_ddT,      nh3_h2s_lewis, nh3_h2s_lewis_ddT,
+    h2s_ideal,       h2s_ideal_ddT,      h2s_antoine,   h2s_antoine_ddT,
+    ch4_ideal,       ch4_ideal_ddT,      so2_antoine,   so2_antoine_ddT,
+    co2_antoine,     co2_antoine_ddT,    kcl_lodders,   kcl_lodders_ddT,
+    na_h2s_visscher, na_h2s_visscher_ddT};
 
 std::vector<std::string> func1_names = {
-    "h2o_ideal",       "h2o_ideal_ddT",      "nh3_ideal",
-    "nh3_ideal_ddT",   "nh3_h2s_lewis",      "nh3_h2s_lewis_ddT",
-    "h2s_ideal",       "h2s_ideal_ddT",      "h2s_antoine",
-    "h2s_antoine_ddT", "ch4_ideal",          "ch4_ideal_ddT",
-    "so2_antoine",     "so2_antoine_ddT",    "co2_antoine",
-    "co2_antoine_ddT", "kcl_lodders",        "kcl_lodders_ddT",
-    "na_h2s_visscher", "na_h2s_visscher_ddT"};
+    "h2o_ideal",          "h2o_ideal_ddT",     "h2o_bryan",
+    "h2o_bryan_ddT",      "nh3_ideal",         "nh3_ideal_ddT",
+    "nh3_h2s_lewis",      "nh3_h2s_lewis_ddT", "h2s_ideal",
+    "h2s_ideal_ddT",      "h2s_antoine",       "h2s_antoine_ddT",
+    "ch4_ideal",          "ch4_ideal_ddT",     "so2_antoine",
+    "so2_antoine_ddT",    "co2_antoine",       "co2_antoine_ddT",
+    "kcl_lodders",        "kcl_lodders_ddT",   "na_h2s_visscher",
+    "na_h2s_visscher_ddT"};
 
 /////  func2 registry  /////
 

src/func_table.cu

@@ -10,6 +10,8 @@ __device__ user_func1 func1_table_cuda[] = {
     nullptr,
     h2o_ideal,
     h2o_ideal_ddT,
+    h2o_bryan,
+    h2o_bryan_ddT,
     nh3_ideal,
     nh3_ideal_ddT,
     nh3_h2s_lewis,

src/vapors/vapor_functions.h

@@ -46,6 +46,18 @@ inline double h2o_ideal_ddT(double T) {
          tr;
 }
 
+DISPATCH_MACRO
+inline double h2o_bryan(double T) {
+  double beta = 24.845, delta = 4.986009, tr = 273.16, pr = 611.7;
+  return logsvp_ideal(T / tr, beta, delta) + log(pr);
+}
+
+DISPATCH_MACRO
+inline double h2o_bryan_ddT(double T) {
+  double beta = 24.845, delta = 4.986009, tr = 273.16;
+  return logsvp_ideal_ddT(T / tr, beta, delta) / tr;
+}
+
 DISPATCH_MACRO
 inline double nh3_ideal(double T) {
   double betal = 20.08, gammal = 5.62, betas = 20.64, gammas = 1.43, tr = 195.4,

tests/CMakeLists.txt

@@ -10,6 +10,7 @@ enable_testing()
 setup_test(test_reaction)
 setup_test(test_species)
 setup_test(test_thermo)
+setup_test(test_vapor_functions)
 setup_test(test_kinetics)
 setup_test(test_photolysis_options)
 setup_test(test_photolysis_kinetics)

tests/test_vapor_functions.cpp

@@ -0,0 +1,63 @@
+// external
+#include <gtest/gtest.h>
+
+// C/C++
+#include <cmath>
+
+// torch
+#include <torch/torch.h>
+
+// kintera
+#include <kintera/vapors/vapor_functions.h>
+
+#include <kintera/thermo/log_svp.hpp>
+
+using namespace kintera;
+
+namespace {
+
+double h2o_bryan_expected(double T) {
+  double beta = 24.845;
+  double delta = 4.986009;
+  double tr = 273.16;
+  double pr = 611.7;
+  return (1. - tr / T) * beta - delta * std::log(T / tr) + std::log(pr);
+}
+
+double h2o_bryan_ddT_expected(double T) {
+  double beta = 24.845;
+  double delta = 4.986009;
+  double tr = 273.16;
+  double t = T / tr;
+  return (beta / (t * t) - delta / t) / tr;
+}
+
+}  // namespace
+
+TEST(VaporFunctions, h2o_bryan_matches_athena_liquid_branch) {
+  for (double temp : {250.0, 273.16, 289.85, 300.0}) {
+    EXPECT_NEAR(h2o_bryan(temp), h2o_bryan_expected(temp), 1.e-12);
+    EXPECT_NEAR(h2o_bryan_ddT(temp), h2o_bryan_ddT_expected(temp), 1.e-12);
+  }
+}
+
+TEST(VaporFunctions, h2o_bryan_keeps_liquid_branch_below_triple_point) {
+  double temp = 250.0;
+  EXPECT_NEAR(h2o_bryan(temp), h2o_bryan_expected(temp), 1.e-12);
+  EXPECT_GT(std::abs(h2o_bryan(temp) - h2o_ideal(temp)), 1.e-3);
+}
+
+TEST(VaporFunctions, h2o_bryan_dispatches_through_log_svp) {
+  auto nucleation = NucleationOptionsImpl::create();
+  nucleation->logsvp({"h2o_bryan"});
+  LogSVPFunc::init(nucleation);
+
+  auto temp = torch::tensor({250.0, 289.85}, torch::kFloat64);
+  auto logsvp = LogSVPFunc::call(temp).squeeze(-1);
+  auto grad = LogSVPFunc::grad(temp).squeeze(-1);
+
+  EXPECT_NEAR(logsvp[0].item<double>(), h2o_bryan_expected(250.0), 1.e-12);
+  EXPECT_NEAR(logsvp[1].item<double>(), h2o_bryan_expected(289.85), 1.e-12);
+  EXPECT_NEAR(grad[0].item<double>(), h2o_bryan_ddT_expected(250.0), 1.e-12);
+  EXPECT_NEAR(grad[1].item<double>(), h2o_bryan_ddT_expected(289.85), 1.e-12);
+}