added comments, fixed incorrect bc, other minor adjustments to two-point method

include/cantera/oneD/Sim1D.h

@@ -196,11 +196,11 @@ public:
     double fixedTemperatureLocation();
 
     //! ------ One and Two-Point flame control methods
-    //! Set the left control point location. This is used for one or two 
+    //! Set the left control point location. This is used for two
     //! point flame control.
     void setLeftControlPoint(double temperature);
 
-    //! Set the right  control point location. This is used for one or two 
+    //! Set the right  control point location. This is used for two
     //! point flame control.
     void setRightControlPoint(double temperature);
     //! -------------------

include/cantera/oneD/StFlow.h

@@ -260,14 +260,21 @@ public:
     //! the value of the solution at these points is fixed. The values of the control
     //! points are dictated and thus serve as a boundary condition that affects the
     //! solution of the governing equations in the 1D domain. The imposition of fixed
-    //! points in the domain means that the original set of governing equations' boundary 
+    //! points in the domain means that the original set of governing equations' boundary
     //! conditions would over-specify the problem. Thus, the boundary conditions are changed
     //! to reflect the fact that the control points are serving as internal boundary conditions.
     //!
+    //! In this method, the imposition of the two internal boundary conditions requires that two other
+    //! boundary conditions be changed. The first is the boundary condition for the continuity equation
+    //! at the left boundary, which is changed to be a value that is derived from the solution at the
+    //! left boundary. The second is the continuity boundary condition at the right boundary, which is
+    //! also determined from the flow solution by using the oxidizer axial velocity equation variable to
+    //! compute the mass flux at the right boundary.
+    //!
     //! This method is based on the work of M. Nishioka, C.K. Law, and T. Takeno (1996) titled
-    //! "A Flame-Controlling Continuation Method for Generating S-Curve Responses with 
+    //! "A Flame-Controlling Continuation Method for Generating S-Curve Responses with
     //! Detailed Chemistry"
-    
+
     //! The current left control point temperature
     double leftControlPointTemperature() const {
         if (m_twoPointControl && (m_zLeft != Undef)) {
@@ -801,7 +808,7 @@ public:
 
     //! -------------
 
-    
+
 private:
     vector<double> m_ybar;
 };

src/oneD/Boundary1D.cpp

@@ -209,6 +209,11 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
             m_mdot = m_flow->density(0) * xb[c_offset_U];
         } else if (m_flow->isStrained()) { //axisymmetric flow
             if (m_flow->twoPointControlEnabled()) {
+                // When using two-point control, the mass flow rate at the left inlet is not
+                // specified. Instead, the mass flow rate is dictated by the velocity at the
+                // left inlet, which comes from the U variable. The default boundary condition specified
+                // in the StFlow.cpp file already handles this case. We only need to update the stored
+                // value of m_mdot so that other equations that use the quantity are consistent.
                 m_mdot = m_flow->density(0)*xb[c_offset_U];
             } else {
                 // The flow domain sets this to -rho*u. Add mdot to specify the mass
@@ -246,12 +251,13 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
 
         if (m_flow->twoPointControlEnabled()) {// For point control adjustments
             // At the right boundary, the mdot is dictated by the velocity at the
-            // right boundary, which comes from the Uo variable. 
+            // right boundary, which comes from the Uo variable. The variable Uo is
+            // the left-moving velocity and has a negative value, so the mass flow has
+            // to be negated to give a positive value when using Uo.
             m_mdot = -(m_flow->density(last_index) * xb[c_offset_Uo]);
-            rb[c_offset_U] += m_mdot; 
+            rb[c_offset_U] += m_mdot;
         } else {
             rb[c_offset_U] += m_mdot;
-            rb[c_offset_Uo] += m_mdot/m_flow->density(last_index);
         }
 
         for (size_t k = 0; k < m_nsp; k++) {

src/oneD/StFlow.cpp

@@ -423,7 +423,7 @@ void StFlow::evalContinuity(double* x, double* rsd, int* diag,
         rsd[index(c_offset_U,jmin)] = -(rho_u(x,jmin + 1) - rho_u(x,jmin))/m_dz[jmin]
                                       -(density(jmin + 1)*V(x,jmin + 1)
                                       + density(jmin)*V(x,jmin));
-        
+
         diag[index(c_offset_U,jmin)] = 0; // Algebraic constraint
     }
 
@@ -447,7 +447,7 @@ void StFlow::evalContinuity(double* x, double* rsd, int* diag,
             // 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));
-            
+
             diag[index(c_offset_U, j)] = 0; // Algebraic constraint
         }
     } else if (m_isFree) { // "free-flow"
@@ -516,7 +516,7 @@ void StFlow::evalLambda(double* x, double* rsd, int* diag,
         }
         return;
     }
- 
+
     if (jmin == 0) { // left boundary
         if (m_twoPointControl) {
             rsd[index(c_offset_L, jmin)] = lambda(x,jmin+1) - lambda(x,jmin);
@@ -533,7 +533,7 @@ void StFlow::evalLambda(double* x, double* rsd, int* diag,
     // j0 and j1 are constrained to only interior points
     size_t j0 = std::max<size_t>(jmin, 1);
     size_t j1 = std::min(jmax, m_points - 2);
-    double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate
+    double epsilon = 1e-8; // Precision threshold for being 'equal' to a coordinate
     for (size_t j = j0; j <= j1; j++) { // interior points
         if (m_twoPointControl) {
             if (std::abs(grid(j) - m_zLeft) < epsilon ) {
@@ -601,15 +601,13 @@ void StFlow::evalUo(double* x, double* rsd, int* diag,
     if (jmin == 0) { // left boundary
         // Because the Uo equation is used for two-point control, the boundary
         // for Uo is located in the domain interior at the right control point,
-        // thus at the boundary, the 
+        // thus at the boundary, the
         rsd[index(c_offset_Uo,jmin)] = Uo(x,jmin+1) - Uo(x,jmin);
     }
 
     if (jmax == m_points - 1) { // right boundary
         if(m_twoPointControl) {
             rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax) - Uo(x,jmax-1);
-        } else {
-            rsd[index(c_offset_Uo, jmax)] = Uo(x,jmax);
         }
         diag[index(c_offset_Uo, jmax)] = 0;
     }
@@ -617,7 +615,7 @@ void StFlow::evalUo(double* x, double* rsd, int* diag,
     // j0 and j1 are constrained to only interior points
     size_t j0 = std::max<size_t>(jmin, 1);
     size_t j1 = std::min(jmax, m_points - 2);
-    double epsilon = 1e-5; // Precision threshold for being 'equal' to a coordinate
+    double epsilon = 1e-8; // Precision threshold for being 'equal' to a coordinate
     for (size_t j = j0; j <= j1; j++) { // interior points
         if (m_twoPointControl) {
             if (std::abs(grid(j) - m_zRight) < epsilon) {