Bump version for 2.0.2 maintenance release

Partial cherry-pick of r2181 that doesn't change the interface for existing
derived classes.

README

@@ -3,7 +3,7 @@
 CANTERA
 *******
 
-Version 2.0.1
+Version 2.0.2
 
 License Information
 ===================

SConstruct

@@ -526,7 +526,7 @@ opts.AddVariables(
         name recognized by the 'dot' program. On linux systems, this
         should be lowercase 'helvetica'.""",
      'Helvetica'),
-    ('cantera_version', '', '2.0.1')
+    ('cantera_version', '', '2.0.2')
     )
 
 opts.Update(env)
@@ -726,6 +726,7 @@ int main(int argc, char** argv) {
 
 env = conf.Finish()
 
+env['python_array_include'] = ''
 if env['python_package'] in ('full','default'):
     # Test to see if we can import the specified array module
     warnNoPython = False
@@ -737,7 +738,6 @@ if env['python_package'] in ('full','default'):
             env['python_array_include'] = np.get_include()
         except AttributeError:
             print """WARNING: Couldn't find include directory for Python array package"""
-            env['python_array_include'] = ''
 
         print """INFO: Building the full Python package using %s.""" % env['python_array']
         env['python_package'] = 'full'
@@ -753,7 +753,6 @@ if env['python_package'] in ('full','default'):
             env['python_package'] = 'minimal'
 else:
     warnNoPython = False
-    env['python_array_include'] = ''
 
 
 # Matlab Toolbox settings

cantera.pc.in

@@ -12,7 +12,7 @@ Name: cantera
 Description: An object-oriented software toolkit for chemical kinetics, thermodynamics, and transport processes.
 URL: http://code.google.com/p/cantera/
 Requires:
-Version: 2.0.1
+Version: 2.0.2
 
 LINK  = -lctcxx -luser -loneD -lzeroD -lequil -lkinetics -ltransport -lthermo  -lctnumerics -lctmath -ltpx -lctspectra -lconverters -lctbase -lcvode -lctlapack -lctblas -lctf2c -lm -lstdc++
 

configure.ac

@@ -2,7 +2,7 @@
 # Process this file with autoconf to produce a configure script.
 
 AC_PREREQ(2.61)
-AC_INIT([cantera], [2.0.1], [nick@ices.utexas.edu])
+AC_INIT([cantera], [2.0.2], [nick@ices.utexas.edu])
 AC_CONFIG_HEADER(config.h)
 AC_CONFIG_AUX_DIR([build-aux])
 #AX_ENABLE_BUILDDIR

doc/sphinx/conf.py

@@ -70,7 +70,7 @@ copyright = u'2012, Cantera Developers'
 # The short X.Y version.
 version = '2.0'
 # The full version, including alpha/beta/rc tags.
-release = '2.0.1'
+release = '2.0.2'
 
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.

doc/sphinx/cti/input-files.rst

@@ -368,7 +368,7 @@ If you are interested in seeing the internals of how the preprocessing works,
 take a look at file ``ctml_writer.py`` in the Cantera Python package. Or simply
 start Python, and type::
 
-    >>> from Cantera import ctml_writer
+    >>> import ctml_writer
     >>> help(ctml_writer)
 
 The ``ctml_writer.py`` module can also be run as a script to convert input .cti

include/cantera/kinetics/GasKinetics.h

@@ -347,14 +347,20 @@ public:
         return m_pgroups[i];
     }
 
+    virtual void update_rates_T();
+    virtual void update_rates_C();
 
-    void _update_rates_T();
+    void _update_rates_T() {
+        update_rates_T();
+    }
 
     //! Update properties that depend on concentrations.
     //! Currently the enhanced collision partner concentrations are updated
     //! here, as well as the pressure-dependent portion of P-log and Chebyshev
     //! reactions.
-    void _update_rates_C();
+    void _update_rates_C() {
+        update_rates_C();
+    }
 
     //@}
 

include/cantera/kinetics/RxnRates.h

@@ -49,7 +49,7 @@ public:
         if (m_A  <= 0.0) {
             m_logA = -1.0E300;
         } else {
-            m_logA = log(m_A);
+            m_logA = std::log(m_A);
         }
     }
 
@@ -98,7 +98,7 @@ public:
      * factor.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return m_A * exp(m_b*logT - m_E*recipT);
+        return m_A * std::exp(m_b*logT - m_E*recipT);
     }
 
 
@@ -231,7 +231,7 @@ public:
         if (m_A <= 0.0) {
             m_logA = -1.0E300;
         } else {
-            m_logA = log(m_A);
+            m_logA = std::log(m_A);
         }
 
         const vector_fp& data = rdata.rateCoeffParameters;
