Move non-trivial functions out of header files

include/cantera/oneD/Domain1D.h

@@ -43,21 +43,7 @@ public:
      * @param points  Number of grid points.
      * @param time    (unused)
      */
-    Domain1D(size_t nv=1, size_t points=1,
+    Domain1D(size_t nv=1, size_t points=1, double time=0.0);
-             doublereal time = 0.0) :
-        m_rdt(0.0),
-        m_nv(0),
-        m_time(time),
-        m_container(0),
-        m_index(npos),
-        m_type(0),
-        m_iloc(0),
-        m_jstart(0),
-        m_left(0),
-        m_right(0),
-        m_bw(-1) {
-        resize(nv, points);
-    }
 
     virtual ~Domain1D() {}
 
@@ -134,28 +120,7 @@ public:
      * is virtual so that subclasses can perform other actions required to
      * resize the domain.
      */
-    virtual void resize(size_t nv, size_t np) {
+    virtual void resize(size_t nv, size_t np);
-        // if the number of components is being changed, then a
-        // new grid refiner is required.
-        if (nv != m_nv || !m_refiner) {
-            m_nv = nv;
-            m_refiner.reset(new Refiner(*this));
-        }
-        m_nv = nv;
-        m_td.resize(m_nv, 1);
-        m_name.resize(m_nv,"");
-        m_max.resize(m_nv, 0.0);
-        m_min.resize(m_nv, 0.0);
-        // Default error tolerances for all domains
-        m_rtol_ss.resize(m_nv, 1.0e-4);
-        m_atol_ss.resize(m_nv, 1.0e-9);
-        m_rtol_ts.resize(m_nv, 1.0e-4);
-        m_atol_ts.resize(m_nv, 1.0e-11);
-        m_points = np;
-        m_z.resize(np, 0.0);
-        m_slast.resize(m_nv * m_points, 0.0);
-        locate();
-    }
 
     //! Return a reference to the grid refiner.
     Refiner& refiner() {
@@ -207,13 +172,7 @@ public:
     }
 
     //! Name of the nth component. May be overloaded.
-    virtual std::string componentName(size_t n) const {
+    virtual std::string componentName(size_t n) const;
-        if (m_name[n] != "") {
-            return m_name[n];
-        } else {
-            return fmt::format("component {}", n);
-        }
-    }
 
     void setComponentName(size_t n, const std::string& name) {
         m_name[n] = name;
@@ -229,16 +188,7 @@ public:
     }
 
     //! index of component with name \a name.
