diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index d70a9500..0cbe6824 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -77,6 +77,10 @@ jobs:
         run: |
             cd ${GITHUB_WORKSPACE}/build/test
             ./steady_ns_dd_mms
+      - name: Test UnsteadyNS DD solver
+        run: |
+            cd ${GITHUB_WORKSPACE}/build/test
+            ./unsteady_ns_dd_mms
       - name: Test LinearElastic DD solver
         run: |
             cd ${GITHUB_WORKSPACE}/build/test
diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp
index b3b72113..13c25fe8 100644
--- a/src/unsteady_ns_solver.cpp
+++ b/src/unsteady_ns_solver.cpp
@@ -96,11 +96,13 @@ void UnsteadyNSSolver::InitializeTimeIntegration()
 {
    assert(dt > 0.0);
 
-   uu = massMat->CreateMonolithic();
+   uu = new SparseMatrix(vblock_offsets[1]);
+   SparseMatrix *tmp = massMat->CreateMonolithic();
+   (*uu) += *tmp;
    (*uu) *= bd0 / dt;
    (*uu) += (*M); // add viscous flux operator
    uu->Finalize();
-   // mass->Finalize();
+   delete tmp;
 
    // This should be run after AssembleOperator
    assert(systemOp);
@@ -153,5 +155,13 @@ void UnsteadyNSSolver::Step(double &time, int step)
    /* Solve for the next step */
    mumps->Mult(*RHS_step, *U_step);
 
+   /* remove pressure scalar if all dirichlet bc */
+   if (!pres_dbc)
+   {
+      double p_const = U_stepview->GetBlock(1).Sum() / U_stepview->GetBlock(1).Size();
+
+      U_stepview->GetBlock(1) -= p_const;
+   }
+
    time += dt;
 }
