Merge branch 'main' into hpcftools/flowSolver-gpu

.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

@@ -11,14 +11,9 @@
 
 if (MFEM_USE_MPI)
   list(APPEND NAVIER_COMMON_SOURCES
-      navier_solver.cpp
-      incompressible_navier_solver.cpp
-      stokes_solver.cpp)
-
+      navier_solver.cpp)
   list(APPEND NAVIER_COMMON_HEADERS
-      navier_solver.hpp
-      incompressible_navier_solver.hpp
-      stokes_solver.hpp)
+      navier_solver.hpp)
 
   convert_filenames_to_full_paths(NAVIER_COMMON_SOURCES)
   convert_filenames_to_full_paths(NAVIER_COMMON_HEADERS)
@@ -57,9 +52,18 @@ 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
-    ${NAVIER_COMMON_FILES}
+    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
@@ -94,4 +98,4 @@ if (MFEM_USE_MPI)
         ${MPIEXEC_POSTFLAGS})
     endforeach()
 endif()
-endif ()
+endif()

miniapps/navier/incompNS_2Dtest.cpp

@@ -12,50 +12,93 @@
 // 3D flow over a cylinder benchmark example
 
 #include "incompressible_navier_solver.hpp"
-#include <fstream>
+
+#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);
-   real_t yi = x(1);
-
+   // 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);
-    }
+   u = 0.0;
+   if (x(0) < 0.001) { u(0) = -0.001 * (std::pow(x(1) - 0.5, 2.0) - 0.25); }
 }
 
-int main(int argc, char *argv[])
+MFEM_EXPORT int navier(int argc, char *argv[], double &u, double &p, double &Ψ)
 {
-   Mpi::Init(argc, argv);
+   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 vOrder = 2;
-   int pOrder = 1;
-   int tOrder = 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;
 
-   //Mesh *mesh = new Mesh("box-cylinder.mesh");
-   Mesh mesh = Mesh::MakeCartesian2D(90, 30, mfem::Element::QUADRILATERAL, true, 3, 1);
+   constexpr int precision = 8;
+   std::cout.precision(precision);
 
-   for (int i = 0; i < serial_refinements; ++i)
+   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())
    {
-      mesh.UniformRefinement();
+      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())
    {
@@ -65,8 +108,9 @@ int main(int argc, char *argv[])
    auto *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
 
    // Create the flow solver.
-   IncompressibleNavierSolver flowsolver(pmesh, vOrder, pOrder, tOrder, kin_vis);
-   flowsolver.EnablePA(false);
+   IncompressibleNavierSolver flowsolver(pmesh, v_order, p_order, t_order,
+                                         kin_vis);
+   flowsolver.EnablePA(pa);
 
    // // Set the initial condition.
    // ParGridFunction *u_ic = flowsolver.GetCurrentVelocity();
@@ -75,8 +119,10 @@ int main(int argc, char *argv[])
 
    // 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;
+   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.
@@ -93,42 +139,67 @@ int main(int argc, char *argv[])
    ParGridFunction *psi_gf = flowsolver.GetCurrentPsi();
 
    ParaViewDataCollection pvdc("3dfoc", pmesh);
-   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();
+   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 (t + dt >= t_final - dt / 2)
-      {
-         last_step = true;
-      }
+      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);
+      flowsolver.Step(t, dt, step, vis_step);
 