@@ -272,7 +272,7 @@ public:
             // changed n to k, dgg 1/22/04
             th = std::max(theta[k], Tiny);
             //                th = fmaxx(theta[n], Tiny);
-            m_mcov += m_mc[n]*log(th);
+            m_mcov += m_mc[n]*std::log(th);
         }
     }
 
@@ -295,7 +295,7 @@ public:
      * factor.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return m_A * exp(m_acov + m_b*logT - (m_E + m_ecov)*recipT + m_mcov);
+        return m_A * std::exp(m_acov + m_b*logT - (m_E + m_ecov)*recipT + m_mcov);
     }
 
     doublereal activationEnergy_R() const {
@@ -394,7 +394,7 @@ public:
         if (m_A  <= 0.0) {
             m_logA = -1.0E300;
         } else {
-            m_logA = log(m_A);
+            m_logA = std::log(m_A);
         }
     }
 
@@ -411,7 +411,7 @@ public:
         if (m_A  <= 0.0) {
             m_logA = -1.0E300;
         } else {
-            m_logA = log(m_A);
+            m_logA = std::log(m_A);
         }
     }
 
@@ -443,7 +443,7 @@ public:
      * factor.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return m_A * exp(m_b*logT - m_E*recipT);
+        return m_A * std::exp(m_b*logT - m_E*recipT);
     }
 
     void writeUpdateRHS(std::ostream& s) const {
@@ -496,7 +496,7 @@ public:
         for (iter_t iter = rdata.plogParameters.begin();
              iter != rdata.plogParameters.end();
              iter++) {
-            double logp = log(iter->first);
+            double logp = std::log(iter->first);
             if (pressures_.empty() || pressures_.rbegin()->first != logp) {
                 // starting a new group
                 pressures_[logp] = std::make_pair(j, j+1);
@@ -520,7 +520,7 @@ public:
              iter != pressures_.end();
              iter++) {
             if (iter->second.first == iter->second.second - 1) {
-                A_[iter->second.first] = log(A_[iter->second.first]);
+                A_[iter->second.first] = std::log(A_[iter->second.first]);
             }
         }
 
@@ -590,9 +590,9 @@ public:
         } else {
             double k = 1e-300; // non-zero to make log(k) finite
             for (size_t m = 0; m < m1_; m++) {
-                k += A1_[m] * exp(n1_[m] * logT - Ea1_[m] * recipT);
+                k += A1_[m] * std::exp(n1_[m] * logT - Ea1_[m] * recipT);
             }
-            log_k1 = log(k);
+            log_k1 = std::log(k);
         }
 
         if (m2_ == 1) {
@@ -600,9 +600,9 @@ public:
         } else {
             double k = 1e-300; // non-zero to make log(k) finite
             for (size_t m = 0; m < m2_; m++) {
-                k += A2_[m] * exp(n2_[m] * logT - Ea2_[m] * recipT);
+                k += A2_[m] * std::exp(n2_[m] * logT - Ea2_[m] * recipT);
             }
-            log_k2 = log(k);
+            log_k2 = std::log(k);
         }
 
         return log_k1 + (log_k2 - log_k1) * (logP_ - logP1_) * rDeltaP_;
@@ -614,7 +614,7 @@ public:
      * This function returns the actual value of the rate constant.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return exp(update(logT, recipT));
+        return std::exp(update(logT, recipT));
     }
 
     doublereal activationEnergy_R() const {
@@ -645,7 +645,7 @@ public:
                     throw CanteraError("Plog::validate",
                         "Invalid rate coefficient for reaction #" +
                         int2str(rdata.number) + ":\n" + rdata.equation + "\n" +
-                        "at P = " + fp2str(exp((++iter)->first)) +
+                        "at P = " + fp2str(std::exp((++iter)->first)) +
                         ", T = " + fp2str(T[i]));
                 }
             }
