diff --git a/src/M2ulPhyS.cpp b/src/M2ulPhyS.cpp
index 2ad80cfe..1576be49 100644
--- a/src/M2ulPhyS.cpp
+++ b/src/M2ulPhyS.cpp
@@ -237,7 +237,7 @@ void M2ulPhyS::initMixtureAndTransportModels() {
                                               TPSCommWorld, trans_data, trans_tables);
         if (!success) exit(ERROR);
 
-          // Instantiate LteTransport class
+        // Instantiate LteTransport class
 #if defined(_CUDA_) || defined(_HIP_)
         // Tables from above have host pointers.  Must get device
         // pointers here before instantiating device class
@@ -1784,7 +1784,7 @@ void M2ulPhyS::initSolutionAndVisualizationVectors() {
         // visualizationNames_.push_back(std::string("rxn_rate: " + config.reactionEquations[r]));
       }
     }  // if (config.workFluid != DRY_AIR)
-  }    // if tpsP->isVisualizationMode()
+  }  // if tpsP->isVisualizationMode()
 
   // If mms, add conserved and exact solution.
 #ifdef HAVE_MASA
@@ -2094,7 +2094,7 @@ void M2ulPhyS::solveStep() {
       exit(ERROR);
 #endif
     }  // plane dump
-  }    // step check
+  }  // step check
 
   average->addSample(iter, d_mixture);
 }
@@ -2465,7 +2465,9 @@ void M2ulPhyS::Check_NAN() {
   int dof = vfes->GetNDofs();
 
 #ifdef _GPU_
-  { local_print = M2ulPhyS::Check_NaN_GPU(U, dof * num_equation, loc_print); }
+  {
+    local_print = M2ulPhyS::Check_NaN_GPU(U, dof * num_equation, loc_print);
+  }
   if (local_print > 0) {
     cout << "Found a NaN!" << endl;
   }
@@ -4234,7 +4236,7 @@ void M2ulPhyS::updateVisualizationVariables() {
         dataVis[visualIdxs.rxn + r][n] = progressRates[r];
       }
     }  // if (!isDryAir)
-  }    // for (int n = 0; n < ndofs; n++)
+  }  // for (int n = 0; n < ndofs; n++)
 }
 
 void M2ulPhyS::evaluatePlasmaConductivityGF() {
diff --git a/src/calorically_perfect.cpp b/src/calorically_perfect.cpp
index c182c14e..1c122b7f 100644
--- a/src/calorically_perfect.cpp
+++ b/src/calorically_perfect.cpp
@@ -58,10 +58,12 @@ MFEM_HOST_DEVICE double Sutherland(const double T, const double mu_star, const d
 }
 
 CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts,
-                                                           temporalSchemeCoefficients &time_coeff, TPS::Tps *tps)
+                                                           temporalSchemeCoefficients &time_coeff,
+                                                           ParGridFunction *gridScale, TPS::Tps *tps)
     : tpsP_(tps), pmesh_(pmesh), time_coeff_(time_coeff) {
   rank0_ = (pmesh_->GetMyRank() == 0);
   order_ = loMach_opts->order;
+  gridScale_gf_ = gridScale;
 
   std::string visc_model;
   tpsP_->getInput("loMach/calperfect/viscosity-model", visc_model, std::string("sutherland"));
@@ -132,6 +134,9 @@ CaloricallyPerfectThermoChem::CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh,
   tpsP_->getInput("loMach/calperfect/msolve-atol", mass_inverse_atol_, default_atol_);
   tpsP_->getInput("loMach/calperfect/msolve-max-iter", mass_inverse_max_iter_, max_iter_);
   tpsP_->getInput("loMach/calperfect/msolve-verbosity", mass_inverse_pl_, pl_solve_);
+
+  // artificial diffusion (SUPG)
+  tpsP_->getInput("loMach/calperfect/streamwise-stabilization", sw_stab_, false);
 }
 
 CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
@@ -164,6 +169,20 @@ CaloricallyPerfectThermoChem::~CaloricallyPerfectThermoChem() {
   delete rho_over_dt_coeff_;
   delete rho_coeff_;
 
+  delete umag_coeff_;
+  delete gscale_coeff_;
+  delete visc_coeff_;
+  delete visc_inv_coeff_;
+  delete reh1_coeff_;
+  delete reh2_coeff_;
+  delete Reh_coeff_;
+  delete csupg_coeff_;
+  delete uw1_coeff_;
+  delete uw2_coeff_;
+  delete upwind_coeff_;
+  delete swdiff_coeff_;
+  delete supg_coeff_;
+
   // allocated in initializeSelf
   delete sfes_;
   delete sfec_;
@@ -448,6 +467,32 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   rhon_next_coeff_ = new GridFunctionCoefficient(&rn_gf_);
   rhou_coeff_ = new ScalarVectorProductCoefficient(*rhon_next_coeff_, *un_next_coeff_);
 
+  // artifical diffusion coefficients
+  if (sw_stab_) {
+    umag_coeff_ = new VectorMagnitudeCoefficient(*un_next_coeff_);
+    gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_);
+    visc_coeff_ = new GridFunctionCoefficient(&visc_gf_);
+    visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0);
+
+    // compute Reh
+    reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_);
+    reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_);
+    Reh_coeff_ = new ProductCoefficient(*reh1_coeff_, *umag_coeff_);
+
+    // Csupg
+    csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor);
+
+    // compute upwind magnitude
+    uw1_coeff_ = new ProductCoefficient(*rho_coeff_, *csupg_coeff_);
+    uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+    upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+    // streamwise diffusion direction
+    swdiff_coeff_ = new TransformedMatrixVectorCoefficient(un_next_coeff_, &streamwiseTensor);
+
+    supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+  }
+
   At_form_ = new ParBilinearForm(sfes_);
   auto *at_blfi = new ConvectionIntegrator(*rhou_coeff_);
   if (numerical_integ_) {
@@ -483,6 +528,7 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
   MsRho_form_->FormSystemMatrix(empty, MsRho_);
   if (rank0_) std::cout << "CaloricallyPerfectThermoChem MsRho operator set" << endl;
 
+  // Helmholtz
   Ht_form_ = new ParBilinearForm(sfes_);
   auto *hmt_blfi = new MassIntegrator(*rho_over_dt_coeff_);
   auto *hdt_blfi = new DiffusionIntegrator(*thermal_diff_total_coeff_);
@@ -491,6 +537,15 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
     hmt_blfi->SetIntRule(&ir_di);
     hdt_blfi->SetIntRule(&ir_di);
   }
+  // SUPG
+  if (sw_stab_) {
+    auto *sdt_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   sdt_blfi->SetIntRule(&ir_di);
+    // }
+    Ht_form_->AddDomainIntegrator(sdt_blfi);
+  }
   Ht_form_->AddDomainIntegrator(hmt_blfi);
   Ht_form_->AddDomainIntegrator(hdt_blfi);
   Ht_form_->Assemble();
@@ -574,6 +629,18 @@ void CaloricallyPerfectThermoChem::initializeOperators() {
     lqd_blfi->SetIntRule(&ir_di);
   }
   LQ_form_->AddDomainIntegrator(lqd_blfi);
+
+  // DiffusionIntegrator *slqd_blfi;
+  if (sw_stab_) {
+    // slqd_blfi = new DiffusionIntegrator(*supg_coeff_);
+    auto *slqd_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   slqd_blfi->SetIntRule(&ir_di);
+    // }
+    LQ_form_->AddDomainIntegrator(slqd_blfi);
+  }
+
   if (partial_assembly_) {
     LQ_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
   }
diff --git a/src/calorically_perfect.hpp b/src/calorically_perfect.hpp
index 118cd3dd..dcd6b69b 100644
--- a/src/calorically_perfect.hpp
+++ b/src/calorically_perfect.hpp
@@ -48,6 +48,7 @@ class Tps;
 #include "io.hpp"
 #include "thermo_chem_base.hpp"
 #include "tps_mfem_wrap.hpp"
+#include "utils.hpp"
 
 using VecFuncT = void(const Vector &x, double t, Vector &u);
 using ScalarFuncT = double(const Vector &x, double t);
@@ -136,6 +137,9 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   // Initial temperature value (if constant IC)
   double T_ic_;
 
+  // streamwise-stabilization
+  bool sw_stab_;
+
   // FEM related fields and objects
 
   // Scalar \f$H^1\f$ finite element collection.
@@ -155,6 +159,8 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   ParGridFunction R0PM0_gf_;
   ParGridFunction Qt_gf_;
 
+  ParGridFunction *gridScale_gf_ = nullptr;
+
   // ParGridFunction *buffer_tInlet_ = nullptr;
   GridFunctionCoefficient *temperature_bc_field_ = nullptr;
 
@@ -171,6 +177,20 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
   GridFunctionCoefficient *rho_over_dt_coeff_ = nullptr;
   GridFunctionCoefficient *rho_coeff_ = nullptr;
 