\ No newline at end of file
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index ae0b0131..5722d39d 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -13,6 +13,7 @@ file(COPY inputs/test.parser.yml DESTINATION ${CMAKE_BINARY_DIR}/test/inputs)
 add_executable(poisson_dd_mms poisson_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
 add_executable(stokes_dd_mms stokes_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
 add_executable(steady_ns_dd_mms steady_ns_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
+add_executable(unsteady_ns_dd_mms unsteady_ns_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
 add_executable(linelast_dd_mms linelast_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
 add_executable(advdiff_dd_mms advdiff_dd_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
 add_executable(dg_integ_mms dg_integ_mms.cpp $<TARGET_OBJECTS:scaleupROMObj> $<TARGET_OBJECTS:mmsSuiteObj>)
diff --git a/test/mms_suite.cpp b/test/mms_suite.cpp
index d498e5d0..eabcf31e 100644
--- a/test/mms_suite.cpp
+++ b/test/mms_suite.cpp
@@ -488,6 +488,150 @@ void CheckConvergence(const double &threshold)
 
 } // namespace steady_ns
 
+namespace unsteady_ns
+{
+
+UnsteadyNSSolver *SolveWithRefinement(const int num_refinement)
+{
+   config.dict_["mesh"]["uniform_refinement"] = num_refinement;
+   UnsteadyNSSolver *test = new UnsteadyNSSolver();
+
+   test->InitVariables();
+   test->InitVisualization();
+
+   test->AddBCFunction(mms::steady_ns::uFun_ex);
+   test->SetBdrType(BoundaryType::DIRICHLET);
+   test->AddRHSFunction(mms::steady_ns::fFun);
+   
+   const int dim = test->GetDim();
+   VectorFunctionCoefficient ucoeff(dim, mms::steady_ns::uFun_ex);
+   FunctionCoefficient pcoeff(mms::steady_ns::pFun_ex);
+   for (int m = 0; m < test->GetNumSubdomains(); m++)
+   {
+      test->GetVelGridFunction(m)->ProjectCoefficient(ucoeff);
+      test->GetPresGridFunction(m)->ProjectCoefficient(pcoeff);
+   }
+
+   test->BuildOperators();
+
+   test->SetupBCOperators();
+
+   test->Assemble();
+
+   test->Solve();
+
+   test->SaveVisualization();
+
+   return test;
+}
+
+void CheckConvergence(const double &threshold)
+{
+   mms::steady_ns::nu = config.GetOption<double>("stokes/nu", 1.0);
+   mms::steady_ns::zeta = config.GetOption<double>("navier-stokes/zeta", 1.0);
+
+   int num_refine = config.GetOption<int>("manufactured_solution/number_of_refinement", 3);
+   int base_refine = config.GetOption<int>("manufactured_solution/baseline_refinement", 0);
+
+   //printf("Num. Elem.\tRel. v err.\tConv Rate\tNorm\tRel. p err.\tConv Rate\tNorm\n");
+   printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n",
+          "Num. Elem.", "Rel v err", "Conv Rate", "Norm", "Rel p err", "Conv Rate", "Norm");
+
+   Vector uconv_rate(num_refine), pconv_rate(num_refine);
+   uconv_rate = 0.0;
+   pconv_rate = 0.0;
+   double uerror1 = 0.0, perror1 = 0.0;
+   for (int r = base_refine; r < num_refine; r++)
+   {
+      UnsteadyNSSolver *test = SolveWithRefinement(r);
+
+      // Compare with exact solution
+      int dim = test->GetDim();
+      VectorFunctionCoefficient exact_usol(dim, mms::steady_ns::uFun_ex);
+      FunctionCoefficient exact_psol(mms::steady_ns::pFun_ex);
+
+      // For all velocity dirichlet bc, pressure does not have the absolute value.
+      // specify the constant scalar for the reference value.
+      double p_const = 0.0;
+      int ps = 0;
+      for (int k = 0; k < test->GetNumSubdomains(); k++)
+      {
+         GridFunction *pk = test->GetPresGridFunction(k);
+         GridFunction p_ex(*pk);
+         p_ex.ProjectCoefficient(exact_psol);
+         ps += p_ex.Size();
+         p_const += p_ex.Sum();
+         // If p_ex is the view vector of pk, then this will prevent false negative test result.
+         p_ex += 1.0;
+      }
+      p_const /= static_cast<double>(ps);
+
+      for (int k = 0; k < test->GetNumSubdomains(); k++)
+      {
+         GridFunction *pk = test->GetPresGridFunction(k);
+         (*pk) += p_const;
+      }
+
+      int uorder = test->GetVelFEOrder();
+      int porder = test->GetPresFEOrder();
+      int order_quad = max(2, 2*uorder+1);
+      const IntegrationRule *irs[Geometry::NumGeom];
+      for (int i=0; i < Geometry::NumGeom; ++i)
+      {
+         irs[i] = &(IntRules.Get(i, order_quad));
+      }
+
+      int numEl = 0;
+      double unorm = 0.0, pnorm = 0.0;
+      for (int k = 0; k < test->GetNumSubdomains(); k++)
+      {
+         Mesh *mk = test->GetMesh(k);
+         unorm += pow(ComputeLpNorm(2.0, exact_usol, *mk, irs), 2);
+         pnorm += pow(ComputeLpNorm(2.0, exact_psol, *mk, irs), 2);
+         numEl += mk->GetNE();
+      }
+      unorm = sqrt(unorm);
+      pnorm = sqrt(pnorm);
+
+      double uerror = 0.0, perror = 0.0;
+      for (int k = 0; k < test->GetNumSubdomains(); k++)
+      {
+         GridFunction *uk = test->GetVelGridFunction(k);
+         GridFunction *pk = test->GetPresGridFunction(k);
+         uerror += pow(uk->ComputeLpError(2, exact_usol), 2);
+         perror += pow(pk->ComputeLpError(2, exact_psol), 2);
+      }
+      uerror = sqrt(uerror);
+      perror = sqrt(perror);
+
+      uerror /= unorm;
+      perror /= pnorm;
+      
+      if (r > base_refine)
+      {
+         uconv_rate(r) = uerror1 / uerror;
+         pconv_rate(r) = perror1 / perror;
+      }
+      printf("%10d\t%10.5E\t%10.5E\t%10.5E\t%10.5E\t%10.5E\t%10.5E\n", numEl, uerror, uconv_rate(r), unorm, perror, pconv_rate(r), pnorm);
+
+      // reported convergence rate
+      if (r > base_refine)
+      {
+         EXPECT_TRUE(uconv_rate(r) > pow(2.0, uorder+1) - threshold);
+         EXPECT_TRUE(pconv_rate(r) > pow(2.0, porder+1) - threshold);
+      }
+
+      uerror1 = uerror;
+      perror1 = perror;
+
+      delete test;
+   }
+
+   return;
+}
+
+}
+
 namespace linelast
 {
    void ExactSolution(const Vector &x, Vector &u)
diff --git a/test/mms_suite.hpp b/test/mms_suite.hpp
index 07e77202..cba2e6ae 100644
--- a/test/mms_suite.hpp
+++ b/test/mms_suite.hpp
@@ -9,6 +9,7 @@
 #include "poisson_solver.hpp"
 #include "stokes_solver.hpp"
 #include "steady_ns_solver.hpp"
+#include "unsteady_ns_solver.hpp"
 #include "linelast_solver.hpp"
 #include "advdiff_solver.hpp"
 #include <fstream>
@@ -77,6 +78,14 @@ void CheckConvergence(const double &threshold = 1.0);
 
 }   // namespace steady_ns
 
+namespace unsteady_ns
+{
+
+UnsteadyNSSolver *SolveWithRefinement(const int num_refinement);
+void CheckConvergence(const double &threshold = 1.0);
+
+}   // namespace steady_ns
+
 namespace linelast
 {
 
diff --git a/test/unsteady_ns_dd_mms.cpp b/test/unsteady_ns_dd_mms.cpp
new file mode 100644
index 00000000..4c99a377
--- /dev/null
+++ b/test/unsteady_ns_dd_mms.cpp
@@ -0,0 +1,39 @@
+// Copyright 2023 Lawrence Livermore National Security, LLC. See the top-level LICENSE file for details.
+//
+// SPDX-License-Identifier: MIT
+
+#include<gtest/gtest.h>
+#include "mms_suite.hpp"
+
+using namespace std;
+using namespace mfem;
+using namespace mms::unsteady_ns;
+
+/**
+ * Simple smoke test to make sure Google Test is properly linked
+ */
+TEST(GoogleTestFramework, GoogleTestFrameworkFound) {
+   SUCCEED();
+}
+
+TEST(DDSerialTest, Test_convergence)
+{
+   config = InputParser("inputs/dd_mms.yml");
+   config.dict_["navier-stokes"]["operator-type"] = "lf";
+   config.dict_["discretization"]["order"] = 1;
+   config.dict_["manufactured_solution"]["number_of_refinement"] = 3;
+   config.dict_["navier-stokes"]["timestep_size"] = 0.01;
+   config.dict_["navier-stokes"]["number_of_timesteps"] = 50;
+   CheckConvergence();
+
+   return;
+}
+
+int main(int argc, char* argv[])
+{
+   MPI_Init(&argc, &argv);
+   ::testing::InitGoogleTest(&argc, argv);
+   int result = RUN_ALL_TESTS();
+   MPI_Finalize();
+   return result;
+}
\ No newline at end of file