-      if (step % 1 == 0)
+      if (vis_step)
       {
-         pvdc.SetCycle(step);
-         pvdc.SetTime(t);
-         pvdc.Save();
-      }
-
-      if (Mpi::Root())
-      {
-         printf("%11s %11s\n", "Time", "dt");
-         printf("%.5E %.5E\n", t, dt);
-         fflush(stdout);
+         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 0;
+   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

@@ -10,24 +10,28 @@
 // CONTRIBUTING.md for details.
 
 #include "incompressible_navier_solver.hpp"
-#include "../../general/forall.hpp"
-#include <fstream>
-#include <iomanip>
+
+#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), velGF(torder_+1,nullptr), pGF(torder_+1,nullptr)
-{
-   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);
+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))
    {
@@ -38,13 +42,12 @@ IncompressibleNavierSolver::IncompressibleNavierSolver(ParMesh *mesh, int velord
       pres_ess_attr = 0;
    }
 
-   int vfes_truevsize = vfes->GetTrueVSize();
-   int pfes_truevsize = pfes->GetTrueVSize();
-
-   for( int i = 0; i<torder+1; i++)
+   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;
+      velGF[i] = new ParGridFunction(vfes);
+      *velGF[i] = 0.0;
+      pGF[i] = new ParGridFunction(pfes);
+      *pGF[i] = 0.0;
    }
 
    psiGF.SetSpace(psifes);
@@ -62,69 +65,63 @@ void IncompressibleNavierSolver::Setup(real_t dt)
       {
          mfem::out << "Using Partial Assembly" << std::endl;
       }
-      else
-      {
-         mfem::out << "Using Full Assembly" << std::endl;
-      }
+      else { mfem::out << "Using Full Assembly" << std::endl; }
    }
 
-   this->Setup_velocity( dt );
+   this->Setup_velocity(dt);
 
-   this->Setup_auxiliary( dt );
+   this->Setup_auxiliary(dt);
 
-   this->Setup_pressure( 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);
+   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 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 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);
+      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);
-   }
+   if (partial_assembly) { velBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
 
    velBForm->Assemble();
    velBForm->FormSystemMatrix(vel_ess_tdof, vOp);
 
    //-------------------------------------------------------------------------
-   
-   //Setup of coefficient for Eq(18)
+
+   // Setup of coefficient for Eq(18)
    pUnitVectorCoeff = new UnitVectorGridFunctionCoeff(pmesh->Dimension());
    auto *pvel_lfi = new VectorDomainLFGradIntegrator(*pUnitVectorCoeff);
 
-   //Setup of coefficient for Eq(20)
+   // Setup of coefficient for Eq(20)
    nonlinTermCoeff = new NonLinTermVectorGridFunctionCoeff(pmesh->Dimension());
    auto *p_nonlintermlfi = new VectorDomainLFIntegrator(*nonlinTermCoeff);
 
-   //Setup of coefficient for Eq(21)
+   // Setup of coefficient for Eq(21)
    prevVelLoadCoeff = new PrevVelVectorGridFunctionCoeff(pmesh->Dimension());
    auto *prevVelLoadLFi = new VectorDomainLFIntegrator(*prevVelLoadCoeff);
 
-   //Setup of linear form of Eq(13)
+   // Setup of linear form of Eq(13)
    velLForm = new ParLinearForm(vfes);
    if (numerical_integ)
    {
@@ -147,7 +144,8 @@ void IncompressibleNavierSolver::Setup_velocity(real_t dt)
    else
    {
       velInvPC = new HypreSmoother(*vOp.As<HypreParMatrix>());
-      dynamic_cast<HypreSmoother *>(velInvPC)->SetType(HypreSmoother::Jacobi, 1);
+      dynamic_cast<HypreSmoother *>(velInvPC)->SetType(HypreSmoother::Jacobi,
+                                                       1);
    }
 
    velInv = new CGSolver(vfes->GetComm());
@@ -156,54 +154,48 @@ void IncompressibleNavierSolver::Setup_velocity(real_t dt)
    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);
+   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); 
-   }
+   if (numerical_integ) { psidiff_blfi->SetIntRule(&ir_ni); }
    psiBForm->AddDomainIntegrator(psidiff_blfi);
