[1D] Move functions from FreeFlame and AxiStagnFlow into StFlow

This makes it possible to implement alternative constitutive relations
(e.g. ionized or non-ideal gases) as a derived class from StFlow and have them
support all of the standard flame configurations (freely propagating, burner
stabilized, counterflow).

include/cantera/oneD/Domain1D.h

@@ -16,6 +16,8 @@ namespace Cantera
 
 // domain types
 const int cFlowType = 50;
+const int cFreeFlow = 51;
+const int cAxisymmetricStagnationFlow = 52;
 const int cConnectorType = 100;
 const int cSurfType = 102;
 const int cInletType = 104;

include/cantera/oneD/StFlow.h

@@ -147,9 +147,14 @@ public:
     virtual void restore(const XML_Node& dom, doublereal* soln,
                          int loglevel);
 
-    // overloaded in subclasses
     virtual std::string flowType() {
-        return "<none>";
+        if (m_type == cFreeFlow) {
+            return "Free Flame";
+        } else if (m_type == cAxisymmetricStagnationFlow) {
+            return "Axisymmetric Stagnation";
+        } else {
+            throw CanteraError("StFlow::flowType", "Unknown value for 'm_type'");
+        }
     }
 
     void solveEnergyEqn(size_t j=npos);
@@ -201,7 +206,7 @@ public:
     }
 
     virtual bool fixed_mdot() {
-        return true;
+        return (domainType() != cFreeFlow);
     }
     void setViscosityFlag(bool dovisc) {
         m_dovisc = dovisc;
@@ -218,13 +223,13 @@ public:
                       integer* mask, doublereal rdt);
 
     //! Evaluate all residual components at the right boundary.
-    virtual void evalRightBoundary(doublereal* x, doublereal* res,
+    virtual void evalRightBoundary(double* x, double* res, int* diag,
-                                   integer* diag, doublereal rdt) = 0;
+                                   double rdt);
 
     //! Evaluate the residual corresponding to the continuity equation at all
     //! interior grid points.
-    virtual void evalContinuity(size_t j, doublereal* x, doublereal* r,
+    virtual void evalContinuity(size_t j, double* x, double* r,
-                                integer* diag, doublereal rdt) = 0;
+                                int* diag, double rdt);
 
     //! Index of the species on the left boundary with the largest mass fraction
     size_t leftExcessSpecies() const {
@@ -432,6 +437,13 @@ protected:
     //! to `j1`, based on solution `x`.
     virtual void updateTransport(doublereal* x, size_t j0, size_t j1);
 
+public:
+    //! Location of the point where temperature is fixed
+    double m_zfixed;
+
+    //! Temperature at the point used to fix the flame location
+    double m_tfixed;
+
 private:
     vector_fp m_ybar;
 };
@@ -446,15 +458,7 @@ public:
     AxiStagnFlow(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1) :
         StFlow(ph, nsp, points) {
         m_dovisc = true;
-    }
+        m_type = cAxisymmetricStagnationFlow;
-
-    virtual void evalRightBoundary(doublereal* x, doublereal* res,
-                                   integer* diag, doublereal rdt);
-    virtual void evalContinuity(size_t j, doublereal* x, doublereal* r,
-                                integer* diag, doublereal rdt);
-
-    virtual std::string flowType() {
-        return "Axisymmetric Stagnation";
     }
 };
 
@@ -465,28 +469,11 @@ public:
 class FreeFlame : public StFlow
 {
 public:
-    FreeFlame(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1);
+    FreeFlame(IdealGasPhase* ph = 0, size_t nsp = 1, size_t points = 1) :
-    virtual void evalRightBoundary(doublereal* x, doublereal* res,
+        StFlow(ph, nsp, points) {
-                                   integer* diag, doublereal rdt);
+        m_dovisc = false;
-    virtual void evalContinuity(size_t j, doublereal* x, doublereal* r,
+        m_type = cFreeFlow;
-                                integer* diag, doublereal rdt);
-
-    virtual std::string flowType() {
-        return "Free Flame";
     }
-    virtual bool fixed_mdot() {
-        return false;
-    }
-    virtual void _finalize(const doublereal* x);
-    virtual void restore(const XML_Node& dom, doublereal* soln, int loglevel);
-
-    virtual XML_Node& save(XML_Node& o, const doublereal* const sol);
-
-    //! Location of the point where temperature is fixed
-    doublereal m_zfixed;
-
-    //! Temperature at the point used to fix the flame location
-    doublereal m_tfixed;
 };
 
 }

src/oneD/StFlow.cpp

@@ -26,7 +26,9 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
     m_do_multicomponent(false),
     m_do_radiation(false),
     m_kExcessLeft(0),
-    m_kExcessRight(0)
+    m_kExcessRight(0),
+    m_zfixed(Undef),
+    m_tfixed(Undef)
 {
     m_type = cFlowType;
     m_points = points;
@@ -85,7 +87,6 @@ StFlow::StFlow(IdealGasPhase* ph, size_t nsp, size_t points) :
         gr.push_back(1.0*ng/m_points);
     }
     setupGrid(m_points, gr.data());
