Improve test coverage of transport calculations and fix related bugs

include/cantera/transport/GasTransport.h

@@ -154,6 +154,8 @@ class GasTransport : public Transport
         return m_mode == CK_Mode;
     }
 
+    void invalidateCache() override;
+
 protected:
     GasTransport();
 

include/cantera/transport/MultiTransport.h

@@ -102,6 +102,8 @@ class MultiTransport : public GasTransport
 
     void init(ThermoPhase* thermo, int mode=0, int log_level=-7) override;
 
+    void invalidateCache() override;
+
 protected:
     //! Update basic temperature-dependent quantities if the temperature has
     //! changed.

include/cantera/transport/Transport.h

@@ -427,6 +427,11 @@ class Transport
 
     //! @}
 
+    //! Invalidate any cached values which are normally updated only when a
+    //! change in state is detected
+    //! @since New in %Cantera 3.1.
+    virtual void invalidateCache() {}
+
 protected:
     //! pointer to the object representing the phase
     ThermoPhase* m_thermo;

interfaces/cython/cantera/transport.pyx

@@ -308,7 +308,7 @@ cdef class Transport(_SolutionBase):
     def get_collision_integral_polynomials(self, i, j):
         """
         Get the coefficients of the polynomial fit of the collision integral for
-        species ``i`` and ``j`` as a function of temperature. See
+        species ``i`` and ``j`` as a function of the reduced temperature. See
         :ct:`GasTransport::fitCollisionIntegrals`.
         """
         n_values = 7 if self.transport.CKMode() else 9
@@ -365,11 +365,13 @@ cdef class Transport(_SolutionBase):
         tran_setBinDiffusivityPolynomial(self.transport, i, j, &data[0])
 
     def set_collision_integral_polynomial(self, i, j, avalues, bvalues, cvalues,
-                                          actualT=True):
-        """
+                                          actualT=False):
+        r"""
         Set the coefficients of the polynomial fit of the collision integral for
         species ``i`` and ``j`` as a function of temperature. See