-   if (partial_assembly)
-   {
-      psiBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   }
+   if (partial_assembly) { psiBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
 
    psiBForm->Assemble();
    psiBForm->FormSystemMatrix(empty, psiOp);
 
    //-------------------------------------------------------------------------
 
-   //Setup of coefficient for linear form in Eq(14)
+   // 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)
+
+   // Setup of linear form of Eq(14)
    psiLForm = new ParLinearForm(psifes);
-   
-   if (numerical_integ)
-   {
-      Dvel_lfi->SetIntRule(&ir_ni);
-   }
+
+   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;
+      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());
@@ -226,14 +218,15 @@ void IncompressibleNavierSolver::Setup_auxiliary(real_t dt)
    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);
+   const IntegrationRule &ir_ni =
+      gll_rules.Get(vfes->GetFE(0)->GetGeomType(), 2 * velorder - 1);
    Array<int> empty;
 
    //-------------------------------------------------------------------------
@@ -242,34 +235,25 @@ void IncompressibleNavierSolver::Setup_pressure(real_t dt)
    pBForm = new ParBilinearForm(pfes);
    auto *pmass_blfi = new MassIntegrator;
 
-   if (numerical_integ)
-   {
-       pmass_blfi->SetIntRule(&ir_ni); 
-   }
+   if (numerical_integ) { pmass_blfi->SetIntRule(&ir_ni); }
    pBForm->AddDomainIntegrator(pmass_blfi);
-   if (partial_assembly)
-   {
-      pBForm->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   }
+   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)
+   // 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)
+   // 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);
-   }
+   // Setup of linear form of Eq(15)
+   pLForm = new ParLinearForm(pfes);
+   if (numerical_integ) { p_lfi->SetIntRule(&ir_ni); }
    pLForm->AddDomainIntegrator(p_lfi);
 
    //-------------------------------------------------------------------------
@@ -291,15 +275,14 @@ void IncompressibleNavierSolver::Setup_pressure(real_t dt)
    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::UpdateTimestepHistory(real_t dt) {}
 
-}
-
-void IncompressibleNavierSolver::Step(real_t &time, real_t dt, int current_step)
+void IncompressibleNavierSolver::Step(real_t &time, real_t dt, int current_step,
+                                      const bool vis_step)
 {
    this->Step_velocity(time, dt, current_step);
 
@@ -308,16 +291,24 @@ void IncompressibleNavierSolver::Step(real_t &time, real_t dt, int current_step)
    this->Step_pressure(time, dt, current_step);
 
    *velGF[1] = *velGF[0];
-   *pGF[1]   = *pGF[0];
-   
-   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";
+   *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 ++;
+   iter++;
 }
 
-void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt, int current_step)
+void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt,
+                                               int current_step)
 {
    for (auto &vel_dbc : vel_dbcs)
    {
@@ -325,14 +316,14 @@ void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt, int curr
       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(18)
+   pUnitVectorCoeff->SetGridFunction(pGF[1]);
 
-   //Update state in coefficient for Eq(21)
-   prevVelLoadCoeff ->SetGridFunction( velGF[1], dt ); 
+   // 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);
@@ -346,7 +337,8 @@ void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt, int curr
    }
    else
    {
-      velBForm->FormLinearSystem(vel_ess_tdof, *velGF[0], velLF, vOp  , X1, B1, 1);
+      velBForm->FormLinearSystem(vel_ess_tdof, *velGF[0], velLF, vOp, X1, B1,
+                                 1);
    }
 
    velInv->Mult(B1, X1);
@@ -355,26 +347,24 @@ void IncompressibleNavierSolver::Step_velocity(real_t &time, real_t dt, int curr
    velBForm->RecoverFEMSolution(X1, velLF, *velGF[0]);
 }
 
-void IncompressibleNavierSolver::Step_auxiliary(real_t &time, real_t dt, int current_step)
+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 );
+   subtract(1.0 / dt, *velGF[0], *velGF[1], DvGF);
+   DvelCoeff->SetGridFunction(&DvGF);
 
    psiLForm->Assemble();
    psiLForm->ParallelAssemble(psiLF);
 
