Cleanup

Serial DFEM diffusion tests

examples/dfem/laghos.cpp

@@ -486,8 +486,8 @@ public:
          Vector wvc, yvc;
          for (int c = 0; c < hydro.H1.GetMesh()->Dimension(); c++)
          {
-            wvc.MakeRef(wv, c*hydro.H1c.GetVSize(), hydro.H1c.GetVSize());
-            yvc.MakeRef(yv, c*hydro.H1c.GetVSize(), hydro.H1c.GetVSize());
+            wvc.MakeRef(wv, c*hydro.H1c.GetTrueVSize(), hydro.H1c.GetTrueVSize());
+            yvc.MakeRef(yv, c*hydro.H1c.GetTrueVSize(), hydro.H1c.GetTrueVSize());
             hydro.Mv->FullAddMult(wvc, yvc);
             yvc.SyncAliasMemory(yv);
          }
@@ -556,6 +556,9 @@ public:
       u(x.Size()),
       H1tsize(hydro.H1.GetTrueVSize()),
       L2tsize(hydro.L2.GetTrueVSize()),
+      ux_l(hydro.H1.GetVSize()),
+      uv_l(hydro.H1.GetVSize()),
+      ue_l(hydro.L2.GetVSize()),
       fd_gradient(fd_gradient) {}
 
    void Mult(const Vector &k, Vector &R) const override
@@ -584,18 +587,22 @@ public:
       Rv.MakeRef(R, H1tsize, H1tsize);
       Re.MakeRef(R, 2*H1tsize, L2tsize);
 
+      hydro.H1.GetProlongationMatrix()->Mult(ux, ux_l);
+      hydro.H1.GetProlongationMatrix()->Mult(uv, uv_l);
+      hydro.L2.GetProlongationMatrix()->Mult(ue, ue_l);
+
       Rx = kx;
       Rx -= uv;
 
-      hydro.momentum_mf->SetParameters({&hydro.rho0, &hydro.x0, &ux, &hydro.material, &ue, &hydro.qdata->h0, &hydro.qdata->order_v});
+      hydro.momentum_mf->SetParameters({&hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &ue_l, &hydro.qdata->h0, &hydro.qdata->order_v});
       hydro.momentum_mf->Mult(uv, Rv);
 
       // hydro.Mv.TrueAddMult(kv, Rv);
       Vector kvc, Rvc;
       for (int c = 0; c < hydro.H1.GetMesh()->Dimension(); c++)
       {
-         kvc.MakeRef(kv, c*hydro.H1c.GetVSize(), hydro.H1c.GetVSize());
-         Rvc.MakeRef(Rv, c*hydro.H1c.GetVSize(), hydro.H1c.GetVSize());
+         kvc.MakeRef(kv, c*hydro.H1c.GetTrueVSize(), hydro.H1c.GetTrueVSize());
+         Rvc.MakeRef(Rv, c*hydro.H1c.GetTrueVSize(), hydro.H1c.GetTrueVSize());
          hydro.Mv->FullAddMult(kvc, Rvc);
          Rvc.SyncAliasMemory(Rv);
       }
@@ -604,7 +611,7 @@ public:
       Rv.SetSubVector(hydro.ess_tdof, 0.0);
       // Rv = 0.0;
 
