[1D] Split StFlow evaluation into separate functions for each equation

- pull radiation computation into a separate method
- Various comment formatting/cleaning up
- eliminate evalResidual and evalRightBoundary methods
- moved wdot eval to update_thermo() method
- removed setGas() call as it is redundant

include/cantera/oneD/IonFlow.h

@@ -68,11 +68,17 @@ public:
 
 protected:
     /*!
-     * This function overloads the original function. The residual function
+     * This function overloads the original electric field function.
-     * of electric field is added.
+     * The residual function of a non-zero electric field is added.
      */
-    void evalResidual(double* x, double* rsd, int* diag,
+    void evalElectricField(double* x, double* rsd, int* diag,
-                      double rdt, size_t jmin, size_t jmax) override;
+                           double rdt, size_t jmin, size_t jmax) override;
+    /*!
+     * This function overloads the original species function.
+     * A Neumann boundary for the charged species at the left boundary is added.
+     */
+    void evalSpecies(double* x, double* rsd, int* diag,
+                     double rdt, size_t jmin, size_t jmax) override;
     void updateTransport(double* x, size_t j0, size_t j1) override;
     void updateDiffFluxes(const double* x, size_t j0, size_t j1) override;
     //! Solving phase one: the fluxes of charged species are turned off

include/cantera/oneD/StFlow.h

@@ -26,7 +26,7 @@ enum offset
     , c_offset_V //! strain rate
     , c_offset_T //! temperature
     , c_offset_L //! (1/r)dP/dr
-    , c_offset_E //! electric poisson's equation
+    , c_offset_E //! electric field equation
     , c_offset_Y //! mass fractions
 };
 
@@ -232,6 +232,10 @@ public:
         return m_qdotRadiation[j];
     }
 
+    //! Computes the radiative heat loss vector over points jmin to jmax and stores
+    //! the data in the qdotRadiation variable.
+    void computeRadiation(double* x, size_t jmin, size_t jmax);
+
     //! Set the emissivities for the boundary values
     /*!
      * Reads the emissivities for the left and right boundary values in the
@@ -300,13 +304,6 @@ public:
      */
     void eval(size_t j, double* x, double* r, integer* mask, double rdt) override;
 
