diff --git a/CMakeLists.txt b/CMakeLists.txt index ce611fa8cb127548cbc725b22bd58f22341b1774..73fca8824705e2f49e7c0693f2d049f28f841877 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -8,10 +8,9 @@ include_directories("${CMAKE_CURRENT_SOURCE_DIR}/externals/mcx/pymcx/include") include_directories("${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen") include_directories("${CMAKE_CURRENT_SOURCE_DIR}/externals/nlohmann_json") -add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen") - -set(Eigen3_DIR "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen" CACHE INTERNAL "Path to Eigen") -add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/autodiff") +#add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen") +#set(Eigen3_DIR "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen" CACHE INTERNAL "Path to Eigen") +#add_subdirectory("${CMAKE_CURRENT_SOURCE_DIR}/externals/autodiff") ### TARGETS add_executable(main "${CMAKE_CURRENT_SOURCE_DIR}/src/main.cpp") @@ -25,4 +24,4 @@ add_executable(multifluid "${CMAKE_CURRENT_SOURCE_DIR}/src/multifluid.cpp") if (MSVC) target_sources(multifluid PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/externals/Eigen/debug/msvc/eigen.natvis") endif() -target_link_libraries(multifluid autodiff) \ No newline at end of file +#target_link_libraries(multifluid autodiff) \ No newline at end of file diff --git a/src/multifluid.cpp b/src/multifluid.cpp index 9981b14b9a510f96a77d950c07f10e1a34cbe44a..8c9e61e0475a29385e9307469a6e73c37d98fa50 100644 --- a/src/multifluid.cpp +++ b/src/multifluid.cpp @@ -2,8 +2,8 @@ #include "teqp/models/multifluid.hpp" #include "teqp/critical_tracing.hpp" -#include "autodiff/forward.hpp" -#include "autodiff/reverse.hpp" +//#include "autodiff/forward.hpp" +//#include "autodiff/reverse.hpp" auto build_multifluid_model(const std::vector<std::string>& components) { using namespace nlohmann; @@ -26,10 +26,10 @@ auto build_multifluid_model(const std::vector<std::string>& components) { } void trace() { - auto model = build_multifluid_model({ "methane", "ethane" }); + auto model = build_multifluid_model({ "Methane", "Ethane" }); auto rhoc0 = 1.0/model.redfunc.vc[0]; auto T = model.redfunc.Tc[0]; - const auto dT = 1; + const auto dT = -1; std::valarray<double> rhovec = { rhoc0, 0.0 }; for (auto iter = 0; iter < 1000; ++iter) { auto drhovecdT = get_drhovec_dT_crit(model, T, rhovec); @@ -44,9 +44,49 @@ void trace() { } } +void trace_arclength() { + auto model = build_multifluid_model({ "R32", "R1234yf" }); + auto rhoc0 = 1.0 / model.redfunc.vc[1]; + auto T = model.redfunc.Tc[1]; + double t = 0.0, dt = 200; + std::valarray<double> last_drhodt; + std::valarray<double> rhovec = { 0.0, rhoc0 }; + auto norm = [](const auto &v){ return sqrt((v*v).sum()); }; + auto dot = [](const auto& v1, const auto& v2) { return (v1*v2).sum(); }; + for (auto iter = 0; iter < 1000; ++iter) { + + auto drhodT = get_drhovec_dT_crit(model, T, rhovec); + auto dTdt = 1.0 / norm(drhodT); + auto drhodt = drhodT * dTdt; + + // Flip the sign if the tracing wants to go backwards, or if the first step would take you to negative concentrations + if (iter > 0 && dot(drhodt, last_drhodt) < 0){ + drhodt *= -1; + dTdt *= -1; + } + else if (iter == 0 && any(rhovec + drhodt * dt < 0)){ + drhodt *= -1; + dTdt *= -1; + } + + rhovec += drhodt*dt; + T += dTdt*dt; + + auto rhotot = rhovec.sum(); + auto z0 = rhovec[0] / rhotot; + + std::cout << z0 << " ," << rhovec[0] << "," << T << "," << rhotot*model.R*T + get_pr(model, T, rhovec) << std::endl; + if (z0 < 0 || z0 > 1) { + break; + } + last_drhodt = drhodt; + } +} + int main(){ //test_dummy(); - trace(); + //trace(); + trace_arclength(); auto model = build_multifluid_model({ "methane", "ethane" }); std::valarray<double> rhovec = { 1.0, 2.0 }; double T = 300;