From 7f0658d85eb457f0b79308a4c4fbc8003b4d9390 Mon Sep 17 00:00:00 2001
From: Ian Bell <ian.bell@nist.gov>
Date: Fri, 9 Apr 2021 12:31:37 -0400
Subject: [PATCH] Fix one genre of problem with critical tracing pure fluid
init; still not quite right...
---
include/teqp/algorithms/critical_tracing.hpp | 6 +++---
notebooks/time_teqp.py | 19 ++++++++++---------
notebooks/trace_critical.py | 3 +++
3 files changed, 16 insertions(+), 12 deletions(-)
diff --git a/include/teqp/algorithms/critical_tracing.hpp b/include/teqp/algorithms/critical_tracing.hpp
index f2777b8..8fd199e 100644
--- a/include/teqp/algorithms/critical_tracing.hpp
+++ b/include/teqp/algorithms/critical_tracing.hpp
@@ -185,7 +185,7 @@ bool any(const Iterable& foo) {
template<typename Model, typename TType, typename RhoType>
auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhovec) {
- // The derivatives of total Psi w.r.t.sigma_1(mcx for residual, analytic for ideal)
+ // The derivatives of total Psi w.r.t.sigma_1 (numerical for residual, analytic for ideal)
// Returns a tuple, with residual, ideal, total dicts with of number of derivatives, value of derivative
auto all_derivs = get_derivs(model, T, rhovec);
auto derivs = all_derivs.tot;
@@ -213,7 +213,7 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
deriv_sigma2 = (plus_sigma2.tot - minus_sigma2.tot) / (2.0 * sigma2);
stepping_desc = "conventional centered";
}
- else if (any(eval(rhovec_minus < 0))) {
+ else if (all(eval(rhovec_plus > 0))) {
// Forward derivative in the direction of v1
auto plus_sigma2 = get_derivs(model, T, rhovec_plus);
auto rhovec_2plus = (rhovec + 2 * ei.v1 * sigma2).eval();
@@ -221,7 +221,7 @@ auto get_drhovec_dT_crit(const Model& model, const TType T, const RhoType& rhove
deriv_sigma2 = (-3 * derivs + 4 * plus_sigma2.tot - plus2_sigma2.tot) / (2.0 * sigma2);
stepping_desc = "forward";
}
- else if (any(eval(rhovec_minus > 0))) {
+ else if (all(eval(rhovec_minus > 0))) {
// Negative derivative in the direction of v1
auto minus_sigma2 = get_derivs(model, T, rhovec_minus);
auto rhovec_2minus = (rhovec - 2 * ei.v1 * sigma2).eval();
diff --git a/notebooks/time_teqp.py b/notebooks/time_teqp.py
index 006ddba..ca0c6dc 100644
--- a/notebooks/time_teqp.py
+++ b/notebooks/time_teqp.py
@@ -1,6 +1,6 @@
import timeit
import sys
-sys.path.append('bld/Release')
+sys.path.append('../bld/Release')
import teqp
import numpy as np
@@ -11,7 +11,8 @@ import scipy.optimize
def build_models():
return [
teqp.PCSAFTEOS(['Methane']),
- teqp.vdWEOS([150.687], [4863000.0])
+ teqp.vdWEOS([150.687], [4863000.0]),
+ teqp.build_multifluid_model(["Methane"], '../mycp', '../mycp/dev/mixtures/mixture_binary_pairs.json')
]
def time(*, model, n, Nrep):
@@ -19,10 +20,10 @@ def time(*, model, n, Nrep):
rho = 3.0
T = 300
- f = getattr(teqp, f"get_Ar0{n}n") if n > 2 else getattr(teqp, f"get_Ar0{n}")
+ f = getattr(model, f"get_Ar0{n}n")# if n > 2 else getattr(teqp, f"get_Ar0{n}")
tic = timeit.default_timer()
for i in range(Nrep):
- f(model, T, rho, molefrac)
+ f(T, rho, molefrac)
toc = timeit.default_timer()
elap = (toc-tic)/Nrep
return elap
@@ -32,16 +33,16 @@ def timeall(*, models, Nrep):
for model in models:
for n in [1,2,3,4,5,6]:
t = time(model=model, n=n, Nrep=Nrep)
- o.append({'model': str(model), 'n': n, 't / s': t})
- df = pandas.DataFrame(o)
+ o.append({'model': str(model), 'n': n, 't / s': t, 't / us': t*1e6})
+ df = pandas.DataFrame(o)
for model,gp in df.groupby('model'):
- plt.plot(gp['n'], gp['t / s']*1e6, label=model)
+ modname = str(model).split(' object')[0].rsplit('teqp.',1)[-1]
+ plt.plot(gp['n'], gp['t / us'], label=modname)
plt.gca().set(xlabel='n', ylabel=r't / $\mu$s')
plt.legend(loc='best')
- # plt.xscale('log')
plt.yscale('log')
plt.title(r'Timing of $A^{\rm r}_{0n}$')
plt.show()
if __name__ == '__main__':
- timeall(models=build_models(), Nrep= 10000)
\ No newline at end of file
+ timeall(models=build_models(), Nrep= 1000)
\ No newline at end of file
diff --git a/notebooks/trace_critical.py b/notebooks/trace_critical.py
index 4dbfaf4..564d7bb 100644
--- a/notebooks/trace_critical.py
+++ b/notebooks/trace_critical.py
@@ -26,7 +26,10 @@ def trace():
T0 = Tcvec[i]
z = np.array([0, 0]); z[i] = 1
rho0 = z*rhocvec[i]
+ tic = timeit.default_timer()
curveJSON = teqp.trace_critical_arclength_binary(model, T0, rho0, "")
+ toc = timeit.default_timer()
+ print(toc-tic)
df = pandas.DataFrame(curveJSON)
df['z0 / mole frac.'] = df['rho0 / mol/m^3']/(df['rho0 / mol/m^3']+df['rho1 / mol/m^3'])
print(df['z0 / mole frac.'])
--
GitLab