-    //! Evaluate all residual components at the right boundary.
-    virtual void evalRightBoundary(double* x, double* res, int* diag, double rdt);
-
-    //! Evaluate the residual corresponding to the continuity equation at all
-    //! interior grid points.
-    virtual void evalContinuity(size_t j, double* x, double* r, int* diag, double rdt);
-
     //! Index of the species on the left boundary with the largest mass fraction
     size_t leftExcessSpecies() const {
         return m_kExcessLeft;
@@ -336,11 +333,30 @@ protected:
     //! to be updated are defined.
     virtual void updateProperties(size_t jg, double* x, size_t jmin, size_t jmax);
 
-    //! Evaluate the residual function. This function is called in eval
+    //! Evaluate the residual function for the continuity equation.
-    //! after updateProperties is called.
+    virtual void evalContinuity(double* x, double* rsd, int* diag,
-    virtual void evalResidual(double* x, double* rsd, int* diag,
+                                double rdt, size_t jmin, size_t jmax);
+
+    //! Evaluate the residual function for the momentum equation.
+    virtual void evalMomentum(double* x, double* rsd, int* diag,
                               double rdt, size_t jmin, size_t jmax);
 
+    //! Evaluate the residual function for the energy equation.
+    virtual void evalEnergy(double* x, double* rsd, int* diag,
+                            double rdt, size_t jmin, size_t jmax);
+
+    //! Evaluate the residual function for the species equation.
+    virtual void evalSpecies(double* x, double* rsd, int* diag,
+                             double rdt, size_t jmin, size_t jmax);
+
+    //! Evaluate the residual function for the lambda equation.
+    virtual void evalLambda(double* x, double* rsd, int* diag,
+                            double rdt, size_t jmin, size_t jmax);
+
+    //! Evaluate the residual function for the lambda equation.
+    virtual void evalElectricField(double* x, double* rsd, int* diag,
+                                   double rdt, size_t jmin, size_t jmax);
+
     /**
      * Update the thermodynamic properties from point j0 to point j1
      * (inclusive), based on solution x.
@@ -352,6 +368,7 @@ protected:
             m_wtm[j] = m_thermo->meanMolecularWeight();
             m_cp[j] = m_thermo->cp_mass();
             m_thermo->getPartialMolarEnthalpies(&m_hk(0, j));
+            m_kin->getNetProductionRates(&m_wdot(0, j));
         }
     }
 

src/oneD/IonFlow.cpp

@@ -199,41 +199,61 @@ void IonFlow::setSolvingStage(const size_t stage)
     }
 }
 
-void IonFlow::evalResidual(double* x, double* rsd, int* diag,
+void IonFlow::evalElectricField(double* x, double* rsd, int* diag,
-                           double rdt, size_t jmin, size_t jmax)
+                                double rdt, size_t jmin, size_t jmax)
 {
-    StFlow::evalResidual(x, rsd, diag, rdt, jmin, jmax);
+    //-----------------------------------------------
+    //    Electric field by Gauss's law
+    //
+    //    dE/dz = e/eps_0 * sum(q_k*n_k)
+    //
+    //    E = -dV/dz
+    //-----------------------------------------------
+    StFlow::evalElectricField(x, rsd, diag, rdt, jmin, jmax);
     if (m_stage != 2) {
         return;
     }
 
     for (size_t j = jmin; j <= jmax; j++) {
         if (j == 0) {
-            // enforcing the flux for charged species is difficult
-            // since charged species are also affected by electric
-            // force, so Neumann boundary condition is used.
-            for (size_t k : m_kCharge) {
-                rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1);
-            }
             rsd[index(c_offset_E, j)] = E(x,0);
             diag[index(c_offset_E, j)] = 0;
         } else if (j == m_points - 1) {
             rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
             diag[index(c_offset_E, j)] = 0;
         } else {
-            //-----------------------------------------------
-            //    Electric field by Gauss's law
-            //
-            //    dE/dz = e/eps_0 * sum(q_k*n_k)
-            //
-            //    E = -dV/dz
-            //-----------------------------------------------
             rsd[index(c_offset_E, j)] = dEdz(x,j) - rho_e(x,j) / epsilon_0;
             diag[index(c_offset_E, j)] = 0;
         }
     }
 }
 
+void IonFlow::evalSpecies(double* x, double* rsd, int* diag,
+                          double rdt, size_t jmin, size_t jmax)
+{
+    //-----------------------------------------------
+    //    Species equations
+    //-----------------------------------------------
+    StFlow::evalSpecies(x, rsd, diag, rdt, jmin, jmax);
+    if (m_stage != 2) {
+        return;
+    }
+
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            // enforcing the flux for charged species is difficult
+            // since charged species are also affected by electric
+            // force, so Neumann boundary condition is used.
+            for (size_t k : m_kCharge) {
+                rsd[index(c_offset_Y + k, 0)] = Y(x,k,0) - Y(x,k,1);
+            }
+        } else if (j == m_points - 1) { // right boundary
+        } else { // interior points
+        }
+    }
+}
+
+
 void IonFlow::solveElectricField(size_t j)
 {
     bool changed = false;

src/oneD/StFlow.cpp

@@ -317,7 +317,22 @@ void StFlow::eval(size_t jg, double* xg, double* rg, integer* diagg, double rdt)
     }
 
     updateProperties(jg, x, jmin, jmax);
-    evalResidual(x, rsd, diag, rdt, jmin, jmax);
+
+    //----------------------------------------------------
+    // Evaluate the residual equations at all required
+    // grid points.
+    // These residuals may be modified by a boundary object.
+    // The boundary object connected will modify
+    // these equations by subtracting its values for V, T, mdot, etc.
+    // As a result, these residual equations will force the solution
+    // variables to the values for the boundary object.
+    //----------------------------------------------------
+    evalContinuity(x, rsd, diag, rdt, jmin, jmax);
+    evalMomentum(x, rsd, diag, rdt, jmin, jmax);
+    evalEnergy(x, rsd, diag, rdt, jmin, jmax);
+    evalSpecies(x, rsd, diag, rdt, jmin, jmax);
+    evalLambda(x, rsd, diag, rdt, jmin, jmax);
+    evalElectricField(x, rsd, diag, rdt, jmin, jmax);
 }
 
 void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
@@ -344,190 +359,6 @@ void StFlow::updateProperties(size_t jg, double* x, size_t jmin, size_t jmax)
     updateDiffFluxes(x, j0, j1);
 }
 
-void StFlow::evalResidual(double* x, double* rsd, int* diag,
-                          double rdt, size_t jmin, size_t jmax)
-{
-    //----------------------------------------------------
-    // evaluate the residual equations at all required
-    // grid points
-    //----------------------------------------------------
-
-    // calculation of qdotRadiation (see docstring of enableRadiation)
-    if (m_do_radiation) {
-        // variable definitions for the Planck absorption coefficient and the
-        // radiation calculation:
-        double k_P_ref = 1.0*OneAtm;
-
-        // polynomial coefficients:
-        const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
-                                     0.51382, -1.86840e-5};
-        const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
-                                     56.310, -5.8169};
-
-        // calculation of the two boundary values
-        double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
-        double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
-
-        // loop over all grid points
-        for (size_t j = jmin; j < jmax; j++) {
-            // helping variable for the calculation
-            double radiative_heat_loss = 0;
-
-            // calculation of the mean Planck absorption coefficient
-            double k_P = 0;
-            // absorption coefficient for H2O
-            if (m_kRadiating[1] != npos) {
-                double k_P_H2O = 0;
-                for (size_t n = 0; n <= 5; n++) {
-                    k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
-                }
-                k_P_H2O /= k_P_ref;
-                k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
-            }
-            // absorption coefficient for CO2
-            if (m_kRadiating[0] != npos) {
-                double k_P_CO2 = 0;
-                for (size_t n = 0; n <= 5; n++) {
-                    k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
-                }
-                k_P_CO2 /= k_P_ref;
-                k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
-            }
-
-            // calculation of the radiative heat loss term
-            radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
-            - boundary_Rad_left - boundary_Rad_right);
-
-            // set the radiative heat loss vector
-            m_qdotRadiation[j] = radiative_heat_loss;
-        }
-    }
-
-    for (size_t j = jmin; j <= jmax; j++) {
-        //----------------------------------------------
-        //         left boundary
-        //----------------------------------------------
-
-        if (j == 0) {
-            // these may be modified by a boundary object
-
-            // Continuity. This propagates information right-to-left, since
-            // rho_u at point 0 is dependent on rho_u at point 1, but not on
-            // mdot from the inlet.
-            rsd[index(c_offset_U,0)] =
-                -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
-                -(density(1)*V(x,1) + density(0)*V(x,0));
-
-            // the inlet (or other) object connected to this one will modify
-            // these equations by subtracting its values for V, T, and mdot. As
-            // a result, these residual equations will force the solution
-            // variables to the values for the boundary object
-            rsd[index(c_offset_V,0)] = V(x,0);
-            rsd[index(c_offset_T,0)] = T(x,0);
-            if (m_usesLambda) {
-                rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
-            } else {
-                rsd[index(c_offset_L, 0)] = lambda(x, 0);
-                diag[index(c_offset_L, 0)] = 0;
-            }
-
-            // The default boundary condition for species is zero flux. However,
-            // the boundary object may modify this.
-            double sum = 0.0;
-            for (size_t k = 0; k < m_nsp; k++) {
-                sum += Y(x,k,0);
-                rsd[index(c_offset_Y + k, 0)] =
-                    -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
-            }
-            rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
-
-            // set residual of poisson's equ to zero
-            rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)];
-        } else if (j == m_points - 1) {
-            evalRightBoundary(x, rsd, diag, rdt);
-        } else { // interior points
-            evalContinuity(j, x, rsd, diag, rdt);
-            // set residual of poisson's equ to zero
-            rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
-
-            //------------------------------------------------
-            //    Radial momentum equation
-            //
-            //    \rho dV/dt + \rho u dV/dz + \rho V^2
-            //       = d(\mu dV/dz)/dz - lambda
-            //-------------------------------------------------
-            if (m_usesLambda) {
-                rsd[index(c_offset_V,j)] =
-                    (shear(x, j) - lambda(x, j) - rho_u(x, j) * dVdz(x, j)
-                    - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j]
-                    - rdt * (V(x, j) - V_prev(j));
-                diag[index(c_offset_V, j)] = 1;
-            } else {
-                rsd[index(c_offset_V, j)] = V(x, j);
-                diag[index(c_offset_V, j)] = 0;
-            }
-
-            //-------------------------------------------------
-            //    Species equations
-            //
-            //   \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz
-            //   = M_k\omega_k
-            //-------------------------------------------------
-            getWdot(x,j);
-            for (size_t k = 0; k < m_nsp; k++) {
-                double convec = rho_u(x,j)*dYdz(x,k,j);
-                double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
-                                / (z(j+1) - z(j-1));
-                rsd[index(c_offset_Y + k, j)]
-                = (m_wt[k]*(wdot(k,j))
-                   - convec - diffus)/m_rho[j]
-                  - rdt*(Y(x,k,j) - Y_prev(k,j));
-                diag[index(c_offset_Y + k, j)] = 1;
-            }
-
-            //-----------------------------------------------
-            //    energy equation
-            //
-            //    \rho c_p dT/dt + \rho c_p u dT/dz
-            //    = d(k dT/dz)/dz
-            //      - sum_k(\omega_k h_k_ref)
-            //      - sum_k(J_k c_p_k / M_k) dT/dz
-            //-----------------------------------------------
-            if (m_do_energy[j]) {
-
-                setGas(x,j);
-                double dtdzj = dTdz(x,j);
-                double sum = 0.0;
-
-                grad_hk(x, j);
-                for (size_t k = 0; k < m_nsp; k++) {
-                    double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
-                    sum += wdot(k,j)*m_hk(k,j);
-                    sum += flxk * m_dhk_dz(k,j) / m_wt[k];
-                }
-
-                rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
-                                            - divHeatFlux(x,j) - sum;
-                rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
-                rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
-                rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j]));
-                diag[index(c_offset_T, j)] = 1;
-            } else {
-                // residual equations if the energy equation is disabled
-                rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
-                diag[index(c_offset_T, j)] = 0;
-            }
-
-            if (m_usesLambda) {
-                rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
-            } else {
-                rsd[index(c_offset_L, j)] = lambda(x, j);
-            }
-            diag[index(c_offset_L, j)] = 0;
-        }
-    }
-}
-
 void StFlow::updateTransport(double* x, size_t j0, size_t j1)
 {
      if (m_do_multicomponent) {
@@ -558,23 +389,6 @@ void StFlow::updateTransport(double* x, size_t j0, size_t j1)
     }
 }
 
-void StFlow::show(const double* x)
-{
-    writelog("    Pressure:  {:10.4g} Pa\n", m_press);
-
-    Domain1D::show(x);
-
-    if (m_do_radiation) {
-        writeline('-', 79, false, true);
-        writelog("\n          z      radiative heat loss");
-        writeline('-', 79, false, true);
-        for (size_t j = 0; j < m_points; j++) {
-            writelog("\n {:10.4g}        {:10.4g}", m_z[j], m_qdotRadiation[j]);
-        }
-        writelog("\n");
-    }
-}
-
 void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
 {
     if (m_do_multicomponent) {
@@ -617,6 +431,312 @@ void StFlow::updateDiffFluxes(const double* x, size_t j0, size_t j1)
     }
 }
 
+void StFlow::evalContinuity(double* x, double* rsd, int* diag,
+                          double rdt, size_t jmin, size_t jmax)
+{
+    //----------------------------------------------
+    //    Continuity equation
+    //
+    //    d(\rho u)/dz + 2\rho V = 0
+    //----------------------------------------------
+    // This propagates information right-to-left, since rho_u at point 0 is dependent
+    // on rho_u at point 1, but not on mdot from the inlet.
+    // The default value for continuity is zero u at the boundary.
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            rsd[index(c_offset_U,0)] =
+                -(rho_u(x,1) - rho_u(x,0))/m_dz[0]
+                -(density(1)*V(x,1) + density(0)*V(x,0));
+        } else if (j == m_points - 1) { // right boundary
+            if (m_usesLambda) {
+                rsd[index(c_offset_U, j)] = rho_u(x, j);
+            } else {
+                rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1);
+            }
+        } else { // interior points
+            if (m_usesLambda) {
+                // 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.
+                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));
+            } else if (m_isFree) {
+                // terms involving V are zero as V=0 by definition
+                if (grid(j) > m_zfixed) {
+                    rsd[index(c_offset_U,j)] =
+                        - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1];
+                } else if (grid(j) == m_zfixed) {
+                    if (m_do_energy[j]) {
+                        rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
+                    } else {
+                        rsd[index(c_offset_U,j)] = (rho_u(x,j)
+                                                    - m_rho[0]*0.3); // why 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];
+                }
+            } else {
+                // unstrained with fixed mass flow rate
+                rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
+            }
+            diag[index(c_offset_U, j)] = 0; // Algebraic constraint
+        }
+    }
+}
+
+void StFlow::evalMomentum(double* x, double* rsd, int* diag,
+                          double rdt, size_t jmin, size_t jmax)
+{
+    //------------------------------------------------
+    //    Radial momentum equation
+    //
+    //    \rho dV/dt + \rho u dV/dz + \rho V^2
+    //       = d(\mu dV/dz)/dz - lambda
+    //-------------------------------------------------
+    // The default boundary condition is zero V.
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            rsd[index(c_offset_V,0)] = V(x,0);
+        } else if (j == m_points - 1) { // right boundary
+            rsd[index(c_offset_V,j)] = V(x,j);
+            diag[index(c_offset_V,j)] = 0;
+        } else { // interior points
+            if (m_usesLambda) {
+                rsd[index(c_offset_V,j)] =
+                    (shear(x, j) - lambda(x, j) - rho_u(x, j) * dVdz(x, j)
+                    - m_rho[j] * V(x, j) * V(x, j)) / m_rho[j]
+                    - rdt * (V(x, j) - V_prev(j));
+                diag[index(c_offset_V, j)] = 1;
+            } else {
+                rsd[index(c_offset_V, j)] = V(x, j);
+                diag[index(c_offset_V, j)] = 0;
+            }
+        }
+    }
+}
+
+void StFlow::evalEnergy(double* x, double* rsd, int* diag,
+                        double rdt, size_t jmin, size_t jmax)
+{
+    //-----------------------------------------------
+    //    Energy equation
+    //
+    //    \rho c_p dT/dt + \rho c_p u dT/dz
+    //    = d(k dT/dz)/dz
+    //      - sum_k(\omega_k h_k_ref)
+    //      - sum_k(J_k c_p_k / M_k) dT/dz
+    //-----------------------------------------------
+    // The default boundary condition for energy is a zero T.
+
+    // Calculation of qdotRadiation (see docstring of enableRadiation)
+    if (m_do_radiation) {
+        computeRadiation(x, jmin, jmax);
+    }
+
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            rsd[index(c_offset_T,0)] = T(x,0);
+        } else if (j == m_points - 1) { // right boundary
+            rsd[index(c_offset_T, j)] = T(x, j);
+        } else { // interior points
+            if (m_do_energy[j]) {
+                double dtdzj = dTdz(x,j);
+                double sum = 0.0;
+
+                grad_hk(x, j);
+                for (size_t k = 0; k < m_nsp; k++) {
+                    double flxk = 0.5*(m_flux(k,j-1) + m_flux(k,j));
+                    sum += wdot(k,j)*m_hk(k,j);
+                    sum += flxk * m_dhk_dz(k,j) / m_wt[k];
+                }
+
+                rsd[index(c_offset_T, j)] = - m_cp[j]*rho_u(x,j)*dtdzj
+                                            - divHeatFlux(x,j) - sum;
+                rsd[index(c_offset_T, j)] /= (m_rho[j]*m_cp[j]);
+                rsd[index(c_offset_T, j)] -= rdt*(T(x,j) - T_prev(j));
+                rsd[index(c_offset_T, j)] -= (m_qdotRadiation[j] / (m_rho[j] * m_cp[j]));
+                diag[index(c_offset_T, j)] = 1;
+            } else {
+                // residual equations if the energy equation is disabled
+                rsd[index(c_offset_T, j)] = T(x,j) - T_fixed(j);
+                diag[index(c_offset_T, j)] = 0;
+            }
+        }
+    }
+}
+
+void StFlow::evalSpecies(double* x, double* rsd, int* diag,
+                        double rdt, size_t jmin, size_t jmax)
+{
+    //-------------------------------------------------
+    //    Species equations
+    //
+    //   \rho dY_k/dt + \rho u dY_k/dz + dJ_k/dz
+    //   = M_k\omega_k
+    //-------------------------------------------------
+    // The default boundary condition for species is zero flux.
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            double sum = 0.0;
+            for (size_t k = 0; k < m_nsp; k++) {
+                sum += Y(x,k,0);
+                rsd[index(c_offset_Y + k, 0)] =
+                    -(m_flux(k,0) + rho_u(x,0)* Y(x,k,0));
+            }
+            rsd[index(c_offset_Y + leftExcessSpecies(), 0)] = 1.0 - sum;
+        } else if (j == m_points - 1) { // right boundary
+            double 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;
+        } else { // interior points
+            for (size_t k = 0; k < m_nsp; k++) {
+                double convec = rho_u(x,j)*dYdz(x,k,j);
+                double diffus = 2.0*(m_flux(k,j) - m_flux(k,j-1))
+                                / (z(j+1) - z(j-1));
+                rsd[index(c_offset_Y + k, j)]
+                = (m_wt[k]*(wdot(k,j))
+                   - convec - diffus)/m_rho[j]
+                  - rdt*(Y(x,k,j) - Y_prev(k,j));
+                diag[index(c_offset_Y + k, j)] = 1;
+            }
+        }
+    }
+}
+
+void StFlow::evalLambda(double* x, double* rsd, int* diag,
+                        double rdt, size_t jmin, size_t jmax)
+{
+    //-------------------------------------------------
+    //    Lambda equation
+    //-------------------------------------------------
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            if (m_usesLambda) {
+                rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
+            } else {
+                rsd[index(c_offset_L, 0)] = lambda(x, 0);
+                diag[index(c_offset_L, 0)] = 0;
+            }
+        } else if (j == m_points - 1) { // right boundary
+            rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1);
+            diag[index(c_offset_L, j)] = 0;
+        } else { // interior points
+            if (m_usesLambda) {
+                rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j - 1);
+            } else {
+                rsd[index(c_offset_L, j)] = lambda(x, j);
+            }
+            diag[index(c_offset_L, j)] = 0;
+        }
+    }
+}
+
+void StFlow::evalElectricField(double* x, double* rsd, int* diag,
+                               double rdt, size_t jmin, size_t jmax)
+{
+    //-----------------------------------------------
+    //    Electric field (by Gauss's law)
+    //
+    //    dE/dz = 0
+    //
+    //    E = -dV/dz
+    //-----------------------------------------------
+    // The default value for E at the boundary is zero, and E is zero
+    // in the domain as well.
+    for (size_t j = jmin; j <= jmax; j++) {
+        if (j == 0) { // left boundary
+            rsd[index(c_offset_E, 0)] = x[index(c_offset_E, j)];
+        } else if (j == m_points - 1) { // right boundary
+            rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
+        } else { // interior points
+            rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
+        }
+    }
+}
+
+void StFlow::computeRadiation(double* x, size_t jmin, size_t jmax)
+{
+    // Variable definitions for the Planck absorption coefficient and the
+    // radiation calculation:
+    double k_P_ref = 1.0*OneAtm;
+
+    // Polynomial coefficients:
+    const double c_H2O[6] = {-0.23093, -1.12390, 9.41530, -2.99880,
+                                    0.51382, -1.86840e-5};
+    const double c_CO2[6] = {18.741, -121.310, 273.500, -194.050,
+                                    56.310, -5.8169};
+
+    // Calculation of the two boundary values
+    double boundary_Rad_left = m_epsilon_left * StefanBoltz * pow(T(x, 0), 4);
+    double boundary_Rad_right = m_epsilon_right * StefanBoltz * pow(T(x, m_points - 1), 4);
+
+    for (size_t j = jmin; j < jmax; j++) {
+        // calculation of the mean Planck absorption coefficient
+        double k_P = 0;
+        // Absorption coefficient for H2O
+        if (m_kRadiating[1] != npos) {
+            double k_P_H2O = 0;
+            for (size_t n = 0; n <= 5; n++) {
+                k_P_H2O += c_H2O[n] * pow(1000 / T(x, j), (double) n);
+            }
+            k_P_H2O /= k_P_ref;
+            k_P += m_press * X(x, m_kRadiating[1], j) * k_P_H2O;
+        }
+        // Absorption coefficient for CO2
+        if (m_kRadiating[0] != npos) {
+            double k_P_CO2 = 0;
+            for (size_t n = 0; n <= 5; n++) {
+                k_P_CO2 += c_CO2[n] * pow(1000 / T(x, j), (double) n);
+            }
+            k_P_CO2 /= k_P_ref;
+            k_P += m_press * X(x, m_kRadiating[0], j) * k_P_CO2;
+        }
+
+        // Calculation of the radiative heat loss term
+        double radiative_heat_loss = 0;
+        radiative_heat_loss = 2 * k_P *(2 * StefanBoltz * pow(T(x, j), 4)
+                              - boundary_Rad_left - boundary_Rad_right);
+
+        // set the radiative heat loss vector
+        m_qdotRadiation[j] = radiative_heat_loss;
+    }
+}
+
+void StFlow::grad_hk(const double* x, size_t j)
+{
+    for(size_t k = 0; k < m_nsp; k++) {
+        if (u(x, j) > 0.0) {
+            m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1];
+        }
+        else {
+            m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j];
+        }
+    }
+}
+
+void StFlow::show(const double* x)
+{
+    writelog("    Pressure:  {:10.4g} Pa\n", m_press);
+
+    Domain1D::show(x);
+
+    if (m_do_radiation) {
+        writeline('-', 79, false, true);
+        writelog("\n          z      radiative heat loss");
+        writeline('-', 79, false, true);
+        for (size_t j = 0; j < m_points; j++) {
+            writelog("\n {:10.4g}        {:10.4g}", m_z[j], m_qdotRadiation[j]);
+        }
+        writelog("\n");
+    }
+}
+
 string StFlow::componentName(size_t n) const
 {
     switch (n) {
@@ -937,84 +1057,4 @@ void StFlow::fixTemperature(size_t j)
     }
 }
 
-void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double 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_V,j)] = V(x,j);
-    diag[index(c_offset_V,j)] = 0;
-    double sum = 0.0;
-    // set residual of poisson's equ to zero
-    rsd[index(c_offset_E, j)] = x[index(c_offset_E, j)];
-    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;
-    if (m_usesLambda) {
-        rsd[index(c_offset_U, j)] = rho_u(x, j);
-    } else {
-        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j-1);
-    }
-
-    rsd[index(c_offset_L, j)] = lambda(x, j) - lambda(x, j-1);
-    diag[index(c_offset_L, j)] = 0;
-    rsd[index(c_offset_T, j)] = T(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
-    //
-    //    d(\rho u)/dz + 2\rho V = 0
-    //----------------------------------------------
-    if (m_usesLambda) {
-        // 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.
-        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));
-    } else if (m_isFree) {
-        // terms involving V are zero as V=0 by definition
-        if (grid(j) > m_zfixed) {
-            rsd[index(c_offset_U,j)] =
-                - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1];
-        } else if (grid(j) == m_zfixed) {
-            if (m_do_energy[j]) {
-                rsd[index(c_offset_U,j)] = (T(x,j) - m_tfixed);
-            } else {
-                rsd[index(c_offset_U,j)] = (rho_u(x,j)
-                                            - m_rho[0]*0.3); // why 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];
-        }
-    } else {
-        // unstrained with fixed mass flow rate
-        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
-    }
-}
-
-void StFlow::grad_hk(const double* x, size_t j)
-{
-    for(size_t k = 0; k < m_nsp; k++) {
-        if (u(x, j) > 0.0) {
-            m_dhk_dz(k,j) = (m_hk(k,j) - m_hk(k,j-1)) / m_dz[j - 1];
-        }
-        else {
-            m_dhk_dz(k,j) = (m_hk(k,j+1) - m_hk(k,j)) / m_dz[j];
-        }
-    }
-}
-
 } // namespace