-    setID("stagnation flow");
 
     // Find indices for radiating species
     m_kRadiating.resize(2, npos);
@@ -204,6 +205,29 @@ void StFlow::_finalize(const doublereal* x)
     if (e) {
         solveEnergyEqn();
     }
+
+    if (domainType() == cFreeFlow) {
+        // If the domain contains the temperature fixed point, make sure that it
+        // is correctly set. This may be necessary when the grid has been modified
+        // externally.
+        if (m_tfixed != Undef) {
+            for (size_t j = 0; j < m_points; j++) {
+                if (z(j) == m_zfixed) {
+                    return; // fixed point is already set correctly
+                }
+            }
+
+            for (size_t j = 0; j < m_points - 1; j++) {
+                // Find where the temperature profile crosses the current
+                // fixed temperature.
+                if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) {
+                    m_tfixed = T(x, j+1);
+                    m_zfixed = z(j+1);
+                    return;
+                }
+            }
+        }
+    }
 }
 
 void StFlow::eval(size_t jg, doublereal* xg,
@@ -750,6 +774,11 @@ void StFlow::restore(const XML_Node& dom, doublereal* soln, int loglevel)
                               getFloat(ref, "curve"), getFloat(ref, "prune"));
         refiner().setGridMin(getFloat(ref, "grid_min"));
     }
+
+    if (domainType() == cFreeFlow) {
+        getOptionalFloat(dom, "t_fixed", m_tfixed);
+        getOptionalFloat(dom, "z_fixed", m_zfixed);
+    }
 }
 
 XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
@@ -807,6 +836,10 @@ XML_Node& StFlow::save(XML_Node& o, const doublereal* const sol)
     addFloat(ref, "curve", refiner().maxSlope());
     addFloat(ref, "prune", refiner().prune());
     addFloat(ref, "grid_min", refiner().gridMin());
+    if (m_zfixed != Undef) {
+        addFloat(flow, "z_fixed", m_zfixed, "m");
+        addFloat(flow, "t_fixed", m_tfixed, "K");
+    }
     return flow;
 }
 
@@ -872,62 +905,7 @@ void StFlow::fixTemperature(size_t j)
     }
 }
 