+  VectorMagnitudeCoefficient *umag_coeff_ = nullptr;
+  GridFunctionCoefficient *gscale_coeff_ = nullptr;
+  GridFunctionCoefficient *visc_coeff_ = nullptr;
+  PowerCoefficient *visc_inv_coeff_ = nullptr;
+  ProductCoefficient *reh1_coeff_ = nullptr;
+  ProductCoefficient *reh2_coeff_ = nullptr;
+  ProductCoefficient *Reh_coeff_ = nullptr;
+  TransformedCoefficient *csupg_coeff_ = nullptr;
+  ProductCoefficient *uw1_coeff_ = nullptr;
+  ProductCoefficient *uw2_coeff_ = nullptr;
+  ProductCoefficient *upwind_coeff_ = nullptr;
+  TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr;
+  ScalarMatrixProductCoefficient *supg_coeff_ = nullptr;
+
   // operators and solvers
   ParBilinearForm *At_form_ = nullptr;
   ParBilinearForm *Ms_form_ = nullptr;
@@ -233,7 +253,7 @@ class CaloricallyPerfectThermoChem : public ThermoChemModelBase {
 
  public:
   CaloricallyPerfectThermoChem(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff,
-                               TPS::Tps *tps);
+                               ParGridFunction *gridScale, TPS::Tps *tps);
   virtual ~CaloricallyPerfectThermoChem();
 
   // Functions overriden from base class