-        :ct:`GasTransport::fitCollisionIntegrals`.
+        :ct:`GasTransport::fitCollisionIntegrals`. The flag ``actualT`` determines
+        whether the fit is done in terms of the reduced temperature
+        (:math:`T^*_{ij} = T k_B / \epsilon_{ij}`; default) or the actual temperature.
         """
         n_values = 7 if self.transport.CKMode() else 9
         if len(avalues) != n_values:

src/transport/DustyGasTransport.cpp

@@ -162,11 +162,7 @@ void DustyGasTransport::updateMultiDiffCoeffs()
     eval_H_matrix();
 
     // invert H
-    int ierr = invert(m_multidiff);
-    if (ierr != 0) {
-        throw CanteraError("DustyGasTransport::updateMultiDiffCoeffs",
-                           "invert returned ierr = {}", ierr);
-    }
+    invert(m_multidiff);
 }
 
 void DustyGasTransport::getMultiDiffCoeffs(const size_t ld, double* const d)

src/transport/GasTransport.cpp

@@ -357,6 +357,16 @@ void GasTransport::setupCollisionParameters()
     }
 }
 
+void GasTransport::invalidateCache()
+{
+    Transport::invalidateCache();
+    m_temp = Undef;
+    m_visc_ok = false;
+    m_spvisc_ok = false;
+    m_viscwt_ok = false;
+    m_bindiff_ok = false;
+}
+
 void GasTransport::setupCollisionIntegral()
 {
     double tstar_min = 1.e8, tstar_max = 0.0;
@@ -818,20 +828,24 @@ void GasTransport::getBinDiffCorrection(double t, MMCollisionInt& integrals,
 
 void GasTransport::getViscosityPolynomial(size_t i, double* coeffs) const
 {
+    checkSpeciesIndex(i);
     for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
         coeffs[k] = m_visccoeffs[i][k];
     }
 }
 
 void GasTransport::getConductivityPolynomial(size_t i, double* coeffs) const
 {
+    checkSpeciesIndex(i);
     for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
         coeffs[k] = m_condcoeffs[i][k];
     }
 }
 
 void GasTransport::getBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs) const
 {
+    checkSpeciesIndex(i);
+    checkSpeciesIndex(j);
     size_t mi = (j >= i? i : j);
     size_t mj = (j >= i? j : i);
     size_t ic = 0;
@@ -850,6 +864,8 @@ void GasTransport::getCollisionIntegralPolynomial(size_t i, size_t j,
                                                  double* bstar_coeffs,
                                                  double* cstar_coeffs) const
 {
+    checkSpeciesIndex(i);
+    checkSpeciesIndex(j);
     for (int k = 0; k < (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE) + 1; k++) {
         astar_coeffs[k] = m_astar_poly[m_poly[i][j]][k];
         bstar_coeffs[k] = m_bstar_poly[m_poly[i][j]][k];
@@ -859,32 +875,26 @@ void GasTransport::getCollisionIntegralPolynomial(size_t i, size_t j,
 
 void GasTransport::setViscosityPolynomial(size_t i, double* coeffs)
 {
+    checkSpeciesIndex(i);
     for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
         m_visccoeffs[i][k] = coeffs[k];
     }
-
-    m_visc_ok = false;
-    m_spvisc_ok = false;
-    m_viscwt_ok = false;
-    m_bindiff_ok = false;
-    m_temp = -1;
+    invalidateCache();
 }
 
 void GasTransport::setConductivityPolynomial(size_t i, double* coeffs)
 {
+    checkSpeciesIndex(i);
     for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
         m_condcoeffs[i][k] = coeffs[k];
     }
-
-    m_visc_ok = false;
-    m_spvisc_ok = false;
-    m_viscwt_ok = false;
-    m_bindiff_ok = false;
-    m_temp = -1;
+    invalidateCache();
 }
 
 void GasTransport::setBinDiffusivityPolynomial(size_t i, size_t j, double* coeffs)
 {
+    checkSpeciesIndex(i);
+    checkSpeciesIndex(j);
     size_t mi = (j >= i? i : j);
     size_t mj = (j >= i? j : i);
     size_t ic = 0;
@@ -896,19 +906,16 @@ void GasTransport::setBinDiffusivityPolynomial(size_t i, size_t j, double* coeff
     for (int k = 0; k < (m_mode == CK_Mode ? 4 : 5); k++) {
         m_diffcoeffs[ic][k] = coeffs[k];
     }
-
-    m_visc_ok = false;
-    m_spvisc_ok = false;
-    m_viscwt_ok = false;
-    m_bindiff_ok = false;
-    m_temp = -1;
+    invalidateCache();
 }
 
 void GasTransport::setCollisionIntegralPolynomial(size_t i, size_t j,
                                                   double* astar_coeffs,
                                                   double* bstar_coeffs,
                                                   double* cstar_coeffs, bool actualT)
 {
+    checkSpeciesIndex(i);
+    checkSpeciesIndex(j);
     size_t degree = (m_mode == CK_Mode ? 6 : COLL_INT_POLY_DEGREE);
     vector<double> ca(degree+1), cb(degree+1), cc(degree+1);
 
@@ -921,18 +928,9 @@ void GasTransport::setCollisionIntegralPolynomial(size_t i, size_t j,
     m_astar_poly.push_back(ca);
     m_bstar_poly.push_back(cb);
     m_cstar_poly.push_back(cc);
-    m_poly[i][j] = static_cast<int>(m_astar_poly.size()) - 1;
-    m_poly[j][i] = m_poly[i][j];
-    if (actualT) {
-        m_star_poly_uses_actualT[i][j] = 1;
-        m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j];
-    }
-
-    m_visc_ok = false;
-    m_spvisc_ok = false;
-    m_viscwt_ok = false;
-    m_bindiff_ok = false;
-    m_temp = -1;
+    m_poly[j][i] = m_poly[i][j] = static_cast<int>(m_astar_poly.size()) - 1;
+    m_star_poly_uses_actualT[j][i] = m_star_poly_uses_actualT[i][j] = (actualT) ? 1 : 0;
+    invalidateCache();
 }
 
 }

src/transport/HighPressureGasTransport.cpp

@@ -242,11 +242,7 @@ void HighPressureGasTransport::getMultiDiffCoeffs(const size_t ld, double* const
     }
 
     // invert L00,00
-    int ierr = invert(m_Lmatrix, m_nsp);
-    if (ierr != 0) {
-        throw CanteraError("HighPressureGasTransport::getMultiDiffCoeffs",
-                           "invert returned ierr = {}", ierr);
-    }
+    invert(m_Lmatrix, m_nsp);
     m_l0000_ok = false; // matrix is overwritten by inverse
     m_lmatrix_soln_ok = false;
 

src/transport/MixTransport.cpp

@@ -86,10 +86,6 @@ void MixTransport::update_T()
     if (t == m_temp && m_nsp == m_thermo->nSpecies()) {
         return;
     }
-    if (t < 0.0) {
-        throw CanteraError("MixTransport::update_T",
-                           "negative temperature {}", t);
-    }
     GasTransport::update_T();
     // temperature has changed, so polynomial fits will need to be redone.
     m_spcond_ok = false;

src/transport/MultiTransport.cpp

@@ -54,7 +54,7 @@ void MultiTransport::init(ThermoPhase* thermo, int mode, int log_level)
     // set flags all false
     m_l0000_ok = false;
     m_lmatrix_soln_ok = false;
-    m_thermal_tlast = 0.0;
+    m_thermal_tlast = Undef;
 
     // some work space
     m_spwork1.resize(m_nsp);
@@ -81,6 +81,15 @@ void MultiTransport::init(ThermoPhase* thermo, int mode, int log_level)
     }
 }
 
+void MultiTransport::invalidateCache()
+{
+    GasTransport::invalidateCache();
+    m_thermal_tlast = Undef;
+    m_l0000_ok = false;
+    m_lmatrix_soln_ok = false;
+    m_molefracs_last[0] += 1.23e-4;
+}
+
 double MultiTransport::thermalConductivity()
 {
     solveLMatrixEquation();
@@ -304,7 +313,7 @@ void MultiTransport::getMassFluxes(const double* state1, const double* state2,
     double gradmax = -1.0;
     for (size_t j = 0; j < m_nsp; j++) {
         if (fabs(x2[j] - x1[j]) > gradmax) {
-            gradmax = fabs(x1[j] - x2[j]);
+            gradmax = fabs(x1[j] - x2[j]) / delta;
             jmax = j;
         }
     }
@@ -313,7 +322,7 @@ void MultiTransport::getMassFluxes(const double* state1, const double* state2,
     // and set the entry in gradx to zero
     for (size_t j = 0; j < m_nsp; j++) {
         m_aa(jmax,j) = y[j];
-        fluxes[j] = x2[j] - x1[j];
+        fluxes[j] = (x2[j] - x1[j]) / delta;
     }
     fluxes[jmax] = 0.0;
 
@@ -329,7 +338,7 @@ void MultiTransport::getMassFluxes(const double* state1, const double* state2,
 
     // thermal diffusion
     if (addThermalDiffusion) {
-        double grad_logt = (t2 - t1)/m_temp;
+        double grad_logt = (t2 - t1) / m_temp / delta;
         for (size_t i = 0; i < m_nsp; i++) {
             fluxes[i] -= m_spwork[i]*grad_logt;
         }
@@ -365,11 +374,7 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, double* const d)
     }
 
     // invert L00,00
-    int ierr = invert(m_Lmatrix, m_nsp);
-    if (ierr != 0) {
-        throw CanteraError("MultiTransport::getMultiDiffCoeffs",
-                           "invert returned ierr = {}", ierr);
-    }
+    invert(m_Lmatrix, m_nsp);
     m_l0000_ok = false; // matrix is overwritten by inverse
     m_lmatrix_soln_ok = false;
 

test/data/ET_test.yaml

@@ -6,9 +6,6 @@ units: {length: cm, quantity: mole, activation-energy: K}
 phases:
 - name: gas
   thermo: ideal-gas
-  elements: [O, N, E]
-  species:
-  - species: [O2, E, O2^-]
   skip-undeclared-third-bodies: true
   kinetics: gas
   reactions:
@@ -69,6 +66,23 @@ species:
     well-depth: 136.5
     polarizability: 1.424
   note: L4/89
+- name: C10H8
+  composition: {C: 10, H: 8}
+  thermo:
+    model: NASA7
+    temperature-ranges: [300.0, 1370.0, 3500.0]
+    data:
+    - [-8.71832426, 0.108564074, -9.40254318e-05, 4.10014902e-08, -7.10574258e-12,
+      1.66434413e+04, 61.5987877]
+    - [15.1184828, 0.0389675576, -1.78248658e-05, 3.92092279e-09, -3.39215689e-13,
+      1.01121562e+04, -60.9041102]
+  transport:
+    model: gas
+    geometry: nonlinear
+    well-depth: 630.4
+    diameter: 6.18
+    polarizability: 16.5
+    rotational-relaxation: 1.0
 
 reactions:
 - equation: E + O2 + O2 => O2^- + O2  # Kossyi (1992) Eq.45