-   Vector X2, B2;   
+   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);
-   }
+   else { psiBForm->FormLinearSystem(empty, psiGF, psiLF, psiOp, X2, B2, 1); }
 
    psiInv->Mult(B2, X2);
    iter_psisolve = psiInv->GetNumIterations();
@@ -382,17 +372,18 @@ void IncompressibleNavierSolver::Step_auxiliary(real_t &time, real_t dt, int cur
    psiBForm->RecoverFEMSolution(X2, psiLF, psiGF);
 }
 
-void IncompressibleNavierSolver::Step_pressure(real_t &time, real_t dt, int current_step)
+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 );
+   divVelCoeff->SetGridFunction(velGF[0]);
+   divVelGF.ProjectCoefficient(*divVelCoeff);
 
-   add( *pGF[1], psiGF, pRHS);
-   add( pRHS, -1.0*kin_vis, divVelGF, pRHS);
-   pRHSCoeff->SetGridFunction( &pRHS );
+   add(*pGF[1], psiGF, pRHS);
+   add(pRHS, -1.0 * kin_vis, divVelGF, pRHS);
+   pRHSCoeff->SetGridFunction(&pRHS);
 
    pLForm->Assemble();
    pLForm->ParallelAssemble(pLF);
@@ -404,10 +395,7 @@ void IncompressibleNavierSolver::Step_pressure(real_t &time, real_t dt, int curr
       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);
-   }
+   else { pBForm->FormLinearSystem(empty, *pGF[0], pLF, pOp, X3, B3, 1); }
 
    pInv->Mult(B3, X3);
    iter_psolve = pInv->GetNumIterations();
@@ -415,35 +403,27 @@ void IncompressibleNavierSolver::Step_pressure(real_t &time, real_t dt, int curr
    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)