@@ -697,8 +697,8 @@ public:
         chebCoeffs_(rdata.chebCoeffs),
         dotProd_(rdata.chebDegreeT)
     {
-        double logPmin = log10(rdata.chebPmin);
-        double logPmax = log10(rdata.chebPmax);
+        double logPmin = std::log10(rdata.chebPmin);
+        double logPmax = std::log10(rdata.chebPmax);
         double TminInv = 1.0 / rdata.chebTmin;
         double TmaxInv = 1.0 / rdata.chebTmax;
 
@@ -754,7 +754,7 @@ public:
      * This function returns the actual value of the rate constant.
      */
     doublereal updateRC(doublereal logT, doublereal recipT) const {
-        return pow(10, update(logT, recipT));
+        return std::pow(10, update(logT, recipT));
     }
 
     doublereal activationEnergy_R() const {

include/cantera/oneD/Domain1D.h

@@ -48,6 +48,7 @@ public:
     Domain1D(size_t nv=1, size_t points=1,
              doublereal time = 0.0) :
         m_rdt(0.0),
+        m_nv(0),
         m_time(time),
         m_container(0),
         m_index(npos),

include/cantera/oneD/Sim1D.h

@@ -101,8 +101,16 @@ public:
     void setRefineCriteria(int dom = -1, doublereal ratio = 10.0,
                            doublereal slope = 0.8, doublereal curve = 0.8, doublereal prune = -0.1);
     void setMaxGridPoints(int dom = -1, int npoints = 300);
-
+ 
     void restore(std::string fname, std::string id);
+    //! Set the minimum grid spacing in the specified domain(s).
+    /*!
+     *  @param dom Domain index. If dom == -1, the specified spacing
+                   is applied to all domains.
+        @param gridmin The minimum allowable grid spacing [m]
+    */
+    void setGridMin(int dom, double gridmin);
+
     void getInitialSoln();
 
     void setSolution(const doublereal* soln) {

include/cantera/oneD/StFlow.h

@@ -492,7 +492,6 @@ protected:
     vector_fp m_zfix;
     vector_fp m_tfix;
 
-    doublereal m_efctr;
     bool m_dovisc;
     void updateTransport(doublereal* x, size_t j0, size_t j1);
 

include/cantera/oneD/Surf1D.h

@@ -1,320 +0,0 @@
-dep
-
-#ifndef CT_SURF1D_H
-#define CT_SURF1D_H
-
-#include "Domain1D.h"
-#include "cantera/thermo/SurfPhase.h"
-#include "cantera/kinetics/InterfaceKinetics.h"
-#include "StFlow.h"
-#include "OneDim.h"
-#include "cantera/base/ctml.h"
-
-namespace Cantera
-{
-
-// A class for surface domains in one-dimensional simulations, The
-// surface is zero-dimensional, and defined by a set of surface
-// species coverages.
-
-class Surf1D : public Domain1D
-{
-
-public:
-
-    Surf1D(InterfaceKinetics* skin = 0) : Domain1D(1, 1, 0.0) {
-        m_type = cSurfType;
-        m_flow_left = 0;
-        m_flow_right = 0;
-        m_kin = 0;
-        m_sphase = 0;
-        if (skin) {
-            setKinetics(skin);
-        }
-    }
-    virtual ~Surf1D() {}
-
-    // Set the kinetics manager for the surface.
-    void setKinetics(InterfaceKinetics* kin) {
-        m_kin = kin;
-        int np = kin->nPhases();
-        m_sphase = 0;
-        for (int n = 0; n < np; n++) {
-            if (kin->phase(n).eosType() == cSurf) {
-                m_sphase = (SurfPhase*)&m_kin->phase(n);
-                m_nsurf = n;
-            } else {
-                m_bulk.push_back(&kin->phase(n));
-                m_nbulk.push_back(n);
-            }
-        }
-        if (!m_sphase) {
-            throw CanteraError("setKinetics","no surface phase defined");
-        }
-
-        m_nsp = m_sphase->nSpecies();
-        resize(m_nsp,1);
-        if (m_bulk.size() == 1) {
-            m_bulk.push_back(0);
-        }
-    }
-
-    void fixSpecies(int k, doublereal c) {
-        if (c >= 0.0) {
-            m_fixed_cov[k] = c;
-        }
-        m_do_surf_species[k] = false;
-        needJacUpdate();
-    }
-
-    void solveSpecies(int k) {
-        m_do_surf_species[k] = true;
-        needJacUpdate();
-    }
-
-    /// Set the surface temperature
-    void setTemperature(doublereal t) {
-        m_sphase->setTemperature(t);
-        needJacUpdate();
-    }
-
-    /// Temperature [K].
-    doublereal temperature() {
-        return m_sphase->temperature();
-    }
-
-    void setCoverages(doublereal* c) {
-        m_sphase->setCoverages(c);
-        copy(c, c + m_nsp, m_fixed_cov.begin());
-    }
-
-    void setMultiplier(int k, doublereal f) {
-        m_mult[k] = f;
-        needJacUpdate();
-    }
-
-    doublereal multiplier(int k) {
-        return m_mult[k];
-    }
-
-    virtual std::string componentName(int n) const {
-        return m_sphase->speciesName(n);
-    }
-
-    virtual void init() {
-        if (m_index < 0) {
-            throw CanteraError("Surf1D",
-                               "install in container before calling init.");
-        }
-        m_nsp = m_sphase->nSpecies();
-        resize(m_nsp,1);
-        m_mult.resize(m_nsp, 1.0);
-        m_do_surf_species.resize(m_nsp, true);
-        m_fixed_cov.resize(m_nsp, 1.0/m_nsp);
-
-        // set bounds
-        vector_fp lower(m_nsp, -1.e-3);
-        vector_fp upper(m_nsp, 1.0);
-        setBounds(m_nsp, lower.begin(), m_nsp, upper.begin());
-
-        // set tolerances
-        vector_fp rtol(m_nsp, 1e-4);
-        vector_fp atol(m_nsp, 1.e-10);
-        setTolerances(m_nsp, rtol.begin(), m_nsp, atol.begin());
-
-        m_left_nsp = 0;
-        m_right_nsp = 0;
-
-        // check for left and right flow objects
-        if (m_index > 0) {
-            Domain1D& r = container().domain(m_index-1);
-            if (r.domainType() == cFlowType) {
-                m_flow_left = (StFlow*)&r;
-                m_left_nv = m_flow_left->nComponents();
-                m_left_points = m_flow_left->nPoints();
-                m_left_loc = container().start(m_index-1);
-                m_left_nsp = m_left_nv - 4;
-                m_phase_left = &m_flow_left->phase();
-                m_molwt_left = m_phase_left->molecularWeights().begin();
-                if (m_phase_left == m_bulk[0]) {
-                    m_start_left = m_kin->start(m_nbulk[0]);
-                } else if (m_phase_left == m_bulk[1]) {
-                    m_start_left = m_kin->start(m_nbulk[1]);
-                } else
-                    throw CanteraError("Surf1D::init",
-                                       "left gas does not match one in surface mechanism");
-            } else
-                throw CanteraError("Surf1D::init",
-                                   "Surface domains can only be "
-                                   "connected to flow domains.");
-        }
-
-        if (m_index < container().nDomains() - 1) {
-            Domain1D& r = container().domain(m_index+1);
-            if (r.domainType() == cFlowType) {
-                m_flow_right = (StFlow*)&r;
-                m_right_nv = m_flow_right->nComponents();
-                m_right_loc = container().start(m_index+1);
-                m_right_nsp = m_right_nv - 4;
-                m_phase_right = &m_flow_right->phase();
-                m_molwt_right = m_phase_right->molecularWeights().begin();
-                if (m_phase_right == m_bulk[0]) {
-                    m_start_right = m_kin->start(m_nbulk[0]);
-                } else if (m_phase_right == m_bulk[1]) {
-                    m_start_right = m_kin->start(m_nbulk[1]);
-                } else
-                    throw CanteraError("Surf1D::init",
-                                       "right gas does not match one in surface mechanism");
-            } else
-                throw CanteraError("Surf1D::init",
-                                   "Surface domains can only be "
-                                   "connected to flow domains.");
-        }
-        m_work.resize(m_kin->nSpecies());
-    }
-
-
-    virtual void eval(int jg, doublereal* xg, doublereal* rg,
-                      integer* diagg, doublereal rdt) {
-        int k;
-
-        if (jg >= 0 && (jg < firstPoint() - 2
-                        || jg > lastPoint() + 2)) {
-            return;
-        }
-
-        // start of local part of global arrays
-        doublereal* x = xg + loc();
-        doublereal* r = rg + loc();
-        integer* diag = diagg + loc();
-
-        // set the coverages
-        doublereal sum = 0.0;
-        for (k = 0; k < m_nsp; k++) {
-            m_work[k] = x[k];
-            sum += x[k];
-        }
-        m_sphase->setCoverages(m_work.begin());
-
-        // set the left gas state to the adjacent point
-
-        int leftloc = 0, rightloc = 0;
-        int pnt = 0;
-
-        if (m_flow_left) {
-            leftloc = m_flow_left->loc();
-            pnt = m_flow_left->nPoints() - 1;
-            m_flow_left->setGas(xg + leftloc, pnt);
-        }
-
-        if (m_flow_right) {
-            rightloc = m_flow_right->loc();
-            m_flow_right->setGas(xg + rightloc, 0);
-        }
-
-        m_kin->getNetProductionRates(m_work.begin());
-        doublereal rs0 = 1.0/m_sphase->siteDensity();
-
-        scale(m_work.begin(), m_work.end(), m_work.begin(), m_mult[0]);
-
-        bool enabled = true;
-        int ioffset = m_kin->start(m_nsurf); // m_left_nsp + m_right_nsp;
-        doublereal maxx = -1.0;
-        int imx = -1;
-        for (k = 0; k < m_nsp; k++) {
-            r[k] = m_work[k + ioffset] * m_sphase->size(k) * rs0;
-            r[k] -= rdt*(x[k] - prevSoln(k,0));
-            diag[k] = 1;
-            if (x[k] > maxx) {
-                maxx = x[k];
-                imx = k;
-            }
-            if (!m_do_surf_species[k]) {
-                r[k] = x[k] - m_fixed_cov[k];
-                diag[k] = 0;
-                enabled = false;
-            }
-        }
-        if (enabled) {
-            r[imx] = 1.0 - sum;
-            diag[imx] = 0;
-        }
-
-        // gas-phase residuals
-        doublereal rho;
-        if (m_flow_left) {
-            rho = m_phase_left->density();
-            doublereal rdz = 2.0/
-                             (m_flow_left->z(m_left_points-1) -
-                              m_flow_left->z(m_left_points - 2));
-
-            for (k = 0; k < m_left_nsp; k++) {
-                m_work[k + m_start_left] *= m_molwt_left[k];
-            }
-
-            int ileft = loc() - m_left_nv;
-
-            // if the energy equation is enabled at this point,
-            // set the gas temperature to the surface temperature
-            if (m_flow_left->doEnergy(pnt)) {
-                rg[ileft + 2] = xg[ileft + 2] - m_sphase->temperature();
-            }
-
-            for (k = 1; k < m_left_nsp; k++) {
-                if (enabled && m_flow_left->doSpecies(k)) {
-                    rg[ileft + 4 + k]  += m_work[k + m_start_left];
-                    //+= rdz*m_work[k + m_sp_left]/rho;
-
-                }
-            }
-        }
-
-        if (m_flow_right) {
-            for (k = 0; k < m_right_nsp; k++) {
-                m_work[k + m_start_right] *= m_molwt_right[k];
-            }
-
-            int iright = loc() + m_nsp;
-            rg[iright + 2] -= m_sphase->temperature();
-            //r[iright + 3] = x[iright];
-            for (k = 0; k < m_right_nsp; k++) {
-                rg[iright + 4 + k] -= m_work[k + m_start_right];
-            }
-        }
-    }
-
-    virtual void save(XML_Node& o, const doublereal* const soln) {
-        doublereal* s = soln + loc();
-        XML_Node& surf = o.addChild("surface");
-        for (int k = 0; k < m_nsp; k++) {
-            ctml::addFloat(surf, componentName(k), s[k], "", "coverage",
-                           0.0, 1.0);
-        }
-    }
-
-
-protected:
-
-    InterfaceKinetics* m_kin;
-    SurfPhase* m_sphase;
-    StFlow* m_flow_left, *m_flow_right;
-    int m_left_nv, m_right_nv;
-    int m_left_loc, m_right_loc;
-    int m_left_points;
-    int m_nsp, m_left_nsp, m_right_nsp;
-    vector_fp m_work;
-    const doublereal* m_molwt_right, *m_molwt_left;
-    int m_sp_left, m_sp_right;
-    int m_start_left, m_start_right, m_start_surf;
-    ThermoPhase* m_phase_left, *m_phase_right;
-    std::vector<ThermoPhase*> m_bulk;
-    std::vector<int> m_nbulk;
-    int m_nsurf;
-    vector_fp m_mult;
-    std::vector<bool> m_do_surf_species;
-    vector_fp m_fixed_cov;
-};
-
-}
-
-#endif

include/cantera/oneD/refine.h

@@ -29,6 +29,18 @@ public:
     void setMaxPoints(int npmax) {
         m_npmax = npmax;
     }
+
+    //! Set the minimum allowable spacing between adjacent grid points [m].
+    void setGridMin(double gridmin) {
+        m_gridmin = gridmin;
+    }
+
+    //! Returns the the minimum allowable spacing between adjacent
+    //! grid points [m].
+    double gridMin() const {
+        return m_gridmin;
+    }
+
     int analyze(size_t n, const doublereal* z, const doublereal* x);
     int getNewGrid(int n, const doublereal* z, int nn, doublereal* znew);
     //int getNewSoln(int n, const doublereal* x, doublereal* xnew);
@@ -67,6 +79,7 @@ protected:
     Domain1D*            m_domain;
     size_t m_nv, m_npmax;
     doublereal          m_thresh;
+    doublereal m_gridmin; //!< minimum grid spacing [m]
 
 };
 

include/cantera/onedim.h

@@ -5,7 +5,6 @@
 #include "oneD/OneDim.h"
 #include "oneD/Domain1D.h"
 #include "oneD/Inlet1D.h"
-#include "oneD/Surf1D.h"
 #include "oneD/MultiNewton.h"
 #include "oneD/MultiJac.h"
 #include "oneD/StFlow.h"

interfaces/python/Cantera/Reactor.py

@@ -972,7 +972,7 @@ class Wall:
         """
         if qfunc:
             self._qfunc = qfunc # hold on to a reference so it doesn't get deleted
-            n = self.qfunc.func_id()
+            n = qfunc.func_id()
         else:
             n = 0
         return _cantera.wall_setHeatFlux(self.__wall_id, n)

interfaces/python/ck2cti.py

@@ -982,6 +982,7 @@ def readKineticsEntry(entry, speciesDict, energyUnits, moleculeUnits):
         pdepArrhenius = None
         efficiencies = {}
         chebyshevCoeffs = []
+        revReaction = None
 
         # Note that the subsequent lines could be in any order
         for line in lines[1:]:
@@ -1097,6 +1098,9 @@ def readKineticsEntry(entry, speciesDict, energyUnits, moleculeUnits):
                 for collider, efficiency in zip(tokens[0::2], tokens[1::2]):
                     efficiencies[collider.strip()] = float(efficiency.strip())
 
+        if revReaction:
+            revReaction.duplicate = reaction.duplicate
+
         # Decide which kinetics to keep and store them on the reaction object
         # Only one of these should be true at a time!
         if chebyshev is not None:

site_scons/wxsgen.py

@@ -71,7 +71,7 @@ class WxsGenerator(object):
                                      UpgradeCode='2340BEE1-279D-11E1-A4AA-001FBC085391',
                                      Language='1033',
                                      Codepage='1252',
-                                     Version='2.0.1',
+                                     Version='2.0.2',
                                      Manufacturer='Cantera Developers'))
 
         fields = {'Platform': 'x64'} if self.x64 else {}

src/clib/ct.cpp

@@ -302,7 +302,7 @@ extern "C" {
     {
         try {
             ThermoPhase& p = ThermoCabinet::item(n);
-            p.checkElementArraySize(lenm);
+            p.checkSpeciesArraySize(lenm);
             const vector_fp& wt = p.molecularWeights();
             copy(wt.begin(), wt.end(), mw);
             return 0;

src/kinetics/AqueousKinetics.cpp

@@ -594,6 +594,16 @@ void AqueousKinetics::finalize()
 {
     if (!m_finalized) {
         m_finalized = true;
+
+        // Guarantee that these arrays can be converted to double* even in the
+        // special case where there are no reactions defined.
+        if (!m_ii) {
+            m_perturb.resize(1, 1.0);
+            m_ropf.resize(1, 0.0);
+            m_ropr.resize(1, 0.0);
+            m_ropnet.resize(1, 0.0);
+            m_rkcn.resize(1, 0.0);
+        }
     }
 }
 

src/kinetics/GasKinetics.cpp

@@ -156,7 +156,7 @@ void GasKinetics::
 update_C() {}
 //====================================================================================================================
 void GasKinetics::
-_update_rates_T()
+update_rates_T()
 {
     doublereal T = thermo().temperature();
     m_logStandConc = log(thermo().standardConcentration());
@@ -188,7 +188,7 @@ _update_rates_T()
 //====================================================================================================================
 
 void GasKinetics::
-_update_rates_C()
+update_rates_C()
 {
     thermo().getActivityConcentrations(&m_conc[0]);
     doublereal ctot = thermo().molarDensity();
@@ -232,7 +232,8 @@ void GasKinetics::updateKc()
     doublereal rrt = 1.0/(GasConstant * thermo().temperature());
     for (size_t i = 0; i < m_nrev; i++) {
         size_t irxn = m_revindex[i];
-        m_rkcn[irxn] = exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc);
+        m_rkcn[irxn] = std::min(exp(m_rkcn[irxn]*rrt - m_dn[irxn]*m_logStandConc),
+                                BigNumber);
     }
 
     for (size_t i = 0; i != m_nirrev; ++i) {
@@ -246,7 +247,7 @@ void GasKinetics::updateKc()
  */
 void GasKinetics::getEquilibriumConstants(doublereal* kc)
 {
-    _update_rates_T();
+    update_rates_T();
     thermo().getStandardChemPotentials(&m_grt[0]);
     fill(m_rkcn.begin(), m_rkcn.end(), 0.0);
 
@@ -500,8 +501,8 @@ void GasKinetics::processFalloffReactions()
 //====================================================================================================================
 void GasKinetics::updateROP()
 {
-    _update_rates_C();
-    _update_rates_T();
+    update_rates_C();
+    update_rates_T();
 
     if (m_ROP_ok) {
         return;
@@ -557,8 +558,8 @@ void GasKinetics::updateROP()
 void GasKinetics::
 getFwdRateConstants(doublereal* kfwd)
 {
-    _update_rates_C();
-    _update_rates_T();
+    update_rates_C();
+    update_rates_T();
 
     // copy rate coefficients into ropf
     copy(m_rfn.begin(), m_rfn.end(), m_ropf.begin());
@@ -838,6 +839,16 @@ void GasKinetics::finalize()
         concm_3b_values.resize(m_3b_concm.workSize());
         concm_falloff_values.resize(m_falloff_concm.workSize());
         m_finalized = true;
+
+        // Guarantee that these arrays can be converted to double* even in the
+        // special case where there are no reactions defined.
+        if (!m_ii) {
+            m_perturb.resize(1, 1.0);
+            m_ropf.resize(1, 0.0);
+            m_ropr.resize(1, 0.0);
+            m_ropnet.resize(1, 0.0);
+            m_rkcn.resize(1, 0.0);
+        }
     }
 }
 //====================================================================================================================

src/kinetics/InterfaceKinetics.cpp

@@ -403,10 +403,10 @@ void InterfaceKinetics::updateKc()
 void InterfaceKinetics::checkPartialEquil()
 {
     vector_fp dmu(nTotalSpecies(), 0.0);
-    vector_fp rmu(nReactions(), 0.0);
-    vector_fp frop(nReactions(), 0.0);
-    vector_fp rrop(nReactions(), 0.0);
-    vector_fp netrop(nReactions(), 0.0);
+    vector_fp frop(std::max<size_t>(nReactions(), 1), 0.0);
+    vector_fp rrop(std::max<size_t>(nReactions(), 1), 0.0);
+    vector_fp netrop(std::max<size_t>(nReactions(), 1), 0.0);
+    vector_fp rmu(std::max<size_t>(nReactions(), 1), 0.0);
     if (m_nrev > 0) {
         doublereal rt = GasConstant*thermo(0).temperature();
         cout << "T = " << thermo(0).temperature() << " " << rt << endl;
@@ -1332,7 +1332,8 @@ void InterfaceKinetics::init()
 void InterfaceKinetics::finalize()
 {
     Kinetics::finalize();
-    m_rwork.resize(nReactions());
+    size_t safe_reaction_size = std::max<size_t>(nReactions(), 1);
+    m_rwork.resize(safe_reaction_size);
     size_t ks = reactionPhaseIndex();
     if (ks == npos) throw CanteraError("InterfaceKinetics::finalize",
                                            "no surface phase is present.");
@@ -1343,13 +1344,19 @@ void InterfaceKinetics::finalize()
                            +int2str(m_surf->nDim()));
 
     m_StandardConc.resize(m_kk, 0.0);
-    m_deltaG0.resize(m_ii, 0.0);
-    m_ProdStanConcReac.resize(m_ii, 0.0);
+    m_deltaG0.resize(safe_reaction_size, 0.0);
+    m_ProdStanConcReac.resize(safe_reaction_size, 0.0);
 
     if (m_thermo.size() != m_phaseExists.size()) {
         throw CanteraError("InterfaceKinetics::finalize", "internal error");
     }
 
+    // Guarantee that these arrays can be converted to double* even in the
+    // special case where there are no reactions defined.
+    if (!m_ii) {
+        m_perturb.resize(1, 1.0);
+    }
+
     m_finalized = true;
 }
 
@@ -1485,7 +1492,7 @@ void InterfaceKinetics::setPhaseStability(const int iphase, const int isStable)
 //================================================================================================
 void EdgeKinetics::finalize()
 {
-    m_rwork.resize(nReactions());
+    m_rwork.resize(std::max<size_t>(nReactions(), 1));
     size_t ks = reactionPhaseIndex();
     if (ks == npos) throw CanteraError("EdgeKinetics::finalize",
                                            "no edge phase is present.");
@@ -1494,6 +1501,13 @@ void EdgeKinetics::finalize()
         throw CanteraError("EdgeKinetics::finalize",
                            "expected interface dimension = 1, but got dimension = "
                            +int2str(m_surf->nDim()));
+
+    // Guarantee that these arrays can be converted to double* even in the
+    // special case where there are no reactions defined.
+    if (!m_ii) {
+        m_perturb.resize(1, 1.0);
+    }
+
     m_finalized = true;
 }
 //================================================================================================

src/kinetics/importKinetics.cpp

@@ -902,6 +902,7 @@ bool installReactionArrays(const XML_Node& p, Kinetics& kin,
     p.getChildren("reactionArray",rarrays);
     int na = static_cast<int>(rarrays.size());
     if (na == 0) {
+        kin.finalize();
         return false;
     }
     for (int n = 0; n < na; n++) {

src/matlab/kineticsmethods.cpp

@@ -24,10 +24,10 @@ void kineticsmethods(int nlhs, mxArray* plhs[],
         int in2 = getInt(prhs[5]);
         int in3 = getInt(prhs[6]);
         int in4 = getInt(prhs[7]);
-        vv = (double) newKineticsFromXML(root, iph, in1, in2, in3, in4);
+        int p = (int) newKineticsFromXML(root, iph, in1, in2, in3, in4);
         plhs[0] = mxCreateNumericMatrix(1,1,mxDOUBLE_CLASS,mxREAL);
         double* h = mxGetPr(plhs[0]);
-        *h = vv;
+        *h = (double) p;
         return;
     }
 

src/oneD/Domain1D.cpp

@@ -102,7 +102,7 @@ eval(size_t jg, doublereal* xg, doublereal* rg,
         jmin = 0;
         jmax = m_points - 1;
     } else {          // evaluate points for Jacobian
-        jmin = std::max<size_t>(jpt-1, 0);
+        jmin = std::max<size_t>(jpt, 1) - 1;
         jmax = std::min(jpt+1,m_points-1);
     }
 

src/oneD/MultiNewton.cpp

@@ -312,6 +312,7 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
     bool frst = true;
     doublereal rdt = r.rdt();
     int j0 = jac.nEvals();
+    int nJacReeval = 0;
 
     while (1 > 0) {
 
@@ -370,6 +371,10 @@ int MultiNewton::solve(doublereal* x0, doublereal* x1,
         else if (m < 0) {
             if (jac.age() > 1) {
                 forceNewJac = true;
+                if (nJacReeval > 3) {
+                    goto done;
+                }
+                nJacReeval++;
                 if (loglevel > 0)
                     writelog("\nRe-evaluating Jacobian, since no damping "
                              "coefficient\ncould be found with this Jacobian.\n");

src/oneD/Sim1D.cpp

@@ -578,6 +578,20 @@ void Sim1D::setRefineCriteria(int dom, doublereal ratio,
     }
 }
 
+void Sim1D::setGridMin(int dom, double gridmin)
+{
+    if (dom >= 0) {
+        Refiner& r = domain(dom).refiner();
+        r.setGridMin(gridmin);
+    } else {
+        for (size_t n = 0; n < m_nd; n++) {
+            Refiner& r = domain(n).refiner();
+            r.setGridMin(gridmin);
+        }
+    }
+}
+
+
 void Sim1D::setMaxGridPoints(int dom, int npoints)
 {
     if (dom >= 0) {

src/oneD/StFlow.cpp

@@ -102,8 +102,7 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     m_jac(0),
     m_ok(false),
     m_do_soret(false),
-    m_transport_option(-1),
-    m_efctr(0.0)
+    m_transport_option(-1)
 {
     m_type = cFlowType;
 
@@ -565,9 +564,6 @@ void AxiStagnFlow::eval(size_t jg, doublereal* xg,
                     - divHeatFlux(x,j) - sum - sum2;
                 rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
 
-                rsd[index(c_offset_T, j)] =
-                    rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j));
-
                 rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
                 diag[index(c_offset_T, j)] = 1;
             }
@@ -869,9 +865,6 @@ void FreeFlame::eval(size_t jg, doublereal* xg,
                     - divHeatFlux(x,j) - sum - sum2;
                 rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
 
-                rsd[index(c_offset_T, j)] =
-                    rsd[index(c_offset_T, j)] + m_efctr*(T_fixed(j) - T(x,j));
-
                 rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
                 diag[index(c_offset_T, j)] = 1;
             }

src/oneD/refine.cpp

@@ -40,7 +40,8 @@ static doublereal eps()
 
 Refiner::Refiner(Domain1D& domain) :
     m_ratio(10.0), m_slope(0.8), m_curve(0.8), m_prune(-0.001),
-    m_min_range(0.01), m_domain(&domain), m_npmax(3000)
+    m_min_range(0.01), m_domain(&domain), m_npmax(3000),
+    m_gridmin(5e-6)
 {
     m_nv = m_domain->nComponents();
     m_active.resize(m_nv, true);
@@ -129,7 +130,7 @@ int Refiner::analyze(size_t n, const doublereal* z,
                 dmax = m_slope*(vmax - vmin) + m_thresh;
                 for (j = 0; j < n-1; j++) {
                     r = fabs(v[j+1] - v[j])/dmax;
-                    if (r > 1.0) {
+                    if (r > 1.0 && dz[j] >= 2 * m_gridmin) {
                         m_loc[j] = 1;
                         m_c[name] = 1;
                         //if (int(m_loc.size()) + n > m_npmax) goto done;
@@ -161,7 +162,8 @@ int Refiner::analyze(size_t n, const doublereal* z,
                 dmax = m_curve*(smax - smin); // + 0.5*m_curve*(smax + smin);
                 for (j = 0; j < n-2; j++) {
                     r = fabs(s[j+1] - s[j]) / (dmax + m_thresh/dz[j]);
-                    if (r > 1.0) {
+                    if (r > 1.0 && dz[j] >= 2 * m_gridmin &&
+                        dz[j+1] >= 2 * m_gridmin) {
                         m_c[name] = 1;
                         m_loc[j] = 1;
                         m_loc[j+1] = 1;

src/transport/MultiTransport.cpp

@@ -579,7 +579,7 @@ void MultiTransport::getMultiDiffCoeffs(const size_t ld, doublereal* const d)
 
     // update the binary diffusion coefficients
     update_T();
-    updateDiff_T();
+    updateThermal_T();
 
     // evaluate L0000 if the temperature or concentrations have
     // changed since it was last evaluated.