Merge branch 'main' into hpcftools/flowSolver-gpu

- added comments to Blf and Lf contributions
- split setup and step into vel, auxiliary and pressure part

.gitignore

@@ -19,6 +19,9 @@ CMakeFiles/
 # Clangd server cache
 *.cache*
 
+# VSCode configuration
+/.vscode/
+
 # Backup files
 *~
 

fem/dfem/parameterspace.hpp

@@ -101,10 +101,11 @@ public:
       // Setup DofToQuad information
       dtq.nqpt = (int)floor(std::pow(ir.GetNPoints(), 1.0 / mesh.Dimension()) + 0.5);
       dtq.ndof = dtq.nqpt;
+      dtq.mode = used_in_tensor_product ? DofToQuad::TENSOR : DofToQuad::FULL;
 
       // Calculate sizes
       const int num_qp = used_in_tensor_product ?
-                         std::pow(dtq.nqpt, mesh.Dimension()) :
+                         static_cast<int>(std::pow(dtq.nqpt, mesh.Dimension())) :
                          ir.GetNPoints();
 
       tsize = vdim * num_qp * mesh.GetNE();

fem/dfem/util.hpp

@@ -1702,12 +1702,12 @@ std::array<DofToQuadMap, N> load_dtq_mem(
    std::array<DofToQuadMap, N> f;
    for (std::size_t i = 0; i < N; i++)
    {
-      const auto [nqp_b, dim_b, ndof_b] = dtq[i].B.GetShape();
-      const auto B = Reshape(&dtq[i].B[0], nqp_b, dim_b, ndof_b);
-      auto mem_Bi = Reshape(reinterpret_cast<real_t *>(mem) + offset, nqp_b, dim_b,
-                            ndof_b);
       if (dtq[i].which_input != -1)
       {
+         const auto [nqp_b, dim_b, ndof_b] = dtq[i].B.GetShape();
+         const auto B = Reshape(&dtq[i].B[0], nqp_b, dim_b, ndof_b);
+         auto mem_Bi = Reshape(reinterpret_cast<real_t *>(mem) + offset, nqp_b, dim_b,
+                               ndof_b);
          MFEM_FOREACH_THREAD(q, x, nqp_b)
          {
             MFEM_FOREACH_THREAD(d, y, ndof_b)

miniapps/navier/CMakeLists.txt

@@ -52,6 +52,20 @@ if (MFEM_USE_MPI)
     ${NAVIER_COMMON_FILES}
     LIBRARIES mfem)
 
+  add_mfem_miniapp(incompressible_navier_dfem
+    MAIN incompressible_navier_dfem.cpp
+    EXTRA_HEADERS incompressible_navier_nvtx.hpp
+    LIBRARIES mfem)
+
+  add_mfem_miniapp(incompNS_2Dtest
+    MAIN incompNS_2Dtest.cpp
+    EXTRA_SOURCES incompressible_navier_solver.cpp
+                  incompressible_navier_tests.cpp
+    EXTRA_HEADERS incompressible_navier_solver.hpp
+                  incompressible_navier_nvtx.hpp
+    EXTRA_DEFINES MFEM_USE_CMAKE_TESTS
+    LIBRARIES mfem)
+
   add_mfem_miniapp(navier_turbchan
     MAIN navier_turbchan.cpp
     ${NAVIER_COMMON_FILES}
@@ -84,4 +98,4 @@ if (MFEM_USE_MPI)
         ${MPIEXEC_POSTFLAGS})
     endforeach()
 endif()
-endif ()
+endif()

miniapps/navier/incompNS_2Dtest.cpp

@@ -0,0 +1,205 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// 3D flow over a cylinder benchmark example
+
+#include "incompressible_navier_solver.hpp"
+
+#define NVTX_COLOR ::gpu::nvtx::kLawnGreen
+#include "incompressible_navier_nvtx.hpp"
+
+using namespace mfem;
+using namespace incompressible_navier;
+
+void vel(const Vector &x, real_t t, Vector &u)
+{
+   // real_t xi = x(0), yi = x(1);
+   u = 0.0;
+}
+
+void vel_inlet(const Vector &x, real_t t, Vector &u)
+{
+   u = 0.0;
+   if (x(0) < 0.001) { u(0) = -0.001 * (std::pow(x(1) - 0.5, 2.0) - 0.25); }
+}
+
+MFEM_EXPORT int navier(int argc, char *argv[], double &u, double &p, double &Ψ)
+{
+   dbg();
+   static mfem::MPI_Session mpi(argc, argv);
+   const int myid = mpi.WorldRank();
+   Hypre::Init();
+
+   const char *device_config = "cpu";
+   int serial_refinements = 1;
+   int nx = 90, ny = 30;
+   int v_order = 2;
+   int p_order = 1;
+   int t_order = 1;
+   real_t kin_vis = 20.0;
+   real_t dt = 1e-2;
+   real_t t = 0.0;
+   real_t t_final = 1.0;
+   bool last_step = false;
+   bool visualization = true;
+   bool use_paraview = false;
+   bool pa = false;
+   int vis_steps = 100;
+   int max_tsteps = -1;
+
+   constexpr int precision = 8;
+   std::cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&serial_refinements, "-sr", "--serial-refinements",
+                  "Number serial refinements.");
+   args.AddOption(&nx, "-nx", "--nx", "Number of elements in X.");
+   args.AddOption(&ny, "-ny", "--ny", "Number of elements in Y.");
+   args.AddOption(&v_order, "-vo", "--vo", "Order.");
+   args.AddOption(&p_order, "-po", "--po", "Order.");
+   args.AddOption(&t_order, "-to", "--to", "Order.");
+   args.AddOption(&kin_vis, "-kv", "--kin-vis",
+                  "Kineic viscosity coefficient.");
+   args.AddOption(&dt, "-dt", "--time-step", "Initial time step size.");
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly",
+                  "Enable or disable partial assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&use_paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Use ParaView.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&max_tsteps, "-ms", "--max-steps",
+                  "Maximum number of steps (negative means no restriction).");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0) { args.PrintUsage(mfem::out); }
+      return EXIT_FAILURE;
+   }
+   if (myid == 0) { args.PrintOptions(mfem::out); }
+
+   // Mesh *mesh = new Mesh("box-cylinder.mesh");
+   const real_t sx = 3.0, sy = 1.0;
+   const bool generate_edges = true;
+   const auto QUAD = Element::QUADRILATERAL;
+   Mesh mesh = Mesh::MakeCartesian2D(nx, ny, QUAD, generate_edges, sx, sy);
+
+   for (int i = 0; i < serial_refinements; ++i) { mesh.UniformRefinement(); }
+
+   if (Mpi::Root())
+   {
+      std::cout << "Number of elements: " << mesh.GetNE() << std::endl;
+   }
+
+   auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
+
+   // Create the flow solver.
+   IncompressibleNavierSolver flowsolver(pmesh, v_order, p_order, t_order,
+                                         kin_vis);
+   flowsolver.EnablePA(pa);
+
+   // // Set the initial condition.
+   // ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
+   // VectorFunctionCoefficient u_excoeff(pmesh->Dimension(), vel);
+   // u_ic->ProjectCoefficient(u_excoeff);
+
+   // Add Dirichlet boundary conditions to velocity space restricted to
+   // selected attributes on the mesh.
+   Array<int> attr(pmesh->bdr_attributes.Max());
+   attr = 0;
+   Array<int> attr_inlet(pmesh->bdr_attributes.Max());
+   attr_inlet = 0;
+   // Inlet is attribute 1.
+   attr[0] = 1;
+   // Walls is attribute 3.
+   attr[2] = 1;
+   flowsolver.AddVelDirichletBC(vel, attr);
+
+   attr_inlet[3] = 1;
+   flowsolver.AddVelDirichletBC(vel_inlet, attr_inlet);
+
+   flowsolver.Setup(dt);
+
+   ParGridFunction *u_gf = flowsolver.GetCurrentVelocity();
+   ParGridFunction *p_gf = flowsolver.GetCurrentPressure();
+   ParGridFunction *psi_gf = flowsolver.GetCurrentPsi();
+
+   ParaViewDataCollection pvdc("3dfoc", pmesh);
+   if (use_paraview)
+   {
+      pvdc.SetDataFormat(VTKFormat::BINARY32);
+      // pvdc.SetHighOrderOutput(true);
+      pvdc.SetCycle(0);
+      pvdc.SetTime(t);
+      pvdc.RegisterField("velocity", u_gf);
+      pvdc.RegisterField("pressure", p_gf);
+      pvdc.RegisterField("psi", psi_gf);
+      pvdc.Save();
+   }
+
+   for (int step = 0; !last_step; ++step)
+   {
+      if (step == max_tsteps) { last_step = true; }
+      if (t + dt >= t_final - dt / 2) { last_step = true; }
+      const bool vis_step = last_step || (step % vis_steps) == 0;
+
+      flowsolver.Step(t, dt, step, vis_step);
+
+      if (vis_step)
+      {
+         if (Mpi::Root() && vis_steps)
+         {
+            printf("%11s %11s\n", "Time", "dt");
+            printf("%.5E %.5E\n", t, dt);
+         }
+         if (use_paraview)
+         {
+            pvdc.SetCycle(step);
+            pvdc.SetTime(t);
+            pvdc.Save();
+         }
+      }
+   }
+
+   // flowsolver.PrintTimingData();
+
+   auto reduce = [](ParGridFunction *gf) -> real_t { return (*gf) * (*gf); };
+   // auto reduce = [](ParGridFunction *gf) -> real_t { return gf->Norml2(); };
+   u = reduce(u_gf), p = reduce(p_gf), Ψ = reduce(psi_gf);
+
+   fflush(stdout);
+   delete pmesh;
+
+   return EXIT_SUCCESS;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+#ifndef MFEM_USE_CMAKE_TESTS
+int main(int argc, char *argv[])
+try
+{
+   dbg();
+   double u, p, Ψ; // unused
+   return navier(argc, argv, u, p, Ψ);
+}
+catch (std::exception &e)
+{
+   std::cerr << "\033[31m..xxxXXX[ERROR]XXXxxx.." << std::endl;
+   std::cerr << "\033[31m{}" << e.what() << std::endl;
+   return EXIT_FAILURE;
+}
+#endif // MFEM_USE_CMAKE_TESTS

miniapps/navier/incompressible_navier_dfem.cpp

