From 278203c51b6088b4d2229ee7f2fc16d1a226b85c Mon Sep 17 00:00:00 2001 From: Ingmar Schoegl Date: Wed, 1 Mar 2023 23:27:30 -0600 Subject: [PATCH] [unittest] Add ctreactor tests --- test/clib/clib_interface.cpp | 60 +++++++++++--------- test/clib/test_ctreactor.cpp | 107 +++++++++++++++++++++++++---------- 2 files changed, 112 insertions(+), 55 deletions(-) diff --git a/test/clib/clib_interface.cpp b/test/clib/clib_interface.cpp index 33d10f60367..a97bd129574 100644 --- a/test/clib/clib_interface.cpp +++ b/test/clib/clib_interface.cpp @@ -76,7 +76,7 @@ TEST(ct, cabinet_exceptions) solution_name(999, 0, 0); string err = reportError(); - EXPECT_THAT(err, HasSubstr("Invalid index 999.")); + EXPECT_THAT(err, HasSubstr("Index 999 out of range.")); solution_del(999); err = reportError(); @@ -94,7 +94,7 @@ TEST(ct, new_solution) int ref = solution_newFromFile("h2o2.yaml", name.c_str(), ""); int buflen = solution_name(ref, 0, 0) + 1; // include \0 - ASSERT_EQ(buflen, name.size() + 1); + ASSERT_EQ(buflen, int(name.size() + 1)); char* buf = new char[buflen]; solution_name(ref, buflen, buf); @@ -143,7 +143,7 @@ TEST(ct, thermo) int ret; int thermo = thermo_newFromFile("gri30.yaml", "gri30"); ASSERT_GE(thermo, 0); - size_t nsp = thermo_nSpecies(thermo); + int nsp = thermo_nSpecies(thermo); ASSERT_EQ(nsp, 53); auto& phase = ThermoCabinet::item(thermo); @@ -170,27 +170,32 @@ TEST(ct, thermo) TEST(ct, kinetics) { int thermo = thermo_newFromFile("gri30.yaml", "gri30"); + int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, -1, -1, -1, -1); + ASSERT_GE(kin, 0); + + int nr = kin_nReactions(kin); + ASSERT_EQ(nr, 325); + thermo_equilibrate(thermo, "HP", 0, 1e-9, 50000, 1000, 0); double T = thermo_temperature(thermo); + thermo_setTemperature(thermo, T - 200); - int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, 0, 0, 0, 0); - ASSERT_GE(kin, 0); + auto sol = newSolution("gri30.yaml", "gri30", "none"); + auto phase = sol->thermo(); + auto kinetics = sol->kinetics(); - size_t nr = kin_nReactions(kin); - ASSERT_EQ(nr, 325); + phase->equilibrate("HP"); + ASSERT_NEAR(T, phase->temperature(), 1e-2); + phase->setTemperature(T - 200); - int ret = thermo_setTemperature(thermo, T - 200); - ASSERT_EQ(ret, 0); + vector c_ropf(nr); + kin_getFwdRatesOfProgress(kin, 325, c_ropf.data()); + vector cpp_ropf(nr); + kinetics->getFwdRatesOfProgress(cpp_ropf.data()); - // char buf [1000]; - // double ropf[325]; - // printf("\n Reaction Forward ROP\n"); - // kin_getFwdRatesOfProgress(kin, 325, ropf); - // int n; // declare this here for C89 compatibility - // for (n = 0; n < nr; n++) { - // kin_getReactionString(kin, n, 1000, buf); - // printf("%35s %8.6e\n", buf, ropf[n]); - // } + for (int n = 0; n < nr; n++) { + ASSERT_NEAR(cpp_ropf[n], c_ropf[n], 1e-6); + } } TEST(ct, transport) @@ -199,16 +204,19 @@ TEST(ct, transport) int tran = trans_newDefault(thermo, 0); ASSERT_GE(tran, 0); - double dkm[53]; - int ret = trans_getMixDiffCoeffs(tran, 53, dkm); + int nsp = thermo_nSpecies(thermo); + vector c_dkm(nsp); + int ret = trans_getMixDiffCoeffs(tran, 53, c_dkm.data()); ASSERT_EQ(ret, 0); - // // printf("\n Species Mix diff coeff\n"); - // int k; // declare this here for C89 compatibility - // for (k = 0; k < nsp; k++) { - // thermo_getSpeciesName(thermo, k, 1000, buf); - // printf("%10s %8.6e\n", buf, dkm[k]); - // } + vector cpp_dkm(nsp); + auto sol = newSolution("gri30.yaml", "gri30"); + auto transport = sol->transport(); + transport->getMixDiffCoeffs(cpp_dkm.data()); + + for (int n = 0; n < nsp; n++) { + ASSERT_NEAR(cpp_dkm[n], c_dkm[n], 1e-6); + } } diff --git a/test/clib/test_ctreactor.cpp b/test/clib/test_ctreactor.cpp index 6d0e4f78b80..9d0d6ba2acb 100644 --- a/test/clib/test_ctreactor.cpp +++ b/test/clib/test_ctreactor.cpp @@ -15,8 +15,10 @@ template<> NetworkCabinet* NetworkCabinet::s_storage; typedef SharedCabinet ThermoCabinet; typedef SharedCabinet KineticsCabinet; +typedef SharedCabinet SolutionCabinet; template<> ThermoCabinet* ThermoCabinet::s_storage; template<> KineticsCabinet* KineticsCabinet::s_storage; +template<> SolutionCabinet* SolutionCabinet::s_storage; TEST(ctreactor, reactor_objects) { @@ -24,7 +26,7 @@ TEST(ctreactor, reactor_objects) KineticsCabinet::reset(); int thermo = thermo_newFromFile("gri30.yaml", "gri30"); - int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, 0, 0, 0, 0); + int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, -1, -1, -1, -1); int reactor = reactor_new("IdealGasReactor"); ASSERT_GT(reactor, 0); @@ -37,48 +39,95 @@ TEST(ctreactor, reactor_objects) ASSERT_EQ(&igr.contents(), ThermoCabinet::at(thermo).get()); } -TEST(ctreactor, reactor) +vector T_blessed = { + 1050.000, 1050.064, 1050.197, 1050.369, 1050.593, 1050.881, 1051.253, 1051.736, + 1052.370, 1053.216, 1054.372, 1056.007, 1058.448, 1062.431, 1070.141, 1094.331, + 2894.921, 2894.921, 2894.921, 2894.921, 2894.921}; + +TEST(ctreactor, reactor_cpp) { - int thermo = thermo_newFromFile("gri30.yaml", "gri30"); - int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, 0, 0, 0, 0); + double T = 1050; + double P = 5 * 101325; + string X = "CH4:1.0, O2:2.0, N2:7.52"; + + auto sol = newSolution("gri30.yaml", "gri30", "none"); + sol->thermo()->setState_TPX(T, P, X); + IdealGasReactor cppReactor; + cppReactor.insert(sol); + cppReactor.initialize(); + ReactorNet network; + network.addReactor(cppReactor); + network.initialize(); + + double t = 0.0; + int count = 0; + while (t < 0.1) { + double Tref = T_blessed[count]; + ASSERT_NEAR(cppReactor.temperature(), Tref, 1e-2); + t = network.time() + 5e-3; + network.advance(t); + count++; + } +} + +TEST(ctreactor, reactor_insert) +{ + double T = 1050; + double P = 5 * 101325; + string X = "CH4:1.0, O2:2.0, N2:7.52"; - auto& phase = ThermoCabinet::item(thermo); - thermo_setTemperature(thermo, 1050); - thermo_setPressure(thermo, 5 * 101325); - ASSERT_NEAR(phase.pressure(), 5 * 101325, 1e-2); - thermo_setMoleFractionsByName(thermo, "CH4:1.0, O2:2.0, N2:7.52"); - thermo_setPressure(thermo, 5 * 101325); - ASSERT_NEAR(phase.pressure(), 5 * 101325, 1e-2); + int sol = solution_newFromFile("gri30.yaml", "gri30", "none"); + int thermo = solution_thermo(sol); + thermo_setMoleFractionsByName(thermo, X.c_str()); + thermo_setTemperature(thermo, T); + thermo_setPressure(thermo, P); int reactor = reactor_new("IdealGasReactor"); - ASSERT_GT(reactor, 0); - reactor_setThermoMgr(reactor, thermo); - reactor_setKineticsMgr(reactor, kin); + int ret = reactor_insert(reactor, sol); + ASSERT_EQ(ret, 0); int net = reactornet_new(); - int ret = reactornet_addreactor(net, reactor); + ret = reactornet_addreactor(net, reactor); ASSERT_EQ(ret, 0); - ret = thermo_print(thermo, 1, 1.e-6); + double t = 0.0; + int count = 0; + while (t < 0.1) { + double Tref = T_blessed[count]; + ASSERT_NEAR(reactor_temperature(reactor), Tref, 1e-2); + t = reactornet_time(net) + 5e-3; + ret = reactornet_advance(net, t); + ASSERT_EQ(ret, 0); + count++; + } +} - auto& igr = ReactorCabinet::item(reactor); +TEST(ctreactor, reactor_from_parts) +{ + double T = 1050; + double P = 5 * 101325; + string X = "CH4:1.0, O2:2.0, N2:7.52"; - ASSERT_NEAR(igr.temperature(), 1050, 1e-2); - ASSERT_NEAR(igr.pressure(), 5 * 101325, 1e-2); + int thermo = thermo_newFromFile("gri30.yaml", "gri30"); + int kin = kin_newFromFile("gri30.yaml", "gri30", thermo, -1, -1, -1, -1); + thermo_setMoleFractionsByName(thermo, X.c_str()); + thermo_setTemperature(thermo, T); + thermo_setPressure(thermo, P); - vector Tvec = { - 1050.000, 1050.063, 1050.192, 1050.361, 1050.578, 1050.858, 1051.216, 1051.680, - 1052.286, 1053.088, 1054.174, 1055.692, 1057.919, 1061.450, 1067.905, 1084.644, - 2894.008, 2894.004, 2894.004, 2894.004, 2894.004}; + int reactor = reactor_new("IdealGasReactor"); + reactor_setThermoMgr(reactor, thermo); + reactor_setKineticsMgr(reactor, kin); + int net = reactornet_new(); + int ret = reactornet_addreactor(net, reactor); + ASSERT_EQ(ret, 0); double t = 0.0; int count = 0; while (t < 0.1) { - double T = reactor_temperature(reactor); - T = thermo_temperature(thermo); - double Tref = Tvec[count]; - ASSERT_NEAR(T, Tref, 1e-2); - t = reactornet_time(net); - ret = reactornet_advance(net, t + 5e-3); + T = reactor_temperature(reactor); + double Tref = T_blessed[count]; + ASSERT_NEAR(reactor_temperature(reactor), Tref, 1e-2); + t = reactornet_time(net) + 5e-3; + ret = reactornet_advance(net, t); ASSERT_EQ(ret, 0); count++; }