This commit adds some simple reactor jacobian tests and fixes bugs

The indexing has been fixed and works with multiple surfaces correctly avoiding
non-interacting phases. Tests were added to compare simple surface jacobians, gas
phase jacobians, and single reactions. Removed surface contribution to energy
equation terms as these contributions are not currently included in the system
of equations and caus spurious entries.

data/simple_surface.yaml

@@ -0,0 +1,118 @@
+units: {length: cm, quantity: mol, activation-energy: J/mol}
+
+phases:
+- name: gas
+  thermo: ideal-gas
+  species: [A, B, C, D]
+  kinetics: gas
+  reactions: [gas-reactions]
+  state:
+    T: 300.0
+    P: 1.01325e+05
+- name: surf
+  thermo: ideal-surface
+  adjacent-phases: [gas]
+  species: [A(S), B(S), C(S), D(S), (S)]
+  kinetics: surface
+  reactions: [surf-reactions]
+  state:
+    T: 900.0
+    coverages: {(S): 1.0}
+  site-density: 2.7063e-09
+
+species:
+- name: A
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: B
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: C
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: D
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: A(S)
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: B(S)
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: C(S)
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: D(S)
+  composition: {C: 1}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+- name: (S)
+  composition: {}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1000.0, 3000.0]
+    data:
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+    - [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
+
+gas-reactions:
+- equation: A + B <=> C + D
+  rate-constant: {A: 1200, b: 0, Ea: 0}
+- equation: A + C <=> B + D
+  rate-constant: {A: 800, b: 0, Ea: 0}
+
+surf-reactions:
+- equation: A + (S) => A(S)
+  rate-constant: {A: 5e3, b: 0, Ea: 0}
+- equation: B + (S) => B(S)
+  rate-constant: {A: 2e4, b: 0, Ea: 0}
+- equation: A(S) => A + (S)
+  rate-constant: {A: 500, b: 0, Ea: 0}
+- equation: B(S) => B + (S)
+  rate-constant: {A: 800, b: 0, Ea: 0}
+- equation: C(S) => C + (S)
+  rate-constant: {A: 500, b: 0, Ea: 0}
+- equation: D(S) => D + (S)
+  rate-constant: {A: 400, b: 0, Ea: 0}
+
+- equation: A(S) + B(S) => C(S) + D(S)
+  rate-constant: {A: 1200, b: 0, Ea: 0}

src/zeroD/IdealGasConstPressureMoleReactor.cpp

@@ -205,26 +205,15 @@ Eigen::SparseMatrix<double> IdealGasConstPressureMoleReactor::jacobian()
         }
         // d T_dot/dnj
         // allocate vectors for whole system
-        Eigen::VectorXd enthalpy(ssize);
-        Eigen::VectorXd specificHeat(ssize);
+        Eigen::VectorXd enthalpy = Eigen::VectorXd::Zero(ssize);
+        Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
         // gas phase
         m_thermo->getPartialMolarCp(specificHeat.data());
         m_thermo->getPartialMolarEnthalpies(enthalpy.data());
         // scale production rates by the volume for gas species
         for (size_t i = 0; i < m_nsp; i++) {
             netProductionRates[i] *= m_vol;
         }
-        // surface phases
-        size_t shift = m_nsp;
-        for (auto S : m_surfaces) {
-            S->thermo()->getPartialMolarCp(specificHeat.data() + shift);
-            S->thermo()->getPartialMolarEnthalpies(enthalpy.data() + shift);
-            // scale production rates by areas
-            for (size_t i = shift; i < shift + S->thermo()->nSpecies(); i++) {
-                netProductionRates[i] *= S->area();
-            }
-            shift += S->thermo()->nSpecies();
-        }
         double qdot = enthalpy.dot(netProductionRates);
         // find denominator ahead of time
         double NCp = 0.0;

src/zeroD/IdealGasMoleReactor.cpp