+                                              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);
-   }
+   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)
+void IncompressibleNavierSolver::AddVelDirichletBC(VectorCoefficient *coeff,
+                                                   Array<int> &attr)
 {
    vel_dbcs.emplace_back(attr, coeff);
 
@@ -452,10 +432,7 @@ void IncompressibleNavierSolver::AddVelDirichletBC(VectorCoefficient *coeff, Arr
       mfem::out << "Adding Velocity Dirichlet BC to attributes ";
       for (int i = 0; i < attr.Size(); ++i)
       {
-         if (attr[i] == 1)
-         {
-            mfem::out << i << " ";
-         }
+         if (attr[i] == 1) { mfem::out << i << " "; }
       }
       mfem::out << std::endl;
    }
@@ -464,17 +441,15 @@ void IncompressibleNavierSolver::AddVelDirichletBC(VectorCoefficient *coeff, Arr
    {
       MFEM_ASSERT((vel_ess_attr[i] && attr[i]) == 0,
                   "Duplicate boundary definition deteceted.");
-      if (attr[i] == 1)
-      {
-         vel_ess_attr[i] = 1;
-      }
+      if (attr[i] == 1) { vel_ess_attr[i] = 1; }
    }
-
 }
 
-void IncompressibleNavierSolver::AddVelDirichletBC(VecFuncT *f, Array<int> &attr)
+void IncompressibleNavierSolver::AddVelDirichletBC(VecFuncT *f,
+                                                   Array<int> &attr)
 {
-   AddVelDirichletBC(new VectorFunctionCoefficient(pmesh->Dimension(), f), attr);
+   AddVelDirichletBC(new VectorFunctionCoefficient(pmesh->Dimension(), f),
+                     attr);
 }
 
 IncompressibleNavierSolver::~IncompressibleNavierSolver()
@@ -486,7 +461,7 @@ IncompressibleNavierSolver::~IncompressibleNavierSolver()
    delete kinvisCoeff;
    delete dtCoeff;
 
-   for( int i = 0; i<torder+1; i++)
+   for (int i = 0; i < torder + 1; i++)
    {
       delete velGF[i];
       delete pGF[i];

miniapps/navier/incompressible_navier_solver.hpp

@@ -9,67 +9,60 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#ifndef MFEM_INCOMP_NAVIER_SOLVER_HPP
-#define MFEM_INCOMP_NAVIER_SOLVER_HPP
+#pragma once
 
 #define INCOMP_NAVIER_VERSION 0.1
 
 #include "mfem.hpp"
 
-namespace mfem
-{
-namespace incompressible_navier
+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)
+// Coefficient which computed contribution of Eq(18)
 class UnitVectorGridFunctionCoeff : public VectorCoefficient
 {
 public:
-   UnitVectorGridFunctionCoeff( int dim)
-   : VectorCoefficient(dim*dim)
-   { }
+   UnitVectorGridFunctionCoeff(int dim): VectorCoefficient(dim * dim) {}
 
-   void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
+   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.SetSize(vdim);
+      V = 0.0; // FIXME
       V[0] = coeffVal;
       V[3] = coeffVal;
    }
 
-   void SetGridFunction( GridFunction * gridfunc )
-   {
-      gridfunc_ = gridfunc;
-   }
+   void SetGridFunction(GridFunction *gridfunc) { gridfunc_ = gridfunc; }
 
    GridFunction *gridfunc_ = nullptr;
 };
 
-//Coefficient which computed contribution of Eq(21)
+// Coefficient which computed contribution of Eq(21)
 class PrevVelVectorGridFunctionCoeff : public VectorCoefficient
 {
 public:
-   PrevVelVectorGridFunctionCoeff( int dim)
-   : VectorCoefficient(dim)
-   { }
+   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_;
+      V *= 1.0 / dt_;
    }
 
-   void SetGridFunction( GridFunction * gridfunc, real_t dt )
+   void SetGridFunction(GridFunction *gridfunc, real_t dt)
    {
       gridfunc_ = gridfunc;
       dt_ = dt;
       delete gridFuncCoeff;
-      gridFuncCoeff = new VectorGridFunctionCoefficient( gridfunc );
+      gridFuncCoeff = new VectorGridFunctionCoefficient(gridfunc);
    }
 
    GridFunction *gridfunc_ = nullptr;
@@ -77,13 +70,11 @@ public:
    real_t dt_;
 };
 
-//Coefficient which computed contribution of Eq(20)
+// Coefficient which computed contribution of Eq(20)
 class NonLinTermVectorGridFunctionCoeff : public VectorCoefficient
 {
 public:
-   NonLinTermVectorGridFunctionCoeff( int dim)
-   : VectorCoefficient(dim)
-   { }
+   NonLinTermVectorGridFunctionCoeff(int dim): VectorCoefficient(dim) {}
 
    void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)
    {
@@ -95,16 +86,16 @@ public:
 
       gridfunc_->GetVectorGradient(T, vecGrad);
 
-      vecGrad.MultTranspose( val, V );
+      vecGrad.MultTranspose(val, V);
 
       V *= -1.0;
    }
 
-   void SetGridFunction( ParGridFunction * gridfunc )
+   void SetGridFunction(ParGridFunction *gridfunc)
    {
       delete gridFuncCoeff;
       gridfunc_ = gridfunc;
-      gridFuncCoeff = new VectorGridFunctionCoefficient( gridfunc );
+      gridFuncCoeff = new VectorGridFunctionCoefficient(gridfunc);
    }
 
    VectorGridFunctionCoefficient *gridFuncCoeff = nullptr;
