[oneD] Propagate mdot for unstrained flows left-to-right

src/oneD/Boundary1D.cpp

@@ -198,10 +198,13 @@ void Inlet1D::eval(size_t jg, double* xg, double* rg,
             // Set mdot equal to rho*u, and also set lambda to zero.
             m_mdot = m_flow->density(0) * xb[c_offset_U];
             rb[c_offset_L] = xb[c_offset_L];
-        } else {
+        } else if (m_flow->usesLambda()) {
             // The flow domain sets this to -rho*u. Add mdot to specify the mass
-            // flow rate (both axisymmetric and unstrained).
+            // flow rate
             rb[c_offset_L] += m_mdot;
+        } else {
+            rb[c_offset_U] = m_flow->density(0) * xb[c_offset_U] - m_mdot;
+            rb[c_offset_L] = xb[c_offset_L];
         }
 
         // add the convective term to the species residual equations
@@ -374,28 +377,18 @@ void Outlet1D::eval(size_t jg, double* xg, double* rg, integer* diagg,
         size_t nc = m_flow_right->nComponents();
         double* xb = x;
         double* rb = r;
-        if (!m_flow_right->usesLambda()) {
-            rb[c_offset_L] = xb[c_offset_L];
-        }
-        if (m_flow_right->doEnergy(0)) {
-            rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T + nc];
-        }
         for (size_t k = c_offset_Y; k < nc; k++) {
             rb[k] = xb[k] - xb[k + nc];
         }
     }
 
     if (m_flow_left) {
+        // flow is left-to-right
         size_t nc = m_flow_left->nComponents();
         double* xb = x - nc;
         double* rb = r - nc;
         int* db = diag - nc;
 
-        // zero Lambda
-        if (!m_flow_left->isFree()) {
-            rb[c_offset_U] = xb[c_offset_L];
-        }
-
         if (m_flow_left->doEnergy(m_flow_left->nPoints()-1)) {
             rb[c_offset_T] = xb[c_offset_T] - xb[c_offset_T - nc]; // zero T gradient
         }

src/oneD/StFlow.cpp

@@ -488,7 +488,7 @@ void StFlow::evalResidual(double* x, double* rsd, int* diag,
             if (m_isFree) {
                 rsd[index(c_offset_L, 0)] = lambda(x, 0);
             } else {
-                // used to set mass flow rate
+                // used to set mass flow rate (both axisymmetric and unstrained)
                 rsd[index(c_offset_L, 0)] = -rho_u(x, 0);
             }
 
@@ -1025,16 +1025,23 @@ void StFlow::evalRightBoundary(double* x, double* rsd, int* diag, double rdt)
     }
     rsd[index(c_offset_Y + rightExcessSpecies(), j)] = 1.0 - sum;
     diag[index(c_offset_Y + rightExcessSpecies(), j)] = 0;
-    if (!m_isFree) {
+    if (m_isFree) {
+        rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
+        rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
+    } else if (m_usesLambda) {
         rsd[index(c_offset_U,j)] = rho_u(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);
         }
     } else {
-        rsd[index(c_offset_U,j)] = rho_u(x,j) - rho_u(x,j-1);
-        rsd[index(c_offset_T,j)] = T(x,j) - T(x,j-1);
+        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
+        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);
+        }
     }
 }
 
@@ -1047,14 +1054,14 @@ void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double
     //
     //    d(\rho u)/dz + 2\rho V = 0
     //----------------------------------------------
-    if (!m_isFree) {
+    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 {
+    } else if (m_isFree) {
         if (grid(j) > m_zfixed) {
             rsd[index(c_offset_U,j)] =
                 - (rho_u(x,j) - rho_u(x,j-1))/m_dz[j-1]
@@ -1071,6 +1078,8 @@ void StFlow::evalContinuity(size_t j, double* x, double* rsd, int* diag, double
                 - (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 {
+        rsd[index(c_offset_U, j)] = rho_u(x, j) - rho_u(x, j - 1);
     }
 }