[1D] Add clarification

src/oneD/Flow1D.cpp

@@ -732,21 +732,25 @@ void Flow1D::evalSpecies(double* x, double* rsd, int* diag,
         for (size_t k = 0; k < m_nsp; k++) {
             sum += Y(x,k,jmin);
             rsd[index(c_offset_Y+k, jmin)] = -(m_flux(k, jmin) +
-                                                rho_u(x, jmin) * Y(x, k, jmin));
+                                               rho_u(x, jmin) * Y(x, k, jmin));
+            // diag[index(c_offset_Y+k, jmin)] are zero and do not need to be reset
         }
+        // Residuals may be modified by left boundary (for k != leftExcessSpecies);
+        // excess species residual is handled as algebraic constraint.
         rsd[index(c_offset_Y + leftExcessSpecies(), jmin)] = 1.0 - sum;
-        diag[index(c_offset_Y + leftExcessSpecies(), jmin)] = 0;
     }
 
     if (jmax == m_points - 1) { // right boundary
         double sum = 0.0;
         for (size_t k = 0; k < m_nsp; k++) {
             sum += Y(x,k,jmax);
-            rsd[index(k+c_offset_Y, jmax)] = m_flux(k, jmax-1) +
+            rsd[index(c_offset_Y+k, jmax)] = m_flux(k, jmax-1) +
                                              rho_u(x, jmax)*Y(x, k, jmax);
+            // diag[index(c_offset_Y+k, jmin)] are zero and do not need to be reset
         }
+        // Residuals may be modified by right boundary (for k != rightExcessSpecies);
+        // excess species residual is handled as algebraic constraint.
         rsd[index(c_offset_Y + rightExcessSpecies(), jmax)] = 1.0 - sum;
-        diag[index(c_offset_Y + rightExcessSpecies(), jmax)] = 0;
     }
 
     // j0 and j1 are constrained to only interior points
@@ -757,7 +761,7 @@ void Flow1D::evalSpecies(double* x, double* rsd, int* diag,
             double convec = rho_u(x, j)*dYdz(x, k, j);
             double diffus = 2*(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]*m_wdot(k, j)
-                                              - convec - diffus) / m_rho[j]
+                                             - convec - diffus) / m_rho[j]
                                             - rdt*(Y(x, k, j) - Y_prev(k, j));
             diag[index(c_offset_Y + k, j)] = 1;
         }