diff --git a/src/em_options.hpp b/src/em_options.hpp
index e78cf70b..d192c88f 100644
--- a/src/em_options.hpp
+++ b/src/em_options.hpp
@@ -52,8 +52,8 @@ class ElectromagneticOptions {
   double atol;                            /**< Linear solver absolute tolerance */
   double preconditioner_background_sigma; /**< Uniform conductivity to use to preconditioner (ignored if <= 0) */
 
-  bool top_only; /**< Flag to specify current in top rings only */
-  bool bot_only; /**< Flag to specify current in bottom rings only */
+  bool top_only;         /**< Flag to specify current in top rings only */
+  bool bot_only;         /**< Flag to specify current in bottom rings only */
   bool variable_current; /**< Flag to specify variable current in each ring */
 
   bool evaluate_magnetic_field;
diff --git a/src/loMach.cpp b/src/loMach.cpp
index 9f95926f..455c33a3 100644
--- a/src/loMach.cpp
+++ b/src/loMach.cpp
@@ -166,11 +166,12 @@ void LoMachSolver::initialize() {
   if (loMach_opts_.thermo_solver == "constant-property") {
     thermo_ = new ConstantPropertyThermoChem(pmesh_, loMach_opts_.order, tpsP_);
   } else if (loMach_opts_.thermo_solver == "calorically-perfect") {
-    thermo_ = new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
+    thermo_ =
+        new CaloricallyPerfectThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_);
   } else if (loMach_opts_.thermo_solver == "lte-thermo-chem") {
     thermo_ = new LteThermoChem(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
   } else if (loMach_opts_.thermo_solver == "reacting-flow") {
-    thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, tpsP_);
+    thermo_ = new ReactingFlow(pmesh_, &loMach_opts_, temporal_coeff_, (meshData_->getGridScale()), tpsP_);
   } else {
     // Unknown choice... die
     if (rank0_) {
@@ -182,10 +183,11 @@ void LoMachSolver::initialize() {
   // Instantiate flow solver
   if (loMach_opts_.flow_solver == "zero-flow") {
     // No flow---set u = 0.  Primarily useful for testing thermochem models in isolation
-    flow_ = new ZeroFlow(pmesh_, 1);
+    flow_ = new ZeroFlow(pmesh_, 1, tpsP_);
   } else if (loMach_opts_.flow_solver == "tomboulides") {
     // Tomboulides flow solver
-    flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_, tpsP_);
+    flow_ = new Tomboulides(pmesh_, loMach_opts_.order, loMach_opts_.order, temporal_coeff_,
+                            (meshData_->getGridScale()), tpsP_);
   } else {
     // Unknown choice... die
     if (rank0_) {
diff --git a/src/lte_thermo_chem.cpp b/src/lte_thermo_chem.cpp
index c41a48eb..19d0c161 100644
--- a/src/lte_thermo_chem.cpp
+++ b/src/lte_thermo_chem.cpp
@@ -453,7 +453,7 @@ void LteThermoChem::initializeSelf() {
 
   // outlet bc
   {
-      // Assumed homogeneous Nuemann on T, so nothing to do
+    // Assumed homogeneous Nuemann on T, so nothing to do
   }
 
   // Wall BCs
diff --git a/src/quasimagnetostatic.cpp b/src/quasimagnetostatic.cpp
index 24d906f0..30288a9b 100644
--- a/src/quasimagnetostatic.cpp
+++ b/src/quasimagnetostatic.cpp
@@ -272,9 +272,7 @@ void QuasiMagnetostaticSolver3D::InitializeCurrent() {
   }
 
   if (rank0_) {
-    std::cout << "J0 = " << J0(0) << ", " << J0(1)
-              << ", " << J0(2) << ", " << J0(3)
-              << ", " << J0(4) << endl;
+    std::cout << "J0 = " << J0(0) << ", " << J0(1) << ", " << J0(2) << ", " << J0(3) << ", " << J0(4) << endl;
   }
 
   PWConstCoefficient J0coef(J0);
@@ -924,9 +922,7 @@ void QuasiMagnetostaticSolverAxiSym::InitializeCurrent() {
   }
 
   if (rank0_) {
-    std::cout << "J0 = " << J0(0) << ", " << J0(1)
-              << ", " << J0(2) << ", " << J0(3)
-              << ", " << J0(4) << endl;
+    std::cout << "J0 = " << J0(0) << ", " << J0(1) << ", " << J0(2) << ", " << J0(3) << ", " << J0(4) << endl;
   }
 
   FunctionCoefficient radius_coeff(radius);
diff --git a/src/reactingFlow.cpp b/src/reactingFlow.cpp
index a5ea063a..bceacda4 100644
--- a/src/reactingFlow.cpp
+++ b/src/reactingFlow.cpp
@@ -55,10 +55,11 @@ static double radius(const Vector &pos) { return pos[0]; }
 static FunctionCoefficient radius_coeff(radius);
 
 ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &time_coeff,
-                           TPS::Tps *tps)
+                           ParGridFunction *gridScale, TPS::Tps *tps)
     : tpsP_(tps), pmesh_(pmesh), dim_(pmesh->Dimension()), time_coeff_(time_coeff) {
   rank0_ = (pmesh_->GetMyRank() == 0);
   order_ = loMach_opts->order;
+  gridScale_gf_ = gridScale;
 
   /// Basic input information
   tpsP_->getInput("loMach/ambientPressure", ambient_pressure_, 101325.0);
@@ -538,6 +539,11 @@ ReactingFlow::ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, tem
     operator_split_ = true;
     assert(nSub_ == 1);
   }
+
+  // artificial diffusion (SUPG)
+  tpsP_->getInput("loMach/reactingFlow/streamwise-stabilization", sw_stab_, false);
+  // NOTE: clipping Qt for debugging
+  // tpsP_->getInput("loMach/reactingFlow/clip-qt", clip_qt_, false);
 }
 
 ReactingFlow::~ReactingFlow() {
@@ -618,6 +624,21 @@ ReactingFlow::~ReactingFlow() {
   delete transport_;
   delete chemistry_;
 
+  delete umag_coeff_;
+  delete gscale_coeff_;
+  delete visc_coeff_;
+  delete visc_inv_coeff_;
+  delete reh1_coeff_;
+  delete reh2_coeff_;
+  delete Reh_coeff_;
+  delete csupg_coeff_;
+  delete uw1_coeff_;
+  delete uw2_coeff_;
+  delete upwind_coeff_;
+  delete swdiff_coeff_;
+  delete supg_coeff_;
+  delete supg_cp_coeff_;
+
   // allocated in initializeSelf
   delete vfes_;
   delete vfec_;
@@ -826,6 +847,22 @@ void ReactingFlow::initializeSelf() {
   Mmix_gf_.SetSpace(sfes_);
   Rmix_gf_.SetSpace(sfes_);
 
+  // qt rhs visualization
+  // rhsqt_bd_.SetSpace(sfes_);
+  // rhsqt_fo_.SetSpace(sfes_);
+  // rhsqt_jh_.SetSpace(sfes_);
+  // rhsqt_hf_.SetSpace(sfes_);
+  // rhsqt_sd_.SetSpace(sfes_);
+  // rhsqt_total_.SetSpace(sfes_);
+  // Xqt_gf_.SetSpace(sfes_);
+  // rhsqt_bd_ = 0.0;
+  // rhsqt_fo_ = 0.0;
+  // rhsqt_jh_ = 0.0;
+  // rhsqt_hf_ = 0.0;
+  // rhsqt_sd_ = 0.0;
+  // rhsqt_total_ = 0.0;
+  // Xqt_gf_ = 0.0;
+
   // exports
   toFlow_interface_.density = &rn_gf_;
   toFlow_interface_.viscosity = &visc_gf_;
@@ -1112,6 +1149,50 @@ void ReactingFlow::initializeOperators() {
     rad_kap_gradT_coeff_ = new ScalarVectorProductCoefficient(radius_coeff, *kap_gradT_coeff_);
   }
 
+  // artifical diffusion coefficients
+  if (sw_stab_) {
+    // run once to avoid division by zero
+    // updateMixture();
+    // updateDensity(1.0);
+    updateDiffusivity();
+    umag_coeff_ = new VectorMagnitudeCoefficient(*un_next_coeff_);
+    gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_);
+    visc_coeff_ = new GridFunctionCoefficient(&visc_gf_);
+    visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.);
+
+    // compute Reh
+    reh1_coeff_ = new ProductCoefficient(*rhon_next_coeff_, *visc_inv_coeff_);
+    reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_);
+    Reh_coeff_ = new ProductCoefficient(*reh2_coeff_, *umag_coeff_);
+
+    // Csupg
+    csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor);
+
+    if (axisym_) {
+      // compute upwind magnitude
+      uw1_coeff_ = new ProductCoefficient(*rad_rho_coeff_, *csupg_coeff_);
+      uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+      upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+      // streamwise diffusion direction
+      swdiff_coeff_ = new TransformedMatrixVectorCoefficient(un_next_coeff_, &streamwiseTensor);
+
+      supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+      supg_cp_coeff_ = new ScalarMatrixProductCoefficient(*cpMix_coeff_, *supg_coeff_);
+    } else {
+      // compute upwind magnitude
+      uw1_coeff_ = new ProductCoefficient(*rhon_next_coeff_, *csupg_coeff_);
+      uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+      upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+      // streamwise diffusion direction
+      swdiff_coeff_ = new TransformedMatrixVectorCoefficient(un_next_coeff_, &streamwiseTensor);
+
+      supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+      supg_cp_coeff_ = new ScalarMatrixProductCoefficient(*cpMix_coeff_, *supg_coeff_);
+    }
+  }
+
   At_form_ = new ParBilinearForm(sfes_);
   ConvectionIntegrator *at_blfi;
   if (axisym_) {
@@ -1218,6 +1299,17 @@ void ReactingFlow::initializeOperators() {
     hmt_blfi->SetIntRule(&ir_di);
     hdt_blfi->SetIntRule(&ir_di);
   }
+  // SUPG
+  DiffusionIntegrator *sdt_blfi;
+  if (sw_stab_) {
+    // auto *sdt_blfi = new DiffusionIntegrator(*supg_coeff_);
+    sdt_blfi = new DiffusionIntegrator(*supg_cp_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   sdt_blfi->SetIntRule(&ir_di);
+    // }
+    Ht_form_->AddDomainIntegrator(sdt_blfi);
+  }
   Ht_form_->AddDomainIntegrator(hmt_blfi);
   Ht_form_->AddDomainIntegrator(hdt_blfi);
   if (partial_assembly_) {
@@ -1244,6 +1336,17 @@ void ReactingFlow::initializeOperators() {
     hmy_blfi->SetIntRule(&ir_di);
     hdy_blfi->SetIntRule(&ir_di);
   }
+  // SUPG
+  DiffusionIntegrator *syt_blfi;
+  if (sw_stab_) {
+    // auto *syt_blfi = new DiffusionIntegrator(*supg_coeff_);
+    syt_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // SUPG diffusion
+    // if (numerical_integ_) {
+    //   sdt_blfi->SetIntRule(&ir_di);
+    // }
+    Hy_form_->AddDomainIntegrator(syt_blfi);
+  }
   Hy_form_->AddDomainIntegrator(hmy_blfi);
   Hy_form_->AddDomainIntegrator(hdy_blfi);
   if (partial_assembly_) {
@@ -1401,6 +1504,15 @@ void ReactingFlow::initializeOperators() {
     lqd_blfi->SetIntRule(&ir_di);
   }
   LQ_form_->AddDomainIntegrator(lqd_blfi);
+
+  DiffusionIntegrator *slqd_blfi;
+  if (sw_stab_) {
+    slqd_blfi = new DiffusionIntegrator(*supg_coeff_);
+    // if (numerical_integ_) {
+    //   slqd_blfi->SetIntRule(&ir_di);
+    // }
+    LQ_form_->AddDomainIntegrator(slqd_blfi);
+  }
   if (partial_assembly_) {
     LQ_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
   }
@@ -1427,6 +1539,15 @@ void ReactingFlow::initializeOperators() {
     lyd_blfi->SetIntRule(&ir_di);
   }
   LY_form_->AddDomainIntegrator(lyd_blfi);
+
+  // auto *slyd_blfi = new DiffusionIntegrator(*supg_coeff_);
+  // if (sw_stab_)
+  //   // SUPG diffusion
+  //   // if (numerical_integ_) {
+  //   //   slyd_blfi->SetIntRule(&ir_di);
+  //   // }
+  //   LY_form_->AddDomainIntegrator(slqd_blfi);
+
   if (partial_assembly_) {
     LY_form_->SetAssemblyLevel(AssemblyLevel::PARTIAL);
   }
@@ -2269,6 +2390,15 @@ void ReactingFlow::initializeViz(ParaViewDataCollection &pvdc) {
   pvdc.RegisterField("weff", &weff_gf_);
   pvdc.RegisterField("emission", &emission_gf_);
 
+  // diagnose Qt issues, rhs contributions
+  // pvdc.RegisterField("rhsqt_bd", &rhsqt_bd_); //boundary terms
+  // pvdc.RegisterField("rhsqt_fo", &rhsqt_fo_); //bilinear form
+  // pvdc.RegisterField("rhsqt_jh", &rhsqt_jh_); //joule heating
+  // pvdc.RegisterField("rhsqt_hf", &rhsqt_hf_); //heat of formation
+  // pvdc.RegisterField("rhsqt_sd", &rhsqt_sd_); ///species-temp diff
+  // pvdc.RegisterField("rhsqt_total", &rhsqt_total_);
+  // pvdc.RegisterField("Xqt", &Xqt_gf_);
+
   vizSpecFields_.clear();
   vizSpecNames_.clear();
   for (int sp = 0; sp < nSpecies_; sp++) {
@@ -2716,43 +2846,81 @@ void ReactingFlow::AddQtDirichletBC(ScalarFuncT *f, Array<int> &attr) {
 
 void ReactingFlow::computeQtTO() {
   tmpR0_ = 0.0;
+  // rhsqt_bd_ = 0.0;
+  // rhsqt_fo_ = 0.0;
+  // rhsqt_jh_ = 0.0;
+  // rhsqt_hf_ = 0.0;
+  // rhsqt_sd_ = 0.0;
+  // rhsqt_total_ = 0.0;
+  // Xqt_gf_ = 0.0;
   LQ_bdry_->Update();
   LQ_bdry_->Assemble();
   LQ_bdry_->ParallelAssemble(tmpR0_);
   tmpR0_.Neg();
+  // LQ_bdry_->ParallelAssemble(rhsqt_bd_);
+  // rhsqt_bd_.Neg();
+  // tmpR0_ += rhsqt_bd_;
+  // printf("%f\n", tmpR0_.Norml2());
 
   Array<int> empty;
   LQ_form_->Update();
   LQ_form_->Assemble();
   LQ_form_->FormSystemMatrix(empty, LQ_);
   LQ_->AddMult(Tn_next_, tmpR0_);  // tmpR0_ += LQ{Tn_next}
+  // LQ_->Mult(Tn_next_, rhsqt_fo_);  // tmpR0_ += LQ{Tn_next}
+  // tmpR0_ += rhsqt_fo_;
+  // printf("%f\n", tmpR0_.Norml2());
 
   // Joule heating (and radiation sink)
   jh_form_->Update();
   jh_form_->Assemble();
   jh_form_->ParallelAssemble(jh_);
   tmpR0_ -= jh_;
+  // rhsqt_jh_.SetFromTrueDofs(jh_);
+  // rhsqt_jh_.Neg();
+  // printf("%f\n", tmpR0_.Norml2());
 
   // heat of formation
   Ms_->AddMult(hw_, tmpR0_, -1.0);
+  // Ms_->Mult(hw_, rhsqt_hf_);
+  // rhsqt_hf_.Neg();
+  // tmpR0_ += rhsqt_hf_;
+  // printf("%f\n", tmpR0_.Norml2());
 
   // species-temp diffusion term, already in int-weak form
   tmpR0_.Add(-1.0, crossDiff_);
+  // rhsqt_sd_.SetFromTrueDofs(crossDiff_);
+  // rhsqt_sd_.Neg();
+  // printf("%f\n", tmpR0_.Norml2());
 
   sfes_->GetRestrictionMatrix()->MultTranspose(tmpR0_, resT_gf_);
+  // rhsqt_total_.SetFromTrueDofs(resT_gf_);
 
   Qt_ = 0.0;
   Qt_gf_.SetFromTrueDofs(Qt_);
 
   Vector Xqt, Bqt;
   Mq_form_->FormLinearSystem(Qt_ess_tdof_, Qt_gf_, resT_gf_, Mq_, Xqt, Bqt, 1);
-
   MqInv_->Mult(Bqt, Xqt);
   Mq_form_->RecoverFEMSolution(Xqt, resT_gf_, Qt_gf_);
+  // Xqt_gf_.SetFromTrueDofs(Xqt);
 
   Qt_gf_ *= Rmix_gf_;
   Qt_gf_ /= CpMix_gf_;
   Qt_gf_ /= thermo_pressure_;
+
+  // NOTE: clipping for diagnosis, goes both directions for now
+  // double clip_val;
+  // if (clip_qt_) {
+  //   tpsP_->getInput("loMach/reactingFlow/clip-qt-value", clip_val, 50000.);
+
+  //   for (int n = 0; n < sfes_->GetNDofs(); n++) {
+  //     // Qt_gf_[n] = std::clamp(Qt_gf_[n], -clip_val, clip_val)
+  //     Qt_gf_[n] = std::min(clip_val, std::max(-clip_val, Qt_gf_[n]));
+  //   }
+
+  // }
+
   Qt_gf_.Neg();
 }
 
diff --git a/src/reactingFlow.hpp b/src/reactingFlow.hpp
index 3ea4fa4c..e9e26d7f 100644
--- a/src/reactingFlow.hpp
+++ b/src/reactingFlow.hpp
@@ -153,6 +153,12 @@ class ReactingFlow : public ThermoChemModelBase {
   // Initial temperature value (if constant IC)
   double T_ic_;
 
+  // streamwise-stabilization
+  bool sw_stab_;
+
+  // clip qt for diagnosis
+  // bool clip_qt_;
+
   // FEM related fields and objects
 
   // Scalar \f$H^1\f$ finite element collection.
@@ -173,6 +179,8 @@ class ReactingFlow : public ThermoChemModelBase {
   // Vector \f$H^1\f$ finite element space.
   ParFiniteElementSpace *vfes_ = nullptr;
 
+  ParGridFunction *gridScale_gf_ = nullptr;
+
   // Fields
   ParGridFunction Tnm1_gf_, Tnm2_gf_;
   ParGridFunction Tn_gf_, Tn_next_gf_, Text_gf_, resT_gf_;
@@ -203,6 +211,15 @@ class ReactingFlow : public ThermoChemModelBase {
   ParGridFunction sigma_gf_;
   ParGridFunction jh_gf_;
 
+  // viz for qt rhs
+  // ParGridFunction rhsqt_bd_;
+  // ParGridFunction rhsqt_fo_;
+  // ParGridFunction rhsqt_jh_;
+  // ParGridFunction rhsqt_hf_;
+  // ParGridFunction rhsqt_sd_;
+  // ParGridFunction rhsqt_total_;
+  // ParGridFunction Xqt_gf_;
+
   // ParGridFunction *buffer_tInlet_ = nullptr;
   GridFunctionCoefficient *temperature_bc_field_ = nullptr;
 
@@ -245,6 +262,21 @@ class ReactingFlow : public ThermoChemModelBase {
   ProductCoefficient *rad_radiation_sink_coeff_ = nullptr;
   ScalarVectorProductCoefficient *rad_kap_gradT_coeff_ = nullptr;
 
+  VectorMagnitudeCoefficient *umag_coeff_ = nullptr;
+  GridFunctionCoefficient *gscale_coeff_ = nullptr;
+  GridFunctionCoefficient *visc_coeff_ = nullptr;
+  PowerCoefficient *visc_inv_coeff_ = nullptr;
+  ProductCoefficient *reh1_coeff_ = nullptr;
+  ProductCoefficient *reh2_coeff_ = nullptr;
+  ProductCoefficient *Reh_coeff_ = nullptr;
+  TransformedCoefficient *csupg_coeff_ = nullptr;
+  ProductCoefficient *uw1_coeff_ = nullptr;
+  ProductCoefficient *uw2_coeff_ = nullptr;
+  ProductCoefficient *upwind_coeff_ = nullptr;
+  TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr;
+  ScalarMatrixProductCoefficient *supg_coeff_ = nullptr;
+  ScalarMatrixProductCoefficient *supg_cp_coeff_ = nullptr;
+
   // operators and solvers
   ParBilinearForm *At_form_ = nullptr;
   ParBilinearForm *Ay_form_ = nullptr;
@@ -356,7 +388,8 @@ class ReactingFlow : public ThermoChemModelBase {
   std::vector<std::string> vizSpecNames_;
 
  public:
-  ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff, TPS::Tps *tps);
+  ReactingFlow(mfem::ParMesh *pmesh, LoMachOptions *loMach_opts, temporalSchemeCoefficients &timeCoeff,
+               ParGridFunction *gridScale, TPS::Tps *tps);
   virtual ~ReactingFlow();
 
   // Functions overriden from base class
diff --git a/src/split_flow_base.cpp b/src/split_flow_base.cpp
index 0094f922..e39c5edf 100644
--- a/src/split_flow_base.cpp
+++ b/src/split_flow_base.cpp
@@ -32,12 +32,19 @@
 
 #include "split_flow_base.hpp"
 
+#include "tps.hpp"
+
 using namespace mfem;
 
-ZeroFlow::ZeroFlow(mfem::ParMesh *pmesh, int vorder) : pmesh_(pmesh), vorder_(vorder), dim_(pmesh->Dimension()) {}
+ZeroFlow::ZeroFlow(mfem::ParMesh *pmesh, int vorder, TPS::Tps *tps)
+    : pmesh_(pmesh), vorder_(vorder), dim_(pmesh->Dimension()), tpsP_(tps) {
+  // nonzero flow option
+  tpsP_->getInput("loMach/zeroflow/nonzero-flow", nonzero_flow_, false);
+}
 
 ZeroFlow::~ZeroFlow() {
   delete velocity_;
+  delete zero_;
   delete fes_;
   delete fec_;
 }
@@ -46,13 +53,24 @@ void ZeroFlow::initializeSelf() {
   fec_ = new H1_FECollection(vorder_, dim_);
   fes_ = new ParFiniteElementSpace(pmesh_, fec_, dim_);
   velocity_ = new ParGridFunction(fes_);
+  zero_ = new ParGridFunction(fes_);
   *velocity_ = 0.0;
+  *zero_ = 0.0;
+  // set background velocity if nonzero
+  if (nonzero_flow_) {
+    Vector zero(dim_);
+    Vector velocity_value(dim_);
+    zero = 0.0;
+    tpsP_->getVec("loMach/zeroflow/nonzero-vel", velocity_value, dim_, zero);
+    VectorConstantCoefficient ub_coeff(velocity_value);
+    velocity_->ProjectCoefficient(ub_coeff);
+  }
 
   toThermoChem_interface_.velocity = velocity_;
   toThermoChem_interface_.swirl_supported = false;
 
-  // no need to create additional zero vectors
-  toTurbModel_interface_.gradU = velocity_;
-  toTurbModel_interface_.gradV = velocity_;
-  toTurbModel_interface_.gradW = velocity_;
+  // need to create additional zero vectors, in case of nonzero vel
+  toTurbModel_interface_.gradU = zero_;
+  toTurbModel_interface_.gradV = zero_;
+  toTurbModel_interface_.gradW = zero_;
 }
diff --git a/src/split_flow_base.hpp b/src/split_flow_base.hpp
index 9b57c4d9..e09999ed 100644
--- a/src/split_flow_base.hpp
+++ b/src/split_flow_base.hpp
@@ -35,6 +35,11 @@
  * @brief Contains base class and simplest possible variant of flow model
  */
 
+// forward-declaration for Tps support class (tps.hpp)
+namespace TPS {
+class Tps;
+}
+
 #include "tps_mfem_wrap.hpp"
 
 class IODataOrganizer;
@@ -143,17 +148,24 @@ class FlowBase {
 
 class ZeroFlow final : public FlowBase {
  protected:
+  // Options
+  bool nonzero_flow_;
+
   mfem::ParMesh *pmesh_;
   const int vorder_;
   const int dim_;
 
+  // Options-related structures
+  TPS::Tps *tpsP_ = nullptr;
+
   mfem::FiniteElementCollection *fec_ = nullptr;
   mfem::ParFiniteElementSpace *fes_ = nullptr;
   mfem::ParGridFunction *velocity_ = nullptr;
+  mfem::ParGridFunction *zero_ = nullptr;
 
  public:
   /// Constructor
-  ZeroFlow(mfem::ParMesh *pmesh, int vorder);
+  ZeroFlow(mfem::ParMesh *pmesh, int vorder, TPS::Tps *tps = nullptr);
 
   /// Destructor
   ~ZeroFlow() final;
diff --git a/src/tomboulides.cpp b/src/tomboulides.cpp
index 184df27b..53f9206a 100644
--- a/src/tomboulides.cpp
+++ b/src/tomboulides.cpp
@@ -67,7 +67,8 @@ void Orthogonalize(Vector &v, const ParFiniteElementSpace *pfes) {
   v -= global_sum / static_cast<double>(global_size);
 }
 
-Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, TPS::Tps *tps)
+Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff,
+                         mfem::ParGridFunction *gridScale, TPS::Tps *tps)
     : gll_rules(0, Quadrature1D::GaussLobatto),
       tpsP_(tps),
       pmesh_(pmesh),
@@ -80,6 +81,7 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS
   rank0_ = (pmesh_->GetMyRank() == 0);
   axisym_ = false;
   nvel_ = dim_;
+  gridScale_gf_ = gridScale;
 
   // make sure there is room for BC attributes
   if (!(pmesh_->bdr_attributes.Size() == 0)) {
@@ -128,6 +130,9 @@ Tomboulides::Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalS
     tps->getInput("loMach/tomboulides/psolve-maxIters", pressure_solve_max_iter_, default_max_iter_);
     tps->getInput("loMach/tomboulides/hsolve-maxIters", hsolve_max_iter_, default_max_iter_);
     tps->getInput("loMach/tomboulides/msolve-maxIters", mass_inverse_max_iter_, default_max_iter_);
+
+    // artificial diffusion (SUPG)
+    tpsP_->getInput("loMach/tomboulides/streamwise-stabilization", sw_stab_, false);
   }
 }
 
@@ -169,6 +174,8 @@ Tomboulides::~Tomboulides() {
   delete L_iorho_lor_;
   delete L_iorho_form_;
 
+  delete Mv_stab_form_;
+
   delete swirl_var_viscosity_coeff_;
   delete utheta_vec_coeff_;
   delete rho_ur_ut_coeff_;
@@ -204,6 +211,20 @@ Tomboulides::~Tomboulides() {
   delete iorho_coeff_;
   delete rho_coeff_;
 
+  delete umag_coeff_;
+  delete gscale_coeff_;
+  delete visc_inv_coeff_;
+  delete reh1_coeff_;
+  delete reh2_coeff_;
+  delete Reh_coeff_;
+  delete csupg_coeff_;
+  delete uw1_coeff_;
+  delete uw2_coeff_;
+  delete upwind_coeff_;
+  delete swdiff_coeff_;
+  delete supg_coeff_;
+  delete visc_coeff_;
+
   // objects allocated by initalizeSelf
   if (axisym_) delete gravity_vec_;
   delete utheta_next_gf_;
@@ -343,6 +364,8 @@ void Tomboulides::initializeSelf() {
 
   tmpR0_.SetSize(sfes_truevsize);
   tmpR1_.SetSize(vfes_truevsize);
+  swDiff_vec_.SetSize(vfes_truevsize);
+  tmpR0b_.SetSize(pfes_truevsize);
 
   gradU_.SetSize(vfes_truevsize);
   gradV_.SetSize(vfes_truevsize);
@@ -626,6 +649,8 @@ void Tomboulides::initializeOperators() {
   S_poisson_coeff_ = new VectorSumCoefficient(*twoS_gradmu_coeff_, *gradmu_Qt_coeff_, 1.0, -2. / 3);
   S_mom_coeff_ = new VectorSumCoefficient(*graduT_gradmu_coeff_, *gradmu_Qt_coeff_, 1.0, -1.0);
 
+  u_next_coeff_ = new VectorGridFunctionCoefficient(u_next_gf_);
+
   // Coefficients for axisymmetric
   if (axisym_) {
     rad_rho_coeff_ = new ProductCoefficient(radius_coeff, *rho_coeff_);
@@ -646,7 +671,6 @@ void Tomboulides::initializeOperators() {
     ur_conv_forcing_coeff_ = new VectorArrayCoefficient(2);
     ur_conv_forcing_coeff_->Set(0, utheta2_coeff_);
 
-    u_next_coeff_ = new VectorGridFunctionCoefficient(u_next_gf_);
     rad_rhou_coeff_ = new ScalarVectorProductCoefficient(*rad_rho_coeff_, *u_next_coeff_);
 
     u_next_rad_coeff_ = new GridFunctionCoefficient(u_next_rad_comp_gf_);
@@ -661,6 +685,47 @@ void Tomboulides::initializeOperators() {
     swirl_var_viscosity_coeff_ = new InnerProductCoefficient(*grad_mu_coeff_, *utheta_vec_coeff_);
   }
 
+  // artifical diffusion coefficients
+  if (sw_stab_) {
+    visc_coeff_ = new GridFunctionCoefficient(thermo_interface_->viscosity);
+    umag_coeff_ = new VectorMagnitudeCoefficient(*u_next_coeff_);
+    gscale_coeff_ = new GridFunctionCoefficient(gridScale_gf_);
+
+    visc_inv_coeff_ = new PowerCoefficient(*visc_coeff_, -1.0);
+    // visc_inv_coeff_ = new PowerCoefficient(*mu_coeff_, -1.0);
+
+    // compute Reh
+    reh1_coeff_ = new ProductCoefficient(*rho_coeff_, *visc_inv_coeff_);
+    reh2_coeff_ = new ProductCoefficient(*reh1_coeff_, *gscale_coeff_);
+    Reh_coeff_ = new ProductCoefficient(*reh2_coeff_, *umag_coeff_);
+
+    // Csupg
+    csupg_coeff_ = new TransformedCoefficient(Reh_coeff_, csupgFactor);
+
+    if (axisym_) {
+      // compute upwind magnitude
+      uw1_coeff_ = new ProductCoefficient(*rad_rho_coeff_, *csupg_coeff_);
+      uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+      upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+      // streamwise diffusion direction
+      swdiff_coeff_ = new TransformedMatrixVectorCoefficient(u_next_coeff_, &streamwiseTensor);
+
+      supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+    } else {
+      // compute upwind magnitude
+      // dividing by rho anyway
+      // uw1_coeff_ = new ProductCoefficient(*rho_coeff_, *csupg_coeff_);
+      uw2_coeff_ = new ProductCoefficient(*uw1_coeff_, *gscale_coeff_);
+      upwind_coeff_ = new ProductCoefficient(*uw2_coeff_, *umag_coeff_);
+
+      // streamwise diffusion direction
+      swdiff_coeff_ = new TransformedMatrixVectorCoefficient(u_next_coeff_, &streamwiseTensor);
+
+      supg_coeff_ = new ScalarMatrixProductCoefficient(*upwind_coeff_, *swdiff_coeff_);
+    }
+  }
+
   // Integration rules (only used if numerical_integ_ is true).  When
   // this is the case, the quadrature degree set such that the
   // Gauss-Lobatto quad pts correspond to the Gauss-Lobatto nodes.
@@ -688,6 +753,18 @@ void Tomboulides::initializeOperators() {
     L_iorho_blfi->SetIntRule(&ir_ni_p);
   }
   L_iorho_form_->AddDomainIntegrator(L_iorho_blfi);
+
+  // DiffusionIntegrator *slio_blfi;
+  // if (sw_stab_) {
+  //   // auto *slio_blfi = new DiffusionIntegrator(*supg_coeff_);
+  //   slio_blfi = new DiffusionIntegrator(*supg_coeff_);
+  //   // SUPG diffusion
+  //   // if (numerical_integ_) {
+  //   //   slio_blfi->SetIntRule(&ir_di);
+  //   // }
+  //   L_iorho_form_->AddDomainIntegrator(slio_blfi);
+  // }
+
   L_iorho_form_->Assemble();
   L_iorho_form_->FormSystemMatrix(pres_ess_tdof_, L_iorho_op_);
 
@@ -809,6 +886,19 @@ void Tomboulides::initializeOperators() {
   Mv_rho_form_->Assemble();
   Mv_rho_form_->FormSystemMatrix(empty, Mv_rho_op_);
 
+  // Vector mass matrix for streamwise stability of the velocity gradients
+  if (sw_stab_) {
+    Mv_stab_form_ = new ParBilinearForm(vfes_);
+    VectorMassIntegrator *mvs_blfi;
+    mvs_blfi = new VectorMassIntegrator(*supg_coeff_);
+    if (numerical_integ_) {
+      mvs_blfi->SetIntRule(&ir_ni_v);
+    }
+    Mv_stab_form_->AddDomainIntegrator(mvs_blfi);
+    Mv_stab_form_->Assemble();
+    Mv_stab_form_->FormSystemMatrix(empty, Mv_stab_op_);
+  }
+
   // Inverse (unweighted) mass operator (velocity space)
   if (partial_assembly_) {
     Vector diag_pa(vfes_->GetTrueVSize());
@@ -886,6 +976,7 @@ void Tomboulides::initializeOperators() {
   Hv_form_ = new ParBilinearForm(vfes_);
   VectorMassIntegrator *hmv_blfi;
   VectorDiffusionIntegrator *hdv_blfi;
+
   if (axisym_) {
     hmv_blfi = new VectorMassIntegrator(*rad_rho_over_dt_coeff_);
     hdv_blfi = new VectorDiffusionIntegrator(*rad_mu_coeff_);
@@ -899,6 +990,17 @@ void Tomboulides::initializeOperators() {
   }
   Hv_form_->AddDomainIntegrator(hmv_blfi);
   Hv_form_->AddDomainIntegrator(hdv_blfi);
+
+  VectorDiffusionIntegrator *shdv_blfi;
+  if (sw_stab_) {
+    // auto *shdv_blfi = new VectorDiffusionIntegrator(*supg_coeff_);
+    shdv_blfi = new VectorDiffusionIntegrator(*supg_coeff_);
+    // if (numerical_integ_) {
+    //   shdv_blfi->SetIntRule(&ir_di);
+    // }
+    Hv_form_->AddDomainIntegrator(shdv_blfi);
+  }
+
   if (axisym_) {
     auto *hfv_blfi = new VectorMassIntegrator(*visc_forcing_coeff_);
     Hv_form_->AddDomainIntegrator(hfv_blfi);
@@ -996,6 +1098,13 @@ void Tomboulides::initializeOperators() {
     auto *hds_blfi = new DiffusionIntegrator(*rad_mu_coeff_);
     auto *hfs_blfi = new MassIntegrator(*mu_over_rad_coeff_);
 
+    DiffusionIntegrator *shds_blfi;
+    if (sw_stab_) {
+      // auto *shds_blfi = new DiffusionIntegrator(*supg_coeff_);
+      shds_blfi = new DiffusionIntegrator(*supg_coeff_);
+      Hs_form_->AddDomainIntegrator(shds_blfi);
+    }
+
     Hs_form_->AddDomainIntegrator(hms_blfi);
     Hs_form_->AddDomainIntegrator(hds_blfi);
     Hs_form_->AddDomainIntegrator(hfs_blfi);
@@ -1434,6 +1543,28 @@ void Tomboulides::step() {
     });
   }
 
+  if (sw_stab_) {
+    // Update matrix
+    Array<int> empty;
+    Mv_stab_form_->Update();
+    Mv_stab_form_->Assemble();
+    Mv_stab_form_->FormSystemMatrix(empty, Mv_stab_op_);
+
+    Mv_stab_op_->Mult(gradU_, tmpR1_);
+    D_op_->Mult(tmpR1_, tmpR0b_);
+    setVectorFromScalar(tmpR0b_, 0, &swDiff_vec_);
+    Mv_stab_op_->Mult(gradV_, tmpR1_);
+    D_op_->Mult(tmpR1_, tmpR0b_);
+    setVectorFromScalar(tmpR0b_, 1, &swDiff_vec_);
+    if (dim_ == 3) {
+      Mv_stab_op_->Mult(gradW_, tmpR1_);
+      D_op_->Mult(tmpR1_, tmpR0b_);
+      setVectorFromScalar(tmpR0b_, 2, &swDiff_vec_);
+    }
+    Mv_rho_inv_->Mult(swDiff_vec_, tmpR1_);
+    pp_div_vec_ += tmpR1_;
+  }
+  // printf("%f\n", tmpR1_.Norml2());
   // Add ustar/dt contribution
   pp_div_vec_ += ustar_vec_;
 
diff --git a/src/tomboulides.hpp b/src/tomboulides.hpp
index 1463b426..a508115d 100644
--- a/src/tomboulides.hpp
+++ b/src/tomboulides.hpp
@@ -181,6 +181,9 @@ class Tomboulides final : public FlowBase {
   // reference here.
   const temporalSchemeCoefficients &coeff_;
 
+  // streamwise-stabilization
+  bool sw_stab_;
+
   StopWatch sw_press_, sw_helm_;
 
   std::string ic_string_;
@@ -234,6 +237,8 @@ class Tomboulides final : public FlowBase {
   mfem::ParGridFunction *resp_gf_ = nullptr;
   mfem::ParGridFunction *pp_div_rad_comp_gf_ = nullptr;
 
+  mfem::ParGridFunction *gridScale_gf_ = nullptr;
+
   /// Swirl
   mfem::ParGridFunction *utheta_gf_ = nullptr;
   mfem::ParGridFunction *utheta_next_gf_ = nullptr;
@@ -286,6 +291,20 @@ class Tomboulides final : public FlowBase {
   mfem::VectorArrayCoefficient *utheta_vec_coeff_ = nullptr;
   mfem::InnerProductCoefficient *swirl_var_viscosity_coeff_ = nullptr;
 
+  mfem::VectorMagnitudeCoefficient *umag_coeff_ = nullptr;
+  mfem::GridFunctionCoefficient *gscale_coeff_ = nullptr;
+  mfem::PowerCoefficient *visc_inv_coeff_ = nullptr;
+  mfem::ProductCoefficient *reh1_coeff_ = nullptr;
+  mfem::ProductCoefficient *reh2_coeff_ = nullptr;
+  mfem::ProductCoefficient *Reh_coeff_ = nullptr;
+  mfem::TransformedCoefficient *csupg_coeff_ = nullptr;
+  mfem::ProductCoefficient *uw1_coeff_ = nullptr;
+  mfem::ProductCoefficient *uw2_coeff_ = nullptr;
+  mfem::ProductCoefficient *upwind_coeff_ = nullptr;
+  mfem::TransformedMatrixVectorCoefficient *swdiff_coeff_ = nullptr;
+  mfem::ScalarMatrixProductCoefficient *supg_coeff_ = nullptr;
+  mfem::GridFunctionCoefficient *visc_coeff_ = nullptr;
+
   // mfem "form" objects used to create operators
   mfem::ParBilinearForm *L_iorho_form_ = nullptr;  // \int (1/\rho) \nabla \phi_i \cdot \nabla \phi_j
   mfem::ParLinearForm *forcing_form_ = nullptr;    // \int \phi_i f
@@ -309,6 +328,9 @@ class Tomboulides final : public FlowBase {
   mfem::ParLinearForm *rho_ur_ut_form_ = nullptr;
   mfem::ParLinearForm *swirl_var_viscosity_form_ = nullptr;
 
+  // streamwise stability
+  mfem::ParBilinearForm *Mv_stab_form_ = nullptr;
+
   // mfem operator objects
   mfem::OperatorHandle L_iorho_op_;
   mfem::OperatorHandle Ms_op_;
@@ -321,6 +343,9 @@ class Tomboulides final : public FlowBase {
   mfem::OperatorHandle Hs_op_;
   mfem::OperatorHandle As_op_;
 
+  // streamwise stability
+  mfem::OperatorHandle Mv_stab_op_;
+
   // solver objects
   mfem::ParLORDiscretization *L_iorho_lor_ = nullptr;
   mfem::HypreBoomerAMG *L_iorho_inv_pc_ = nullptr;
@@ -374,6 +399,9 @@ class Tomboulides final : public FlowBase {
   mfem::Vector utheta_m2_vec_;
   mfem::Vector utheta_next_vec_;
 
+  mfem::Vector tmpR0b_;
+  mfem::Vector swDiff_vec_;
+
   // miscellaneous
   double volume_;
   mfem::ParLinearForm *mass_lform_ = nullptr;
@@ -392,7 +420,8 @@ class Tomboulides final : public FlowBase {
 
  public:
   /// Constructor
-  Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff, TPS::Tps *tps = nullptr);
+  Tomboulides(mfem::ParMesh *pmesh, int vorder, int porder, temporalSchemeCoefficients &coeff,
+              mfem::ParGridFunction *gridScale = nullptr, TPS::Tps *tps = nullptr);
 
   /// Destructor
   ~Tomboulides() final;
diff --git a/src/utils.cpp b/src/utils.cpp
index 553181e6..6eb17e29 100644
--- a/src/utils.cpp
+++ b/src/utils.cpp
@@ -1133,23 +1133,117 @@ void scalarGrad3DV(FiniteElementSpace *fes, FiniteElementSpace *vfes, Vector u,
   R1_gf.GetTrueDofs(*gu);
 }
 
-/*
-void vectorGrad3DV(FiniteElementSpace *fes, Vector u, Vector *gu, Vector *gv, Vector *gw) {
-  ParGridFunction R1_gf;
-  R1_gf.SetSpace(fes);
-
-  ParGridFunction R1a_gf, R1b_gf, R1c_gf;
-  R1a_gf.SetSpace(fes);
-  R1b_gf.SetSpace(fes);
-  R1c_gf.SetSpace(fes);
-
-  R1_gf.SetFromTrueDofs(u);
-  vectorGrad3D(R1_gf, R1a_gf, R1b_gf, R1c_gf);
-  R1a_gf.GetTrueDofs(*gu);
-  R1b_gf.GetTrueDofs(*gv);
-  R1c_gf.GetTrueDofs(*gw);
+void streamwiseTensor(const Vector &vel, DenseMatrix &swMgbl) {
+  int dim = vel.Size();
+
+  // streamwise coordinate system
+  Vector unitNorm;
+  Vector unitT1;
+  Vector unitT2;
+  unitNorm.SetSize(dim);
+  unitT1.SetSize(dim);
+  unitT2.SetSize(dim);
+
+  // streamwise direction
+  for (int i = 0; i < dim; i++) unitNorm[i] = vel[i];
+  double mod = 0.0;
+  for (int i = 0; i < dim; i++) mod += unitNorm[i] * unitNorm[i];
+  mod = std::max(mod, 1.0e-18);
+  double Umag = std::sqrt(mod);
+  unitNorm /= Umag;
+
+  // std::cout << Umag << " " << endl ;
+
+  // for zero-flow
+  if (Umag < 1.0e-8) {
+    // for (int i = 0; i < dim; i++) {
+    //   for (int j = 0; j < dim; j++) {
+    //     swMgbl(i, j) = 0.0;
+    //   }
+    // }
+    swMgbl = 0.0;
+    return;
+  }
+
+  // tangent direction (not unique)
+  int maxInd, minusInd, plusInd;
+  if (unitNorm[0] * unitNorm[0] > unitNorm[1] * unitNorm[1]) {
+    maxInd = 0;
+  } else {
+    maxInd = 1;
+  }
+  if (dim == 3 && unitNorm[maxInd] * unitNorm[maxInd] < unitNorm[2] * unitNorm[2]) {
+    maxInd = 2;
+  }
+  minusInd = ((maxInd - 1) % dim + dim) % dim;
+  plusInd = (maxInd + 1) % dim;
+
+  unitT1[minusInd] = -unitNorm[maxInd];
+  unitT1[maxInd] = unitNorm[minusInd];
+  if (dim == 3) {  // DOUBLE CHECK THIS WHEN TESTING 3D
+    unitT1[plusInd] = 0.0;
+  }
+  mod = 0.0;
+  for (int i = 0; i < dim; i++) mod += unitT1[i] * unitT1[i];
+  mod = std::max(mod, 1.0e-18);
+  unitT1 /= std::sqrt(mod);
+
+  // t2 is then orthogonal to both normal & t1
+  if (dim == 3) {
+    unitT2[0] = +(unitNorm[1] * unitT1[2] - unitNorm[2] * unitT1[1]);
+    unitT2[1] = -(unitNorm[0] * unitT1[2] - unitNorm[2] * unitT1[0]);
+    unitT2[2] = +(unitNorm[0] * unitT1[1] - unitNorm[1] * unitT1[0]);
+  }
+
+  // transform from streamwise coords to global
+  DenseMatrix M(dim, dim);
+  for (int d = 0; d < dim; d++) {
+    M(d, 0) = unitNorm[d];
+    M(d, 1) = unitT1[d];
+    if (dim == 3) M(d, 2) = unitT2[d];
+  }
+
+  // streamwise coeff
+  DenseMatrix swM(dim, dim);
+  swM = 0.0;
+  swM(0, 0) = 1.0;
+
+  // std::cout << " " << endl;
+  // for (int i = 0; i < dim; i++) {
+  //   for (int j = 0; j < dim; j++) {
+  //     std::cout << M(i,j) << " " ;
+  //   }
+  //   std::cout << endl;
+  // }
+
+  // M_{im} swM_{mn} M_{jn} or M*"mu"*M^T (with n,t1,t2 in columns of M)
+  // DenseMatrix swMgbl(dim, dim);
+  swMgbl = 0.0;
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      for (int m = 0; m < dim; m++) {
+        for (int n = 0; n < dim; n++) {
+          swMgbl(i, j) += M(i, m) * M(j, n) * swM(m, n);
+        }
+      }
+    }
+  }
+
+  // std::cout << " " << endl;
+  // for (int i = 0; i < dim; i++) {
+  //   for (int j = 0; j < dim; j++) {
+  //     std::cout << swMgbl(i,j) << " " ;
+  //   }
+  //   std::cout << endl;
+  // }
+}
+
+double csupgFactor(double Reh) {
+  // return  0.5 * (tanh(re_factor * Re - re_offset) + 1.0);
+  // printf("%f\n", Reh);
+  // printf("%f\n", 0.5 * (tanh(Reh) + 1.0));
+  return 0.5 * (tanh(Reh) + 1.0);
 }
-*/
 
 void EliminateRHS(Operator &A, ConstrainedOperator &constrainedA, const Array<int> &ess_tdof_list, Vector &x, Vector &b,
                   Vector &X, Vector &B, int copy_interior) {
@@ -1324,4 +1418,44 @@ void GradientVectorGridFunctionCoefficient::Eval(DenseMatrix &G, ElementTransfor
   }
 }
 
+VectorMagnitudeCoefficient::VectorMagnitudeCoefficient(VectorCoefficient &A) : a(&A), va(A.GetVDim()) {}
+
+void VectorMagnitudeCoefficient::SetTime(double t) {
+  if (a) {
+    a->SetTime(t);
+  }
+  this->Coefficient::SetTime(t);
+}
+
+double VectorMagnitudeCoefficient::Eval(ElementTransformation &T, const IntegrationPoint &ip) {
+  a->Eval(va, T, ip);
+  // double res = 0;
+  // for (int i = 0; i < va.size(); i++) { res += va[i] * va[i]}
+  // res = std::sqrt(res)
+  // return res;
+  double mod = std::max(std::sqrt(va * va), 1.0e-18);
+  return mod;
+}
+
+void TransformedMatrixVectorCoefficient::SetTime(double t) {
+  Q1->SetTime(t);
+  this->MatrixCoefficient::SetTime(t);
+}
+
+void TransformedMatrixVectorCoefficient::Eval(DenseMatrix &G, ElementTransformation &T, const IntegrationPoint &ip) {
+  Vector buf;
+  buf.SetSize(Q1->GetVDim());
+  Q1->Eval(buf, T, ip);
+  Function(buf, G);
+
+  // int dim = Q1->GetVDim();
+  // std::cout << " " << endl;
+  // for (int i = 0; i < dim; i++) {
+  //   for (int j = 0; j < dim; j++) {
+  //     std::cout << G(i,j) << " " ;
+  //   }
+  //   std::cout << endl;
+  // }
+}
+
 }  // namespace mfem
diff --git a/src/utils.hpp b/src/utils.hpp
index 4f23c5d9..a5071a8d 100644
--- a/src/utils.hpp
+++ b/src/utils.hpp
@@ -188,6 +188,11 @@ void makeContinuous(ParGridFunction &u);
 
 bool copyFile(const char *SRC, const char *DEST);
 
+/// upwind diffusion
+// double upwindDiffMag(Vector *state, double gridspace);
+void streamwiseTensor(const Vector &vel, DenseMatrix &swMgbl);
+double csupgFactor(double Reh);
+
 /// Eliminate essential BCs in an Operator and apply to RHS.
 /// rename this to something sensible "ApplyEssentialBC" or something
 void EliminateRHS(Operator &A, ConstrainedOperator &constrainedA, const Array<int> &ess_tdof_list, Vector &x, Vector &b,
@@ -224,6 +229,48 @@ class GradientVectorGridFunctionCoefficient : public MatrixCoefficient {
 
   virtual ~GradientVectorGridFunctionCoefficient() {}
 };
+
+/// Scalar coefficient defined as the magnitude of a vector coefficient
+class VectorMagnitudeCoefficient : public Coefficient {
+ private:
+  VectorCoefficient *a;
+
+  mutable Vector va;
+
+ public:
+  /// Construct with the vector coefficient.  Result is \sqrt{(\f$ A \cdot a \f$}.
+  VectorMagnitudeCoefficient(VectorCoefficient &A);
+
+  /// Set the time for internally stored coefficients
+  void SetTime(double t);
+
+  /// Reset the vector
+  void SetACoef(VectorCoefficient &A) { a = &A; }
+  /// Return the vector coefficient
+  VectorCoefficient *GetACoef() const { return a; }
+
+  /// Evaluate the coefficient at @a ip.
+  virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip);
+};
+
+/// Matrix coefficient computed from a function F(v(x)) of a dim-sized vector coefficient, v(x)
+class TransformedMatrixVectorCoefficient : public MatrixCoefficient {
+ protected:
+  VectorCoefficient *Q1;
+  std::function<void(const Vector &, DenseMatrix &)> Function;
+
+ public:
+  TransformedMatrixVectorCoefficient(VectorCoefficient *vc, std::function<void(const Vector &, DenseMatrix &)> F)
+      : MatrixCoefficient(vc->GetVDim()), Q1(vc), Function(std::move(F)) {}
+
+  /// Set the time for internally stored coefficients
+  void SetTime(double t);
+
+  // using MatrixCoefficient::Eval;
+  virtual void Eval(DenseMatrix &G, ElementTransformation &T, const IntegrationPoint &ip);
+
+  virtual ~TransformedMatrixVectorCoefficient() {}
+};
 }  // namespace mfem
 
 #endif  // UTILS_HPP_
diff --git a/src/wallBC.cpp b/src/wallBC.cpp
index 73bd193a..0c6f77ad 100644
--- a/src/wallBC.cpp
+++ b/src/wallBC.cpp
@@ -814,8 +814,8 @@ void WallBC::interpWalls_gpu(const mfem::Vector &x, const elementIndexingData &e
           d_flux[eq + q * num_equation + n * maxIntPoints * num_equation] = Rflux[eq];
         }
       }  // end quadrature point loop
-    }    // end face loop
-  });    // end element loop
+    }  // end face loop
+  });  // end element loop
 #endif
 }
 
diff --git a/test/.gitattributes b/test/.gitattributes
index 83d6aca2..c68e89fe 100644
--- a/test/.gitattributes
+++ b/test/.gitattributes
@@ -54,3 +54,4 @@ meshes/spongeBox.msh filter=lfs diff=lfs merge=lfs -text
 ref_solns/lequere-varmu/restart_output-lequere-varmu.sol.h5 filter=lfs diff=lfs merge=lfs -text
 ref_solns/lequere-varmu/reference-lequere-varmu.sol.h5 filter=lfs diff=lfs merge=lfs -text
 ref_solns/collisions/Oliver.fit.h5 filter=lfs diff=lfs merge=lfs -text
+ref_solns/stabChan/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
diff --git a/test/Makefile.am b/test/Makefile.am
index 56e794a1..e639d796 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -36,6 +36,7 @@ EXTRA_DIST  = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte-
               ref_solns/reactSingleRx/*.h5 \
               ref_solns/reactTable/*.h5 \
               ref_solns/radDecay/*.h5 \
+              ref_solns/stabChan/*.h5 \
 		      vpath.sh die.sh count_gpus.sh sniff_mpirun.sh \
 		      cyl3d.gpu.test cyl3d.mflow.gpu.test wedge.gpu.test \
 	          averaging.gpu.test cyl3d.test cyl3d.gpu.python.test cyl3d.mflow.test cyl3d.dtconst.test \
@@ -50,8 +51,9 @@ EXTRA_DIST  = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte-
 		      tabulated.test lte_mixture.test distance_fcn.test \
               sgsSmag.test sgsSigma.test heatEq.test sponge.test plate.test pipe.mix.test lte2noneq-restart.test \
               coupled-3d.interface.test plasma.axisym.test plasma.axisym.lte1d.test \
-	      lomach-flow.test lomach-lequere.test interpInlet.test sgsLoMach.test autoPeriodic.test aveLoMach.test \
-              reactFlow-binDiff.test reactFlow-singleRx.test reactFlow-table.test radiativeDecay.test
+	      	  lomach-flow.test lomach-lequere.test interpInlet.test sgsLoMach.test autoPeriodic.test aveLoMach.test \
+              reactFlow-binDiff.test reactFlow-singleRx.test reactFlow-table.test radiativeDecay.test \
+			  lomach-chan-stab.test
 
 
 TESTS = vpath.sh
@@ -125,7 +127,8 @@ TESTS += cyl3d.test \
          reactFlow-binDiff.test \
          reactFlow-singleRx.test \
          reactFlow-table.test \
-         radiativeDecay.test
+         radiativeDecay.test \
+		 lomach-chan-stab.test
 
 if PYTHON_ENABLED
 TESTS += cyl3d.python.test \
diff --git a/test/inputs/input.stabChan.ini b/test/inputs/input.stabChan.ini
new file mode 100644
index 00000000..cdfab8b2
--- /dev/null
+++ b/test/inputs/input.stabChan.ini
@@ -0,0 +1,77 @@
+[solver]
+type = loMach
+
+[loMach]
+mesh = meshes/channel182p4_25x9p4_THIN.msh
+order = 1
+nFilter = 0
+filterWeight = 1.0
+maxIters = 100
+outputFreq = 100
+fluid = dry_air
+refLength = 1.0
+equation_system = navier-stokes
+enablePressureForcing = True
+pressureGrad = '0.00005372 0.0 0.0'
+enableGravity = False
+gravity = '0.0 0.0 0.0'
+openSystem = False
+sgsModel = none
+flow-solver = tomboulides
+thermo-solver = calorically-perfect
+
+[loMach/calperfect]
+viscosity-model = sutherland
+sutherland/mu0 = 1.68e-5
+sutherland/T0 = 273.0
+sutherland/S0 = 110.4
+numerical-integ = false
+Prandtl = 0.72
+#ic = channel
+streamwise-stabilization = true
+
+[loMach/tomboulides]
+ic = channel
+numerical-integ = false
+streamwise-stabilization = true
+
+[io]
+outdirBase = output
+#enableRestart = True
+#restartMode = variableP
+
+[time]
+integrator = curlcurl
+cfl = 0.4
+dt_initial = 1.0e-4
+dtFactor = 0.01
+#dt_fixed = 1.0e-4
+bdfOrder = 2
+
+[spongeMultiplier]
+uniform = true
+uniformMult = 1.0
+
+[initialConditions]
+rho = 1.0
+rhoU = 1.0
+rhoV = 0.0
+rhoW = 0.0
+temperature = 300.0
+pressure = 101325.0
+
+[boundaryConditions/wall1]
+patch = 4
+type = viscous_isothermal
+temperature = 300
+velocity = '0.0 0.0 0.0'
+
+[boundaryConditions]
+numWalls = 1
+numInlets = 0
+numOutlets = 0
+
+[periodicity]
+enablePeriodic = True
+periodicX = True
+periodicZ = True
\ No newline at end of file
diff --git a/test/inputs/input.supgBox.ini b/test/inputs/input.supgBox.ini
new file mode 100644
index 00000000..3548e8f1
--- /dev/null
+++ b/test/inputs/input.supgBox.ini
@@ -0,0 +1,109 @@
+[solver]
+type = loMach
+
+[loMach]
+mesh = meshes/boxStab.msh
+order = 1
+nFilter = 0
+filterWeight = 1.0
+maxIters = 10000
+outputFreq = 500
+fluid = dry_air
+refLength = 1.0
+equation_system = navier-stokes
+enablePressureForcing = True
+# pressureGrad = '0.00005372 0.0 0.0'
+enableGravity = False
+gravity = '0.0 0.0 0.0'
+openSystem = False
+sgsModel = none
+# flow-solver = tomboulides
+flow-solver = zero-flow
+thermo-solver = calorically-perfect
+
+[loMach/calperfect]
+viscosity-model = sutherland
+sutherland/mu0 = 1.68e-5
+sutherland/T0 = 273.0
+sutherland/S0 = 110.4
+numerical-integ = false
+Prandtl = 0.72
+#ic = channel
+#hsolve-atol = 1.0e-12
+#msolve-atol = 1.0e-12
+#hsolve-rtol = 1.0e-10
+#msolve-rtol = 1.0e-10
+#hsolve-maxIters = 2000
+#msolve-maxIters = 2000
+streamwise-stabilization = true
+#Reh_offset = 1.0
+#Reh_factor = 0.01
+
+[loMach/zeroflow]
+nonzero-flow = true
+nonzero-vel = '1.0 0.0 0.0'
+
+[io]
+outdirBase = output_wstab1T
+#enableRestart = True
+#restartMode = variableP
+
+[time]
+enableConstantTimestep = True
+integrator = curlcurl
+#cfl = 0.4
+dt_initial = 1.0e-4
+#dtFactor = 0.01
+dt_fixed = 1.0e-4
+bdfOrder = 1
+
+[spongeMultiplier]
+uniform = true
+uniformMult = 1.0
+
+[initialConditions]
+rho = 1.0
+rhoU = 1.0
+rhoV = 0.0
+rhoW = 0.0
+temperature = 300.0
+pressure = 101325.0
+
+[boundaryConditions]
+numWalls = 0
+numInlets = 2
+numOutlets = 0
+
+#[boundaryConditions/wall1]
+#patch = 1
+#type = viscous_isothermal
+#temperature = 300
+#velocity = '0.0 0.0 0.0'
+#
+#[boundaryConditions/wall2]
+#patch = 3
+#type = viscous_isothermal
+#temperature = 300
+#velocity = '0.0 0.0 0.0'
+
+[boundaryConditions/inlet1]
+patch = 4
+type = uniform
+#swirl = 40.0
+#velocity = '1.0 0.0 0.0'
+temperature = 300.0
+
+[boundaryConditions/inlet2]
+patch = 2
+type = uniform
+temperature = 270
+
+[periodicity]
+enablePeriodic = True
+#periodicX = True
+periodicY = True
+
+# background velocity set by new ZeroFlow option
+
+# how to set up the specific BCs in config? inlets walls and outlets are hard coded
+# mainly need neumann walls, or is that periodicity?
\ No newline at end of file
diff --git a/test/lomach-chan-stab.test b/test/lomach-chan-stab.test
new file mode 100755
index 00000000..40798288
--- /dev/null
+++ b/test/lomach-chan-stab.test
@@ -0,0 +1,30 @@
+#!./bats
+# -*- mode: sh -*-
+
+TEST="stabChan"
+RUNFILE="inputs/input.stabChan.ini"
+EXE="../src/tps"
+RESTART="ref_solns/stabChan/restart_output.sol.h5"
+
+setup() {
+    SOLN_FILE="restart_output.sol.h5"
+    REF_FILE="ref_solns/stabChan/restart_output.sol.h5"
+    OUT_FILE="output_solns/restart_output_stabChan.sol.h5"    
+}
+
+@test "[$TEST] check for input file $RUNFILE" {
+    test -s $RUNFILE
+}
+
+@test "[$TEST] run tps with input -> $RUNFILE" {
+    rm -rf output/*
+    rm -f $SOLN_FILE
+    $EXE --runFile $RUNFILE
+    test -s $SOLN_FILE    
+}
+
+@test "[$TEST] verify tps output with input -> $RUNFILE" {
+    test -s $SOLN_FILE
+    test -s $REF_FILE
+    h5diff -r --delta=1e-10 $SOLN_FILE $REF_FILE /velocity
+}
\ No newline at end of file
diff --git a/test/meshes/boxStab.msh b/test/meshes/boxStab.msh
new file mode 100644
index 00000000..74d166c1
--- /dev/null
+++ b/test/meshes/boxStab.msh
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:cf889720de9350d00f33208efc95855b676a2a69d82ed0be121ee69454d53a26
+size 5004
diff --git a/test/meshes/channel182p4_25x9p4_THIN.msh b/test/meshes/channel182p4_25x9p4_THIN.msh
new file mode 100644
index 00000000..ec72d1aa
--- /dev/null
+++ b/test/meshes/channel182p4_25x9p4_THIN.msh
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:b8bba9f4101291f8cfce812f02d669bdbe7178b74674a08eb745bb96531cd9bb
+size 673400
diff --git a/test/ref_solns/stabChan/restart_output.sol.h5 b/test/ref_solns/stabChan/restart_output.sol.h5
new file mode 100644
index 00000000..8255c890
Binary files /dev/null and b/test/ref_solns/stabChan/restart_output.sol.h5 differ
diff --git a/test/standalone_chemistry.cpp b/test/standalone_chemistry.cpp
index 5af6bc82..6538581f 100644
--- a/test/standalone_chemistry.cpp
+++ b/test/standalone_chemistry.cpp
@@ -32,7 +32,7 @@ int main(int argc, char *argv[]) {
   temporalSchemeCoefficients temporal_coeff;
   temporal_coeff.order = 1;
 
-  ReactingFlow *thermo = new ReactingFlow(pmesh, &loMach_opts, temporal_coeff, &tps);
+  ReactingFlow *thermo = new ReactingFlow(pmesh, &loMach_opts, temporal_coeff, nullptr, &tps);
   thermo->initializeSelf();
 
   int nSpecies;