@@ -0,0 +1,397 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+using namespace mfem;
+using namespace mfem::future;
+
+#include "linalg/tensor.hpp"
+using mfem::future::tensor;
+
+#define NVTX_COLOR ::gpu::nvtx::kOrchid
+#include "incompressible_navier_nvtx.hpp"
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void DiffusionSetup(const ParFiniteElementSpace &sfes,
+                    const IntegrationRule &ir,
+                    const Array<int> &domain_attributes,
+                    ParameterFunction &qdata)
+{
+   NVTX_MARK_FUNCTION;
+   auto pmesh = sfes.GetParMesh();
+   auto nodes = static_cast<ParGridFunction *>(pmesh->GetNodes());
+   auto mfes = nodes->ParFESpace();
+
+   constexpr int U = 0, Ξ = 1, Δ = 2;
+
+   DifferentiableOperator dop(
+   {{ U, &sfes }},
+   {
+      {  { Ξ, mfes },
+         { Δ, &qdata.GetParameterSpace() }
+      }
+   },
+   *pmesh);
+   const auto qfunc =
+      [] MFEM_HOST_DEVICE(const tensor<real_t, DIM, DIM> &J,
+                          const real_t &w)
+   {
+      auto invJ = inv(J);
+      tensor<real_t, DIM, DIM> C{};
+      C(0, 0) = M_PI;
+      C(0, 1) = 0, C(1, 0) = 0;
+      assert(C(0, 1) == 0 && C(1, 0) == 0); // diff otherwise
+      C(1, 1) = 1.0 / M_PI;
+      return tuple{ C * invJ * transpose(invJ) * det(J) * w };
+   };
+   dop.AddDomainIntegrator(qfunc,
+                           tuple{ Gradient<Ξ>{}, Weight{} }, // inputs
+                           tuple{ Identity<Δ>{} },           // outputs
+                           ir, domain_attributes);
+   dop.SetParameters({ nodes, &qdata });
+
+   Vector unused(sfes.GetTrueVSize());
+   dop.Mult(unused, qdata);
+   qdata.HostRead();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void DiffusionApply(const ParFiniteElementSpace &sfes,
+                    const IntegrationRule &ir,
+                    const Array<int> &domain_attributes,
+                    ParameterFunction &qdata,
+                    const Vector &x, Vector &y)
+{
+   NVTX_MARK_FUNCTION;
+   auto pmesh = sfes.GetParMesh();
+   constexpr int U = 0, Q = 1;
+   auto qd_ps = &qdata.GetParameterSpace();
+   DifferentiableOperator dop({ { U, &sfes } }, { { Q, qd_ps } }, *pmesh);
+   const auto qfunc =[] MFEM_HOST_DEVICE(const tensor<real_t, DIM> &∇u,
+                                         const tensor<real_t, DIM, DIM> &Q)
+   {
+      return tuple{ Q * ∇u };
+   };
+   dop.AddDomainIntegrator(
+      qfunc,
+      tuple{ Gradient<U>{}, Identity<Q>{} },
+      tuple{ Gradient<U>{} },
+      ir, domain_attributes);
+   dop.SetParameters({ &qdata });
+   dop.Mult(x, y);
+   y.HostRead();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+int DiffVerification(ParFiniteElementSpace &h1fes,
+                     const IntegrationRule &ir,
+                     const Vector &qdata,
+                     const Vector &x, const Vector &y)
+{
+   NVTX_MARK_FUNCTION;
+   constexpr real_t ϵ = 1e-12;
+
+   MatrixFunctionCoefficient matrix_coeff(
+      DIM, [](const Vector &, DenseMatrix &C)
+   {
+      C.SetSize(DIM);
+      C(0, 0) = M_PI;
+      C(0, 1) = 0, C(1, 0) = 0;
+      assert(C(0, 1) == 0 && C(1, 0) == 0); // diff otherwise
+      C(1, 1) = 1.0 / M_PI;
+   });
+
+   ParBilinearForm a(&h1fes);
+
+   auto diff_integ = new DiffusionIntegrator(matrix_coeff);
+   diff_integ->SetIntRule(&ir);
+   a.AddDomainIntegrator(diff_integ);
+   a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   OperatorPtr A;
+   a.Assemble(), a.Finalize();
+   a.FormSystemMatrix(Array<int> {}, A);
+
+   Vector y2(h1fes.TrueVSize());
+   y2 = 0.0;
+   A->Mult(x, y2);
+   y2.HostRead();
+
+   Vector diff(y2);
+   diff -= y;
+   const auto diff_norm = diff.Norml2();
+   if (diff_norm > ϵ)
+   {
+      dbg("\x1B[31m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+      return EXIT_FAILURE;
+   }
+
+   dbg("\x1B[32m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+   return EXIT_SUCCESS;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void MassApply(const ParFiniteElementSpace &sfes,
+               const IntegrationRule &ir,
+               const Array<int> &domain_attributes,
+               const Vector &x, Vector &y)
+{
+   dbg();
+   auto &pmesh = *sfes.GetParMesh();
+   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
+   auto *mfes = nodes->ParFESpace();
+   constexpr int U = 0, Coords = 1;
+   DifferentiableOperator dop({{ U, &sfes }}, {{ Coords, mfes }}, pmesh);
+   const auto mf_mass_qf =
+      [](const real_t &dudxi,
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{ dudxi * w * det(J) };
+   };
+   dop.AddDomainIntegrator(
+      mf_mass_qf,
+      tuple{ Value<U>{}, Gradient<Coords>{}, Weight{} },
+      tuple{ Value<U>{} },
+      ir, domain_attributes);
+   dop.SetParameters({ nodes });
+   // Vector X(sfes.GetTrueVSize()), Y(sfes.GetTrueVSize());
+   // sfes.GetRestrictionMatrix()->Mult(x, X);
+   // dop.Mult(X, Y);
+   dop.Mult(x, y);
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+int MassVerification(ParFiniteElementSpace &h1fes,
+                     const IntegrationRule &ir,
+                     const Vector &x, Vector &y)
+{
+   ParBilinearForm a(&h1fes);
+   auto mass_integ = new MassIntegrator;
+   mass_integ->SetIntRule(&ir);
+   a.AddDomainIntegrator(mass_integ);
+   a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   a.Assemble(), a.Finalize();
+   Vector y2(h1fes.TrueVSize());
+   a.Mult(x, y2);
+   y2.HostRead();
+
+   Vector diff(y2);
+   diff -= y;
+
+   const auto diff_norm = diff.Norml2();
+   constexpr real_t ϵ = 1e-12;
+   if (diff_norm > ϵ)
+   {
+      dbg("\x1B[31m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+      return EXIT_FAILURE;
+   }
+   dbg("\x1B[32m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+   return EXIT_SUCCESS;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void VectorDiffApply(const ParFiniteElementSpace &vfes,
+                     const IntegrationRule &ir,
+                     const Array<int> &domain_attributes,
+                     const Vector &x,
+                     Vector &y)
+{
+   NVTX_MARK_FUNCTION;
+   auto pmesh = vfes.GetParMesh();
+   auto nodes = static_cast<ParGridFunction *>(pmesh->GetNodes());
+   auto mfes = nodes->ParFESpace();
+   constexpr int U = 0, Coords = 1;
+   DifferentiableOperator dop({{ U, &vfes }}, {{ Coords, mfes }}, *pmesh);
+   const auto qfunc =
+      [] MFEM_HOST_DEVICE(const tensor<real_t, DIM, DIM> &∇u,
+                          const tensor<real_t, DIM, DIM> &J, const real_t &w)
+   {
+      return tuple{ ∇u * inv(J) * det(J) * w * transpose(inv(J)) };
+   };
+   dop.AddDomainIntegrator(qfunc,
+                           tuple{ Gradient<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Gradient<U>{} },
+                           ir, domain_attributes);
+   dop.SetParameters({ nodes });
+   dop.Mult(x, y);
+   y.HostRead();
+}
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM, int VDIM = DIM>
+int VectorDiffVerif(ParFiniteElementSpace &h1fes,
+                    const IntegrationRule &ir, const Vector &x,
+                    Vector &y)
+{
+   NVTX_MARK_FUNCTION;
+   ParBilinearForm a(&h1fes);
+   auto A_integ = new VectorDiffusionIntegrator(VDIM);
+   A_integ->SetIntRule(&ir);
+   a.AddDomainIntegrator(A_integ);
+   a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   a.Assemble(), a.Finalize();
+   Vector y2(h1fes.TrueVSize());
+   a.Mult(x, y2);
+   y2.HostRead();
+
+   Vector diff(y2);
+   diff -= y;
+
+   const auto diff_norm = diff.Norml2();
+   constexpr real_t ϵ = 1e-12;
+   if (diff_norm > ϵ)
+   {
+      dbg("\x1B[31m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+      return EXIT_FAILURE;
+   }
+   dbg("\x1B[32m||dFdu_FD u^* - ex||_l2 = {}", diff_norm);
+   return EXIT_SUCCESS;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[]) try
+{
+   NVTX_MARK_FUNCTION;
+   constexpr int DIM = 2, VDIM = DIM;
+
+   static mfem::MPI_Session mpi(argc, argv);
+   const int myid = mpi.WorldRank();
+   Hypre::Init();
+
+   const char *device_config = "cpu";
+   const char *mesh_file = "none";
+   int serial_refinements = 0;
+   int nx = 1, ny = 1;
+   int p = 1;
+   bool visualization = false;
+   bool pa = false;
+
+   std::cout.precision(8);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&p, "-o", "--order", "Finite element order.");
+   args.AddOption(&serial_refinements, "-sr", "--serial-refinements",
+                  "Number serial refinements.");
+   args.AddOption(&nx, "-nx", "--nx", "Number of elements in X.");
+   args.AddOption(&ny, "-ny", "--ny", "Number of elements in Y.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly",
+                  "Enable or disable partial assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0) { args.PrintUsage(mfem::out); }
+      return EXIT_FAILURE;
+   }
+   if (myid == 0) { args.PrintOptions(mfem::out); }
+
+   Mesh smesh;
+   if (std::string(mesh_file) != "none")
+   {
+      smesh = Mesh(mesh_file);
+   }
+   else
+   {
+      const real_t sx = 3.0, sy = 1.0;
+      const bool generate_edges = true;
+      const auto QUAD = Element::QUADRILATERAL;
+      smesh = Mesh::MakeCartesian2D(nx, ny, QUAD, generate_edges, sx, sy);
+   }
+   MFEM_ASSERT(smesh.Dimension() == 2, "2D mesh required!");
+
+   for (int i = 0; i < serial_refinements; ++i) { smesh.UniformRefinement(); }
+
+   dbg("Number of elements: {}", smesh.GetNE());
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   smesh.Clear();
+
+   pmesh.EnsureNodes();
+   pmesh.SetCurvature(p);
+   assert(DIM == pmesh.Dimension());
+
+   Array<int> domain_attributes;
+   if (pmesh.attributes.Size() > 0)
+   {
+      domain_attributes.SetSize(pmesh.attributes.Max());
+      domain_attributes = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace fes(&pmesh, &fec), vfes(&pmesh, &fec, VDIM);
+   dbg("#dofs:{} ", fes.GetTrueVSize());
+
+   const auto &fe = *fes.GetFE(0);
+   const auto &ir =
+      IntRules.Get(fe.GetGeomType(), fe.GetOrder() + fe.GetOrder() + fe.GetDim() - 1);
+
+   dbg("#ndof per el = {}", fe.GetDof());
+   dbg("#nqp = {}", ir.GetNPoints());
+   dbg("#q1d = {}", (int)floor(pow(ir.GetNPoints(), 1.0 / DIM) + 0.5));
+
+   ParGridFunction f1_gf(&fes);
+
+   auto f1 = [](const Vector &coords)
+   {
+      assert(DIM == 2);
+      const double x = coords(0), y = coords(1);
+      return M_PI + x + x * x + x * y + y;
+   };
+   FunctionCoefficient f1_c(f1);
+   f1_gf.ProjectCoefficient(f1_c);
+
+   Vector x(f1_gf), y(fes.GetTrueVSize());
+
+   UniformParameterSpace qd_ps(pmesh, ir, DIM * DIM);
+   ParameterFunction qdata(qd_ps);
+
+   dbg("Diffusion setup, apply & verification");
+   DiffusionSetup<DIM>(fes, ir, domain_attributes, qdata);
+   DiffusionApply<DIM>(fes, ir, domain_attributes, qdata, x, y);
+   if (DiffVerification<DIM>(fes, ir, qdata, x, y) != EXIT_SUCCESS) { return EXIT_FAILURE; }
+
+   dbg("Mass apply & verification");
+   MassApply<DIM>(fes, ir, domain_attributes, x, y);
+   if (MassVerification<DIM>(fes, ir, x, y) != EXIT_SUCCESS) { return EXIT_FAILURE; }
+
+   dbg("Vector diffusion apply");
+   VectorFunctionCoefficient vf1_c(VDIM, [](const Vector &coords, Vector &u)
+   {
+      assert(DIM == 2);
+      const double x = coords(0), y = coords(1);
+      u(0) = M_PI + 0.25 * x * x * y + y * y * x;
+      u(1) = M_PI - 0.25 * x * y * y + y * x * x;
+   });
+   ParGridFunction vf1_gf(&vfes);
+   vf1_gf.ProjectCoefficient(vf1_c);
+   Vector vx(vf1_gf), vy(vfes.GetTrueVSize());
+   VectorDiffApply<DIM>(vfes, ir, domain_attributes, vx, vy);
+   if (VectorDiffVerif<DIM>(vfes, ir, vx, vy) != EXIT_SUCCESS) { return EXIT_FAILURE; }
+
+   return EXIT_SUCCESS;
+}
+catch (std::exception &e)
+{
+   std::cerr << "\033[31m..xxxXXX[ERROR]XXXxxx.." << std::endl;
+   std::cerr << "\033[31m{}" << e.what() << std::endl;
+   return EXIT_FAILURE;
+}

miniapps/navier/incompressible_navier_nvtx.hpp

@@ -0,0 +1,478 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+#pragma once
+
+#define FMT_HEADER_ONLY
+#include <fmt/format.h>
+
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <cstdint>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <mutex>
+#include <stack>
+#include <string>
+
+#ifdef MFEM_USE_CALIPER
+#include <caliper/cali.h>
+#endif
+
+#ifdef MFEM_USE_CUDA
+#include <cudaProfiler.h>
+#include <cuda_runtime_api.h>
+#include <nvToolsExt.h>
+#else
+struct nvtxEventAttributes_t
+{
+   int version;
+   int size;
+   int category;
+   int colorType;
+   uint32_t color;
+   int payloadType;
+   uint64_t payload;
+   int messageType;
+   struct
+   {
+      std::string ascii;
+   } message;
+};
+#define NVTX_VERSION 1
+#define NVTX_EVENT_ATTRIB_STRUCT_SIZE 256
+#define NVTX_COLOR_ARGB 0
+#define NVTX_MESSAGE_TYPE_ASCII 0
+#define nvtxRangePushEx(...)
+#define nvtxRangePop(...)
+#define cudaStreamSynchronize(...)
+#endif
+
+namespace gpu::nvtx
+{
+
+///////////////////////////////////////////////////////////////////////////////
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+// http://www.calmar.ws/vim/256-xterm-24bit-rgb-color-chart.html
+// clang-format off
+enum color_names
+{
+   kBlack = 0, kNavyBlue, kDarkBlue, kMediumBlue, kBlue, kDarkGreen, kWebGreen, kTeal,
+   kDarkCyan, kDeepSkyBlue, kDarkTurquoise, kMediumSpringGreen, kGreen, kLime,
+   kSpringGreen, kAqua, kCyan, kMidnightBlue, kDodgerBlue, kLightSeaGreen, kForestGreen,
+   kSeaGreen, kDarkSlateGray, kLimeGreen, kMediumSeaGreen, kTurquoise, kRoyalBlue,
+   kSteelBlue, kDarkSlateBlue, kMediumTurquoise, kIndigo, kDarkOliveGreen, kCadetBlue,
+   kCornflower, kRebeccaPurple, kMediumAquamarine, kDimGray, kSlateBlue, kOliveDrab,
+   kSlateGray, kLightSlateGray, kMediumSlateBlue, kLawnGreen, kWebMaroon, kWebPurple,
+   kChartreuse, kAquamarine, kOlive, kWebGray, kSkyBlue, kLightSkyBlue, kBlueViolet,
+   kDarkRed, kDarkMagenta, kSaddleBrown, kDarkSeaGreen, kLightGreen, kMediumPurple,
+   kDarkViolet, kPaleGreen, kDarkOrchid, kYellowGreen, kPurple, kSienna, kBrown,
+   kDarkGray, kLightBlue, kGreenYellow, kPaleTurquoise, kMaroon, kLightSteelBlue,
+   kPowderBlue, kFirebrick, kDarkGoldenrod, kMediumOrchid, kRosyBrown, kDarkKhaki,
+   kGray, kSilver, kMediumVioletRed, kIndianRed, kPeru, kChocolate, kTan, kLightGray,
+   kThistle, kOrchid, kGoldenrod, kPaleVioletRed, kCrimson, kGainsboro, kPlum, kBurlywood,
+   kLightCyan, kLavender, kDarkSalmon, kViolet, kPaleGoldenrod, kLightCoral, kKhaki,
+   kAliceBlue, kHoneydew, kAzure, kSandyBrown, kWheat, kBeige, kWhiteSmoke, kMintCream,
+   kGhostWhite, kSalmon, kAntiqueWhite, kLinen, kLightGoldenrod, kOldLace, kRed,
+   kFuchsia, kMagenta, kDeepPink, kOrangeRed, kTomato, kHotPink, kCoral, kDarkOrange,
+   kLightSalmon, kOrange, kLightPink, kPink, kGold, kPeachPuff, kNavajoWhite, kMoccasin,
+   kBisque, kMistyRose, kBlanchedAlmond, kPapayaWhip, kLavenderBlush, kSeashell,
+   kCornsilk, kLemonChiffon, kFloralWhite, kSnow, kYellow, kLightYellow, kIvory, kWhite,
+   kNvidia
+};
+// clang-format on
+
+static constexpr int kNumHexColors = 146;
+static constexpr std::array<uint32_t, kNumHexColors> kHexColors =
+{
+   {
+      0x000000, 0x000080, 0x00008B, 0x0000CD, 0x0000FF, 0x006400, 0x008000,
+      0x008080, 0x008B8B, 0x00BFFF, 0x00CED1, 0x00FA9A, 0x00FF00, 0x00FF00,
+      0x00FF7F, 0x00FFFF, 0x00FFFF, 0x191970, 0x1E90FF, 0x20B2AA, 0x228B22,
+      0x2E8B57, 0x2F4F4F, 0x32CD32, 0x3CB371, 0x40E0D0, 0x4169E1, 0x4682B4,
+      0x483D8B, 0x48D1CC, 0x4B0082, 0x556B2F, 0x5F9EA0, 0x6495ED, 0x663399,
+      0x66CDAA, 0x696969, 0x6A5ACD, 0x6B8E23, 0x708090, 0x778899, 0x7B68EE,
+      0x7CFC00, 0x7F0000, 0x7F007F, 0x7FFF00, 0x7FFFD4, 0x808000, 0x808080,
+      0x87CEEB, 0x87CEFA, 0x8A2BE2, 0x8B0000, 0x8B008B, 0x8B4513, 0x8FBC8F,
+      0x90EE90, 0x9370DB, 0x9400D3, 0x98FB98, 0x9932CC, 0x9ACD32, 0xA020F0,
+      0xA0522D, 0xA52A2A, 0xA9A9A9, 0xADD8E6, 0xADFF2F, 0xAFEEEE, 0xB03060,
+      0xB0C4DE, 0xB0E0E6, 0xB22222, 0xB8860B, 0xBA55D3, 0xBC8F8F, 0xBDB76B,
+      0xBEBEBE, 0xC0C0C0, 0xC71585, 0xCD5C5C, 0xCD853F, 0xD2691E, 0xD2B48C,
+      0xD3D3D3, 0xD8BFD8, 0xDA70D6, 0xDAA520, 0xDB7093, 0xDC143C, 0xDCDCDC,
+      0xDDA0DD, 0xDEB887, 0xE0FFFF, 0xE6E6FA, 0xE9967A, 0xEE82EE, 0xEEE8AA,
+      0xF08080, 0xF0E68C, 0xF0F8FF, 0xF0FFF0, 0xF0FFFF, 0xF4A460, 0xF5DEB3,
+      0xF5F5DC, 0xF5F5F5, 0xF5FFFA, 0xF8F8FF, 0xFA8072, 0xFAEBD7, 0xFAF0E6,
+      0xFAFAD2, 0xFDF5E6, 0xFF0000, 0xFF00FF, 0xFF00FF, 0xFF1493, 0xFF4500,
+      0xFF6347, 0xFF69B4, 0xFF7F50, 0xFF8C00, 0xFFA07A, 0xFFA500, 0xFFB6C1,
+      0xFFC0CB, 0xFFD700, 0xFFDAB9, 0xFFDEAD, 0xFFE4B5, 0xFFE4C4, 0xFFE4E1,
+      0xFFEBCD, 0xFFEFD5, 0xFFF0F5, 0xFFF5EE, 0xFFF8DC, 0xFFFACD, 0xFFFAF0,
+      0xFFFAFA, 0xFFFF00, 0xFFFFE0, 0xFFFFF0, 0xFFFFFF, 0x76B900
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+constexpr size_t static_strlen(const char *str)
+{
+   return *str == '\0' ? 0 : static_strlen(str + 1) + 1;
+}
+
+constexpr uint8_t static_checksum8(const char *bfr)
+{
+   unsigned int chk = 0;
+   size_t len = static_strlen(bfr);
+   for (; len; len--, bfr++) { chk += static_cast<unsigned int>(*bfr); }
+   return static_cast<uint8_t>(chk);
+}
+
+constexpr char *static_strrnchr(const char *str, const char c, int n)
+{
+   size_t len = static_strlen(str);
+   char *p = const_cast<char *>(str) + len - 1;
+   for (; n; n--, p--, len--)
+   {
+      for (; len; p--, len--)
+      {
+         if (*p == c) { break; }
+      }
+      if (!len) { return nullptr; }
+      if (n == 1) { return p; }
+   }
+   return nullptr;
+}
+
+inline uint32_t static_color(const uint8_t COLOR, const int RANK,
+                             const char *FILE)
+{
+   constexpr auto kMpiColorShift = 1;
+   const auto rank_shift = kMpiColorShift * RANK;
+   if (COLOR > 0) { return kHexColors[COLOR + rank_shift]; }
+   const auto file_color = static_checksum8(FILE);
+   return kHexColors[(file_color + rank_shift) % kNumHexColors];
+}
+
+///////////////////////////////////////////////////////////////////////////////
+// Helpers to generate unique variable names
+#define NVTX_FLF __FILE__, __LINE__, __FUNCTION__
+#define NVTX_PRIVATE_NAME(prefix) NVTX_PRIVATE_CONCAT(prefix, __LINE__)
+#define NVTX_PRIVATE_CONCAT(a, b) NVTX_PRIVATE_CONCAT2(a, b)
+#define NVTX_PRIVATE_CONCAT2(a, b) a##b
+
+#ifndef NVTX_COLOR
+#define NVTX_COLOR ::gpu::nvtx::kBlack
+#endif
+
+///////////////////////////////////////////////////////////////////////////////
+struct Debug
+{
+   const bool debug = false, end = true;
+
+   inline Debug() = default;
+
+   inline Debug(const int RANK, const char *FILE, const int LINE,
+                const char *FUNC, uint8_t COLOR, bool ini = true,
+                bool END = true): debug(true), end(END)
+   {
+      const char *base = static_strrnchr(FILE, '/', 2);
+      const char *file = base ? base + 1 : FILE;
+      const uint32_t rgb = static_color(COLOR, RANK, FILE);
+      const uint8_t r = (rgb >> 16) & 0xFF, g = (rgb >> 8) & 0xFF,
+                    b = rgb & 0xFF;
+      std::cout << "\033[38;2;";
+      std::cout << std::to_string(r) << ";";
+      std::cout << std::to_string(g) << ";";
+      std::cout << std::to_string(b) << "m";
+      if (ini)
+      {
+         std::cout << RANK << std::setw(64) << file << ":";
+         std::cout << "\033[2m" << std::setw(4) << std::left << LINE
+                   << "\033[22m: ";
+         if (FUNC) { std::cout << "[" << FUNC << "] "; }
+      }
+      std::cout << std::right << "\033[1m";
+   }
+
+   inline ~Debug()
+   {
+      if (debug) { std::cout << "\033[m" << (end ? "\n" : "") << std::flush; }
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept
+   {
+      if (debug) { std::cout << arg; }
+   }
+
+   template <typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (debug) { this->operator<<(arg); }
+   }
+
+   template <typename... Args>
+   inline void operator()(const char *fmt, Args &&...args) const noexcept
+   {
+      if (debug) { std::cout << fmt::format(fmt, std::forward<Args>(args)...); }
+   }
+
+   inline void operator()() const noexcept {}
+
+   static Debug Set(const char *FILE, const int LINE, const char *FUNC,
+                    uint8_t COLOR, bool INI = true, bool END = true)
+   {
+      static int mpi_rank = 0, dbg_mpi_rank = 0;
+      static bool env_mpi = false, env_dbg = false;
+      if (static bool ini = false; !std::exchange(ini, true))
+      {
+         env_dbg = (getenv("MFEM_DEBUG") != nullptr);
+         env_mpi = getenv("MFEM_DEBUG_MPI") != nullptr;
+         // int mpi_flag = 0;
+         // MPI_Initialized(&mpi_flag);
+         // if (mpi_flag) { MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); }
+         dbg_mpi_rank = atoi(env_mpi ? getenv("MFEM_DEBUG_MPI") : "0");
+      }
+      const bool debug = (env_dbg && (!env_mpi || (dbg_mpi_rank == mpi_rank)));
+      return debug ? Debug(mpi_rank, FILE, LINE, FUNC, COLOR, INI, END)
+             : Debug();
+   }
+};
+
+// Debug console traces, unnamed
+#define NVTX_DEBUG(...)                                                  \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR).operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_INI(...)                                \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, true) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_APPEND(...)                                 \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, false) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_END(...)                                \
+   ::gpu::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, true, false) \
+      .operator()(__VA_ARGS__)
+
+///////////////////////////////////////////////////////////////////////////////
+struct Nvtx
+{
+   const bool nvtx = false, enforce_kernel_sync = false;
+   const char *base, *file;
+   const uint32_t color = kBlack;
+   mutable std::string ascii;
+   mutable nvtxEventAttributes_t event;
+   mutable bool pushed = false;
+
+   inline Nvtx() = default;
+
+   Nvtx(bool enforce_kernel_sync, const char *FILE, const int LINE,
+        const char *FUNC, uint8_t COLOR):
+      nvtx(true), enforce_kernel_sync(enforce_kernel_sync),
+      base(static_strrnchr(FILE, '/', 2)), file(base ? base + 1 : FILE),
+      color(COLOR), ascii(file), event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = static_color(COLOR, 0, FILE);
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+
+      ascii += ":";
+      ascii += std::to_string(LINE);
+      ascii += ":[";
+      ascii += FUNC;
+      ascii += "] ";
+
+      pushed = false;
+   }
+
+   explicit Nvtx(const char *title, uint8_t color = kWheat,
+                 bool enforce_kernel_sync = true):
+      nvtx(true), enforce_kernel_sync(enforce_kernel_sync), color(color),
+      ascii(title), event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = static_color(color, 0, "");
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   inline void operator()() const
+   {
+      if (!nvtx) { return; }
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename T>
+   inline void operator()(const T &arg) const
+   {
+      if (!nvtx) { return; }
+      this->operator<<(arg);
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename... Args>
+   inline void operator()(fmt::format_string<Args...> fmt, Args &&...args) const
+   {
+      if (!nvtx) { return; }
+      ascii += fmt::format(fmt, std::forward<Args>(args)...);
+      event.message.ascii = ascii.c_str();
+      assert(!pushed);
+      nvtxRangePushEx(&event);
+      pushed = true;
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const
+   {
+      if (nvtx) { ascii += arg; }
+   }
+
+   inline ~Nvtx()
+   {
+      if (!nvtx) { return; }
+      if (enforce_kernel_sync)
+      {
+         nvtxEventAttributes_t eks = {};
+         eks.version = NVTX_VERSION;
+         eks.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+         eks.category = 0; // user value
+         eks.colorType = NVTX_COLOR_ARGB;
+         eks.messageType = NVTX_MESSAGE_TYPE_ASCII;
+         eks.message.ascii = "!"; // enforce kernel synchronization
+         eks.color = kHexColors[kYellow];
+         nvtxRangePushEx(&eks);
+         cudaStreamSynchronize(nullptr);
+         nvtxRangePop(/*eks*/);
+      }
+      assert(pushed);
+      nvtxRangePop(/*event*/);
+   }
+
+   using nvtx_ptr = std::unique_ptr<Nvtx>;
+   using nvtx_stack_t = std::stack<nvtx_ptr>;
+
+   static nvtx_ptr Set(const char *FILE, const int LINE, const char *FUNC,
+                       uint8_t COLOR)
+   {
+      static bool nvtx = false, eks = false;
+      if (static bool ini = false; !std::exchange(ini, true))
+      {
+         eks = getenv("MFEM_EKS") != nullptr;
+         nvtx = getenv("MFEM_NVTX") != nullptr;
+         Nvtx force_first_eks("Init EKS", kYellow, true);
+      }
+      return nvtx_ptr(nvtx ? new Nvtx(eks, FILE, LINE, FUNC, COLOR)
+                      : new Nvtx());
+   }
+
+   static nvtx_stack_t &Stack()
+   {
+      auto nvtx_events = []() -> nvtx_stack_t &
+      {
+         static nvtx_stack_t events;
+         return events;
+      };
+      static std::once_flag ready;
+      // one touch to guarantee the object is ready
+      std::call_once(ready, [&] { nvtx_events(); });
+      return nvtx_events();
+   }
+};
+
+// Temporary object only alive for the current statement
+#define NVTX_(COLOR, ...)                                       \
+   NVTX_DEBUG(__VA_ARGS__);                                     \
+   std::unique_ptr<::gpu::nvtx::Nvtx> NVTX_PRIVATE_NAME(nvtx) = \
+      ::gpu::nvtx::Nvtx::Set(NVTX_FLF, COLOR);                  \
+   NVTX_PRIVATE_NAME(nvtx)->operator()(__VA_ARGS__)
+
+// Temporary object only alive for the current statement
+#define NVTX(...) NVTX_(NVTX_COLOR, __VA_ARGS__)
+
+// Begin(with color)/End NVTX event traces
+#define NVTX_BEGIN_(COLOR, ...)                                              \
+   NVTX_DEBUG(__VA_ARGS__);                                                  \
+   ::gpu::nvtx::Nvtx::Stack().push(::gpu::nvtx::Nvtx::Set(NVTX_FLF, COLOR)); \
+   ::gpu::nvtx::Nvtx::Stack().top()->operator()(__VA_ARGS__)
+
+// Begin/End NVTX event traces
+#define NVTX_BEGIN(...) NVTX_BEGIN_(NVTX_COLOR, __VA_ARGS__);
+
+#define NVTX_END(...)                        \
+   ::gpu::nvtx::Nvtx::Stack().top().reset(); \
+   ::gpu::nvtx::Nvtx::Stack().pop()
+
+#ifdef USE_CALIPER
+// CALIPER & NVTX marks
+#define NVTX_MARK_FUNCTION                                            \
+      NVTX();                                                            \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;              \
+      __cali_ann##__func__ = std::make_unique<cali::Function>(__func__);
+
+#define NVTX_MARK(...)                                                   \
+      NVTX(__VA_ARGS__);                                                    \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;                 \
+      __cali_ann##__func__ = std::make_unique<cali::Function>(__VA_ARGS__);
+
+#define NVTX_MARK_FUNCTION_NAME(STR_NAME)                                \
+      NVTX(STR_NAME);                                                       \
+      std::unique_ptr<cali::Function> __cali_ann##__func__;                 \
+      if (g_caliper)                                                        \
+      {                                                                     \
+         __cali_ann##__func__ = std::make_unique<cali::Function>(STR_NAME); \
+      }
+
+#define NVTX_MARK_BEGIN(...)     \
+      CALI_MARK_BEGIN(__VA_ARGS__); \
+      NVTX_BEGIN(__VA_ARGS__);
+
+#define NVTX_MARK_END(...)     \
+      NVTX_END(__VA_ARGS__);      \
+      CALI_MARK_END(__VA_ARGS__);
+#else
+#define NVTX_MARK_FUNCTION NVTX()
+#define NVTX_MARK(...) NVTX(__VA_ARGS__)
+#define NVTX_MARK_FUNCTION_NAME(...) NVTX(__VA_ARGS__)
+#define NVTX_MARK_BEGIN(...) NVTX_BEGIN(__VA_ARGS__)
+#define NVTX_MARK_END(...) NVTX_END(__VA_ARGS__)
+#endif
+
+} // namespace gpu::nvtx
+
+// Debug console traces, unnamed
+#if 1
+#define dbg(...) NVTX_DEBUG(__VA_ARGS__)
+#define dbl(...) NVTX_DEBUG_NO_END(__VA_ARGS__)
+#define dba(...) NVTX_DEBUG_APPEND(__VA_ARGS__)
+#define dbc(...) NVTX_DEBUG_NO_INI(__VA_ARGS__)
+#else
+#define dbg(...)
+#define dbl(...) (void)0
+#define dba(...)
+#define dbc(...)
+#endif
+
+inline bool ClearScreen()
+{
+   dbg("\x1B[2J\x1B[3J\x1B[H");
+   return true;
+}

miniapps/navier/incompressible_navier_solver.cpp

@@ -0,0 +1,490 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "incompressible_navier_solver.hpp"
+
+#define NVTX_COLOR ::gpu::nvtx::kCyan
+#include "incompressible_navier_nvtx.hpp"
+
+using namespace mfem;
+using namespace incompressible_navier;
+
+IncompressibleNavierSolver::IncompressibleNavierSolver(ParMesh *mesh,
+                                                       int velorder, int porder,
+                                                       int torder,
+                                                       real_t kin_vis):
+   pmesh(mesh), velorder(velorder), porder(porder), torder(torder),
+   kin_vis(kin_vis), gll_rules(0, Quadrature1D::GaussLobatto),
+   vfec(new H1_FECollection(velorder, pmesh->Dimension())),
+   psifec(new H1_FECollection(porder)), pfec(new H1_FECollection(porder)),
+   vfes(new ParFiniteElementSpace(pmesh, vfec, pmesh->Dimension())),
+   psifes(new ParFiniteElementSpace(pmesh, pfec)),
+   pfes(new ParFiniteElementSpace(pmesh, pfec)),
+
+   velGF(torder + 1, nullptr), pGF(torder + 1, nullptr)
+{
+   NVTX();
+   // Check if fully periodic mesh
+   if (!(pmesh->bdr_attributes.Size() == 0))
+   {
+      vel_ess_attr.SetSize(pmesh->bdr_attributes.Max());
+      vel_ess_attr = 0;
+
+      pres_ess_attr.SetSize(pmesh->bdr_attributes.Max());
+      pres_ess_attr = 0;
+   }
+
+   for (int i = 0; i < torder + 1; i++)
+   {
+      velGF[i] = new ParGridFunction(vfes);
+      *velGF[i] = 0.0;
+      pGF[i] = new ParGridFunction(pfes);
+      *pGF[i] = 0.0;
+   }
+
+   psiGF.SetSpace(psifes);
+   DvGF.SetSpace(vfes);
+   divVelGF.SetSpace(pfes);
+   pRHS.SetSpace(pfes);
+}
+
+void IncompressibleNavierSolver::Setup(real_t dt)
+{
+   if (verbose && pmesh->GetMyRank() == 0)
+   {
+      mfem::out << "Setup" << std::endl;
+      if (partial_assembly)
+      {
+         mfem::out << "Using Partial Assembly" << std::endl;
+      }
+      else { mfem::out << "Using Full Assembly" << std::endl; }
+   }
+
+   this->Setup_velocity(dt);
+
+   this->Setup_auxiliary(dt);
+
+   this->Setup_pressure(dt);
+}
+
+void IncompressibleNavierSolver::Setup_velocity(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni =
+      gll_rules.Get(vfes->GetFE(0)->GetGeomType(), 2 * velorder - 1);
+
+   vfes->GetEssentialTrueDofs(vel_ess_attr, vel_ess_tdof);
+
+   //-------------------------------------------------------------------------
+
+   // Setup of coefficient for mass term of Eq(13)
+   dtCoeff = new ConstantCoefficient(1.0 / dt);
+   auto *vmass_blfi = new VectorMassIntegrator(*dtCoeff);
+
+   // Setup of coefficient for stiffness term of Eq(13)
+   kinvisCoeff = new ConstantCoefficient(kin_vis);
+   auto *vdiff_blfi = new VectorDiffusionIntegrator(*kinvisCoeff);
+
+   // setup of Bilinear form of Eq(13)
+   velBForm = new ParBilinearForm(vfes);
+   if (numerical_integ)
+   {
+      vmass_blfi->SetIntRule(&ir_ni);
+      vdiff_blfi->SetIntRule(&ir_ni);
+   }
+   velBForm->AddDomainIntegrator(vmass_blfi);
+   velBForm->AddDomainIntegrator(vdiff_blfi);
+   if (partial_assembly) { velBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+
+   velBForm->Assemble();
+   velBForm->FormSystemMatrix(vel_ess_tdof, vOp);
+
+   //-------------------------------------------------------------------------
+
+   // Setup of coefficient for Eq(18)
+   pUnitVectorCoeff = new UnitVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *pvel_lfi = new VectorDomainLFGradIntegrator(*pUnitVectorCoeff);
+
+   // Setup of coefficient for Eq(20)
+   nonlinTermCoeff = new NonLinTermVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *p_nonlintermlfi = new VectorDomainLFIntegrator(*nonlinTermCoeff);
+
+   // Setup of coefficient for Eq(21)
+   prevVelLoadCoeff = new PrevVelVectorGridFunctionCoeff(pmesh->Dimension());
+   auto *prevVelLoadLFi = new VectorDomainLFIntegrator(*prevVelLoadCoeff);
+
+   // Setup of linear form of Eq(13)
+   velLForm = new ParLinearForm(vfes);
+   if (numerical_integ)
+   {
+      prevVelLoadLFi->SetIntRule(&ir_ni);
+      pvel_lfi->SetIntRule(&ir_ni);
+      p_nonlintermlfi->SetIntRule(&ir_ni);
+   }
+   velLForm->AddDomainIntegrator(prevVelLoadLFi);
+   velLForm->AddDomainIntegrator(pvel_lfi);
+   velLForm->AddDomainIntegrator(p_nonlintermlfi);
+
+   //-------------------------------------------------------------------------
+
+   if (partial_assembly)
+   {
+      Vector diag_pa(vfes->GetTrueVSize());
+      velBForm->AssembleDiagonal(diag_pa);
+      velInvPC = new OperatorJacobiSmoother(diag_pa, vel_ess_tdof);
+   }
+   else
+   {
+      velInvPC = new HypreSmoother(*vOp.As<HypreParMatrix>());
+      dynamic_cast<HypreSmoother *>(velInvPC)->SetType(HypreSmoother::Jacobi,
+                                                       1);
+   }
+
+   velInv = new CGSolver(vfes->GetComm());
+   velInv->iterative_mode = true;
+   velInv->SetOperator(*vOp);
+   velInv->SetPreconditioner(*velInvPC);
+   velInv->SetPrintLevel(pl_velsolve);
+   velInv->SetRelTol(rtol_velsolve);
+   velInv->SetAbsTol(0.0);
+   velInv->SetMaxIter(1200);
+}
+
+void IncompressibleNavierSolver::Setup_auxiliary(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni =
+      gll_rules.Get(vfes->GetFE(0)->GetGeomType(), 2 * velorder - 1);
+   Array<int> empty;
+
+   // setup of Bilinear form of Eq(14)
+   psiBForm = new ParBilinearForm(psifes);
+   auto *psidiff_blfi = new DiffusionIntegrator;
+
+   if (numerical_integ) { psidiff_blfi->SetIntRule(&ir_ni); }
+   psiBForm->AddDomainIntegrator(psidiff_blfi);
+   if (partial_assembly) { psiBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+
+   psiBForm->Assemble();
+   psiBForm->FormSystemMatrix(empty, psiOp);
+
+   //-------------------------------------------------------------------------
+
+   // Setup of coefficient for linear form in Eq(14)
+   DvelCoeff = new VectorGridFunctionCoefficient;
+   auto *Dvel_lfi = new DomainLFGradIntegrator(*DvelCoeff);
+
+   // Setup of linear form of Eq(14)
+   psiLForm = new ParLinearForm(psifes);
+
+   if (numerical_integ) { Dvel_lfi->SetIntRule(&ir_ni); }
+   psiLForm->AddDomainIntegrator(Dvel_lfi);
+
+   //-------------------------------------------------------------------------
+   if (partial_assembly)
+   {
+      int psifes_truevsize = psifes->GetTrueVSize();
+      mfem::Vector psin(psifes_truevsize);
+      psin = 0.0;
+      mfem::Vector respsi(psifes_truevsize);
+      respsi = 0.0;
+
+      lor = new ParLORDiscretization(*psiBForm, empty);
+      psiInvPC = new HypreBoomerAMG(lor->GetAssembledMatrix());
+      psiInvPC->SetPrintLevel(0);
+      psiInvPC->Mult(respsi, psin);
+      SpInvOrthoPC = new OrthoSolver(psifes->GetComm());
+      SpInvOrthoPC->SetSolver(*psiInvPC);
+   }
+   else
+   {
+      psiInvPC = new HypreBoomerAMG(*psiOp.As<HypreParMatrix>());
+      psiInvPC->SetPrintLevel(0);
+      SpInvOrthoPC = new OrthoSolver(psifes->GetComm());
+      SpInvOrthoPC->SetSolver(*psiInvPC);
+   }
+
+   psiInv = new CGSolver(psifes->GetComm());
+   psiInv->iterative_mode = true;
+   psiInv->SetOperator(*psiOp);
+   psiInv->SetPreconditioner(*SpInvOrthoPC);
+   psiInv->SetPrintLevel(pl_psisolve);
+   psiInv->SetRelTol(rtol_psisolve);
+   psiInv->SetAbsTol(0.0);
+   psiInv->SetMaxIter(1000);
+}
+
+void IncompressibleNavierSolver::Setup_pressure(real_t dt)
+{
+   // GLL integration rule (Numerical Integration)
+   const IntegrationRule &ir_ni =
+      gll_rules.Get(vfes->GetFE(0)->GetGeomType(), 2 * velorder - 1);
+   Array<int> empty;
+
+   //-------------------------------------------------------------------------
+
+   // setup of Bilinear form of Eq(15)
+   pBForm = new ParBilinearForm(pfes);
+   auto *pmass_blfi = new MassIntegrator;
+
+   if (numerical_integ) { pmass_blfi->SetIntRule(&ir_ni); }
+   pBForm->AddDomainIntegrator(pmass_blfi);
+   if (partial_assembly) { pBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+
+   pBForm->Assemble();
+   pBForm->FormSystemMatrix(empty, pOp);
+
+   //-------------------------------------------------------------------------
+
+   // Setup of divergence of velocity coefficient for linear form in Eq(15)
+   divVelCoeff = new DivergenceGridFunctionCoefficient(velGF[0]);
+
+   // Setup of coefficient for linear form in Eq(14)
+   pRHSCoeff = new GridFunctionCoefficient(&pRHS);
+   auto *p_lfi = new DomainLFIntegrator(*pRHSCoeff);
+
+   // Setup of linear form of Eq(15)
+   pLForm = new ParLinearForm(pfes);
+   if (numerical_integ) { p_lfi->SetIntRule(&ir_ni); }
+   pLForm->AddDomainIntegrator(p_lfi);
+
+   //-------------------------------------------------------------------------
+   if (partial_assembly)
+   {
+      Vector diag_pa(pfes->GetTrueVSize());
+      pBForm->AssembleDiagonal(diag_pa);
+      pInvPC = new OperatorJacobiSmoother(diag_pa, empty);
+   }
+   else
+   {
+      pInvPC = new HypreSmoother(*pOp.As<HypreParMatrix>());
+      dynamic_cast<HypreSmoother *>(pInvPC)->SetType(HypreSmoother::Jacobi, 1);
+   }
+
+   pInv = new CGSolver(pfes->GetComm());
+   pInv->iterative_mode = true;
+   pInv->SetOperator(*pOp);
+   pInv->SetPreconditioner(*pInvPC);
+   pInv->SetPrintLevel(pl_psolve);
+   pInv->SetRelTol(rtol_psolve);
+   pInv->SetAbsTol(0.0);
+   pInv->SetMaxIter(1000);
+}
+
+void IncompressibleNavierSolver::UpdateTimestepHistory(real_t dt) {}
+
+void IncompressibleNavierSolver::Step(real_t &time, real_t dt, int current_step,
+                                      const bool vis_step)
+{
+   this->Step_velocity(time, dt, current_step);
+
+   this->Step_auxiliary(time, dt, current_step);
+
+   this->Step_pressure(time, dt, current_step);
+
+   *velGF[1] = *velGF[0];
+   *pGF[1] = *pGF[0];
+
+   if (vis_step)
+   {
+      mfem::out << "It: " << iter << " | Iter_U: " << iter_vsolve
+                << " | Iter_Psi: " << iter_psisolve
+                << " | Iter_P: " << iter_psolve << "\n";
+      mfem::out << "It: " << iter << " | Resid_U: " << res_vsolve
+                << " | Resid_Psi: " << res_psisolve
+                << " | Resid_P: " << res_psisolve << "\n";
+   }
+
+   time += dt;
+   iter++;
+}
+
+void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt,
+                                               int current_step)
+{
+   for (auto &vel_dbc : vel_dbcs)
+   {
+      velGF[0]->ProjectBdrCoefficient(*vel_dbc.coeff, vel_dbc.attr);
+      velGF[1]->ProjectBdrCoefficient(*vel_dbc.coeff, vel_dbc.attr);
+   }
+
+   // Update state in coefficient for Eq(18)
+   pUnitVectorCoeff->SetGridFunction(pGF[1]);
+
+   // Update state in coefficient for Eq(20)
+   nonlinTermCoeff->SetGridFunction(velGF[1]);
+
+   // Update state in coefficient for Eq(21)
+   prevVelLoadCoeff->SetGridFunction(velGF[1], dt);
+
+   velLForm->Assemble();
+   velLForm->ParallelAssemble(velLF);
+
+   Vector X1, B1;
+
+   if (partial_assembly)
+   {
+      auto *vpC = vOp.As<ConstrainedOperator>();
+      EliminateRHS(*velBForm, *vpC, vel_ess_tdof, *velGF[0], velLF, X1, B1, 1);
+   }
+   else
+   {
+      velBForm->FormLinearSystem(vel_ess_tdof, *velGF[0], velLF, vOp, X1, B1,
+                                 1);
+   }
+
+   velInv->Mult(B1, X1);
+   iter_vsolve = velInv->GetNumIterations();
+   res_vsolve = velInv->GetFinalNorm();
+   velBForm->RecoverFEMSolution(X1, velLF, *velGF[0]);
+}
+
+void IncompressibleNavierSolver::Step_auxiliary(real_t &time, real_t dt,
+                                                int current_step)
+{
+   // Compute new increment GF for LF of Eq(14)  and update state in coefficient
+   subtract(1.0 / dt, *velGF[0], *velGF[1], DvGF);
+   DvelCoeff->SetGridFunction(&DvGF);
+
+   psiLForm->Assemble();
+   psiLForm->ParallelAssemble(psiLF);
+
+   Vector X2, B2;
+   Array<int> empty;
+   if (partial_assembly)
+   {
+      auto *psipC = psiOp.As<ConstrainedOperator>();
+      EliminateRHS(*psiBForm, *psipC, empty, psiGF, psiLF, X2, B2, 1);
+   }
+   else { psiBForm->FormLinearSystem(empty, psiGF, psiLF, psiOp, X2, B2, 1); }
+
+   psiInv->Mult(B2, X2);
+   iter_psisolve = psiInv->GetNumIterations();
+   res_psisolve = psiInv->GetFinalNorm();
+   psiBForm->RecoverFEMSolution(X2, psiLF, psiGF);
+}
+
+void IncompressibleNavierSolver::Step_pressure(real_t &time, real_t dt,
+                                               int current_step)
+{
+   Array<int> empty;
+
+   // Compute new GF for LF of Eq(15) and update state in coefficient
+   divVelCoeff->SetGridFunction(velGF[0]);
+   divVelGF.ProjectCoefficient(*divVelCoeff);
+
+   add(*pGF[1], psiGF, pRHS);
+   add(pRHS, -1.0 * kin_vis, divVelGF, pRHS);
+   pRHSCoeff->SetGridFunction(&pRHS);
+
+   pLForm->Assemble();
+   pLForm->ParallelAssemble(pLF);
+
+   Vector X3, B3;
+
+   if (partial_assembly)
+   {
+      auto *ppC = pOp.As<ConstrainedOperator>();
+      EliminateRHS(*pBForm, *ppC, empty, *pGF[0], pLF, X3, B3, 1);
+   }
+   else { pBForm->FormLinearSystem(empty, *pGF[0], pLF, pOp, X3, B3, 1); }
+
+   pInv->Mult(B3, X3);
+   iter_psolve = pInv->GetNumIterations();
+   res_psisolve = pInv->GetFinalNorm();
+   pBForm->RecoverFEMSolution(X3, pLF, *pGF[0]);
+}
+
+void IncompressibleNavierSolver::EliminateRHS(Operator &A,
+                                              ConstrainedOperator &constrainedA,
+                                              const Array<int> &ess_tdof_list,
+                                              Vector &x, Vector &b, Vector &X,
+                                              Vector &B, int copy_interior)
+{
+   const Operator *Po = A.GetOutputProlongation();
+   const Operator *Pi = A.GetProlongation();
+   const Operator *Ri = A.GetRestriction();
+   A.InitTVectors(Po, Ri, Pi, x, b, X, B);
+   if (!copy_interior) { X.SetSubVectorComplement(ess_tdof_list, 0.0); }
+   constrainedA.EliminateRHS(X, B);
+}
+
+real_t IncompressibleNavierSolver::ComputeCFL(ParGridFunction &u, real_t dt)
+{
+   return 0;
+}
+
+void IncompressibleNavierSolver::AddVelDirichletBC(VectorCoefficient *coeff,
+                                                   Array<int> &attr)
+{
+   vel_dbcs.emplace_back(attr, coeff);
+
+   if (verbose && pmesh->GetMyRank() == 0)
+   {
+      mfem::out << "Adding Velocity Dirichlet BC to attributes ";
+      for (int i = 0; i < attr.Size(); ++i)
+      {
+         if (attr[i] == 1) { mfem::out << i << " "; }
+      }
+      mfem::out << std::endl;
+   }
+
+   for (int i = 0; i < attr.Size(); ++i)
+   {
+      MFEM_ASSERT((vel_ess_attr[i] && attr[i]) == 0,
+                  "Duplicate boundary definition deteceted.");
+      if (attr[i] == 1) { vel_ess_attr[i] = 1; }
+   }
+}
+
+void IncompressibleNavierSolver::AddVelDirichletBC(VecFuncT *f,
+                                                   Array<int> &attr)
+{
+   AddVelDirichletBC(new VectorFunctionCoefficient(pmesh->Dimension(), f),
+                     attr);
+}
+
+IncompressibleNavierSolver::~IncompressibleNavierSolver()
+{
+   delete velBForm;
+   delete psiBForm;
+   delete pBForm;
+
+   delete kinvisCoeff;
+   delete dtCoeff;
+
+   for (int i = 0; i < torder + 1; i++)
+   {
+      delete velGF[i];
+      delete pGF[i];
+   }
+
+   delete DvelCoeff;
+   delete divVelCoeff;
+   delete pRHSCoeff;
+   delete pUnitVectorCoeff;
+
+   delete velInv;
+   delete velInvPC;
+   delete psiInv;
+   delete SpInvOrthoPC;
+   delete psiInvPC;
+   delete lor;
+   delete pInv;
+   delete pInvPC;
+
+   delete vfec;
+   delete psifec;
+   delete pfec;
+   delete vfes;
+   delete psifes;
+   delete pfes;
+}

miniapps/navier/incompressible_navier_solver.hpp

@@ -0,0 +1,338 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#pragma once
+
+#define INCOMP_NAVIER_VERSION 0.1
+
+#include "mfem.hpp"
+
+namespace mfem::incompressible_navier
+{
+
+using VecFuncT = void(const Vector &x, real_t t, Vector &u);
+using ScalarFuncT = real_t(const Vector &x, real_t t);
+
+// Coefficient which computed contribution of Eq(18)
+class UnitVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   UnitVectorGridFunctionCoeff(int dim): VectorCoefficient(dim * dim) {}
+
+   void Eval(Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip) override
+   {
+      real_t coeffVal = gridfunc_->GetValue(T, ip);
+
+      V.SetSize(vdim);
+      V = 0.0; // FIXME
+      V[0] = coeffVal;
+      V[3] = coeffVal;
+   }
+
+   void SetGridFunction(GridFunction *gridfunc) { gridfunc_ = gridfunc; }
+
+   GridFunction *gridfunc_ = nullptr;
+};
+
+// Coefficient which computed contribution of Eq(21)
+class PrevVelVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   PrevVelVectorGridFunctionCoeff(int dim): VectorCoefficient(dim) {}
+
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      V.SetSize(vdim);
+      gridFuncCoeff->Eval(V, T, ip);
+
+      V *= 1.0 / dt_;
+   }
+
+   void SetGridFunction(GridFunction *gridfunc, real_t dt)
+   {
+      gridfunc_ = gridfunc;
+      dt_ = dt;
+      delete gridFuncCoeff;
+      gridFuncCoeff = new VectorGridFunctionCoefficient(gridfunc);
+   }
+
+   GridFunction *gridfunc_ = nullptr;
+   VectorGridFunctionCoefficient *gridFuncCoeff = nullptr;
+   real_t dt_;
+};
+
+// Coefficient which computed contribution of Eq(20)
+class NonLinTermVectorGridFunctionCoeff : public VectorCoefficient
+{
+public:
+   NonLinTermVectorGridFunctionCoeff(int dim): VectorCoefficient(dim) {}
+
+   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   {
+      Vector val(vdim);
+      Vector resultVal(vdim);
+      DenseMatrix vecGrad;
+      V.SetSize(vdim);
+      gridFuncCoeff->Eval(val, T, ip);
+
+      gridfunc_->GetVectorGradient(T, vecGrad);
+
+      vecGrad.MultTranspose(val, V);
+
+      V *= -1.0;
+   }
+
+   void SetGridFunction(ParGridFunction *gridfunc)
+   {
+      delete gridFuncCoeff;
+      gridfunc_ = gridfunc;
+      gridFuncCoeff = new VectorGridFunctionCoefficient(gridfunc);
+   }
+
+   VectorGridFunctionCoefficient *gridFuncCoeff = nullptr;
+   ParGridFunction *gridfunc_ = nullptr;
+};
+
+/// Container for a Dirichlet boundary condition of the velocity field.
+class VelDirichletBC_T
+{
+public:
+   VelDirichletBC_T(Array<int> attr, VectorCoefficient *coeff):
+      attr(attr), coeff(coeff)
+   {
+   }
+
+   VelDirichletBC_T(VelDirichletBC_T &&obj)
+   {
+      // Deep copy the attribute array
+      this->attr = obj.attr;
+
+      // Move the coefficient pointer
+      this->coeff = obj.coeff;
+      obj.coeff = nullptr;
+   }
+
+   ~VelDirichletBC_T() { delete coeff; }
+
+   Array<int> attr;
+   VectorCoefficient *coeff;
+};
+
+/// Transient incompressible Navier Stokes solver in a split scheme formulation.
+/**
+ * This implementation of a transient incompressible Navier Stokes solver uses
+ * the non-dimensionalized formulation. The coupled momentum and
+ * incompressibility equations are decoupled using the split scheme described in
+ * [1]. This leads to three solving steps.
+ *
+ */
+class IncompressibleNavierSolver
+{
+public:
+   /// Initialize data structures, set FE space order and kinematic viscosity.
+   /**
+    * The ParMesh @a mesh can be a linear or curved parallel mesh. The @a order
+    * of the finite element spaces is
+    */
+   IncompressibleNavierSolver(ParMesh *mesh, int velorder, int porder,
+                              int tOrder, real_t kin_vis);
+
+   /// Initialize forms, solvers and preconditioners.
+   void Setup(real_t dt);
+
+   void Setup_velocity(real_t dt);
+
+   void Setup_auxiliary(real_t dt);
+
+   void Setup_pressure(real_t dt);
+
+   /// Compute solution at the next time step t+dt.
+   /**
+    * This method can
+    */
+   void Step(real_t &time, real_t dt, int cur_step, bool vis_step);
+
+   void Step_velocity(real_t &time, real_t dt, int cur_step);
+
+   void Step_auxiliary(real_t &time, real_t dt, int cur_step);
+
+   void Step_pressure(real_t &time, real_t dt, int cur_step);
+
+   /// Return a pointer to the provisional velocity ParGridFunction.
+   ParGridFunction *GetProvisionalVelocity() { return velGF[1]; }
+
+   /// Return a pointer to the current velocity ParGridFunction.
+   ParGridFunction *GetCurrentVelocity() { return velGF[0]; }
+
+   /// Return a pointer to the current pressure ParGridFunction.
+   ParGridFunction *GetCurrentPressure() { return pGF[0]; }
+
+   /// Return a pointer to the current pressure ParGridFunction.
+   ParGridFunction *GetCurrentPsi() { return &psiGF; }
+
+   /// Add a Dirichlet boundary condition to the velocity field.
+   void AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr);
+
+   void AddVelDirichletBC(VecFuncT *f, Array<int> &attr);
+
+   /// Add a Dirichlet boundary condition to the pressure field.
+   // void AddPresDirichletBC(Coefficient *coeff, Array<int> &attr);
+
+   // void AddPresDirichletBC(ScalarFuncT *f, Array<int> &attr);
+
+   /// Enable partial assembly for every operator.
+   void EnablePA(bool pa) { partial_assembly = pa; }
+
+   /// Enable numerical integration rules. This means collocated quadrature at
+   /// the nodal points.
+   void EnableNI(bool ni) { numerical_integ = ni; }
+
+   /// Print timing summary of the solving routine.
+   void PrintTimingData();
+
+   ~IncompressibleNavierSolver();
+
+   /// Rotate entries in the time step and solution history arrays.
+   void UpdateTimestepHistory(real_t dt);
+
+   /// Compute CFL
+   real_t ComputeCFL(ParGridFunction &u, real_t dt);
+
+protected:
+   /// Eliminate essential BCs in an Operator and apply to RHS.
+   void EliminateRHS(Operator &A, ConstrainedOperator &constrainedA,
+                     const Array<int> &ess_tdof_list, Vector &x, Vector &b,
+                     Vector &X, Vector &B, int copy_interior = 0);
+
+   /// Enable/disable debug output.
+   bool debug = false;
+
+   /// Enable/disable verbose output.
+   bool verbose = true;
+
+   /// Enable/disable partial assembly of forms.
+   bool partial_assembly = false;
+
+   /// Enable/disable numerical integration rules of forms.
+   bool numerical_integ = false;
+
+   /// The parallel mesh.
+   ParMesh *pmesh = nullptr;
+
+   /// The order of the velocity and pressure space.
+   const int velorder;
+   const int porder;
+   const int torder;
+
+   /// Kinematic viscosity (dimensionless).
+   const real_t kin_vis;
+   Coefficient *kinvisCoeff = nullptr;
+
+   Coefficient *dtCoeff = nullptr;
+
+   IntegrationRules gll_rules;
+
+   /// Velocity $H^1$ finite element collection.
+   FiniteElementCollection *vfec = nullptr;
+
+   /// Psi $H^1$ finite element collection.
+   FiniteElementCollection *psifec = nullptr;
+
+   /// Pressure $H^1$ finite element collection.
+   FiniteElementCollection *pfec = nullptr;
+
+   /// Velocity $(H^1)^d$ finite element space.
+   ParFiniteElementSpace *vfes = nullptr;
+
+   /// Psi $(H^1)^d$ finite element space.
+   ParFiniteElementSpace *psifes = nullptr;
+
+   /// Pressure $H^1$ finite element space.
+   ParFiniteElementSpace *pfes = nullptr;
+
+   ParBilinearForm *velBForm = nullptr; // vmass + vdiff
+   ParBilinearForm *psiBForm = nullptr; // diffusion
+   ParBilinearForm *pBForm = nullptr;   // mass
+
+   ParLinearForm *velLForm = nullptr; // vLF + vLFGrad + vLF
+   ParLinearForm *psiLForm = nullptr; // LFGrad
+   ParLinearForm *pLForm = nullptr;   // LF
+
+   // current (0) and provisional (1) velocity
+   std::vector<ParGridFunction *> velGF;
+
+   // current (0) pressure ParGridFunction.
+   std::vector<ParGridFunction *> pGF;
+   ParGridFunction psiGF;
+
+   ParGridFunction DvGF, divVelGF, pRHS;
+   VectorGridFunctionCoefficient *DvelCoeff = nullptr;
+   DivergenceGridFunctionCoefficient *divVelCoeff = nullptr;
+   GridFunctionCoefficient *pRHSCoeff = nullptr;
+   UnitVectorGridFunctionCoeff *pUnitVectorCoeff = nullptr;
+   NonLinTermVectorGridFunctionCoeff *nonlinTermCoeff = nullptr;
+   PrevVelVectorGridFunctionCoeff *prevVelLoadCoeff = nullptr;
+
+   OperatorHandle vOp;
+   OperatorHandle psiOp;
+   OperatorHandle pOp;
+
+   Solver *velInvPC = nullptr;
+   CGSolver *velInv = nullptr;
+
+   ParLORDiscretization *lor = nullptr;
+   HypreBoomerAMG *psiInvPC = nullptr;
+   OrthoSolver *SpInvOrthoPC = nullptr;
+   CGSolver *psiInv = nullptr;
+
+   Solver *pInvPC = nullptr;
+   CGSolver *pInv = nullptr;
+
+   Vector velLF, psiLF, pLF;
+
+   // All essential attributes.
+   Array<int> vel_ess_attr;
+   Array<int> pres_ess_attr;
+
+   // All essential true dofs.
+   Array<int> vel_ess_tdof;
+   Array<int> pres_ess_tdof;
+
+   // Bookkeeping for velocity dirichlet bcs.
+   std::vector<VelDirichletBC_T> vel_dbcs;
+
+   // Print levels.
+   int pl_psolve = 0;
+   int pl_psisolve = 0;
+   int pl_velsolve = 0;
+   int pl_amg = 0;
+
+#if defined(MFEM_USE_DOUBLE)
+   real_t rtol_psolve = 1e-12;
+   real_t rtol_psisolve = 1e-12;
+   real_t rtol_velsolve = 1e-12;
+#elif defined(MFEM_USE_SINGLE)
+   real_t rtol_psolve = 1e-9;
+   real_t rtol_psisolve = 1e-5;
+   real_t rtol_velsolve = 1e-7;
+#else
+   #error "Only single and double precision are supported!"
+#endif
+
+   // Iteration counts.
+   int iter = 1, iter_vsolve = 0, iter_psolve = 0, iter_psisolve = 0;
+
+   // Residuals.
+   real_t res_vsolve = 0.0, res_psolve = 0.0, res_psisolve = 0.0;
+};
+
+} // namespace mfem::incompressible_navier

miniapps/navier/incompressible_navier_tests.cpp

@@ -0,0 +1,206 @@
+#include <algorithm>
+#include <iostream>
+#include <memory>
+#include <sstream>
+#include <unistd.h>
+
+#define NVTX_COLOR ::gpu::nvtx::kMagenta
+#include "incompressible_navier_nvtx.hpp"
+
+///////////////////////////////////////////////////////////////////////////////
+int navier(int argc, char *argv[], double &u, double &p, double &Ψ);
+
+///////////////////////////////////////////////////////////////////////////////
+template <class T>
+std::enable_if_t<!std::numeric_limits<T>::is_integer, bool>
+AlmostEq(T x, T y, T tolerance = 100.0 * std::numeric_limits<T>::epsilon())
+{
+   const T neg = std::abs(x - y);
+   constexpr T min = std::numeric_limits<T>::min();
+   constexpr T eps = std::numeric_limits<T>::epsilon();
+   const T min_abs = std::min(std::abs(x), std::abs(y));
+   if (std::abs(min_abs) == 0.0) { return neg < eps; }
+   return (neg / (1.0 + std::max(min, min_abs))) < tolerance;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+using char_uptr = std::unique_ptr<char[]>;
+using args_ptr_t = std::vector<char_uptr>;
+using args_t = std::vector<char *>;
+
+///////////////////////////////////////////////////////////////////////////////
+struct Results
+{
+   double u{}, p{}, Ψ {};
+};
+
+///////////////////////////////////////////////////////////////////////////////
+struct Test
+{
+   static constexpr const char *binary = "incompNS_2Dtest ";
+   static constexpr const char *common = "-no-vis -no-pv";
+   const std::string options;
+   const Results results;
+   Test(const char *args, const Results &res):
+      options(std::string(args) + " " + common), results(res)
+   {
+      dbg("options: {}", options.c_str());
+      dbg("results: U={:.15e}, P={:.15e}, Ψ={:.15e}",
+          results.u, results.p, results.Ψ);
+   }
+   std::string Command() const { return binary + options; }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+#if 1 // dot product reduction (miniapps/navier/incompNS_2Dtest.cpp#L182)
+static const Test gold[] =
+{
+   {
+      "-nx 9 -ny 3 -sr 0",
+      { 3.056430866716070e-06, 1.504027632950462e-01, 7.132601171183242e-08 }
+   },
+   {
+      "-nx 16 -ny 8 -sr 0",
+      { 1.409287729554512e-05, 5.718053801962010e-01, 1.904938419441012e-07 }
+   },
+   // {
+   //    "-nx 9 -ny 3 -sr 1",
+   //    { 1.23258426138828e-05, 5.18207619597952e-01, 1.956175418199867e-07 }
+   // },
+   // {
+   //    "-nx 9 -ny 3 -sr 2",
+   //    { 4.74381190869396e-05, 1.80653923294262e+00, 5.90778654333642e-07 }
+   // },
+};
+#else // Norml2 reduction
+///////////////////////////////////////////////////////////////////////////////
+const Test runs[] =
+{
+   {
+      "-nx 9 -ny 3 -sr 0",
+      { 1.746844586767688e-03, 3.874421956807944e-01, 2.670296315544763e-04 }
+   },
+   // 1.748265101955712e-03, 3.878179512284757e-01, 2.670693013280680e-04 //
+   // Release { "-nx 9 -ny 3 -sr 1",
+   //  { 1.232584261388279e-05, 5.182076195979519e-01, 1.956175418199866e-07 }
+   //  },
+   // { "-nx 9 -ny 3 -sr 2",
+   //  { 4.743811908693962e-05, 1.806539232942622e+00, 5.907786543336419e-07 }
+   //  },
+};
+#endif
+
+///////////////////////////////////////////////////////////////////////////////
+int NavierTest(const int k, const Test &run)
+{
+   dbg();
+   static args_ptr_t args_ptr;
+   args_t args;
+
+   std::istringstream iss(run.Command());
+
+   auto add_arg = [&](std::string token) -> char_uptr
+   {
+      auto arg_ptr = std::make_unique<char[]>(token.size() + 1);
+      std::memcpy(arg_ptr.get(), token.c_str(), token.size() + 1);
+      arg_ptr[token.size()] = '\0';
+      return arg_ptr;
+   };
+
+   std::string token;
+   while (iss >> token)
+   {
+      auto arg_ptr = add_arg(token);
+      args.push_back(arg_ptr.get());
+      args_ptr.emplace_back(std::move(arg_ptr));
+   }
+   args.push_back(nullptr);
+
+   auto launch = [&args, &run, &k]() -> int
+   {
+      // dbg("Launching test #{}: \x1B[33m{}\x1B[m", k, gold.Command().c_str());
+      Results res{};
+      navier(args.size() - 1, args.data(), res.u, res.p, res.Ψ);
+      // dbg("Results: U={:.15e}, P={:.15e}, Ψ={:.15e}", res.u, res.p, res.Ψ);
+
+      const bool u = AlmostEq(res.u, run.results.u);
+      const bool p = AlmostEq(res.p, run.results.p);
+      const bool Ψ = AlmostEq(res.Ψ, run.results.Ψ);
+
+      constexpr auto ok = [](bool ok) -> int { return ok ? 32 : 31; };
+      constexpr auto to_string = [](args_t &args) -> std::string
+      {
+         std::string args_str;
+         for (auto &arg : args)
+         {
+            if (!arg) { break; }
+            args_str += std::string(arg) + " ";
+         }
+         return args_str;
+      };
+
+      dbg("#{} \x1B[33m{}\x1B[m", k, to_string(args).c_str());
+      dbg("U: \x1B[33m{:.15e} \x1B[{}m{:.15e}", run.results.u, ok(u), res.u );
+      dbg("P: \x1B[33m{:.15e} \x1B[{}m{:.15e}", run.results.p, ok(p), res.p);
+      dbg("Ψ: \x1B[33m{:.15e} \x1B[{}m{:.15e}", run.results.Ψ, ok(Ψ), res.Ψ);
+
+      if (u && p && Ψ) { return std::cout << "✅" << std::endl, EXIT_SUCCESS; }
+      else { return std::cout << "❌" << std::endl, EXIT_FAILURE; }
+   };
+
+   // first launch with default arguments
+   if (launch() != EXIT_SUCCESS) { return EXIT_FAILURE; }
+
+   // second launch with the same arguments, but with -pa
+   args.pop_back(); // nullptr
+   auto arg_pa_ptr = add_arg("-pa");
+   args.push_back(arg_pa_ptr.get());
+   args_ptr.emplace_back(std::move(arg_pa_ptr));
+   args.push_back(nullptr);
+   if (launch() != EXIT_SUCCESS) { return EXIT_FAILURE; }
+
+   return EXIT_SUCCESS;
+}
+
+///////////////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+try
+{
+   dbg();
+
+   int opt;
+   int test = -1;
+   auto show_usage = [](const int ret = EXIT_FAILURE)
+   {
+      printf("Usage: program [-a <arg>] [-b <arg>] [-h]\n");
+      printf("  -t <test>  Optional test number \n");
+      printf("  -h         Show this help message\n");
+      exit(ret);
+   };
+
+   while ((opt = getopt(argc, argv, "t:h")) != -1)
+   {
+      switch (opt)
+      {
+         case 't': test = std::atoi(optarg); break;
+         case 'h': show_usage(EXIT_SUCCESS);
+         default:  show_usage(EXIT_FAILURE);
+      }
+   }
+
+   constexpr int N_TESTS = sizeof(gold) / sizeof(Test);
+   if (test >= 0 && test < N_TESTS) { return NavierTest(test, gold[test]); }
+
+   int k = 0;
+   for (auto &run : gold)
+   {
+      if (NavierTest(k++, run) != EXIT_SUCCESS) { return EXIT_FAILURE; }
+   }
+   return EXIT_SUCCESS;
+}
+catch (std::exception &e)
+{
+   std::cerr << "\033[31m..xxxXXX[ERROR]XXXxxx.." << std::endl;
+   std::cerr << "\033[31m{}" << e.what() << std::endl;
+   return EXIT_FAILURE;
+}