diff --git a/include/teqp/algorithms/VLE.hpp b/include/teqp/algorithms/VLE.hpp index 2514d7d138c8f9450fd98dd7b303a9c4f2fa5737..903582459f7bee6cf7340de3907158bd24cf6e1c 100644 --- a/include/teqp/algorithms/VLE.hpp +++ b/include/teqp/algorithms/VLE.hpp @@ -63,11 +63,11 @@ public: }; template<typename Residual, typename Scalar> -auto do_pure_VLE(Residual &resid, Scalar T, Scalar rhoL, Scalar rhoV) { +auto do_pure_VLE_T(Residual &resid, Scalar T, Scalar rhoL, Scalar rhoV, int maxiter) { auto rhovec = (Eigen::ArrayXd(2) << rhoL, rhoV).finished(); auto r0 = resid.call(rhovec); auto J = resid.Jacobian(rhovec); - for (int iter = 0; iter < 100; ++iter){ + for (int iter = 0; iter < maxiter; ++iter){ if (iter > 0) { r0 = resid.call(rhovec); J = resid.Jacobian(rhovec); @@ -76,11 +76,12 @@ auto do_pure_VLE(Residual &resid, Scalar T, Scalar rhoL, Scalar rhoV) { auto rhovecnew = (rhovec + v).eval(); // If the solution has stopped improving, stop. The change in rhovec is equal to v in infinite precision, but - // not when finite precision is involved, so use the minimum non-denormal float as the determination of whether + // not when finite precision is involved, use the minimum non-denormal float as the determination of whether // the values are done changing if (((rhovecnew - rhovec).cwiseAbs() < std::numeric_limits<Scalar>::min()).all()) { break; } rhovec += v; } + return std::make_tuple(rhovec[0], rhovec[1]); } \ No newline at end of file diff --git a/src/tests/catch_tests.cxx b/src/tests/catch_tests.cxx index e6999ab04064669fa7f56273261ef8d7602455e4..ee9b009d6c02dfb01784e9fcba19d748b7184461 100644 --- a/src/tests/catch_tests.cxx +++ b/src/tests/catch_tests.cxx @@ -243,5 +243,5 @@ TEST_CASE("Test pure VLE", "") { auto r1 = resid.call(rhovec1); CHECK((r0.cwiseAbs() > r1.cwiseAbs()).eval().all()); - do_pure_VLE(resid, T, 22834.056386882046, 1025.106554560764); + do_pure_VLE_T(resid, T, 22834.056386882046, 1025.106554560764, 20); } \ No newline at end of file