@@ -115,9 +106,10 @@ public:
 class VelDirichletBC_T
 {
 public:
-   VelDirichletBC_T(Array<int> attr, VectorCoefficient *coeff)
-      : attr(attr), coeff(coeff)
-   {}
+   VelDirichletBC_T(Array<int> attr, VectorCoefficient *coeff):
+      attr(attr), coeff(coeff)
+   {
+   }
 
    VelDirichletBC_T(VelDirichletBC_T &&obj)
    {
@@ -151,7 +143,8 @@ public:
     * 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);
+   IncompressibleNavierSolver(ParMesh *mesh, int velorder, int porder,
+                              int tOrder, real_t kin_vis);
 
    /// Initialize forms, solvers and preconditioners.
    void Setup(real_t dt);
@@ -164,9 +157,9 @@ public:
 
    /// Compute solution at the next time step t+dt.
    /**
-    * This method can 
+    * This method can
     */
-   void Step(real_t &time, real_t dt, int cur_step);
+   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);
 
@@ -184,8 +177,7 @@ public:
    ParGridFunction *GetCurrentPressure() { return pGF[0]; }
 
    /// Return a pointer to the current pressure ParGridFunction.
-   ParGridFunction *GetCurrentPsi() { return &psiGF ; }
-
+   ParGridFunction *GetCurrentPsi() { return &psiGF; }
 
    /// Add a Dirichlet boundary condition to the velocity field.
    void AddVelDirichletBC(VectorCoefficient *coeff, Array<int> &attr);
@@ -205,10 +197,6 @@ public:
    void EnableNI(bool ni) { numerical_integ = ni; }
 
    /// Print timing summary of the solving routine.
-   /**
-    * The summary shows the timing in seconds in the first row of
-    *
-    */
    void PrintTimingData();
 
    ~IncompressibleNavierSolver();
@@ -216,21 +204,14 @@ public:
    /// 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);
+   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;
@@ -248,13 +229,13 @@ protected:
    ParMesh *pmesh = nullptr;
 
    /// The order of the velocity and pressure space.
-   int velorder;
-   int porder;
-   int torder;
+   const int velorder;
+   const int porder;
+   const int torder;
 
    /// Kinematic viscosity (dimensionless).
-   real_t kin_vis;
-   Coefficient * kinvisCoeff = nullptr;
+   const real_t kin_vis;
+   Coefficient *kinvisCoeff = nullptr;
 
    Coefficient *dtCoeff = nullptr;
 
@@ -278,26 +259,29 @@ protected:
    /// Pressure $H^1$ finite element space.
    ParFiniteElementSpace *pfes = nullptr;
 
-   ParBilinearForm *velBForm = nullptr;
-   ParBilinearForm *psiBForm = nullptr;
-   ParBilinearForm *pBForm = nullptr;
+   ParBilinearForm *velBForm = nullptr; // vmass + vdiff
+   ParBilinearForm *psiBForm = nullptr; // diffusion
+   ParBilinearForm *pBForm = nullptr;   // mass
 
-   ParLinearForm *velLForm = nullptr;
-   ParLinearForm *psiLForm = nullptr;
-   ParLinearForm *pLForm = nullptr;
+   ParLinearForm *velLForm = nullptr; // vLF + vLFGrad + vLF
+   ParLinearForm *psiLForm = nullptr; // LFGrad
+   ParLinearForm *pLForm = nullptr;   // LF
 
-   std::vector<ParGridFunction*> velGF;
-   std::vector<ParGridFunction*> pGF;
+   // 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;
-   
+   VectorGridFunctionCoefficient *DvelCoeff = nullptr;
+   DivergenceGridFunctionCoefficient *divVelCoeff = nullptr;
+   GridFunctionCoefficient *pRHSCoeff = nullptr;
+   UnitVectorGridFunctionCoeff *pUnitVectorCoeff = nullptr;
+   NonLinTermVectorGridFunctionCoeff *nonlinTermCoeff = nullptr;
+   PrevVelVectorGridFunctionCoeff *prevVelLoadCoeff = nullptr;
+
    OperatorHandle vOp;
    OperatorHandle psiOp;
    OperatorHandle pOp;
@@ -332,19 +316,16 @@ protected:
    int pl_velsolve = 0;
    int pl_amg = 0;
 