-void AxiStagnFlow::evalRightBoundary(doublereal* x, doublereal* rsd,
+void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
-                                     integer* diag, doublereal rdt)
-{
-    size_t j = m_points - 1;
-    // the boundary object connected to the right of this one may modify or
-    // replace these equations. The default boundary conditions are zero u, V,
-    // and T, and zero diffusive flux for all species.
-    rsd[index(c_offset_U,j)] = rho_u(x,j);
-    rsd[index(c_offset_V,j)] = V(x,j);
-    if (m_do_energy[j]) {
-        rsd[index(c_offset_T,j)] = T(x,j);
-    } else {
-        rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
-    }
-    rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
-    diag[index(c_offset_L, j)] = 0;
-    doublereal sum = 0.0;
-    for (size_t k = 0; k < m_nsp; k++) {
-        sum += Y(x,k,j);
-        rsd[index(k+c_offset_Y,j)] = m_flux(k,j-1) + rho_u(x,j)*Y(x,k,j);
-    }
-    rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
-    diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
-}
-
-void AxiStagnFlow::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
-                                  integer* diag, doublereal rdt)
-{
-    //----------------------------------------------
-    //    Continuity equation
-    //
-    //    Note that this propagates the mass flow rate information to the left
-    //    (j+1 -> j) from the value specified at the right boundary. The
-    //    lambda information propagates in the opposite direction.
-    //
-    //    d(\rho u)/dz + 2\rho V = 0
-    //------------------------------------------------
-    rsd[index(c_offset_U,j)] =
-        -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
-        -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
-
-    //algebraic constraint
-    diag[index(c_offset_U, j)] = 0;
-}
-
-FreeFlame::FreeFlame(IdealGasPhase* ph, size_t nsp, size_t points) :
-    StFlow(ph, nsp, points),
-    m_zfixed(Undef),
-    m_tfixed(Undef)
-{
-    m_dovisc = false;
-    setID("flame");
-}
-
-void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
-                                  integer* diag, doublereal rdt)
 {
     size_t j = m_points - 1;
 
@@ -935,10 +913,7 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
     // replace these equations. The default boundary conditions are zero u, V,
     // and T, and zero diffusive flux for all species.
 
-    // zero gradient
-    rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
     rsd[index(c_offset_V,j)] = V(x,j);
-    rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
     doublereal sum = 0.0;
     rsd[index(c_offset_L, j)] = lambda(x,j) - lambda(x,j-1);
     diag[index(c_offset_L, j)] = 0;
@@ -948,76 +923,53 @@ void FreeFlame::evalRightBoundary(doublereal* x, doublereal* rsd,
     }
     rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
     diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
-}
+    if (domainType() == cAxisymmetricStagnationFlow) {
-
+        rsd[index(c_offset_U,j)] = rho_u(x,j);
-void FreeFlame::evalContinuity(size_t j, doublereal* x, doublereal* rsd,
-                               integer* diag, doublereal rdt)
-{
-    //----------------------------------------------
-    //    Continuity equation
-    //
-    //    d(\rho u)/dz + 2\rho V = 0
-    //----------------------------------------------
-    if (grid(j) > m_zfixed) {
-        rsd[index(c_offset_U,j)] =
-            - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
-            - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
-    } else if (grid(j) == m_zfixed) {
         if (m_do_energy[j]) {
-            rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
+            rsd[index(c_offset_T,j)] = T(x,j);
         } else {
-            rsd[index(c_offset_U,j)] = (rho_u(x,j)
+            rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
-                                        - m_rho[0]*0.3);
         }
-    } else if (grid(j) < m_zfixed) {
+    } else if (domainType() == cFreeFlow) {
-        rsd[index(c_offset_U,j)] =
+        rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
-            - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
+        rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
-            - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
     }
+}
+
+void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double rdt)
+{
     //algebraic constraint
     diag[index(c_offset_U, j)] = 0;
-}
+    //----------------------------------------------
-
+    //    Continuity equation
-void FreeFlame::_finalize(const doublereal* x)
+    //
-{
+    //    d(\rho u)/dz + 2\rho V = 0
-    StFlow::_finalize(x);
+    //----------------------------------------------
-    // If the domain contains the temperature fixed point, make sure that it
+    if (domainType() == cAxisymmetricStagnationFlow) {
-    // is correctly set. This may be necessary when the grid has been modified
+        // Note that this propagates the mass flow rate information to the left
-    // externally.
+        // (j+1 -> j) from the value specified at the right boundary. The
-    if (m_tfixed != Undef) {
+        // lambda information propagates in the opposite direction.
-        for (size_t j = 0; j < m_points; j++) {
+        rsd[index(c_offset_U,j)] =
-            if (z(j) == m_zfixed) {
+            -(rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
-                return; // fixed point is already set correctly
+            -(density(j+1)*V(x,j+1) + density(j)*V(x,j));
-            }
+    } else if (domainType() == cFreeFlow) {
-        }
+        if (grid(j) > m_zfixed) {
-
+            rsd[index(c_offset_U,j)] =
-        for (size_t j = 0; j < m_points - 1; j++) {
+                - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
-            // Find where the temperature profile crosses the current
+                - (density(j-1)*V(x,j-1) + density(j)*V(x,j));
-            // fixed temperature.
+        } else if (grid(j) == m_zfixed) {
-            if ((T(x, j) - m_tfixed) * (T(x, j+1) - m_tfixed) <= 0.0) {
+            if (m_do_energy[j]) {
-                m_tfixed = T(x, j+1);
+                rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
-                m_zfixed = z(j+1);
+            } else {
-                return;
+                rsd[index(c_offset_U,j)] = (rho_u(x,j)
+                                            - m_rho[0]*0.3);
             }
+        } else if (grid(j) < m_zfixed) {
+            rsd[index(c_offset_U,j)] =
+                - (rho_u(x,j+1) - rho_u(x,j))/m_dz[j]
+                - (density(j+1)*V(x,j+1) + density(j)*V(x,j));
         }
     }
 }
 
-void FreeFlame::restore(const XML_Node& dom, doublereal* soln, int loglevel)
-{
-    StFlow::restore(dom, soln, loglevel);
-    getOptionalFloat(dom, "t_fixed", m_tfixed);
-    getOptionalFloat(dom, "z_fixed", m_zfixed);
-}
-
-XML_Node& FreeFlame::save(XML_Node& o, const doublereal* const sol)
-{
-    XML_Node& flow = StFlow::save(o, sol);
-    if (m_zfixed != Undef) {
-        addFloat(flow, "z_fixed", m_zfixed, "m");
-        addFloat(flow, "t_fixed", m_tfixed, "K");
-    }
-    return flow;
-}
-
 } // namespace

src/oneD/boundaries1D.cpp

@@ -42,7 +42,7 @@ void Bdry1D::_init(size_t n)
     // check for left and right flow objects
     if (m_index > 0) {
         Domain1D& r = container().domain(m_index-1);
-        if (r.domainType() == cFlowType) {
+        if (!r.isConnector()) { // flow domain
             m_flow_left = (StFlow*)&r;
             m_left_nv = m_flow_left->nComponents();
             m_left_points = m_flow_left->nPoints();
@@ -59,7 +59,7 @@ void Bdry1D::_init(size_t n)
     // if this is not the last domain, see what is connected on the right
     if (m_index + 1 < container().nDomains()) {
         Domain1D& r = container().domain(m_index+1);
-        if (r.domainType() == cFlowType) {
+        if (!r.isConnector()) { // flow domain
             m_flow_right = (StFlow*)&r;
             m_right_nv = m_flow_right->nComponents();
             m_right_loc = container().start(m_index+1);