-    size_t componentIndex(const std::string& name) const {
+    size_t componentIndex(const std::string& name) const;
-        size_t nc = nComponents();
-        for (size_t n = 0; n < nc; n++) {
-            if (name == componentName(n)) {
-                return n;
-            }
-        }
-        throw CanteraError("Domain1D::componentIndex",
-                           "no component named "+name);
-    }
 
     void setBounds(size_t n, doublereal lower, doublereal upper) {
         m_min[n] = lower;
@@ -445,25 +395,7 @@ public:
      * Find the index of the first grid point in this domain, and
      * the start of its variables in the global solution vector.
      */
-    void locate() {
+    void locate();
-        if (m_left) {
-            // there is a domain on the left, so the first grid point
-            // in this domain is one more than the last one on the left
-            m_jstart = m_left->lastPoint() + 1;
-
-            // the starting location in the solution vector
-            m_iloc = m_left->loc() + m_left->size();
-        } else {
-            // this is the left-most domain
-            m_jstart = 0;
-            m_iloc = 0;
-        }
-        // if there is a domain to the right of this one, then
-        // repeat this for it
-        if (m_right) {
-            m_right->locate();
-        }
-    }
 
     /**
      * Location of the start of the local solution vector in the global
@@ -571,18 +503,7 @@ public:
         return m_z[m_points - 1];
     }
 
-    void setProfile(const std::string& name, doublereal* values, doublereal* soln) {
+    void setProfile(const std::string& name, double* values, double* soln);
-        for (size_t n = 0; n < m_nv; n++) {
-            if (name == componentName(n)) {
-                for (size_t j = 0; j < m_points; j++) {
-                    soln[index(n, j) + m_iloc] = values[j];
-                }
-                return;
-            }
-        }
-        throw CanteraError("Domain1D::setProfile",
-                           "unknown component: "+name);
-    }
 
     vector_fp& grid() {
         return m_z;

include/cantera/oneD/Inlet1D.h

@@ -109,10 +109,7 @@ protected:
 class Inlet1D : public Bdry1D
 {
 public:
-    Inlet1D() : Bdry1D(), m_V0(0.0), m_nsp(0), m_flow(0) {
+    Inlet1D();
-        m_type = cInletType;
-        m_xstr = "";
-    }
 
     /// set spreading rate
     virtual void setSpreadRate(doublereal V0) {
@@ -125,25 +122,9 @@ public:
         return m_V0;
     }
 
-    virtual void showSolution(const doublereal* x) {
+    virtual void showSolution(const double* x);
-        writelog("    Mass Flux:   {:10.4g} kg/m^2/s \n", m_mdot);
-        writelog("    Temperature: {:10.4g} K \n", m_temp);
-        if (m_flow) {
-            writelog("    Mass Fractions: \n");
-            for (size_t k = 0; k < m_flow->phase().nSpecies(); k++) {
-                if (m_yin[k] != 0.0) {
-                    writelog("        {:>16s}  {:10.4g} \n",
-                            m_flow->phase().speciesName(k), m_yin[k]);
-                }
-            }
-        }
-        writelog("\n");
-    }
 
-    virtual void _getInitialSoln(doublereal* x) {
+    virtual void _getInitialSoln(double* x);
-        x[0] = m_mdot;
-        x[1] = m_temp;
-    }
 
     virtual size_t nSpecies() {
         return m_nsp;
@@ -256,10 +237,7 @@ public:
 class OutletRes1D : public Bdry1D
 {
 public:
-    OutletRes1D() : Bdry1D(), m_nsp(0), m_flow(0) {
+    OutletRes1D();
-        m_type = cOutletResType;
-        m_xstr = "";
-    }
 
     virtual void showSolution(const doublereal* x) {}
 
@@ -317,10 +295,7 @@ public:
         x[0] = m_temp;
     }
 
-    virtual void showSolution_s(std::ostream& s, const doublereal* x) {
+    virtual void showSolution_s(std::ostream& s, const double* x);
-        s << "-------------------  Surface " << domainIndex() << " ------------------- " << std::endl;
-        s << "  temperature: " << m_temp << " K" << "    " << x[0] << std::endl;
-    }
 
     virtual void showSolution(const doublereal* x) {
         writelog("    Temperature: {:10.4g} K \n\n", m_temp);
@@ -334,18 +309,9 @@ public:
 class ReactingSurf1D : public Bdry1D
 {
 public:
-    ReactingSurf1D() : Bdry1D(),
+    ReactingSurf1D();
-        m_kin(0), m_surfindex(0), m_nsp(0) {
-        m_type = cSurfType;
-    }
 
-    void setKineticsMgr(InterfaceKinetics* kin) {
+    void setKineticsMgr(InterfaceKinetics* kin);
-        m_kin = kin;
-        m_surfindex = kin->surfacePhaseIndex();
-        m_sphase = (SurfPhase*)&kin->thermo(m_surfindex);
-        m_nsp = m_sphase->nSpecies();
-        m_enabled = true;
-    }
 
     void enableCoverageEquations(bool docov) {
         m_enabled = docov;
@@ -371,14 +337,7 @@ public:
         std::copy(x+1,x+1+m_nsp,m_fixed_cov.begin());
     }
 
-    virtual void showSolution(const doublereal* x) {
+    virtual void showSolution(const doublereal* x);
-        writelog("    Temperature: {:10.4g} K \n", x[0]);
-        writelog("    Coverages: \n");
-        for (size_t k = 0; k < m_nsp; k++) {
-            writelog("    {:>20s} {:10.4g} \n", m_sphase->speciesName(k), x[k+1]);
-        }
-        writelog("\n");
-    }
 
 protected:
     InterfaceKinetics* m_kin;

include/cantera/oneD/OneDim.h

@@ -237,14 +237,7 @@ public:
         return m_nsteps_max;
     }
 
-    void setJacAge(int ss_age, int ts_age=-1) {
+    void setJacAge(int ss_age, int ts_age=-1);
-        m_ss_jac_age = ss_age;
-        if (ts_age > 0) {
-            m_ts_jac_age = ts_age;
-        } else {
-            m_ts_jac_age = m_ss_jac_age;
-        }
-    }
 
     /**
      * Save statistics on function and Jacobian evaluation, and reset the

include/cantera/oneD/StFlow.h

@@ -97,12 +97,7 @@ public:
     }
 
     //! Write the initial solution estimate into array x.
-    virtual void _getInitialSoln(doublereal* x) {
+    virtual void _getInitialSoln(double* x);
-        for (size_t j = 0; j < m_points; j++) {
-            T(x,j) = m_thermo->temperature();
-            m_thermo->getMassFractions(&Y(x, 0, j));
-        }
-    }
 
     virtual void _finalize(const doublereal* x);
 
@@ -154,28 +149,7 @@ public:
         return "<none>";
     }
 
-    void solveEnergyEqn(size_t j=npos) {
+    void solveEnergyEqn(size_t j=npos);
-        bool changed = false;
-        if (j == npos) {
-            for (size_t i = 0; i < m_points; i++) {
-                if (!m_do_energy[i]) {
-                    changed = true;
-                }
-                m_do_energy[i] = true;
-            }
-        } else {
-            if (!m_do_energy[j]) {
-                changed = true;
-            }
-            m_do_energy[j] = true;
-        }
-        m_refiner->setActive(0, true);
-        m_refiner->setActive(1, true);
-        m_refiner->setActive(2, true);
-        if (changed) {
-            needJacUpdate();
-        }
-    }
 
     //! Turn radiation on / off.
     /*!
@@ -199,41 +173,9 @@ public:
      * radiative term and writes them into the variables, which are used for the
      * calculation.
      */
-    void setBoundaryEmissivities(doublereal e_left, doublereal e_right) {
+    void setBoundaryEmissivities(doublereal e_left, doublereal e_right);
-        if (e_left < 0 || e_left > 1) {
-            throw CanteraError("setBoundaryEmissivities",
-                "The left boundary emissivity must be between 0.0 and 1.0!");
-        } else if (e_right < 0 || e_right > 1) {
-            throw CanteraError("setBoundaryEmissivities",
-                "The right boundary emissivity must be between 0.0 and 1.0!");
-        } else {
-            m_epsilon_left = e_left;
-            m_epsilon_right = e_right;
-        }
-    }
 
-    void fixTemperature(size_t j=npos) {
+    void fixTemperature(size_t j=npos);
-        bool changed = false;
-        if (j == npos) {
-            for (size_t i = 0; i < m_points; i++) {
-                if (m_do_energy[i]) {
-                    changed = true;
-                }
-                m_do_energy[i] = false;
-            }
-        } else {
-            if (m_do_energy[j]) {
-                changed = true;
-            }
-            m_do_energy[j] = false;
-        }
-        m_refiner->setActive(0, false);
-        m_refiner->setActive(1, false);
-        m_refiner->setActive(2, false);
-        if (changed) {
-            needJacUpdate();
-        }
-    }
 
     bool doEnergy(size_t j) {
         return m_do_energy[j];

include/cantera/thermo/ThermoPhase.h

@@ -484,13 +484,7 @@ public:
      * @param mu  Output vector of species electrochemical
      *            potentials. Length: m_kk. Units: J/kmol
      */
-    void getElectrochemPotentials(doublereal* mu) const {
+    void getElectrochemPotentials(doublereal* mu) const;
-        getChemPotentials(mu);
-        double ve = Faraday * electricPotential();
-        for (size_t k = 0; k < m_kk; k++) {
-            mu[k] += ve*charge(k);
-        }
-    }
 
     //! Returns an array of partial molar enthalpies for the species
     //! in the mixture. Units (J/kmol)

include/cantera/zeroD/Reactor.h

@@ -52,14 +52,7 @@ public:
         setKineticsMgr(contents);
     }
 
-    void setKineticsMgr(Kinetics& kin) {
+    void setKineticsMgr(Kinetics& kin);
-        m_kin = &kin;
-        if (m_kin->nReactions() == 0) {
-            disableChemistry();
-        } else {
-            enableChemistry();
-        }
-    }
 
     //! Disable changes in reactor composition due to chemical reactions.
     void disableChemistry() {

include/cantera/zeroD/ReactorNet.h

@@ -30,46 +30,21 @@ public:
 
     //! Set initial time. Default = 0.0 s. Restarts integration from this time
     //! using the current mixture state as the initial condition.
-    void setInitialTime(doublereal time) {
+    void setInitialTime(double time);
-        m_time = time;
-        m_integrator_init = false;
-    }
 
     //! Set the maximum time step.
-    void setMaxTimeStep(double maxstep) {
+    void setMaxTimeStep(double maxstep);
-        m_maxstep = maxstep;
-        m_init = false;
-    }
 
     //! Set the maximum number of error test failures permitted by the CVODES
     //! integrator in a single time step.
-    void setMaxErrTestFails(int nmax) {
+    void setMaxErrTestFails(int nmax);
-        m_maxErrTestFails = nmax;
-        m_init = false;
-    }
 
     //! Set the relative and absolute tolerances for the integrator.
-    void setTolerances(doublereal rtol, doublereal atol) {
+    void setTolerances(double rtol, double atol);
-        if (rtol >= 0.0) {
-            m_rtol = rtol;
-        }
-        if (atol >= 0.0) {
-            m_atols = atol;
-        }
-        m_init = false;
-    }
 
     //! Set the relative and absolute tolerances for integrating the
     //! sensitivity equations.
-    void setSensitivityTolerances(doublereal rtol, doublereal atol) {
+    void setSensitivityTolerances(double rtol, double atol);
-        if (rtol >= 0.0) {
-            m_rtolsens = rtol;
-        }
-        if (atol >= 0.0) {
-            m_atolsens = atol;
-        }
-        m_init = false;
-    }
 
     //! Current value of the simulation time.
     doublereal time() {
@@ -153,16 +128,7 @@ public:
      *  rate constant (and implicitly on the reverse rate constant for
      *  reversible reactions).
      */
-    double sensitivity(size_t k, size_t p) {
+    double sensitivity(size_t k, size_t p);
-        if (!m_init) {
-            initialize();
-        }
-        if (p >= m_sensIndex.size()) {
-            throw IndexError("ReactorNet::sensitivity",
-                             "m_sensIndex", p, m_sensIndex.size()-1);
-        }
-        return m_integ->sensitivity(k, m_sensIndex[p])/m_integ->solution(k);
-    }
 
     //! Return the sensitivity of the component named *component* with respect to
     //! the *p*-th sensitivity parameter.

src/oneD/Domain1D.cpp

@@ -11,6 +11,67 @@ using namespace std;
 namespace Cantera
 {
 
+Domain1D::Domain1D(size_t nv, size_t points, double time) :
+    m_rdt(0.0),
+    m_nv(0),
+    m_time(time),
+    m_container(0),
+    m_index(npos),
+    m_type(0),
+    m_iloc(0),
+    m_jstart(0),
+    m_left(0),
+    m_right(0),
+    m_bw(-1)
+{
+    resize(nv, points);
+}
+
+void Domain1D::resize(size_t nv, size_t np)
+{
+    // if the number of components is being changed, then a
+    // new grid refiner is required.
+    if (nv != m_nv || !m_refiner) {
+        m_nv = nv;
+        m_refiner.reset(new Refiner(*this));
+    }
+    m_nv = nv;
+    m_td.resize(m_nv, 1);
+    m_name.resize(m_nv,"");
+    m_max.resize(m_nv, 0.0);
+    m_min.resize(m_nv, 0.0);
+    // Default error tolerances for all domains
+    m_rtol_ss.resize(m_nv, 1.0e-4);
+    m_atol_ss.resize(m_nv, 1.0e-9);
+    m_rtol_ts.resize(m_nv, 1.0e-4);
+    m_atol_ts.resize(m_nv, 1.0e-11);
+    m_points = np;
+    m_z.resize(np, 0.0);
+    m_slast.resize(m_nv * m_points, 0.0);
+    locate();
+}
+
+std::string Domain1D::componentName(size_t n) const
+{
+    if (m_name[n] != "") {
+        return m_name[n];
+    } else {
+        return fmt::format("component {}", n);
+    }
+}
+
+size_t Domain1D::componentIndex(const std::string& name) const
+{
+    size_t nc = nComponents();
+    for (size_t n = 0; n < nc; n++) {
+        if (name == componentName(n)) {
+            return n;
+        }
+    }
+    throw CanteraError("Domain1D::componentIndex",
+                       "no component named "+name);
+}
+
 void Domain1D::setTransientTolerances(doublereal rtol, doublereal atol, size_t n)
 {
     if (n == npos) {
@@ -147,6 +208,26 @@ void Domain1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     }
 }
 
+void Domain1D::locate()
+{
+    if (m_left) {
+        // there is a domain on the left, so the first grid point in this domain
+        // is one more than the last one on the left
+        m_jstart = m_left->lastPoint() + 1;
+
+        // the starting location in the solution vector
+        m_iloc = m_left->loc() + m_left->size();
+    } else {
+        // this is the left-most domain
+        m_jstart = 0;
+        m_iloc = 0;
+    }
+    // if there is a domain to the right of this one, then repeat this for it
+    if (m_right) {
+        m_right->locate();
+    }
+}
+
 void Domain1D::setupGrid(size_t n, const doublereal* z)
 {
     if (n > 1) {
@@ -195,6 +276,19 @@ void Domain1D::showSolution(const doublereal* x)
     writelog("\n");
 }
 
+void Domain1D::setProfile(const std::string& name, double* values, double* soln)
+{
+    for (size_t n = 0; n < m_nv; n++) {
+        if (name == componentName(n)) {
+            for (size_t j = 0; j < m_points; j++) {
+                soln[index(n, j) + m_iloc] = values[j];
+            }
+            return;
+        }
+    }
+    throw CanteraError("Domain1D::setProfile", "unknown component: "+name);
+}
+
 void Domain1D::_getInitialSoln(doublereal* x)
 {
     for (size_t j = 0; j < m_points; j++) {

src/oneD/OneDim.cpp

@@ -101,6 +101,16 @@ MultiNewton& OneDim::newton()
     return *m_newt;
 }
 
+void OneDim::setJacAge(int ss_age, int ts_age)
+{
+    m_ss_jac_age = ss_age;
+    if (ts_age > 0) {
+        m_ts_jac_age = ts_age;
+    } else {
+        m_ts_jac_age = m_ss_jac_age;
+    }
+}
+
 void OneDim::writeStats(int printTime)
 {
     saveStats();

src/oneD/StFlow.cpp

@@ -181,6 +181,14 @@ void StFlow::enableSoret(bool withSoret)
     }
 }
 
+void StFlow::_getInitialSoln(double* x)
+{
+    for (size_t j = 0; j < m_points; j++) {
+        T(x,j) = m_thermo->temperature();
+        m_thermo->getMassFractions(&Y(x, 0, j));
+    }
+}
+
 void StFlow::setGas(const doublereal* x, size_t j)
 {
     m_thermo->setTemperature(T(x,j));
@@ -819,6 +827,68 @@ XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
     return flow;
 }
 
+void StFlow::solveEnergyEqn(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (!m_do_energy[i]) {
+                changed = true;
+            }
+            m_do_energy[i] = true;
+        }
+    } else {
+        if (!m_do_energy[j]) {
+            changed = true;
+        }
+        m_do_energy[j] = true;
+    }
+    m_refiner->setActive(0, true);
+    m_refiner->setActive(1, true);
+    m_refiner->setActive(2, true);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
+void StFlow::setBoundaryEmissivities(doublereal e_left, doublereal e_right)
+{
+    if (e_left < 0 || e_left > 1) {
+        throw CanteraError("setBoundaryEmissivities",
+            "The left boundary emissivity must be between 0.0 and 1.0!");
+    } else if (e_right < 0 || e_right > 1) {
+        throw CanteraError("setBoundaryEmissivities",
+            "The right boundary emissivity must be between 0.0 and 1.0!");
+    } else {
+        m_epsilon_left = e_left;
+        m_epsilon_right = e_right;
+    }
+}
+
+void StFlow::fixTemperature(size_t j)
+{
+    bool changed = false;
+    if (j == npos) {
+        for (size_t i = 0; i < m_points; i++) {
+            if (m_do_energy[i]) {
+                changed = true;
+            }
+            m_do_energy[i] = false;
+        }
+    } else {
+        if (m_do_energy[j]) {
+            changed = true;
+        }
+        m_do_energy[j] = false;
+    }
+    m_refiner->setActive(0, false);
+    m_refiner->setActive(1, false);
+    m_refiner->setActive(2, false);
+    if (changed) {
+        needJacUpdate();
+    }
+}
+
 void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
                                      integer* diag, doublereal rdt)
 {

src/oneD/boundaries1D.cpp

@@ -73,6 +73,37 @@ void Bdry1D::_init(size_t n)
 
 // ---------------- Inlet1D methods ----------------
 
+Inlet1D::Inlet1D()
+    : m_V0(0.0)
+    , m_nsp(0)
+    , m_flow(0)
+{
+    m_type = cInletType;
+    m_xstr = "";
+}
+
+void Inlet1D::showSolution(const double* x)
+{
+    writelog("    Mass Flux:   {:10.4g} kg/m^2/s \n", m_mdot);
+    writelog("    Temperature: {:10.4g} K \n", m_temp);
+    if (m_flow) {
+        writelog("    Mass Fractions: \n");
+        for (size_t k = 0; k < m_flow->phase().nSpecies(); k++) {
+            if (m_yin[k] != 0.0) {
+                writelog("        {:>16s}  {:10.4g} \n",
+                        m_flow->phase().speciesName(k), m_yin[k]);
+            }
+        }
+    }
+    writelog("\n");
+}
+
+void Inlet1D::_getInitialSoln(double* x)
+{
+    x[0] = m_mdot;
+    x[1] = m_temp;
+}
+
 void Inlet1D::setMoleFractions(const std::string& xin)
 {
     m_xstr = xin;
@@ -367,6 +398,14 @@ void Symm1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
 
 // -------- Outlet1D --------
 
+OutletRes1D::OutletRes1D()
+    : m_nsp(0)
+    , m_flow(0)
+{
+    m_type = cOutletResType;
+    m_xstr = "";
+}
+
 string Outlet1D::componentName(size_t n) const
 {
     switch (n) {
@@ -671,8 +710,31 @@ void Surf1D::restore(const XML_Node& dom, doublereal* soln, int loglevel)
     resize(1,1);
 }
 
+void Surf1D::showSolution_s(std::ostream& s, const double* x)
+{
+    s << "-------------------  Surface " << domainIndex() << " ------------------- " << std::endl;
+    s << "  temperature: " << m_temp << " K" << "    " << x[0] << std::endl;
+}
+
 // -------- ReactingSurf1D --------
 
+ReactingSurf1D::ReactingSurf1D()
+    : m_kin(0)
+    , m_surfindex(0)
+    , m_nsp(0)
+{
+    m_type = cSurfType;
+}
+
+void ReactingSurf1D::setKineticsMgr(InterfaceKinetics* kin)
+{
+    m_kin = kin;
+    m_surfindex = kin->surfacePhaseIndex();
+    m_sphase = (SurfPhase*)&kin->thermo(m_surfindex);
+    m_nsp = m_sphase->nSpecies();
+    m_enabled = true;
+}
+
 string ReactingSurf1D::componentName(size_t n) const
 {
     if (n == 0) {
@@ -820,4 +882,14 @@ void ReactingSurf1D::restore(const XML_Node& dom, doublereal* soln,
 
     resize(m_nsp+1,1);
 }
+
+void ReactingSurf1D::showSolution(const double* x)
+{
+    writelog("    Temperature: {:10.4g} K \n", x[0]);
+    writelog("    Coverages: \n");
+    for (size_t k = 0; k < m_nsp; k++) {
+        writelog("    {:>20s} {:10.4g} \n", m_sphase->speciesName(k), x[k+1]);
+    }
+    writelog("\n");
+}
 }

src/thermo/ThermoPhase.cpp

@@ -125,6 +125,15 @@ void ThermoPhase::getLnActivityCoefficients(doublereal* lnac) const
     }
 }
 
+void ThermoPhase::getElectrochemPotentials(doublereal* mu) const
+{
+    getChemPotentials(mu);
+    double ve = Faraday * electricPotential();
+    for (size_t k = 0; k < m_kk; k++) {
+        mu[k] += ve*charge(k);
+    }
+}
+
 void ThermoPhase::setState_TPX(doublereal t, doublereal p, const doublereal* x)
 {
     setMoleFractions(x);

src/zeroD/Reactor.cpp

@@ -25,6 +25,16 @@ Reactor::Reactor() :
     m_nsens(npos)
 {}
 
+void Reactor::setKineticsMgr(Kinetics& kin)
+{
+    m_kin = &kin;
+    if (m_kin->nReactions() == 0) {
+        disableChemistry();
+    } else {
+        enableChemistry();
+    }
+}
+
 void Reactor::getInitialConditions(double t0, size_t leny, double* y)
 {
     warn_deprecated("Reactor::getInitialConditions",

src/zeroD/ReactorNet.cpp

@@ -31,6 +31,46 @@ ReactorNet::~ReactorNet()
     delete m_integ;
 }
 
+void ReactorNet::setInitialTime(double time)
+{
+    m_time = time;
+    m_integrator_init = false;
+}
+
+void ReactorNet::setMaxTimeStep(double maxstep)
+{
+    m_maxstep = maxstep;
+    m_init = false;
+}
+
+void ReactorNet::setMaxErrTestFails(int nmax)
+{
+    m_maxErrTestFails = nmax;
+    m_init = false;
+}
+
+void ReactorNet::setTolerances(double rtol, double atol)
+{
+    if (rtol >= 0.0) {
+        m_rtol = rtol;
+    }
+    if (atol >= 0.0) {
+        m_atols = atol;
+    }
+    m_init = false;
+}
+
+void ReactorNet::setSensitivityTolerances(double rtol, double atol)
+{
+    if (rtol >= 0.0) {
+        m_rtolsens = rtol;
+    }
+    if (atol >= 0.0) {
+        m_atolsens = atol;
+    }
+    m_init = false;
+}
+
 void ReactorNet::initialize()
 {
     size_t n, nv;
@@ -141,6 +181,18 @@ void ReactorNet::eval(doublereal t, doublereal* y,
     checkFinite("ydot", ydot, m_nv);
 }
 
+double ReactorNet::sensitivity(size_t k, size_t p)
+{
+    if (!m_init) {
+        initialize();
+    }
+    if (p >= m_sensIndex.size()) {
+        throw IndexError("ReactorNet::sensitivity",
+                         "m_sensIndex", p, m_sensIndex.size()-1);
+    }
+    return m_integ->sensitivity(k, m_sensIndex[p])/m_integ->solution(k);
+}
+
 void ReactorNet::evalJacobian(doublereal t, doublereal* y,
                               doublereal* ydot, doublereal* p, Array2D* j)
 {