-   #if defined(MFEM_USE_DOUBLE)
-   real_t rtol_psolve = 1e-10;
-   real_t rtol_psisolve = 1e-10;
+#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!"
-   real_t rtol_psolve = 1e-12;
-   real_t rtol_psisolve = 1e-6;
-   real_t rtol_velsolve = 1e-8;
+   #error "Only single and double precision are supported!"
 #endif
 
    // Iteration counts.
@@ -352,11 +333,6 @@ protected:
 
    // Residuals.
    real_t res_vsolve = 0.0, res_psolve = 0.0, res_psisolve = 0.0;
-
 };
 
-} // namespace incompressible_navier
-
-} // namespace mfem
-
-#endif
+} // 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;
+}

miniapps/navier/stokes_solver.cpp

@@ -1,119 +0,0 @@
-#include "stokes_solver.hpp"
-
-namespace mfem {
-
-StokesOperator::StokesOperator(ParFiniteElementSpace &vel_fes,
-                               ParFiniteElementSpace &pres_fes):
-    Operator(vel_fes.GetTrueVSize()+pres_fes.GetTrueVSize()),
-    vfes(vel_fes),
-    pfes(pres_fes),
-    offsets({0, vel_fes.GetTrueVSize(), pres_fes.GetTrueVSize()}),
-    intrules(0, Quadrature1D::GaussLobatto),
-    zero_coeff(0.0)
-{
-    if (vel_fes.GetParMesh()->bdr_attributes.Size() > 0)
-    {
-        vel_ess_bdr.SetSize(vel_fes.GetParMesh()->bdr_attributes.Max());
-        vel_ess_bdr = 0.0;
-        pres_ess_bdr.SetSize(vel_fes.GetParMesh()->bdr_attributes.Max());
-        pres_ess_bdr = 0.0;
-    }
-
-    vfes.GetEssentialTrueDofs(vel_ess_bdr, vel_ess_tdofs);
-    pfes.GetEssentialTrueDofs(pres_ess_bdr, pres_ess_tdofs);
-
-    offsets.PartialSum();
-
-
-    vel_bc_gf.reset(new ParGridFunction(&vfes));
-    *vel_bc_gf = 0.0; //set the velocity grid function to zero
-
-    pres_bc_gf.reset(new ParGridFunction(&pfes));
-    *pres_bc_gf = 0.0; //set the pressure grid function to zero
-
-    // The nonlinear convective integrators use over-integration (dealiasing) as
-    // a stabilization mechanism.
-    ir_nl = intrules.Get(vfes.GetFE(0)->GetGeomType(),
-                         (int)(ceil(1.5 * 2*(vel_fes.GetOrder(0)+1) - 3)));
-
-    ir = intrules.Get(vfes.GetFE(0)->GetGeomType(),
-                      (int)(2*(vel_fes.GetOrder(0)+1) - 3));
-
-    ir_face = intrules.Get(vfes.GetFaceElement(0)->GetGeomType(),
-                           (int)(2*(vel_fes.GetOrder(0)+1) - 3));
-
-    b11_form=nullptr;
-    b22_form=nullptr;
-    b12_form=nullptr;
-    b21_form=nullptr;
-
-}
-
-
-void StokesOperator::SetVelBC(std::vector<VelDirichletBC>& vvbc)
-{
-    for(auto vbc=vvbc.begin();vbc!=vvbc.end();vbc++)
-    {
-        for (int i = 0; i < vbc->second->Size(); i++)
-        {
-            if (*(vbc->second)[i] == 1)
-            {
-                vel_ess_bdr[i] = 1;
-            }
-        }
-    }
-
-    vfes.GetEssentialTrueDofs(vel_ess_bdr, vel_ess_tdofs);
-}
-
-void StokesOperator::SetPressBC(std::vector<PresDirichletBC>& vpbc)
-{
-    for(auto pbc=vpbc.begin();pbc!=vpbc.end();pbc++)
-    {
-        for(int i=0;i<pbc->second->Size();i++){
-            if (*(pbc->second)[i] == 1)
-            {
-                vel_ess_bdr[i] = 1;
-            }
-        }
-    }
-
-    pfes.GetEssentialTrueDofs(pres_ess_bdr, pres_ess_tdofs);
-}
-
-void StokesOperator::Mult(const Vector &x, Vector &y) const
-{
-
-}
-
-void StokesOperator::Setup()
-{
-    BilinearFormIntegrator *integrator;
-
-    delete b11_form;
-    b11_form=new ParBilinearForm(&vfes);
-    integrator=new ElasticityIntegrator(zero_coeff,*viscosity);
-    integrator->SetIntRule(&ir);
-    b11_form->AddDomainIntegrator(integrator);
-
-    delete b12_form;
-    b12_form=new ParMixedBilinearForm(&pfes,&vfes);
-    integrator=new VectorDivergenceIntegrator();
-    integrator->SetIntRule(&ir);
-    b12_form->AddDomainIntegrator(integrator);
-
-    delete b21_form;
-    b21_form=new ParMixedBilinearForm(&vfes,&pfes);
-    integrator=new GradientIntegrator();
-    integrator->SetIntRule(&ir);
-    b21_form->AddDomainIntegrator(integrator);
-
-    if (matrix_free)
-    {
-        b11_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-        b12_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-        b21_form->SetAssemblyLevel(AssemblyLevel::PARTIAL);
-    }
-}
-
-}

