From 66d5fc0c1d0f6b305bea3caae0dc04ba42d4b1b2 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Thu, 21 Oct 2021 10:30:36 -0400
Subject: [PATCH] Restructure storing of points so that the first point also
gets saved
---
include/teqp/algorithms/critical_tracing.hpp | 50 ++++++++++++--------
1 file changed, 29 insertions(+), 21 deletions(-)
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index 99b9f35..42648a8 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -402,9 +402,33 @@ struct CriticalTracing {
// Make T and rhovec references to the contents of x0 vector
// The views are mutable (danger!)
- double &T = x0[0];
+ double& T = x0[0];
auto rhovec = Eigen::Map<Eigen::ArrayXd>(&(x0[0]) + 1, x0.size() - 1);
+ auto store_point = [&]() {
+
+ // Calculate some other parameters, for debugging, or scientific interest
+ auto rhotot = rhovec.sum();
+ using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
+ double p = rhotot * model.R(rhovec / rhovec.sum()) * T + id::get_pr(model, T, rhovec);
+ auto conditions = get_criticality_conditions(model, T, rhovec);
+ double splus = id::get_splus(model, T, rhovec);
+
+ // Store the data in a JSON structure
+ nlohmann::json point = {
+ {"t", t},
+ {"T / K", T},
+ {"rho0 / mol/m^3", static_cast<double>(rhovec[0])},
+ {"rho1 / mol/m^3", static_cast<double>(rhovec[1])},
+ {"c", c},
+ {"s^+", splus},
+ {"p / Pa", p},
+ {"lambda1", conditions[0]},
+ {"dirderiv(lambda1)/dalpha", conditions[1]},
+ };
+ JSONdata.push_back(point);
+ };
+
{
auto dxdt = x0;
xprime(x0, dxdt, -1.0);
@@ -440,6 +464,9 @@ struct CriticalTracing {
write_line();
}
}
+ if (iter == 0 && retry_count == 0) {
+ store_point();
+ }
double dtold = dt;
auto x0_previous = x0;
@@ -523,26 +550,7 @@ struct CriticalTracing {
if (!filename.empty()) {
write_line();
}
-
- // Calculate some other parameters, for debugging, or scientific interest
- using id = IsochoricDerivatives<decltype(model), Scalar, VecType>;
- double p = rhotot * model.R(rhovec / rhovec.sum()) * T + id::get_pr(model, T, rhovec);
- conditions = get_criticality_conditions(model, T, rhovec);
- double splus = id::get_splus(model, T, rhovec);
-
- // Store the data in a JSON structure
- nlohmann::json point = {
- {"t", t},
- {"T / K", T},
- {"rho0 / mol/m^3", static_cast<double>(rhovec[0])},
- {"rho1 / mol/m^3", static_cast<double>(rhovec[1])},
- {"c", c},
- {"s^+", splus},
- {"p / Pa", p},
- {"lambda1", conditions[0]},
- {"dirderiv(lambda1)/dalpha", conditions[1]},
- };
- JSONdata.push_back(point);
+ store_point();
if (counter_T_converged > options.small_T_count) {
break;
--
GitLab