-      hydro.energy_conservation_mf->SetParameters({&uv, &hydro.rho0, &hydro.x0, &ux, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
+      hydro.energy_conservation_mf->SetParameters({&uv_l, &hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
       hydro.energy_conservation_mf->Mult(ue, Re);
 
       Re.Neg();
@@ -649,6 +656,10 @@ public:
       uv.MakeRef(u, H1tsize, H1tsize);
       ue.MakeRef(u, 2*H1tsize, L2tsize);
 
+      hydro.H1.GetProlongationMatrix()->Mult(ux, ux_l);
+      hydro.H1.GetProlongationMatrix()->Mult(uv, uv_l);
+      hydro.L2.GetProlongationMatrix()->Mult(ue, ue_l);
+
       if (fd_gradient)
       {
          fd_jacobian.reset(new FDJacobian(*this, k));
@@ -656,24 +667,24 @@ public:
       }
       else
       {
-         auto dRvdx = hydro.momentum_mf->GetDerivative(COORDINATES, {&uv},
-         {&hydro.rho0, &hydro.x0, &ux, &hydro.material, &ue, &hydro.qdata->h0, &hydro.qdata->order_v});
+         auto dRvdx = hydro.momentum_mf->GetDerivative(COORDINATES, {&uv_l},
+         {&hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &ue_l, &hydro.qdata->h0, &hydro.qdata->order_v});
 
-         auto dRvdv = hydro.momentum_mf->GetDerivative(VELOCITY, {&uv},
-         {&hydro.rho0, &hydro.x0, &ux, &hydro.material, &ue, &hydro.qdata->h0, &hydro.qdata->order_v});
+         auto dRvdv = hydro.momentum_mf->GetDerivative(VELOCITY, {&uv_l},
+         {&hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &ue_l, &hydro.qdata->h0, &hydro.qdata->order_v});
 
-         auto dRvde = hydro.momentum_mf->GetDerivative(SPECIFIC_INTERNAL_ENERGY, {&uv},
-         {&hydro.rho0, &hydro.x0, &ux, &hydro.material, &ue, &hydro.qdata->h0, &hydro.qdata->order_v});
+         auto dRvde = hydro.momentum_mf->GetDerivative(SPECIFIC_INTERNAL_ENERGY, {&uv_l},
+         {&hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &ue_l, &hydro.qdata->h0, &hydro.qdata->order_v});
 
-         auto dRedx = hydro.energy_conservation_mf->GetDerivative(COORDINATES, {&ue},
-         {&uv, &hydro.rho0, &hydro.x0, &ux, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
+         auto dRedx = hydro.energy_conservation_mf->GetDerivative(COORDINATES, {&ue_l},
+         {&uv_l, &hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
 
-         auto dRedv = hydro.energy_conservation_mf->GetDerivative(VELOCITY, {&ue},
-         {&uv, &hydro.rho0, &hydro.x0, &ux, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
+         auto dRedv = hydro.energy_conservation_mf->GetDerivative(VELOCITY, {&ue_l},
+         {&uv_l, &hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
 
          auto dRede = hydro.energy_conservation_mf->GetDerivative(
-         SPECIFIC_INTERNAL_ENERGY, {&ue},
-         {&uv, &hydro.rho0, &hydro.x0, &ux, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
+         SPECIFIC_INTERNAL_ENERGY, {&ue_l},
+         {&uv_l, &hydro.rho0, &hydro.x0, &ux_l, &hydro.material, &hydro.qdata->h0, &hydro.qdata->order_v});
 
          jacobian->Setup(hydro, dRvdx, dRvdv, dRvde, dRedx, dRedv, dRede);
          return *jacobian;
@@ -683,7 +694,7 @@ public:
    hydro_t &hydro;
    const real_t dt;
    const Vector &x;
-   mutable Vector u;
+   mutable Vector u, ux_l, uv_l, ue_l;
    const int H1tsize;
    const int L2tsize;
    mutable std::shared_ptr<FDJacobian> fd_jacobian;

fem/dfem/doperator.hpp

@@ -18,6 +18,10 @@
 #include "qfunction.hpp"
 #include "integrate.hpp"
 
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kOrchid
+#include "general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -75,7 +79,7 @@ public:
 
    void Mult(const Vector &direction_t, Vector &y) const override
    {
-      daction_l.SetSize(height);
+      // daction_l.SetSize(height);
       daction_l = 0.0;
 
       prolongation(direction, direction_t, direction_l);
@@ -88,7 +92,7 @@ public:
 
    void MultTranspose(const Vector &direction_t, Vector &y) const override
    {
-      daction_l.SetSize(width);
+      // daction_l.SetSize(width);
       daction_l = 0.0;
 
       prolongation(transpose_direction, direction_t, direction_l);
@@ -134,6 +138,18 @@ public:
       prolongation_transpose(residual_l, y);
    }
 
+   void Mult(ParGridFunction &x, ParGridFunction &y) const
+   {
+      MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
+      MFEM_VERIFY(y.Size() == residual_l.Size(), "output size mismatch");
+      prolongation(solutions, x.GetTrueVector(), solutions_l);
+      for (auto &action : action_callbacks)
+      {
+         action(solutions_l, parameters_l, residual_l);
+      }
+      y = residual_l;
+   }
+
    template <
       typename func_t,
       typename... input_ts,
@@ -179,7 +195,7 @@ public:
                 residual_l.Size(),
                 daction_transpose_callbacks[derivative_id],
                 fields[test_space_field_idx],
-                GetVSize(fields[test_space_field_idx]),
+                GetTrueVSize(fields[test_space_field_idx]),
                 solutions_l,
                 parameters_l,
                 restriction_callback,
@@ -325,11 +341,21 @@ void DifferentiableOperator::AddDomainIntegrator(
       inputs_vdim[i] = mfem::get<i>(inputs).vdim;
    });
 
-   Array<int> elem_attributes;
-   elem_attributes.SetSize(mesh.GetNE());
-   for (int i = 0; i < mesh.GetNE(); ++i)
+   const int NE = mesh.GetNE();
+   if (NE == 0)
    {
-      elem_attributes[i] = mesh.GetAttribute(i);
+      // use of GetElement(0), GetFE(0) in GetDofToQuad assume that NE > 0
+      MFEM_ABORT("Mesh with no elements is not yet supported!");
+   }
+
+   Array<int> elem_attributes;
+   if (NE > 0)
+   {
+      elem_attributes.SetSize(NE);
+      for (int i = 0; i < NE; ++i)
+      {
+         elem_attributes[i] = mesh.GetAttribute(i);
+      }
    }
 
    const auto output_fop = mfem::get<0>(outputs);
@@ -418,6 +444,8 @@ void DifferentiableOperator::AddDomainIntegrator(
 
    const int test_vdim = output_fop.vdim;
    const int test_op_dim = output_fop.size_on_qp / output_fop.vdim;
+   MFEM_VERIFY(num_entities > 0,
+               "The number of entities must be greater than zero");
    const int num_test_dof = output_e_size / output_fop.vdim /
                             num_entities;
 
@@ -470,12 +498,13 @@ void DifferentiableOperator::AddDomainIntegrator(
                                           action_shmem_info.field_sizes,
                                           num_entities);
 
+      const bool has_attr = domain_attributes.Size() > 0;
       const auto d_domain_attr = domain_attributes.Read();
       const auto d_elem_attr = elem_attributes.Read();
 
       forall([=] MFEM_HOST_DEVICE (int e, void *shmem)
       {
-         if (!d_domain_attr[d_elem_attr[e] - 1]) { return; }
+         if (has_attr && !d_domain_attr[d_elem_attr[e] - 1]) { return; }
 
          auto [input_dtq_shmem, output_dtq_shmem, fields_shmem, input_shmem,
                                 residual_shmem, scratch_shmem] =
@@ -687,7 +716,8 @@ void DifferentiableOperator::AddDomainIntegrator(
                            apply_kernel_fwddiff_enzyme(r, qfunc, qf_args, qf_shadow_args, input_shmem,
                                                        shadow_shmem, q);
 #else
-                           apply_kernel_native_dual(r, qfunc, qf_args, input_shmem, shadow_shmem, q);
+                           MFEM_ABORT("Native dual support is not enabled!");
+                           // apply_kernel_native_dual(r, qfunc, qf_args, input_shmem, shadow_shmem, q);
 #endif
                            d_qp(j, m, q) = 0.0;
 

fem/dfem/parametricspace.hpp

@@ -46,6 +46,18 @@ public:
       dtq.nqpt = dtq.ndof;
    }
 
+   ParametricSpace(int spatial_dim, int local_size, int element_size,
+                   int total_size, int d1d, int q1d) :
+      spatial_dim(spatial_dim),
+      local_size(local_size),
+      element_size(element_size),
+      total_size(total_size),
+      identity(total_size)
+   {
+      dtq.ndof = d1d;
+      dtq.nqpt = q1d;
+   }
+
    int Dimension() const
    {
       return spatial_dim;

fem/dfem/util.hpp

@@ -30,9 +30,20 @@
 
 using std::size_t;
 
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kGold
+#include "general/nvtx.hpp"
+
 namespace mfem
 {
 
+auto print_vec = [](const char *header, const Vector &v)
+{
+   dbl("{}:", header);
+   for (int i=0; i < v.Size(); i++) { dba("{:f} ", v(i)); }
+   dbc();
+};
+
 template<typename... Ts>
 constexpr auto to_array(const std::tuple<Ts...>& tuple)
 {
@@ -776,10 +787,12 @@ const Operator *get_prolongation(const FieldDescriptor &f)
       if constexpr (std::is_same_v<T, const FiniteElementSpace *> ||
                     std::is_same_v<T, const ParFiniteElementSpace *>)
       {
+         dbg("FES GetProlongationMatrix");
          return arg->GetProlongationMatrix();
       }
       else if constexpr (std::is_same_v<T, const ParametricSpace *>)
       {
+         dbg("ParametricSpace GetProlongation");
          return arg->GetProlongation();
       }
       else
@@ -917,6 +930,8 @@ void prolongation(const std::vector<FieldDescriptor> fields,
       const Vector x_i(const_cast<Vector&>(x), data_offset, width);
       fields_l[i].SetSize(P->Height());
       P->Mult(x_i, fields_l[i]);
+      print_vec("[P  ] x", x_i);
+      print_vec("[P  ]lx", fields_l[i]);
       data_offset += width;
    }
 }
@@ -930,6 +945,7 @@ auto get_prolongation_transpose(const FieldDescriptor &f, const fop_t &fop,
    {
       auto PT = [=](const Vector &r_local, Vector &y)
       {
+         dbg("[PT ] 1");
          double local_sum = r_local.Sum();
          MPI_Allreduce(&local_sum, y.GetData(), 1, MPI_DOUBLE, MPI_SUM, mpi_comm);
          MFEM_ASSERT(y.Size() == 1, "output size doesn't match kernel description");
@@ -940,6 +956,7 @@ auto get_prolongation_transpose(const FieldDescriptor &f, const fop_t &fop,
    {
       auto PT = [](Vector &r_local, Vector &y)
       {
+         dbg("[PT ] 0");
          y = r_local;
       };
       return PT;
@@ -949,6 +966,7 @@ auto get_prolongation_transpose(const FieldDescriptor &f, const fop_t &fop,
       const Operator *P = get_prolongation(f);
       auto PT = [=](const Vector &r_local, Vector &y)
       {
+         dbg("[PT ] P");
          P->MultTranspose(r_local, y);
       };
       return PT;

fem/pfespace.cpp

@@ -23,6 +23,10 @@
 #include <limits>
 #include <list>
 
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kLavender
+#include "general/nvtx.hpp"
+
 namespace mfem
 {
 
@@ -1215,6 +1219,7 @@ const Operator *ParFiniteElementSpace::GetProlongationMatrix() const
 
       if (nd_strias) { return Dof_TrueDof_Matrix(); }
 
+      dbg();
       if (NRanks == 1)
       {
          Pconf = new IdentityOperator(GetTrueVSize());
@@ -1234,6 +1239,7 @@ const Operator *ParFiniteElementSpace::GetProlongationMatrix() const
    }
    else
    {
+      assert(false);
       return Dof_TrueDof_Matrix();
    }
 }
@@ -3646,6 +3652,7 @@ ConformingProlongationOperator::ConformingProlongationOperator(
 
 void ConformingProlongationOperator::Mult(const Vector &x, Vector &y) const
 {
+   dbg();
    MFEM_ASSERT(x.Size() == Width(), "");
    MFEM_ASSERT(y.Size() == Height(), "");
 

general/nvtx.hpp

@@ -0,0 +1,471 @@
+#pragma once
+
+#include <fmt/format.h>
+
+#include <array>
+#include <cassert>
+#include <cstddef>
+#include <cstdint>
+#include <iomanip>
+#include <iostream>
+#include <memory>
+#include <mutex>
+#include <stack>
+#include <string>
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+#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 nvtx
+{
+
+///////////////////////////////////////////////////////////////////////////////
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+// 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
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+inline size_t static_strlen(const char *str)
+{
+   return *str == '\0' ? 0 : static_strlen(str + 1) + 1;
+}
+
+inline 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);
+}
+
+inline 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 ::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)...); }
+      if (debug) { std::cout << fmt::format(fmt::runtime(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;
+      static bool ini = false;
+      if (!ini)
+      {
+         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");
+         ini = true;
+      }
+      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(...)                                                  \
+   ::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR).operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_INI(...)                                \
+   ::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, true) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_APPEND(...)                                 \
+   ::nvtx::Debug::Set(NVTX_FLF, NVTX_COLOR, false, false) \
+      .operator()(__VA_ARGS__)
+
+#define NVTX_DEBUG_NO_END(...)                                \
+   ::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;
+      static bool ini = false;
+      if (!ini)
+      {
+         eks = ::getenv("MFEM_EKS") != nullptr;
+         nvtx = ::getenv("MFEM_NVTX") != nullptr;
+         Nvtx force_first_eks("Init EKS", kYellow, true);
+         ini = 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<::nvtx::Nvtx> NVTX_PRIVATE_NAME(nvtx) = \
+      ::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__);                                                  \
+   ::nvtx::Nvtx::Stack().push(::nvtx::Nvtx::Set(NVTX_FLF, COLOR)); \
+   ::nvtx::Nvtx::Stack().top()->operator()(__VA_ARGS__)
+
+// Begin/End NVTX event traces
+#define NVTX_BEGIN(...) NVTX_BEGIN_(NVTX_COLOR, __VA_ARGS__);
+
+#define NVTX_END(...)                        \
+   ::nvtx::Nvtx::Stack().top().reset(); \
+   ::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 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

tests/unit/CMakeLists.txt

@@ -17,6 +17,7 @@ include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 # The following list can be updated using (in bash):
 #    for d in general linalg mesh fem enzyme; do ls -1 $d/*.cpp; done
 set(UNIT_TESTS_SRCS
+  dfem/test_diffusion.cpp
   general/test_array.cpp
   general/test_arrays_by_name.cpp
   general/test_error.cpp

tests/unit/dfem/test_diffusion.cpp

@@ -0,0 +1,240 @@
+// 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.
+
+#include "unit_tests.hpp"
+
+#include "mfem.hpp"
+#include "fem/dfem/doperator.hpp"
+#include "linalg/tensor.hpp"
+
+using namespace mfem;
+using mfem::internal::tensor;
+using DOperator = DifferentiableOperator;
+
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kAquamarine
+#include "general/nvtx.hpp"
+
+namespace dfem_pa_kernels
+{
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM> struct Diffusion
+{
+   using vecd_t = tensor<real_t, DIM>;
+   using matd_t = tensor<real_t, DIM, DIM>;
+
+   struct MFApply
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const vecd_t &dudxi,
+                                              const real_t &rho,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         const auto invJ = inv(J), TinJ = transpose(invJ);
+         return mfem::tuple{ (dudxi * invJ) * TinJ * det(J) * w * rho };
+      }
+   };
+
+   struct PASetup
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const real_t &u,
+                                              const real_t &rho,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         return mfem::tuple{ inv(J) * transpose(inv(J)) * det(J) * w * rho };
+      }
+   };
+
+   struct PAApply
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const vecd_t &dudxi,
+                                              const matd_t &q) const
+      {
+         return mfem::tuple{ q * dudxi };
+      };
+   };
+};
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void DFemDiffusion(const char *filename, int p, const int r)
+{
+   CAPTURE(filename, DIM, p, r);
+
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   MFEM_VERIFY(pmesh.Dimension() == DIM, "Mesh dimension mismatch");
+
+   pmesh.EnsureNodes();
+   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
+   p = std::max(p, pmesh.GetNodalFESpace()->GetMaxElementOrder());
+   smesh.Clear();
+
+   Array<int> all_domain_attr;
+   if (pmesh.bdr_attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.bdr_attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace pfes(&pmesh, &fec);
+   ParFiniteElementSpace *mfes = nodes->ParFESpace();
+
+   const int NE = pfes.GetNE(), d1d(p + 1), q = 2 * p + r;
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), q);
+   const int q1d(IntRules.Get(Geometry::SEGMENT, ir->GetOrder()).GetNPoints());
+   MFEM_VERIFY(d1d <= q1d, "q1d should be >= d1d");
+   MFEM_VERIFY(NE > 0, "Mesh with no elements is not yet supported!");
+
+   ParGridFunction x(&pfes), y(&pfes), z(&pfes);
+
+   x.Randomize(1);
+   x.SetTrueVector();
+   x.SetFromTrueVector();
+
+   auto rho = [](const Vector &xyz)
+   {
+      const real_t x = xyz(0), y = xyz(1), z = DIM == 3 ? xyz(2) : 0.0;
+      real_t r = M_PI * pow(x, 2);
+      if (DIM >= 2) { r += pow(y, 3); }
+      if (DIM >= 3) { r += pow(z, 4); }
+      return r;
+   };
+   FunctionCoefficient rho_coeff(rho);
+
+   ParBilinearForm blf_fa(&pfes);
+   blf_fa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
+   blf_fa.Assemble();
+   blf_fa.Finalize();
+
+   SECTION("Partial assembly")
+   {
+      dbg("Partial assembly");
+      ParBilinearForm blf_pa(&pfes);
+      blf_pa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
+      blf_pa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      blf_pa.Assemble();
+      blf_pa.Mult(x, z);
+
+      blf_fa.Mult(x, y);
+      y -= z;
+      REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   QuadratureSpace qs(pmesh, *ir);
+   CoefficientVector rho_coeff_cv(rho_coeff, qs);
+   MFEM_VERIFY(rho_coeff_cv.GetVDim() == 1, "Coefficient should be scalar");
+   MFEM_VERIFY(rho_coeff_cv.Size() == q1d * q1d * (DIM == 3 ? q1d : 1) * NE, "");
+
+   const int rho_local_size = 1;
+   const int rho_elem_size(rho_local_size * ir->GetNPoints());
+   const int rho_total_size(rho_elem_size * NE);
+   ParametricSpace rho_ps(DIM, rho_local_size, rho_elem_size, rho_total_size,
+                          DIM == 3 ? d1d : d1d * d1d, // 🔥 2D workaround
+                          DIM == 3 ? q1d : q1d * q1d);
+
+   static constexpr int U = 0, Coords = 1, Rho = 3;
+   const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
+
+   SECTION("DFEM Matrix free")
+   {
+      DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
+      typename Diffusion<DIM>::MFApply mf_apply_qf;
+      dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                 mfem::tuple{ Gradient<U>{}, None<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                 mfem::tuple{ Gradient<U>{} }, *ir,
+                                 all_domain_attr);
+      dop_mf.SetParameters({ &rho_coeff_cv, nodes });
+      dop_mf.Mult(x, z);
+      z.SetTrueVector(), z.SetFromTrueVector();
+      blf_fa.Mult(x, y);
+      y.SetTrueVector(), y.SetFromTrueVector();
+      y -= z;
+      REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   SECTION("DFEM Partial assembly")
+   {
+      static constexpr int QData = 2;
+      const int qd_local_size = DIM * DIM;
+      const int qd_elem_size(qd_local_size * ir->GetNPoints());
+      const int qd_total_size(qd_elem_size * NE);
+      ParametricSpace qd_ps(DIM, qd_local_size, qd_elem_size, qd_total_size,
+                            DIM == 3 ? d1d : d1d * d1d, // 🔥 2D workaround
+                            DIM == 3 ? q1d : q1d * q1d);
+      ParametricFunction qdata(qd_ps);
+      qdata.UseDevice(true);
+
+      DOperator dSetup(sol, {{Rho, &rho_ps}, {Coords, mfes}, {QData, &qd_ps}}, pmesh);
+      typename Diffusion<DIM>::PASetup pa_setup_qf;
+      dSetup.AddDomainIntegrator(
+         pa_setup_qf,
+         mfem::tuple{ None<U>{}, None<Rho>{}, Gradient<Coords>{}, Weight{} },
+         mfem::tuple{ None<QData>{} }, *ir, all_domain_attr);
+      dSetup.SetParameters({ &rho_coeff_cv, nodes, &qdata });
+      pfes.GetRestrictionMatrix()->Mult(x, x.GetTrueVector());
+      dSetup.Mult(x.GetTrueVector(), qdata);
+
+      DOperator dop_pa(sol, { { QData, &qd_ps } }, pmesh);
+      typename Diffusion<DIM>::PAApply pa_apply_qf;
+      dop_pa.AddDomainIntegrator(pa_apply_qf,
+                                 mfem::tuple{ Gradient<U>{}, None<QData>{} },
+                                 mfem::tuple{ Gradient<U>{} },
+                                 *ir, all_domain_attr);
+      dop_pa.SetParameters({ &qdata });
+      dop_pa.Mult(x, z);
+      z.SetTrueVector(), z.SetFromTrueVector();
+      blf_fa.Mult(x, y);
+      y.SetTrueVector(), y.SetFromTrueVector();
+      y -= z;
+      REQUIRE(y.Normlinf() == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+}
+
+///////////////////////////////////////////////////////////////////////////////
+TEST_CASE("DFEM Diffusion", "[Parallel][DFEM]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 1 : GENERATE(1, 2, 3);
+   const auto r = !all_tests ? 0 : GENERATE(0, 1, 2, 3);
+
+   SECTION("2D p=" + std::to_string(p) + " r=" + std::to_string(r))
+   {
+      const auto filename =
+         GENERATE("../../data/star.mesh",
+                  "../../data/star-q3.mesh",
+                  "../../data/rt-2d-q3.mesh",
+                  "../../data/inline-quad.mesh",
+                  "../../data/periodic-square.mesh");
+      DFemDiffusion<2>(filename, p, r);
+   }
+
+   SECTION("3D p=" + std::to_string(p) + " r=" + std::to_string(r))
+   {
+      const auto filename =
+         GENERATE("../../data/fichera.mesh",
+                  "../../data/fichera-q3.mesh",
+                  "../../data/inline-hex.mesh",
+                  "../../data/toroid-hex.mesh",
+                  "../../data/periodic-cube.mesh");
+      DFemDiffusion<3>(filename, p, r);
+   }
+}
+
+} // namespace dfem_pa_kernels

tests/unit/punit_test_main.cpp

@@ -18,6 +18,10 @@
 #error "This test should be disabled without MFEM_USE_MPI!"
 #endif
 
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kOrange
+#include "general/nvtx.hpp"
+
 int main(int argc, char *argv[])
 {
 #ifdef MFEM_USE_SINGLE
@@ -32,6 +36,8 @@ int main(int argc, char *argv[])
 #endif
    mfem::Device device("cpu"); // make sure hypre runs on CPU, if possible
 
+   dbg();
+
    // Only run tests that are labeled with Parallel.
    return RunCatchSession(argc, argv, {"[Parallel]"}, Root());
 }