miniapps/navier/stokes_solver.hpp

@@ -1,76 +0,0 @@
-#ifndef STOKESSOLVER_H
-#define STOKESSOLVER_H
-
-#define STOKES_VERSION 0.1
-
-#include "mfem.hpp"
-
-
-namespace mfem {
-
-using VelDirichletBC = std::pair<VectorCoefficient *, Array<int> *>;
-using PresDirichletBC = std::pair<Coefficient *, Array<int> *>;
-
-class StokesOperator:public Operator
-{
-public:
-    StokesOperator(ParFiniteElementSpace &vel_fes,
-                   ParFiniteElementSpace &pres_fes);
-
-    void SetVelBC(std::vector<VelDirichletBC>& vvbc);
-    void SetPressBC(std::vector<PresDirichletBC>& vpbc);
-
-    virtual
-    void Mult(const Vector &x, Vector &y) const override;
-
-    const Array<int>& GetOffsets() const
-    {
-        return offsets;
-    }
-
-    void Setup();
-    void Assemble();
-
-
-private:
-    ParFiniteElementSpace &vfes;
-    ParFiniteElementSpace &pfes;
-
-// ParGridFunction &kinematic_viscosity;
-
-    std::unique_ptr<ParGridFunction> vel_bc_gf;
-    std::unique_ptr<ParGridFunction> pres_bc_gf;
-
-    Array<int> vel_ess_bdr;
-    Array<int> pres_ess_bdr;
-
-    Array<int> vel_ess_tdofs;
-    Array<int> pres_ess_tdofs;
-
-    bool matrix_free;
-
-    Array<int> offsets;
-
-
-    IntegrationRules intrules;
-    IntegrationRule ir; //general integraion rule
-    IntegrationRule ir_nl; //non-linear integration rule
-    IntegrationRule ir_face; //face integration rule
-
-    ConstantCoefficient zero_coeff;
-
-    std::unique_ptr<Coefficient> viscosity;
-
-    ParBilinearForm *b11_form; //velocity
-    ParBilinearForm *b22_form; //pressure
-    ParMixedBilinearForm *b12_form; //mixed (velocity,pressure)
-    ParMixedBilinearForm *b21_form; //mized (pressure,velocity)
-
-    BlockOperator* A;
-
-};
-
-
-}
-
-#endif // STOKESSOLVER_H