diff --git a/include/cantera/zeroD/FlowDevice.h b/include/cantera/zeroD/FlowDevice.h index 5bf946261e7..ffa61f66b71 100644 --- a/include/cantera/zeroD/FlowDevice.h +++ b/include/cantera/zeroD/FlowDevice.h @@ -140,7 +140,9 @@ class FlowDevice //! @since New in %Cantera 3.0. //! virtual void buildNetworkJacobian(vector>& jacVector) { - throw NotImplementedError(type() + "::buildNetworkJacobian"); + if (!m_jac_calculated) { + throw NotImplementedError(type() + "::buildNetworkJacobian"); + } } //! Specify the jacobian terms have been calculated and should not be recalculated. diff --git a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h index 3064a641134..e0aa1231f77 100644 --- a/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h +++ b/include/cantera/zeroD/IdealGasConstPressureMoleReactor.h @@ -49,12 +49,6 @@ class IdealGasConstPressureMoleReactor : public ConstPressureMoleReactor bool preconditionerSupported() const override { return true; }; - double temperatureDerivative() override { return 1; }; - - double temperatureRadiationDerivative() override { - return 4 * std::pow(temperature(), 3); - } - double moleDerivative(size_t index) override; double moleRadiationDerivative(size_t index) override; diff --git a/include/cantera/zeroD/IdealGasMoleReactor.h b/include/cantera/zeroD/IdealGasMoleReactor.h index 82cc9f70fbc..c25248dde90 100644 --- a/include/cantera/zeroD/IdealGasMoleReactor.h +++ b/include/cantera/zeroD/IdealGasMoleReactor.h @@ -45,12 +45,6 @@ class IdealGasMoleReactor : public MoleReactor bool preconditionerSupported() const override {return true;}; - double temperatureDerivative() override { return 1; }; - - double temperatureRadiationDerivative() override { - return 4 * std::pow(temperature(), 3); - } - double moleDerivative(size_t index) override; double moleRadiationDerivative(size_t index) override; diff --git a/include/cantera/zeroD/ReactorBase.h b/include/cantera/zeroD/ReactorBase.h index 29b7b495cd0..9513aad3be1 100644 --- a/include/cantera/zeroD/ReactorBase.h +++ b/include/cantera/zeroD/ReactorBase.h @@ -248,32 +248,6 @@ class ReactorBase //! Set the ReactorNet that this reactor belongs to. void setNetwork(ReactorNet* net); - //! Calculate the derivative of temperature with respect to the temperature in the - //! heat transfer equation based on the reactor specific equation of state. - //! This function should also transform the state of the derivative to that - //! appropriate for the jacobian's state/ - //! @warning This function is an experimental part of the %Cantera API and may be changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - virtual double temperatureDerivative() { - throw NotImplementedError("Reactor::temperatureDerivative"); - } - - //! Calculate the derivative of temperature with respect to the temperature in the - //! heat transfer radiation equation based on the reactor specific equation of - //! state. - //! This function should also transform the state of the derivative to that - //! appropriate for the jacobian's state/ - //! @warning This function is an experimental part of the %Cantera API and may be - //! changed - //! or removed without notice. - //! @since New in %Cantera 3.0. - //! - virtual double temperatureRadiationDerivative() { - throw NotImplementedError("Reactor::temperatureRadiationDerivative"); - } - //! Calculate the derivative of T with respect to the ith species in the heat //! transfer equation based on the reactor specific equation of state. //! @param index index of the species the derivative is with respect too diff --git a/interfaces/cython/cantera/_utils.pxd b/interfaces/cython/cantera/_utils.pxd index ab5808cf67f..bf24486a941 100644 --- a/interfaces/cython/cantera/_utils.pxd +++ b/interfaces/cython/cantera/_utils.pxd @@ -110,7 +110,4 @@ cdef anymap_to_py(CxxAnyMap& m) cdef CxxAnyValue python_to_anyvalue(item, name=*) except * cdef anyvalue_to_python(string name, CxxAnyValue& v) -cdef vector[CxxEigenTriplet] python_to_triplets(triplets) except * -cdef triplets_to_python(vector[CxxEigenTriplet]& triplets) - cdef CxxEigenTriplet get_triplet(row, col, val) except * diff --git a/interfaces/cython/cantera/_utils.pyx b/interfaces/cython/cantera/_utils.pyx index dae7e2f033f..2ce1e18539e 100644 --- a/interfaces/cython/cantera/_utils.pyx +++ b/interfaces/cython/cantera/_utils.pyx @@ -520,41 +520,10 @@ cdef vector[vector[string]] list2_string_to_anyvalue(data): v[i][j] = stringify(jtem) return v -cdef vector[CxxEigenTriplet] python_to_triplets(triplets): - # check that triplet dimensions are met - trip_shape = np.shape(triplets) - assert len(trip_shape) == 2 - assert trip_shape[1] == 3 - cdef vector[CxxEigenTriplet] trips - # c++ variables - cdef size_t row - cdef size_t col - cdef double val - cdef CxxEigenTriplet* trip_ptr - for r, c, v in triplets: - row = r - col = c - val = v - trip_ptr = new CxxEigenTriplet(row, col, val) - trips.push_back(dereference(trip_ptr)) - return trips - cdef CxxEigenTriplet get_triplet(row, col, val): cdef CxxEigenTriplet* trip_ptr = new CxxEigenTriplet(row, col, val) return dereference(trip_ptr) -cdef triplets_to_python(vector[CxxEigenTriplet]& triplets): - values = [] - cdef size_t row - cdef size_t col - cdef double val - for t in triplets: - row = t.row() - col = t.col() - val = t.value() - values.append((row, col, val)) - return values - def _py_to_any_to_py(dd): # used for internal testing purposes only cdef string name = stringify("test") diff --git a/interfaces/cython/cantera/delegator.pyx b/interfaces/cython/cantera/delegator.pyx index 46aa03f9af5..9ecaf95ebe5 100644 --- a/interfaces/cython/cantera/delegator.pyx +++ b/interfaces/cython/cantera/delegator.pyx @@ -9,7 +9,7 @@ from libc.stdlib cimport malloc from libc.string cimport strcpy from ._utils import CanteraError -from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, triplets_to_python, python_to_triplets, get_triplet +from ._utils cimport stringify, pystr, anymap_to_py, py_to_anymap, get_triplet from .units cimport Units, UnitStack # from .reaction import ExtensibleRate, ExtensibleRateData from .reaction cimport (ExtensibleRate, ExtensibleRateData, CxxReaction, diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp index 5dfd31433b9..71b55deaae3 100644 --- a/src/zeroD/Reactor.cpp +++ b/src/zeroD/Reactor.cpp @@ -616,7 +616,7 @@ void Reactor::buildFlowJacobian(vector>& jacVector) m_outlet[i]->buildReactorJacobian(this, jacVector); } - for (size_t i = 0; i buildReactorJacobian(this, jacVector); } } diff --git a/src/zeroD/ReactorNet.cpp b/src/zeroD/ReactorNet.cpp index a037f9aee76..dd4cb18de82 100644 --- a/src/zeroD/ReactorNet.cpp +++ b/src/zeroD/ReactorNet.cpp @@ -642,7 +642,7 @@ void ReactorNet::buildJacobian(vector>& jacVector) } } // flow devices - if (m_jac_skip_flow_devices) { + if (!m_jac_skip_flow_devices) { // outlets for (size_t i = 0; i < r->nOutlets(); i++) { r->outlet(i).buildNetworkJacobian(jacVector); diff --git a/test/python/test_reactor.py b/test/python/test_reactor.py index b8ea200e89b..9f66359d56d 100644 --- a/test/python/test_reactor.py +++ b/test/python/test_reactor.py @@ -188,6 +188,61 @@ def test_finite_difference_jacobian(self): if name in constant: assert all(J[i, species_start:] == 0), (i, name) + def test_network_finite_difference_jacobian(self): + self.make_reactors(T1=900, P1=101325, X1="H2:0.4, O2:0.4, N2:0.2") + k1H2 = self.gas1.species_index("H2") + k2H2 = self.gas1.species_index("H2") + while self.r1.thermo.X[k1H2] > 0.3 or self.r2.thermo.X[k2H2] > 0.3: + self.net.step() + + J = self.net.finite_difference_jacobian + assert J.shape == (self.net.n_vars, self.net.n_vars) + + # state variables that should be constant, depending on reactor type + constant = {"mass", "volume", "int_energy", "enthalpy", "pressure"} + variable = {"temperature"} + for i in range(3): + name = self.r1.component_name(i) + if name in constant: + assert all(J[i,:] == 0), (i, name) + elif name in variable: + assert any(J[i,:] != 0) + # check in second reactor + name = self.r2.component_name(i) + if name in constant: + assert all(J[i + self.r1.n_vars,:] == 0), (i, name) + elif name in variable: + assert any(J[i + self.r1.n_vars,:] != 0) + + # Disabling energy equation should zero these terms + self.r1.energy_enabled = False + self.r2.energy_enabled = False + J = self.net.finite_difference_jacobian + for i in range(3): + name = self.r1.component_name(i) + if name == "temperature": + assert all(J[i,:] == 0) + name = self.r2.component_name(i) + if name == "temperature": + assert all(J[i + self.r1.n_vars,:] == 0) + + # Disabling species equations should zero these terms + self.r1.energy_enabled = True + self.r1.chemistry_enabled = False + self.r2.energy_enabled = True + self.r2.chemistry_enabled = False + J = self.net.finite_difference_jacobian + constant = set(self.gas1.species_names + self.gas2.species_names) + r1_species_start = self.r1.component_index(self.gas1.species_name(0)) + r2_species_start = self.r2.component_index(self.gas2.species_name(0)) + for i in range(self.r1.n_vars): + name = self.r1.component_name(i) + if name in constant: + assert all(J[i, r1_species_start:] == 0), (i, name) + name = self.r2.component_name(i) + if name in constant: + assert all(J[i + self.r1.n_vars, (r2_species_start + self.r1.n_vars):] == 0), (i, name) + def test_timestepping(self): self.make_reactors() @@ -1187,7 +1242,7 @@ def create_reactors(self, **kwargs): self.precon = ct.AdaptivePreconditioner() self.net2.preconditioner = self.precon self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True, - "skip-coverage-dependence":True} + "skip-coverage-dependence":True, "skip-flow-devices": True} def test_get_solver_type(self): self.create_reactors() @@ -1195,6 +1250,18 @@ def test_get_solver_type(self): self.net2.initialize() self.assertEqual(self.net2.linear_solver_type, "GMRES") + def test_mass_flow_jacobian(self): + self.create_reactors(add_mdot=True) + # reset derivative settings + self.net2.derivative_settings = {"skip-third-bodies":True, "skip-falloff":True, + "skip-coverage-dependence":True, "skip-flow-devices": False} + + with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"): + J = self.net2.jacobian + + with pytest.raises(NotImplementedError, match="MassFlowController::buildReactorJacobian"): + J = self.r2.jacobian + def test_heat_transfer_network(self): # create first reactor gas1 = ct.Solution("h2o2.yaml", "ohmech") diff --git a/test/zeroD/test_zeroD.cpp b/test/zeroD/test_zeroD.cpp index bcd8a0823e0..aa5574c3585 100644 --- a/test/zeroD/test_zeroD.cpp +++ b/test/zeroD/test_zeroD.cpp @@ -247,6 +247,12 @@ TEST(JacobianTests, test_wall_jacobian_build) EXPECT_LT(it.col(), reactor2.neq()); } } + // check that switch works + wallJac.clear(); + w.jacobianCalculated(); + w.buildNetworkJacobian(wallJac); + EXPECT_EQ(wallJac.size(), 0); + w.jacobianNotCalculated(); // build jac for network terms wallJac.clear(); w.buildNetworkJacobian(wallJac); @@ -271,6 +277,37 @@ TEST(JacobianTests, test_wall_jacobian_build) } } +TEST(JacobianTests, test_flow_jacobian_build) +{ + // create reservoir reactor + auto sol = newSolution("h2o2.yaml"); + sol->thermo()->setState_TPY(1000.0, OneAtm, "O2:1.0"); + Reservoir res; + res.insert(sol); + // create reactor + IdealGasConstPressureMoleReactor reactor; + reactor.insert(sol); + reactor.setInitialVolume(1.0); + // create the flow device + MassFlowController mfc; + mfc.install(res, reactor); + mfc.setMassFlowCoeff(1.0); + // setup reactor network and integrate + ReactorNet network; + network.addReactor(reactor); + network.initialize(); + // manually get wall jacobian elements + vector> flowJac; + // expect errors from building jacobians + EXPECT_THROW(mfc.buildReactorJacobian(&reactor, flowJac), NotImplementedError); + // check the jacobian calculated flag and throw/catch errors accordingly + mfc.jacobianCalculated(); + mfc.buildNetworkJacobian(flowJac); + EXPECT_EQ(flowJac.size(), 0); + mfc.jacobianNotCalculated(); + EXPECT_THROW(mfc.buildNetworkJacobian(flowJac), NotImplementedError); +} + int main(int argc, char** argv) { printf("Running main() from test_zeroD.cpp\n");