diff --git a/.gitattributes b/.gitattributes
index f98bb4510..84d62c6cd 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -30,3 +30,4 @@ test/ref_solns/aveLoMach/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -te
 test/ref_solns/reactNitrogen/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/supg/restart_output_wstab1T.sol.h5 filter=lfs diff=lfs merge=lfs -text
 test/ref_solns/supg/restart_output_lte_supg.sol.h5 filter=lfs diff=lfs merge=lfs -text
+test/ref_solns/stabChan/restart_output.sol.h5 filter=lfs diff=lfs merge=lfs -text
diff --git a/src/loMach.cpp b/src/loMach.cpp
index 521c8ad3c..befa822cf 100644
--- a/src/loMach.cpp
+++ b/src/loMach.cpp
@@ -186,7 +186,8 @@ void LoMachSolver::initialize() {
     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/tomboulides.cpp b/src/tomboulides.cpp
index 184df27b6..53f9206a0 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 1463b426e..a508115d9 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/test/.gitattributes b/test/.gitattributes
index 83d6aca25..c68e89fe1 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 879ab4596..260bb6609 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -39,6 +39,7 @@ EXTRA_DIST  = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte-
               ref_solns/radDecay/*.h5 \
               ref_solns/reactNitrogen/*.h5 \
               ref_solns/supg/*.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 \
@@ -52,10 +53,11 @@ EXTRA_DIST  = tap-driver.sh test_tps_splitcomm.py soln_differ inputs meshes lte-
 		      gradient.test coupled-3d.test coupled-3d.splitcomm.test coupled-3d.py-loop.test \
 		      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 \
+	      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 \
-              reactFlow-nitrogen.test argon_minimal.test argon_minimal.binary.test supg-scalars.test
+              reactFlow-nitrogen.test argon_minimal.test argon_minimal.binary.test supg-scalars.test \
+	      lomach-chan-stab.test
 
 TESTS = vpath.sh
 XFAIL_TESTS =
@@ -130,7 +132,8 @@ TESTS += cyl3d.test \
          reactFlow-table.test \
          radiativeDecay.test \
          reactFlow-nitrogen.test \
-	 supg-scalars.test
+	 supg-scalars.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 000000000..cdfab8b28
--- /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/lomach-chan-stab.test b/test/lomach-chan-stab.test
new file mode 100755
index 000000000..40798288c
--- /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/channel182p4_25x9p4_THIN.msh b/test/meshes/channel182p4_25x9p4_THIN.msh
new file mode 100644
index 000000000..ec72d1aa7
--- /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 000000000..c8b8feb37
--- /dev/null
+++ b/test/ref_solns/stabChan/restart_output.sol.h5
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:2737ebf44db149c10aab0c9ecc049448df6ae6a1e8ff3d77de5b3314408c0ab3
+size 150168