@@ -169,7 +169,7 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
     // map derivatives from the surface chemistry jacobian
     // to the reactor jacobian
     if (!m_surfaces.empty()) {
-        std::vector<Eigen::Triplet<double>> species_trips(dnk_dnj.nonZeros());
+        std::vector<Eigen::Triplet<double>> species_trips;
         for (int k = 0; k < dnk_dnj.outerSize(); k++) {
             for (Eigen::SparseMatrix<double>::InnerIterator it(dnk_dnj, k); it; ++it) {
                 species_trips.emplace_back(it.row(), it.col(), it.value());
@@ -217,33 +217,12 @@ Eigen::SparseMatrix<double> IdealGasMoleReactor::jacobian()
         }
         // d T_dot/dnj
         Eigen::VectorXd netProductionRates = Eigen::VectorXd::Zero(ssize);
-        Eigen::VectorXd internal_energy(ssize);
-        Eigen::VectorXd specificHeat(ssize);
+        Eigen::VectorXd internal_energy = Eigen::VectorXd::Zero(ssize);
+        Eigen::VectorXd specificHeat = Eigen::VectorXd::Zero(ssize);
         // getting species data
         m_thermo->getPartialMolarIntEnergies(internal_energy.data());
         m_kin->getNetProductionRates(netProductionRates.data());
         m_thermo->getPartialMolarCp(specificHeat.data());
-        // surface phases
-        size_t offset = m_nsp;
-        for (auto S : m_surfaces) {
-            S->thermo()->getPartialMolarCp(specificHeat.data() + offset);
-            S->thermo()->getPartialMolarEnthalpies(internal_energy.data() + offset);
-            // get additional net production rates
-            auto curr_kin = S->kinetics();
-            vector_fp prod_rates(curr_kin->nTotalSpecies());
-            curr_kin->getNetProductionRates(prod_rates.data());
-            for (size_t i = 0; i < curr_kin->nTotalSpecies(); i++) {
-                size_t row = speciesIndex(curr_kin->kineticsSpeciesName(i));
-                if (row != npos) {
-                    netProductionRates[row] += prod_rates[i];
-                }
-            }
-            // scale production rates by areas
-            for (size_t i = offset; i < offset + S->thermo()->nSpecies(); i++) {
-                netProductionRates[i] *= S->area();
-            }
-            offset += S->thermo()->nSpecies();
-        }
         // convert Cp to Cv for ideal gas as Cp - Cv = R
         for (size_t i = 0; i < m_nsp; i++) {
             specificHeat[i] -= GasConstant;

src/zeroD/MoleReactor.cpp

@@ -92,10 +92,10 @@ void MoleReactor::addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets)
         double A = S->area();
         auto kin = S->kinetics();
         size_t nk = S->thermo()->nSpecies();
-        // limit on the current idx
-        size_t limit = offset + nk;
-        // find location of first gas phase species in the kinetics object
-        size_t gasIdx = kin->kineticsSpeciesIndex(m_thermo->speciesName(0));
+        // index of gas and surface phases to check if the species is in gas or surface
+        int spi = kin->reactionPhaseIndex();
+        int gpi = kin->speciesPhaseIndex(kin->kineticsSpeciesIndex(
+            m_thermo->speciesName(0)));
         // get surface jacobian in concentration units
         Eigen::SparseMatrix<double> surfJac = kin->netProductionRates_ddCi();
         // loop through surface specific jacobian and add elements to triplets vector
@@ -104,19 +104,31 @@ void MoleReactor::addSurfaceJacobian(vector<Eigen::Triplet<double>> &triplets)
             for (Eigen::SparseMatrix<double>::InnerIterator it(surfJac, k); it; ++it) {
                 size_t row = it.row();
                 size_t col = it.col();
-                // check gas location and number of species against row and col to
-                // avoid any solid phases which may be included and
-                // map surf row and col to reactor by adding the offset if the index
-                // is less than the species in the SurfPhase or subtracting the
-                // number of SurfPhase species otherwise
-                row = (row >= gasIdx && row < limit) ? row - gasIdx : row + offset;
-                col = (col >= gasIdx && col < limit) ? col - gasIdx : col + offset;
-                if (row < limit && col < limit) {
-                    // determine appropriate scalar by location
+                auto& rowPhase = kin->speciesPhase(row);
+                auto& colPhase = kin->speciesPhase(col);
+                int rpi = kin->phaseIndex(rowPhase.name());
+                int cpi = kin->phaseIndex(colPhase.name());
+                // check if the reactor kinetics object contains both phases to avoid
+                // any solid phases which may be included then use phases to map surf
+                // kinetics indicies to reactor kinetic indices
+                if ((rpi == spi || rpi == gpi) && (cpi == spi || cpi == gpi) ) {
+                    // subtract start of phase
+                    row -= kin->kineticsSpeciesIndex(0, rpi);
+                    col -= kin->kineticsSpeciesIndex(0, cpi);
+                    // since the gas phase is the first phase in the reactor state
+                    // vector add the offset only if it is a surf species index to
+                    // both row and col
+                    row = (rpi != spi) ? row : row + offset;
+                    // determine appropriate scalar to account for dimensionality
                     // gas phase species indices will be less than m_nsp
                     // so use volume if that is the case or area otherwise
                     double scalar = A;
-                    scalar /= (col < m_nsp) ? m_vol : A;
+                    if (cpi == spi) {
+                        col += offset;
+                        scalar /= A;
+                    } else {
+                        scalar /= m_vol;
+                    }
                     // push back scaled value triplet
                     triplets.emplace_back(row, col, scalar * it.value());
                 }

test/python/test_reactor.py

@@ -1208,6 +1208,99 @@ def test_adaptive_precon_integration(self):
 
 class TestReactorJacobians(utilities.CanteraTest):
 
+    def test_multi_surface_simple(self):
+        # conditions for simulation
+        yml = "simple_surface.yaml"
+        fuel = "A:1.0, B:1.0"
+        # gas kinetics
+        gas = ct.Solution(yml, "gas")
+        gas.TPX = 1000, 2e5, fuel
+        gas.set_multiplier(0)
+        # surface kinetics for the simulation
+        surf = ct.Interface(yml, 'surf', [gas])
+        surf2 = ct.Interface(yml, 'surf', [gas])
+        surf.coverages = 'A(S):0.1, B(S):0.2, C(S):0.3, D(S):0.2, (S):0.2'
+        surf2.coverages = 'A(S):0.1, D(S):0.2, (S):0.2'
+        # create reactor
+        r = ct.IdealGasMoleReactor(gas)
+        r.volume = 3
+        # create surfaces
+        rsurf1 = ct.ReactorSurface(surf, r, A=9e-4)
+        rsurf2 = ct.ReactorSurface(surf2, r, A=5e-4)
+        # create network
+        net = ct.ReactorNet([r])
+        net.step()
+        # get jacobians
+        jacobian = r.jacobian
+        fd_jacobian = r.finite_difference_jacobian
+        # the volume row is not considered in comparisons because it is presently
+        # not calculated.
+        # check first row is near, terms which are generally on the order of 1e5 to 1e7
+        self.assertArrayNear(jacobian[0, 2:], fd_jacobian[0, 2:], 1e-1, 1e-2)
+        # check first col is near, these are finite difference terms and should be close
+        self.assertArrayNear(jacobian[2:, 0], fd_jacobian[2:, 0], 1e-3, 1e-4)
+        # check all species are near, these terms are usually ~ 1e2
+        self.assertArrayNear(jacobian[2:, 2:], fd_jacobian[2:, 2:], 1e-3, 1e-4)
+
+    def test_gas_simple(self):
+        # conditions for simulation
+        yml = "simple_surface.yaml"
+        fuel = "A:1.0, B:1.0, C:1.0, D:1.0"
+        # gas kinetics
+        gas = ct.Solution(yml, "gas")
+        gas.TPX = 1000, 2e5, fuel
+        # create reactor
+        r = ct.IdealGasMoleReactor(gas)
+        r.volume = 1
+        # create network
+        net = ct.ReactorNet([r])
+        net.initialize()
+        # compare jacobians
+        self.assertArrayNear(r.jacobian, r.finite_difference_jacobian, rtol=1e-6)
+
+    def test_const_volume_hydrogen_single(self):
+        # conditions for simulation
+        yml = "h2o2.yaml"
+        # gas kinetics
+        gas = ct.Solution(yml, "ohmech")
+        gas.TPX = 1000, ct.one_atm, "H2:1.0, H:2.0, H2O:2.0"
+        gas.set_multiplier(0)
+        gas.set_multiplier(1, 14)
+        # create reactor
+        r = ct.IdealGasMoleReactor(gas)
+        r.volume = 2
+        # create network
+        net = ct.ReactorNet([r])
+        net.step()
+        # get jacobians
+        jacobian = r.jacobian
+        fd_jacobian = r.finite_difference_jacobian
+        # the volume row is not considered in comparisons because it is presently
+        # not calculated.
+        # check first row is near, terms which are generally on the order of 1e5 to 1e7
+        self.assertArrayNear(jacobian[0, 2:], fd_jacobian[0, 2:], 1e-2, 1e-3)
+        # check first col is near, these are finite difference terms and should be close
+        self.assertArrayNear(jacobian[2:, 0], fd_jacobian[2:, 0], 1e-3, 1e-4)
+        # check all species are near, these terms are usually ~ 1e2
+        self.assertArrayNear(jacobian[2:, 2:], fd_jacobian[2:, 2:], 1e-3, 1e-4)
+
+    def test_const_pressure_hydrogen_single(self):
+        # conditions for simulation
+        yml = "h2o2.yaml"
+        # gas kinetics
+        gas = ct.Solution(yml, "ohmech")
+        gas.TPX = 1000, ct.one_atm, "H2:1.0, H:2.0, H2O:2.0"
+        gas.set_multiplier(0)
+        gas.set_multiplier(1, 14)
+        # create reactor
+        r = ct.IdealGasConstPressureMoleReactor(gas)
+        r.volume = 2
+        # create network
+        net = ct.ReactorNet([r])
+        net.step()
+        # the volume row is not considered in comparisons because it is presently
+        self.assertArrayNear(r.jacobian, r.finite_difference_jacobian, 1e-3, 1e-4)
+
     def test_coverage_dependence_flags(self):
         gas = ct.Solution("ptcombust.yaml", "gas")
         surf = ct.Interface("ptcombust.yaml", "Pt_surf", [gas])