Add tests, fix bugs, remove added unused code in cython and cpp

Author: anthony-walker

include/cantera/zeroD/FlowDevice.h

@@ -140,7 +140,9 @@ class FlowDevice
     //! @since New in %Cantera 3.0.
     //!
     virtual void buildNetworkJacobian(vector<Eigen::Triplet<double>>& 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.

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;

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;

include/cantera/zeroD/ReactorBase.h

@@ -283,32 +283,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

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 *

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")

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,

src/zeroD/Reactor.cpp

@@ -619,7 +619,7 @@ void Reactor::buildFlowJacobian(vector<Eigen::Triplet<double>>& jacVector)
             m_outlet[i]->buildReactorJacobian(this, jacVector);
         }
 
-        for (size_t i = 0; i <m_outlet.size(); i++) {
+        for (size_t i = 0; i <m_inlet.size(); i++) {
             m_inlet[i]->buildReactorJacobian(this, jacVector);
         }
     }

src/zeroD/ReactorNet.cpp

@@ -642,7 +642,7 @@ void ReactorNet::buildJacobian(vector<Eigen::Triplet<double>>& 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);

test/python/test_reactor.py

@@ -198,6 +198,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()
 
@@ -1378,14 +1433,26 @@ 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()
         assert self.precon.side == "right"
         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")