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

Author: anthony-walker

include/cantera/zeroD/FlowDevice.h

@@ -189,7 +189,9 @@ class FlowDevice : public ConnectorNode
     //! @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

@@ -43,12 +43,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;

interfaces/cython/cantera/_utils.pxd

@@ -119,7 +119,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

@@ -527,41 +527,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

@@ -684,7 +684,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

@@ -803,7 +803,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

@@ -182,6 +182,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()
 
@@ -1429,14 +1484,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()
         assert 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")

test/zeroD/test_zeroD.cpp

@@ -358,6 +358,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);
@@ -382,6 +388,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<Eigen::Triplet<double>>  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");