Use external PA data

version 1: PA new
version 2: PA new (layout by vdim)
version 3: MF ∂fem
version 4: PA ∂fem
version 5: PA ∂fem new
version 6: Auto PA ∂fem
version 7: Auto PA ∂fem new

.gitignore

@@ -413,8 +413,8 @@ miniapps/diag-smoothers/mg-abs-l1-jacobi
 tests/unit/output_meshes
 tests/unit/unit_tests
 tests/unit/punit_tests
-tests/unit/gpu_unit_tests
-tests/unit/pgpu_unit_tests
+tests/unit/cunit_tests
+tests/unit/pcunit_tests
 tests/unit/sedov_tests_*
 tests/unit/psedov_tests_*
 tests/unit/tmop_pa_tests_*

config/config.hpp

@@ -61,6 +61,15 @@ constexpr real_t operator""_r(unsigned long long v)
 
 } // namespace mfem
 
+// Defines 'MFEM_THROW' depending if we are compiled with -fno-exceptions
+#if defined(__EXCEPTIONS) || defined(__cpp_exceptions) || defined(_CPPUNWIND)
+#define MFEM_HAS_EXCEPTIONS
+#define MFEM_THROW(exception_type, msg) throw exception_type(msg)
+#else
+#undef MFEM_HAS_EXCEPTIONS
+#define MFEM_THROW(...)
+#endif
+
 // Return value for main function in examples that should be skipped by testing
 // in some case. This return value prevents failures in testing.
 #define MFEM_SKIP_RETURN_VALUE 242

fem/bilininteg.hpp

@@ -2176,18 +2176,22 @@ class DiffusionIntegrator: public BilinearFormIntegrator
 {
 public:
 
-   using ApplyKernelType = void(*)(const int, const bool, const Array<real_t>&,
-                                   const Array<real_t>&, const Array<real_t>&,
-                                   const Array<real_t>&,
-                                   const Vector&, const Vector&,
-                                   Vector&, const int, const int);
+   using DiffusionApplyKernelType = void(*)(const int, const bool,
+                                            const Array<real_t>&,
+                                            const Array<real_t>&, const Array<real_t>&,
+                                            const Array<real_t>&,
+                                            const Vector&, const Vector&,
+                                            Vector&, const int, const int);
 
-   using DiagonalKernelType = void(*)(const int, const bool, const Array<real_t>&,
-                                      const Array<real_t>&, const Vector&, Vector&,
-                                      const int, const int);
+   using DiffusionDiagonalKernelType = void(*)(const int, const bool,
+                                               const Array<real_t>&,
+                                               const Array<real_t>&, const Vector&, Vector&,
+                                               const int, const int);
 
-   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
-   MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
+   MFEM_REGISTER_KERNELS(DiffusionApplyPAKernel, DiffusionApplyKernelType,
+                         (int, int, int));
+   MFEM_REGISTER_KERNELS(DiffusionDiagonalPAKernel, DiffusionDiagonalKernelType,
+                         (int, int, int));
    struct Kernels { Kernels(); };
 
 protected:
@@ -2207,6 +2211,7 @@ private:
    const FiniteElementSpace *fespace;
    const DofToQuad *maps;         ///< Not owned
    const GeometricFactors *geom;  ///< Not owned
+public:
    int dim, ne, dofs1D, quad1D;
    Vector pa_data;
    bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
@@ -2348,8 +2353,8 @@ public:
    template <int DIM, int D1D, int Q1D>
    static void AddSpecialization()
    {
-      ApplyPAKernels::Specialization<DIM,D1D,Q1D>::Add();
-      DiagonalPAKernels::Specialization<DIM,D1D,Q1D>::Add();
+      DiffusionApplyPAKernel::Specialization<DIM,D1D,Q1D>::Add();
+      DiffusionDiagonalPAKernel::Specialization<DIM,D1D,Q1D>::Add();
    }
 protected:
    const IntegrationRule* GetDefaultIntegrationRule(

fem/dfem/action.hpp

@@ -0,0 +1,488 @@
+// 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
+
+#include <cstddef>
+#include <utility>
+
+#include "fem/kernels.hpp"
+#include "fem/kernel_dispatch.hpp"
+
+#include "util.hpp"
+
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kOrchid
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
+using restriction_callback_t =
+   std::function<void(std::vector<mfem::Vector> &,
+                      const std::vector<mfem::Vector> &,
+                      std::vector<mfem::Vector> &)>;
+
+template<typename T, typename = void>
+struct GetTensorDim
+{
+   static constexpr int ndim = 0;
+};
+
+template<typename T>
+struct GetTensorDim<T, std::void_t<decltype(T::ndim)>>
+{
+   static constexpr int ndim = T::ndim;
+};
+
+template<typename T>
+using TensorDim = GetTensorDim<std::remove_cv_t<T>>;
+
+namespace mfem::future
+{
+
+template <std::size_t num_fields>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+std::array<DeviceTensor<1>, num_fields>
+load_field_address(const std::array<int, num_fields> &sizes,
+                   const std::array<DeviceTensor<2>, num_fields> &fields_e,
+                   const int &entity_idx)
+{
+   std::array<DeviceTensor<1>, num_fields> f;
+   for_constexpr<num_fields>([&](auto field_idx)
+   {
+      f[field_idx] =
+         DeviceTensor<1>(&fields_e[field_idx](0, entity_idx), sizes[field_idx]);
+   });
+   return f;
+}
+
+template <std::size_t N>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+std::array<real_t*, N>
+load_field_e_ptr(const std::array<DeviceTensor<2>, N> &fields_e,
+                 const int e)
+{
+   std::array<real_t*, N> f;
+   for_constexpr<N>([&](auto i) { f[i] = &fields_e[i](0, e); });
+   return f;
+}
+
+namespace qf
+{
+
+template <typename reg_t, typename T, int n>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void process_qf_result_from_reg(reg_t &r0,
+                                const int qx, const int qy, const int qz,
+                                const tensor<T, n> &v)
+{
+   r0[0][qz][qy][qx] = v[0];
+   r0[1][qz][qy][qx] = v[1];
+   r0[2][qz][qy][qx] = v[2];
+}
+
+template <int T_Q1D,
+          size_t num_args,
+          typename reg_t,
+          typename qfunc_t,
+          typename args_ts>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void apply_kernel(reg_t &r0, reg_t &r1,
+                  real_t *r2, const int Q1D,
+                  const int qx, const int qy, const int qz,
+                  const qfunc_t &qfunc,
+                  args_ts &args)
+{
+   if constexpr (num_args == 2)
+   {
+      // ∇u
+      tensor<real_t, 3> &arg_0 = get<0>(args);
+      arg_0[0] = r1[0][qz][qy][qx];
+      arg_0[1] = r1[1][qz][qy][qx];
+      arg_0[2] = r1[2][qz][qy][qx];
+
+      // D (PA data)
+      alignas(64) tensor<real_t, 3, 3> &arg_1 = get<1>(args);
+
+      if constexpr (T_Q1D > 0)
+      {
+         auto *D = (real_t (*)[T_Q1D][T_Q1D][3][3]) r2;
+         for (int j = 0; j < 3; j++)
+         {
+            for (int k = 0; k < 3; k++)
+            {
+               arg_1[k][j] = D[qx][qy][qz][k][j];
+            }
+         }
+      }
+      else
+      {
+         // dbg("Q1D:{}", Q1D);
+         const auto D = Reshape(r2, 3, 3, Q1D, Q1D, Q1D);
+         for (int j = 0; j < 3; j++)
+         {
+            for (int k = 0; k < 3; k++)
+            {
+               arg_1[k][j] = D(j, k, qz, qy, qx);
+            }
+         }
+      }
+   }
+
+   const auto r = get<0>(apply(qfunc, args));
+
+
+   if constexpr (decltype(r)::ndim == 1)
+      // if constexpr (TensorDim<decltype(r)>::ndim == 1)
+   {
+      // process_qf_result_from_reg(r0, qx, qy, qz, r);
+      r0[0][qz][qy][qx] = r[0];
+      r0[1][qz][qy][qx] = r[1];
+      r0[2][qz][qy][qx] = r[2];
+   }
+   if constexpr (TensorDim<decltype(r)>::ndim == 0)
+   {
+      MFEM_ABORT("qfunc returned a scalar, expected a vector");
+   }
+}
+
+} // namespace qf
+
+template<size_t num_fields,
+         size_t num_inputs,
+         size_t num_outputs,
+         typename qfunc_t,
+         typename input_t,
+         typename output_fop_t>
+class NewActionCallback
+{
+   const restriction_callback_t restriction_cb;
+   const qfunc_t qfunc;
+   input_t &inputs;
+   const std::array<int, num_inputs> input_to_field;
+   const std::array<DofToQuadMap, num_inputs> input_dtq_maps;
+   const std::array<DofToQuadMap, num_outputs> output_dtq_maps;
+   const int num_entities;
+   const int test_vdim;
+   const int num_test_dof;
+   const int dimension;
+   const int q1d_;
+   const ThreadBlocks thread_blocks;
+   SharedMemoryInfo<num_fields, num_inputs, num_outputs> shmem_info;
+   Array<int> &elem_attributes;
+   const output_fop_t &output_fop;
+   const Array<int> &domain_attributes;
+   // refs
+   std::vector<Vector> &fields_e;
+   Vector &residual_e;
+   std::function<void(Vector &, Vector &)> output_restriction_transpose;
+   // args
+   std::vector<Vector> &solutions_l;
+   const std::vector<Vector> &parameters_l;
+   Vector &residual_l;
+
+public:
+   inline MFEM_ALWAYS_INLINE
+   NewActionCallback(const bool use_kernels_specialization,
+                     const restriction_callback_t restriction_cb,
+                     const qfunc_t qfunc,
+                     input_t &inputs,
+                     const std::array<int, num_inputs> input_to_field,
+                     const std::array<DofToQuadMap, num_inputs> input_dtq_maps,
+                     const std::array<DofToQuadMap, num_outputs> output_dtq_maps,
+                     const int num_entities,
+                     const int test_vdim,
+                     const int num_test_dof,
+                     const int dimension,
+                     const int q1d,
+                     const ThreadBlocks thread_blocks,
+                     SharedMemoryInfo<num_fields, num_inputs, num_outputs> shmem_info,
+                     Array<int> &elem_attributes,
+                     const output_fop_t &output_fop,
+                     const Array<int> &domain_attributes,
+                     // refs
+                     std::vector<Vector> &fields_e,
+                     Vector &residual_e,
+                     std::function<void(Vector &, Vector &)> output_restriction_transpose,
+                     // args
+                     std::vector<Vector> &solutions_l,
+                     const std::vector<Vector> &parameters_l,
+                     Vector &residual_l):
+      restriction_cb(restriction_cb),
+      qfunc(qfunc),
+      inputs(inputs),
+      input_to_field(input_to_field),
+      input_dtq_maps(input_dtq_maps),
+      output_dtq_maps(output_dtq_maps),
+      num_entities(num_entities),
+      test_vdim(test_vdim),
+      num_test_dof(num_test_dof),
+      dimension(dimension),
+      q1d_(q1d),
+      thread_blocks(thread_blocks),
+      shmem_info(shmem_info),
+      elem_attributes(elem_attributes),
+      output_fop(output_fop),
+      domain_attributes(domain_attributes),
+      fields_e(fields_e),
+      residual_e(residual_e),
+      output_restriction_transpose(std::move(output_restriction_transpose)),
+      solutions_l(solutions_l),
+      parameters_l(parameters_l),
+      residual_l(residual_l)
+   {
+      if (!use_kernels_specialization) { return; }
+      NewActionCallbackKernels::template Specialization<3,4>::Add();
+      // NewActionCallbackKernels::template Specialization<4,5>::Add();
+      NewActionCallbackKernels::template Specialization<5,6>::Add();
+      // NewActionCallbackKernels::template Specialization<6,7>::Add();
+      NewActionCallbackKernels::template Specialization<7,8>::Add();
+      // NewActionCallbackKernels::template Specialization<8,9>::Add();
+      // NewActionCallbackKernels::template Specialization<9,10>::Add();
+   }
+
+   template<int T_D1D = 0, int T_Q1D = 0>
+   static inline MFEM_ALWAYS_INLINE
+   void action_callback_new(const restriction_callback_t &restriction_cb,
+                            const qfunc_t &qfunc,
+                            input_t &inputs,
+                            const std::array<int, num_inputs> &input_to_field,
+                            const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
+                            const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
+                            const int &ne,
+                            const int &test_vdim,
+                            const int &num_test_dof,
+                            const int &dimension,
+                            //   const int q1d,
+                            const ThreadBlocks &thread_blocks,
+                            SharedMemoryInfo<num_fields, num_inputs, num_outputs> &shmem_info,
+                            Array<int> &elem_attributes,
+                            const output_fop_t &output_fop,
+                            const Array<int> &domain_attributes,
+                            // refs
+                            std::vector<Vector> &fields_e,
+                            Vector &residual_e,
+                            std::function<void(Vector &, Vector &)> output_restriction_transpose,
+                            // args
+                            std::vector<Vector> &solutions_l,
+                            const std::vector<Vector> &parameters_l,
+                            Vector &residual_l,
+                            // fallback arguments
+                            const int d1d, const int q1d)
+   {
+      db1();
+      assert(dimension == 3);
+
+      // constexpr int D1D = T_D1D;// ? T_D1D : d1d; // 🔥
+      // constexpr int Q1D = T_Q1D;// ? T_Q1D : q1d; // 🔥
+
+      constexpr int DIM = 3;
+      constexpr int MQ1 = T_Q1D > 0 ? kernels::internal::SetMaxOf(T_Q1D) : 32;
+      constexpr int MD1 = T_D1D > 0 ? kernels::internal::SetMaxOf(T_D1D) : 32;
+      // constexpr int MQ1 = T_Q1D ? T_Q1D : 32;
+      // constexpr int MD1 = T_D1D ? T_D1D : 32;
+
+      // types
+      using qf_signature =
+         typename create_function_signature<decltype(&qfunc_t::operator())>::type;
+      using qf_param_ts = typename qf_signature::parameter_ts;
+
+      // #warning 0xDEADBEEF
+      // auto solutions_l_pi = solutions_l;
+      // solutions_l_pi[0].Randomize(0xDEADBEEF);
+
+      restriction_cb(solutions_l, parameters_l, fields_e); // 🔥 to inline
+      /*for (size_t i = 0; i < 1; i++)
+      {
+         dbl("fields_e: size:{} dot:{}\n", fields_e[i].Size(), fields_e[i]*fields_e[i]);
+         for (int j = 0; j < fields_e[i].Size(); j++)
+         {
+            dba("{} ", fields_e[i][j]);
+         }
+         dbc();
+      }*/
+
+      residual_e = 0.0;
+      auto ye = Reshape(residual_e.ReadWrite(), num_test_dof, test_vdim, ne);
+
+      // std::array<DeviceTensor<2>, num_fields>: (field_sizes, num_entities)
+      auto wrapped_fields_e = wrap_fields(fields_e,
+                                          shmem_info.field_sizes,
+                                          ne);
+
+      const bool has_attr = domain_attributes.Size() > 0;
+      const auto d_domain_attr = domain_attributes.Read();
+      const auto d_elem_attr = elem_attributes.Read();
+
+      // db1("forall");
+      forall([=] MFEM_HOST_DEVICE (int e, void *)
+      {
+         const int D1D = T_D1D ? T_D1D : d1d;
+         const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+         // this could be optimized out
+         if (has_attr && !d_domain_attr[d_elem_attr[e] - 1]) { return; }
+
+         alignas(64) MFEM_SHARED real_t smem[MQ1][MQ1];
+         alignas(64) MFEM_SHARED real_t sB[MD1][MQ1];
+         alignas(64) MFEM_SHARED real_t sG[MD1][MQ1];
+         alignas(64) kernels::internal::d_regs3d_t<DIM, MQ1> r0, r1;
+         alignas(64) real_t *r2;
+
+         const auto fields_e_ptr = load_field_e_ptr(wrapped_fields_e, e);
+
+         // Interpolate
+         // const auto dummy_field_weight = DeviceTensor<1>(nullptr, 0);
+         for_constexpr<num_inputs>([&](auto i)
+         {
+            const auto input = get<i>(inputs);
+            using field_operator_t = std::decay_t<decltype(input)>;
+
+            if constexpr (is_gradient_fop<field_operator_t>::value) // Grad
+            {
+               // const int vdim = input.vdim;
+               constexpr int VDIM = 1; // 🔥
+               const real_t *field_e_r = fields_e_ptr[input_to_field[i]];
+               const auto field = Reshape(field_e_r, D1D, D1D, D1D, VDIM);
+               const auto dtq = input_dtq_maps[i];
+               const auto B = dtq.B, G = dtq.G;
+               // db1("B:{} {} {} {} {} {}", B[0], B[1], B[2], B[3], B[4], B[5]);
+               // db1("G:{} {} {} {} {} {}", G[0], G[1], G[2], G[3], G[4], G[5]);
+               kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
+               kernels::internal::LoadMatrix(D1D, Q1D, G, sG);
+               // const auto B = reinterpret_cast<const real_t (*)[MQ1]>(Bi[i]);
+               // const auto G = reinterpret_cast<const real_t (*)[MQ1]>(Gi[i]);
+               for (int c = 0; c < VDIM; c++)
+               {
+                  kernels::internal::LoadDofs3d(D1D, c, field, r0);
+                  kernels::internal::Grad3d(D1D, Q1D, smem, sB, sG, r0, r1, c);
+               }
+            }
+
+            if constexpr (is_identity_fop<field_operator_t>::value) // Identity
+            {
+               // db1("Identity");
+               r2 = fields_e_ptr[input_to_field[i]];
+            }
+         }); // for_constexpr<num_inputs>
+
+         // db1("Now calling qfunction");
+         auto qf_args = decay_tuple<qf_param_ts> {};
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            // MFEM_UNROLL(MQ1)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               // MFEM_UNROLL(MQ1)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  // dbg("r1: ({},{},{}) {} {} {}", qx, qy, qz,
+                  //     r1[0][qz][qy][qx], r1[1][qz][qy][qx], r1[2][qz][qy][qx]);
+                  qf::apply_kernel<T_Q1D, num_inputs>(r0, r1, r2, Q1D, qx, qy, qz,
+                                                      qfunc, qf_args);
+                  // dbg("r0: ({},{},{}) {} {} {}", qx, qy, qz,
+                  //     r0[0][qz][qy][qx], r0[1][qz][qy][qx], r0[2][qz][qy][qx]);
+               }
+            }
+         }
+
+         // db1("Integrate");
+         auto y = Reshape(&ye(0, 0, e), num_test_dof, test_vdim);
+         if constexpr (is_gradient_fop<std::decay_t<output_fop_t>>::value) // GradientT
+         {
+            // const int vdim = output_fop.vdim;
+            constexpr int VDIM = 1; // 🔥
+            auto yd = Reshape(&y(0, 0), D1D, D1D, D1D, VDIM);
+            // const auto B = reinterpret_cast<const real_t (*)[MQ1]>(Bo);
+            // const auto G = reinterpret_cast<const real_t (*)[MQ1]>(Go);
+            for (int c = 0; c < VDIM; c++)
+            {
+               kernels::internal::GradTranspose3d(D1D, Q1D, smem, sB, sG, r0, r1, c);
+               kernels::internal::WriteDofs3d(D1D, c, r1, yd);
+            }
+         }
+      },
+      ne,
+      thread_blocks,
+      0,
+      nullptr);
+
+      // dbg("residual_e: size:{} dot:{}",
+      //     residual_e.Size(), residual_e * residual_e);
+      output_restriction_transpose(residual_e, residual_l);  // 🔥 to inline
+   }
+
+   using NewActionCallbackType = decltype(
+                                    &NewActionCallback::action_callback_new<>);
+   MFEM_REGISTER_KERNELS(NewActionCallbackKernels, NewActionCallbackType, (int,
+                                                                           int));
+
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   void Apply(const int d1d, const int q1d)
+   {
+      NewActionCallbackKernels::Run(d1d, q1d,
+                                    // args
+                                    restriction_cb,
+                                    qfunc,
+                                    inputs,
+                                    input_to_field,
+                                    input_dtq_maps,
+                                    output_dtq_maps,
+                                    num_entities,
+                                    test_vdim,
+                                    num_test_dof,
+                                    dimension,
+                                    thread_blocks,
+                                    shmem_info,
+                                    elem_attributes,
+                                    output_fop,
+                                    domain_attributes,
+                                    fields_e,
+                                    residual_e,
+                                    output_restriction_transpose,
+                                    solutions_l,
+                                    parameters_l,
+                                    residual_l,
+                                    // fallback arguments
+                                    d1d, q1d);
+   }
+};
+
+template<size_t num_fields,
+         size_t num_inputs,
+         size_t num_outputs,
+         typename qfunc_t,
+         typename input_t,
+         typename output_fop_t>
+template<int D1D, int Q1D>
+inline MFEM_ALWAYS_INLINE
+typename NewActionCallback<num_fields, num_inputs, num_outputs, qfunc_t, input_t, output_fop_t>::NewActionCallbackType
+NewActionCallback<num_fields, num_inputs, num_outputs, qfunc_t, input_t, output_fop_t>::NewActionCallbackKernels::Kernel()
+{
+   return action_callback_new<D1D, Q1D>;
+}
+
+template<size_t num_fields,
+         size_t num_inputs,
+         size_t num_outputs,
+         typename qfunc_t,
+         typename input_t,
+         typename output_fop_t>
+inline MFEM_ALWAYS_INLINE
+typename NewActionCallback<num_fields, num_inputs, num_outputs, qfunc_t, input_t, output_fop_t>::NewActionCallbackType
+NewActionCallback<num_fields, num_inputs, num_outputs, qfunc_t, input_t, output_fop_t>::NewActionCallbackKernels::Fallback
+(int d1d, int q1d)
+{
+   return action_callback_new<>;
+}
+
+} // namespace mfem::future

fem/dfem/autopa.hpp

@@ -0,0 +1,406 @@
+// 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
+
+#include <cstddef>
+
+#include "fem/kernels.hpp"
+#include "fem/kernel_dispatch.hpp"
+
+#include "util.hpp"
+
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kSandyBrown
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
+namespace mfem::future
+{
+
+template<size_t num_inputs,
+         size_t num_outputs,
+         typename input_t,
+         size_t num_fields,
+         typename output_fop_t>
+class NewAutoActionCallback
+{
+   input_t &inputs;
+   const std::array<DofToQuadMap, num_inputs> &input_dtq_maps;
+   const std::array<DofToQuadMap, num_outputs> &output_dtq_maps;
+   const int num_entities;
+   const int test_vdim;
+   const int num_test_dof;
+   const int dimension;
+   const FieldDescriptor direction;
+   const int test_op_dim;
+   const int trial_vdim;
+   const int total_trial_op_dim;
+   const int num_qp;
+   Vector &dependent_inputs_trial_op_dim;
+   const int num_dependent_inputs;
+   const ElementDofOrdering element_dof_ordering;
+   const int q1d_;
+   const ThreadBlocks &thread_blocks;
+   Vector &shmem_cache;
+   SharedMemoryInfo<num_fields, num_inputs, num_outputs> &shmem_info;
+   Array<int> &elem_attributes;
+   const output_fop_t &output_fop;
+   const Array<int> &domain_attributes;
+   const DeviceTensor<1, const real_t> &ir_weights;
+   const bool use_sum_factorization;
+   const std::array<bool, num_inputs> &input_is_dependent;
+   Vector &direction_e;
+   Vector &derivative_action_e;
+   const real_t *pa_data; // pointer to the PA data
+   // refs
+   Vector &qpdc_mem; // derivative_qp_caches
+   std::function<void(Vector &, Vector &)> &output_restriction_transpose;
+   // args
+   std::vector<Vector> &f_e;
+   const Vector &dir_l;
+   Vector &der_action_l;
+
+public:
+   inline MFEM_ALWAYS_INLINE
+   NewAutoActionCallback(const bool &use_kernels_specialization,
+                         input_t &inputs,
+                         const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
+                         const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
+                         const int dimension,
+                         const int num_entities,
+                         const int test_vdim,
+                         const int num_test_dof,
+                         const FieldDescriptor &direction,
+                         const int test_op_dim,
+                         const int trial_vdim,
+                         const int total_trial_op_dim,
+                         const int num_qp,
+                         Vector &dependent_inputs_trial_op_dim,
+                         const int num_dependent_inputs,
+                         const ElementDofOrdering element_dof_ordering,
+                         const int q1d,
+                         const ThreadBlocks &thread_blocks,
+                         Vector &shmem_cache,
+                         SharedMemoryInfo<num_fields, num_inputs, num_outputs> &shmem_info,
+                         Array<int> &elem_attributes,
+                         const output_fop_t &output_fop,
+                         const Array<int> &domain_attributes,
+                         const DeviceTensor<1, const real_t> &ir_weights,
+                         const bool use_sum_factorization,
+                         const std::array<bool, num_inputs> &input_is_dependent,
+                         Vector &direction_e,
+                         Vector &derivative_action_e,
+                         const real_t *pa_data,
+                         // refs
+                         Vector &qpdc_mem,
+                         std::function<void(Vector &, Vector &)> &output_restriction_transpose,
+                         // args
+                         std::vector<Vector> &f_e,
+                         const Vector &dir_l,
+                         Vector &der_action_l):
+      inputs(inputs),
+      input_dtq_maps(input_dtq_maps),
+      output_dtq_maps(output_dtq_maps),
+      num_entities(num_entities),
+      test_vdim(test_vdim),
+      num_test_dof(num_test_dof),
+      dimension(dimension),
+      direction(direction),
+      test_op_dim(test_op_dim),
+      trial_vdim(trial_vdim),
+      total_trial_op_dim(total_trial_op_dim),
+      num_qp(num_qp),
+      dependent_inputs_trial_op_dim(dependent_inputs_trial_op_dim),
+      num_dependent_inputs(num_dependent_inputs),
+      element_dof_ordering(element_dof_ordering),
+      q1d_(q1d),
+      thread_blocks(thread_blocks),
+      shmem_cache(shmem_cache),
+      shmem_info(shmem_info),
+      elem_attributes(elem_attributes),
+      output_fop(output_fop),
+      domain_attributes(domain_attributes),
+      ir_weights(ir_weights),
+      use_sum_factorization(use_sum_factorization),
+      input_is_dependent(input_is_dependent),
+      direction_e(direction_e),
+      derivative_action_e(derivative_action_e),
+      pa_data(pa_data),
+      qpdc_mem(qpdc_mem),
+      output_restriction_transpose(output_restriction_transpose),
+      f_e(f_e),
+      dir_l(dir_l),
+      der_action_l(der_action_l)
+   {
+      if (!use_kernels_specialization) { return; }
+      NewAutoActionCallbackKernels::template Specialization<3,4>::Add();
+      // NewAutoActionCallbackKernels::template Specialization<4,5>::Add();
+      NewAutoActionCallbackKernels::template Specialization<5,6>::Add();
+      // NewAutoActionCallbackKernels::template Specialization<6,7>::Add();
+      NewAutoActionCallbackKernels::template Specialization<7,8>::Add();
+      // NewAutoActionCallbackKernels::template Specialization<8,9>::Add();
+      // NewAutoActionCallbackKernels::template Specialization<9,10>::Add();
+   }
+
+   template<int T_D1D = 0, int T_Q1D = 0>
+   static inline MFEM_ALWAYS_INLINE
+   void auto_pa_action_callback(const input_t &inputs,
+                                const std::array<DofToQuadMap, num_inputs> input_dtq_maps,
+                                const int ne,
+                                const int test_vdim,
+                                const int num_test_dof,
+                                const int dimension,
+                                const FieldDescriptor &direction,
+                                const int test_op_dim,
+                                const int trial_vdim,
+                                const int total_trial_op_dim,
+                                const int num_qp,
+                                const int num_dependent_inputs,
+                                const ElementDofOrdering ordering,
+                                const ThreadBlocks thread_blocks,
+                                Array<int> elem_attributes,
+                                const output_fop_t output_fop,
+                                const Array<int> domain_attributes,
+                                const bool use_sum_factorization,
+                                Vector &direction_e,
+                                Vector &derivative_action_e,
+                                const real_t *pa_data,
+                                // refs
+                                Vector &qpdc_mem,
+                                std::function<void(Vector &, Vector &)> &or_transpose,
+                                // args
+                                // std::vector<Vector> &f_e, // unused
+                                const Vector &direction_l,
+                                Vector &der_action_l,
+                                // fallback arguments
+                                const int d1d, const int q1d)
+   {
+      db1("d1d:{} q1d:{}", d1d, q1d);
+      db1("T_D1D:{} T_Q1D:{}", T_D1D, T_Q1D);
+      db1("num_qp: {}, ne: {}", num_qp, ne);
+      db1("test_op_dim: {}, test_vdim: {}", test_op_dim, test_vdim);
+      db1("total_trial_op_dim: {}, trial_vdim: {}", total_trial_op_dim, trial_vdim);
+      db1("num_inputs:{} num_fields:{} num_outputs:{}", num_inputs, num_fields,
+          num_outputs);
+      assert(dimension == 3);
+
+      constexpr int DIM = 3;
+      constexpr int MQ1 = T_Q1D > 0 ? T_Q1D : 32;
+      constexpr int MD1 = T_D1D > 0 ? T_D1D : 32;
+
+      assert(ordering == ElementDofOrdering::LEXICOGRAPHIC);
+      // db1("direction_l: size:{} dot:{}", direction_l.Size(),
+      //     direction_l*direction_l);
+      restriction<Entity::Element>(direction, direction_l, direction_e, ordering);
+      // db1("direction_e: size:{} dot:{}", direction_e.Size(),
+      //     direction_e*direction_e);
+      const auto dir_e = Reshape(direction_e.Read(), d1d,d1d,d1d, 1, ne);
+
+      assert(pa_data);
+      // const auto qpdc = Reshape(pa_data ? nullptr :qpdc_mem.Read(),
+      //                           test_vdim, test_op_dim,
+      //                           trial_vdim, total_trial_op_dim,
+      //                           num_qp, ne);
+      const auto DX = Reshape(pa_data, 3, 3, q1d, q1d, q1d, ne);
+
+      // db1("qpdc: size:{} dot:{}", qpdc_mem.Size(), qpdc_mem*qpdc_mem);
+
+      // const auto d_elem_attr = elem_attributes.Read();
+      // const bool has_attr = domain_attributes.Size() > 0;
+      // const auto d_domain_attr = domain_attributes.Read();
+
+      derivative_action_e = 0.0;
+      assert(test_vdim == 1);
+      assert(num_test_dof == d1d*d1d*d1d);
+      auto ye = Reshape(derivative_action_e.ReadWrite(), num_test_dof, test_vdim, ne);
+
+      forall([=] MFEM_HOST_DEVICE (int e, real_t *)
+      {
+         const int D1D = T_D1D ? T_D1D : d1d;
+         const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+         // if (has_attr && !d_domain_attr[d_elem_attr[e] - 1]) { return; }
+
+         MFEM_SHARED real_t smem[MQ1][MQ1];
+         MFEM_SHARED real_t sB[MD1][MQ1];
+         MFEM_SHARED real_t sG[MD1][MQ1];
+         kernels::internal::d_regs3d_t<DIM, MQ1> r0, r1;
+
+         const auto dir_fop = get<0>(inputs);
+         // const int vdim = dir_fop.vdim;
+         const int vd = 0;
+
+         // Interpolate
+         const auto input = get<0>(inputs);
+         using field_operator_t = std::decay_t<decltype(input)>;
+         if constexpr (is_gradient_fop<field_operator_t>::value) // Grad
+         {
+            constexpr int VDIM = 1; // 🔥
+            // const auto dtq = input_dtq_maps[0];
+            const auto dtq = get<0>(input_dtq_maps);
+            const auto B = dtq.B, G = dtq.G;
+            // db1("B:{} {} {} {} {} {}", B[0], B[1], B[2], B[3], B[4], B[5]);
+            // db1("G:{} {} {} {} {} {}", G[0], G[1], G[2], G[3], G[4], G[5]);
+            kernels::internal::LoadMatrix(D1D, Q1D, B, sB);
+            kernels::internal::LoadMatrix(D1D, Q1D, G, sG);
+            for (int c = 0; c < VDIM; c++)
+            {
+               kernels::internal::LoadDofs3d(e, D1D, c, dir_e, r0);
+               kernels::internal::Grad3d(D1D, Q1D, smem, sB, sG, r0, r1, c);
+            }
+         }
+         else if constexpr (is_identity_fop<field_operator_t>::value) // Identity
+         {
+            // db1("Id");
+            assert(false); // not here
+         }
+         else { assert(false); }
+
+
+         // db1("Qfunction");
+         for (int qz = 0; qz < Q1D; ++qz)
+         {
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  // const int q = qx + q1d * (qy + q1d * qz);
+
+                  const real_t u = r1[0][qz][qy][qx];
+                  const real_t v = r1[1][qz][qy][qx];
+                  const real_t w = r1[2][qz][qy][qx];
+
+                  for (int k = 0; k < test_op_dim; k++)
+                  {
+                     // const auto trial_op_dim = dpitod(0, 1);
+
+                     // size_t m_offset = 0;
+                     for (int j = 0; j < trial_vdim; j++)
+                     {
+                        // const real_t val = qpdc(vd, k, j, 0 + m_offset, q, e) * u
+                        //                    + qpdc(vd, k, j, 1 + m_offset, q, e) * v
+                        //                    + qpdc(vd, k, j, 2 + m_offset, q, e) * w;
+                        const real_t val = DX(k, 0, qx, qy, qz, e) * u +
+                                           DX(k, 1, qx, qy, qz, e) * v +
+                                           DX(k, 2, qx, qy, qz, e) * w;
+                        r0[k][qz][qy][qx] = val;
+                     }
+                     // m_offset += 3;//trial_op_dim;
+                  }
+               }
+            }
+         }
+
+         // db1("Integrate");
+         auto y = Reshape(&ye(0, 0, e), num_test_dof, test_vdim);
+         if constexpr (is_gradient_fop<std::decay_t<output_fop_t>>::value) // GradientT
+         {
+            // db1("GradTranspose3d");
+            // const int vdim = output_fop.vdim;
+            constexpr int VDIM = 1; // 🔥
+            auto yd = Reshape(&y(0, 0), D1D, D1D, D1D, VDIM);
+            // const auto B = reinterpret_cast<const real_t (*)[MQ1]>(Bo);
+            // const auto G = reinterpret_cast<const real_t (*)[MQ1]>(Go);
+            for (int c = 0; c < VDIM; c++)
+            {
+               kernels::internal::GradTranspose3d(D1D, Q1D, smem, sB, sG, r0, r1, c);
+               kernels::internal::WriteDofs3d(D1D, c, r1, yd);
+            }
+         }
+         else
+         {
+            assert(false);
+         }
+
+      },
+      ne,
+      thread_blocks,
+      0,
+      nullptr);
+
+      // dbg("derivative_action_e: size:{} dot:{}",
+      //     derivative_action_e.Size(), derivative_action_e * derivative_action_e);
+      or_transpose(derivative_action_e, der_action_l);
+   }
+
+   using NewAutoActionCallbackType =
+      decltype(&NewAutoActionCallback::auto_pa_action_callback<>);
+   MFEM_REGISTER_KERNELS(NewAutoActionCallbackKernels,
+                         NewAutoActionCallbackType, (int, int));
+
+
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   void Apply(const int d1d, const int q1d)
+   {
+      NewAutoActionCallbackKernels::Run(d1d, q1d,
+                                        // args
+                                        inputs,
+                                        input_dtq_maps,
+                                        num_entities,
+                                        test_vdim,
+                                        num_test_dof,
+                                        dimension,
+                                        direction,
+                                        test_op_dim,
+                                        trial_vdim,
+                                        total_trial_op_dim,
+                                        num_qp,
+                                        num_dependent_inputs,
+                                        element_dof_ordering,
+                                        thread_blocks,
+                                        elem_attributes,
+                                        output_fop,
+                                        domain_attributes,
+                                        use_sum_factorization,
+                                        direction_e,
+                                        derivative_action_e,
+                                        pa_data,
+                                        // refs
+                                        qpdc_mem,
+                                        output_restriction_transpose,
+                                        // args
+                                        dir_l,
+                                        der_action_l,
+                                        // fallback arguments
+                                        d1d, q1d);
+   }
+};
+
+template<size_t num_inputs,
+         size_t num_outputs,
+         typename input_t,
+         size_t num_fields,
+         typename output_fop_t>
+template<int D1D, int Q1D>
+inline MFEM_ALWAYS_INLINE
+typename NewAutoActionCallback<num_inputs, num_outputs, input_t, num_fields, output_fop_t>::NewAutoActionCallbackType
+NewAutoActionCallback<num_inputs, num_outputs, input_t, num_fields, output_fop_t>::NewAutoActionCallbackKernels::Kernel()
+{
+   return auto_pa_action_callback<D1D, Q1D>;
+}
+
+template<size_t num_inputs,
+         size_t num_outputs,
+         typename input_t,
+         size_t num_fields,
+         typename output_fop_t>
+inline MFEM_ALWAYS_INLINE
+typename NewAutoActionCallback<num_inputs, num_outputs, input_t, num_fields, output_fop_t>::NewAutoActionCallbackType
+NewAutoActionCallback<num_inputs, num_outputs, input_t, num_fields, output_fop_t>::NewAutoActionCallbackKernels::Fallback
+(int d1d, int q1d)
+{
+   return auto_pa_action_callback<>;
+}
+
+} // namespace mfem::future

fem/dfem/integrate.hpp

@@ -102,7 +102,7 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
    const DeviceTensor<3, real_t> &f,
    const output_t &output,
    const DofToQuadMap &dtq,
-   std::array<DeviceTensor<1>, 6> &scratch_mem)
+   const std::array<DeviceTensor<1>, 6> &scratch_mem)
 {
    [[maybe_unused]] auto B = dtq.B;
    [[maybe_unused]] auto G = dtq.G;
@@ -239,7 +239,7 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
    const DeviceTensor<3, real_t> &f,
    const output_t &output,
    const DofToQuadMap &dtq,
-   std::array<DeviceTensor<1>, 6> &scratch_mem)
+   const std::array<DeviceTensor<1>, 6> &scratch_mem)
 {
    [[maybe_unused]] auto B = dtq.B;
    [[maybe_unused]] auto G = dtq.G;
@@ -258,11 +258,11 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
 
       for (int vd = 0; vd < vdim; vd++)
       {
-         MFEM_FOREACH_THREAD(qy, y, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
          {
-            MFEM_FOREACH_THREAD(dx, x, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
             {
-               MFEM_FOREACH_THREAD(qz, z, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
                {
                   real_t acc = 0.0;
                   for (int qx = 0; qx < q1d; qx++)
@@ -275,11 +275,11 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(dy, y, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
          {
-            MFEM_FOREACH_THREAD(dx, x, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
             {
-               MFEM_FOREACH_THREAD(qz, z, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
                {
                   real_t acc = 0.0;
                   for (int qy = 0; qy < q1d; qy++)
@@ -293,11 +293,11 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
          MFEM_SYNC_THREAD;
 
 
-         MFEM_FOREACH_THREAD(dy, y, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
          {
-            MFEM_FOREACH_THREAD(dx, x, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
             {
-               MFEM_FOREACH_THREAD(dz, z, d1d)
+               MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
                {
                   real_t acc = 0.0;
                   for (int qz = 0; qz < q1d; qz++)
@@ -328,62 +328,69 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
 
       for (int vd = 0; vd < vdim; vd++)
       {
-         MFEM_FOREACH_THREAD(qz, z, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
          {
-            MFEM_FOREACH_THREAD(qy, y, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
-               MFEM_FOREACH_THREAD(dx, x, d1d)
+               MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
                {
-                  real_t uvw[3] = {0.0, 0.0, 0.0};
+                  real_t u = 0.0, v = 0.0, w = 0.0;
                   for (int qx = 0; qx < q1d; qx++)
                   {
-                     uvw[0] += fqp(vd, 0, qx, qy, qz) * G(qx, 0, dx);
-                     uvw[1] += fqp(vd, 1, qx, qy, qz) * B(qx, 0, dx);
-                     uvw[2] += fqp(vd, 2, qx, qy, qz) * B(qx, 0, dx);
+                     const real_t b = B(qx, 0, dx);
+                     const real_t g = G(qx, 0, dx);
+
+                     u += fqp(vd, 0, qx, qy, qz) * g;
+                     v += fqp(vd, 1, qx, qy, qz) * b;
+                     w += fqp(vd, 2, qx, qy, qz) * b;
                   }
-                  s0(qz, qy, dx) = uvw[0];
-                  s1(qz, qy, dx) = uvw[1];
-                  s2(qz, qy, dx) = uvw[2];
+                  s0(qz, qy, dx) = u;
+                  s1(qz, qy, dx) = v;
+                  s2(qz, qy, dx) = w;
                }
             }
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(qz, z, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
          {
-            MFEM_FOREACH_THREAD(dy, y, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
             {
-               MFEM_FOREACH_THREAD(dx, x, d1d)
+               MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
                {
-                  real_t uvw[3] = {0.0, 0.0, 0.0};
+                  real_t u = 0.0, v = 0.0, w = 0.0;
                   for (int qy = 0; qy < q1d; qy++)
                   {
-                     uvw[0] += s0(qz, qy, dx) * B(qy, 0, dy);
-                     uvw[1] += s1(qz, qy, dx) * G(qy, 0, dy);
-                     uvw[2] += s2(qz, qy, dx) * B(qy, 0, dy);
+                     const real_t b = B(qy, 0, dy);
+                     const real_t g = G(qy, 0, dy);
+                     u += s0(qz, qy, dx) * b;
+                     v += s1(qz, qy, dx) * g;
+                     w += s2(qz, qy, dx) * b;
                   }
-                  s3(qz, dy, dx) = uvw[0];
-                  s4(qz, dy, dx) = uvw[1];
-                  s5(qz, dy, dx) = uvw[2];
+                  s3(qz, dy, dx) = u;
+                  s4(qz, dy, dx) = v;
+                  s5(qz, dy, dx) = w;
                }
             }
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(dz, z, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
          {
-            MFEM_FOREACH_THREAD(dy, y, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
             {
-               MFEM_FOREACH_THREAD(dx, x, d1d)
+               MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
                {
-                  real_t uvw[3] = {0.0, 0.0, 0.0};
+                  real_t u = 0.0, v = 0.0, w = 0.0;
                   for (int qz = 0; qz < q1d; qz++)
                   {
-                     uvw[0] += s3(qz, dy, dx) * B(qz, 0, dz);
-                     uvw[1] += s4(qz, dy, dx) * B(qz, 0, dz);
-                     uvw[2] += s5(qz, dy, dx) * G(qz, 0, dz);
+                     const real_t b = B(qz, 0, dz);
+                     const real_t g = G(qz, 0, dz);
+                     u += s3(qz, dy, dx) * b;
+                     v += s4(qz, dy, dx) * b;
+                     w += s5(qz, dy, dx) * g;
                   }
-                  yd(dx, dy, dz, vd) += uvw[0] + uvw[1] + uvw[2];
+                  yd(dx, dy, dz, vd) += u + v + w;
                }
             }
          }
@@ -398,11 +405,11 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
 
       for (int sq = 0; sq < output.size_on_qp; sq++)
       {
-         MFEM_FOREACH_THREAD(qx, x, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
          {
-            MFEM_FOREACH_THREAD(qy, y, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
-               MFEM_FOREACH_THREAD(qz, z, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
                {
                   yqp(sq, qx, qy, qz) = fqp(sq, qx, qy, qz);
                }
@@ -425,7 +432,7 @@ void map_quadrature_data_to_fields(
    const DeviceTensor<3, real_t> &f,
    const output_t &output,
    const DofToQuadMap &dtq,
-   std::array<DeviceTensor<1>, 6> &scratch_mem,
+   const std::array<DeviceTensor<1>, 6> &scratch_mem,
    const int &dimension,
    const bool &use_sum_factorization)
 {

fem/dfem/interpolate.hpp

@@ -11,6 +11,23 @@
 #pragma once
 
 #include "util.hpp"
+#include "fem/kernels.hpp"
+
+///////////////////////////////////////////////////////////////////////////////
+template <class T>
+inline std::enable_if_t<!std::numeric_limits<T>::is_integer, bool>
+AlmostEq(T x, T y, T tolerance = 15.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;
+}
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kPeru
+#include "general/nvtx.hpp"
 
 namespace mfem::future
 {
@@ -30,6 +47,7 @@ void map_field_to_quadrature_data_tensor_product_3d(
 
    if constexpr (is_value_fop<std::decay_t<field_operator_t>>::value)
    {
+      dbg("Value");
       auto [q1d, unused, d1d] = B.GetShape();
       const int vdim = input.vdim;
       const auto field = Reshape(&field_e[0], d1d, d1d, d1d, vdim);
@@ -39,11 +57,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
 
       for (int vd = 0; vd < vdim; vd++)
       {
-         MFEM_FOREACH_THREAD(dz, z, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
          {
-            MFEM_FOREACH_THREAD(dy, y, d1d)
+            MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
             {
-               MFEM_FOREACH_THREAD(qx, x, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
                {
                   real_t acc = 0.0;
                   for (int dx = 0; dx < d1d; dx++)
@@ -56,11 +74,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(dz, z, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
          {
-            MFEM_FOREACH_THREAD(qx, x, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
             {
-               MFEM_FOREACH_THREAD(qy, y, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
                {
                   real_t acc = 0.0;
                   for (int dy = 0; dy < d1d; dy++)
@@ -73,11 +91,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(qz, z, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
          {
-            MFEM_FOREACH_THREAD(qy, y, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
-               MFEM_FOREACH_THREAD(qx, x, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
                {
                   real_t acc = 0.0;
                   for (int dz = 0; dz < d1d; dz++)
@@ -94,10 +112,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
    else if constexpr (
       is_gradient_fop<std::decay_t<field_operator_t>>::value)
    {
-      const auto [q1d, unused, d1d] = B.GetShape();
+      // dbg("Gradient");
+      const auto [q1d, B_dim, d1d] = B.GetShape();
       const int vdim = input.vdim;
       const int dim = input.dim;
-      const auto field = Reshape(&field_e[0], d1d, d1d, d1d, vdim);
+      const auto field = Reshape(&std::as_const(field_e[0]), d1d, d1d, d1d, vdim);
       auto fqp = Reshape(&field_qp[0], vdim, dim, q1d, q1d, q1d);
 
       auto s0 = Reshape(&scratch_mem[0](0), d1d, d1d, q1d);
@@ -106,18 +125,41 @@ void map_field_to_quadrature_data_tensor_product_3d(
       auto s3 = Reshape(&scratch_mem[3](0), d1d, q1d, q1d);
       auto s4 = Reshape(&scratch_mem[4](0), d1d, q1d, q1d);
 
-      for (int vd = 0; vd < vdim; vd++)
+      // constexpr int MQ1 = T_Q1D > 0 ? T_Q1D : 8;
+      // static constexpr int DIM = 3;
+      // MFEM_VERIFY(q1d <= MQ1, "q1d > MQ1");
+      // MFEM_SHARED real_t smem[MQ1][MQ1];
+
+      // kernels::internal::d_regs3d_t<DIM, MQ1> r0, r1;
+      // real_t sB[MQ1][MQ1], sG[MQ1][MQ1];
+
+      /*
       {
-         MFEM_FOREACH_THREAD(dz, z, d1d)
+         assert(B_dim == 1 && "1D B required!");
+         kernels::internal::LoadMatrix(d1d, q1d, B, sB);
+         kernels::internal::LoadMatrix(d1d, q1d, G, sG);
+         for (int qx = 0; qx < q1d; qx++)
          {
-            MFEM_FOREACH_THREAD(dy, y, d1d)
+            for (int dx = 0; dx < d1d; dx++)
             {
-               MFEM_FOREACH_THREAD(qx, x, q1d)
+               assert(AlmostEq(B(qx, 0, dx), sB[dx][qx]));
+               assert(AlmostEq(G(qx, 0, dx), sG[dx][qx]));
+            }
+         }
+      }*/
+
+      for (int c = 0; c < vdim; c++)
+      {
+         MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
+         {
+            MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
+            {
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
                {
                   real_t uv[2] = {0.0, 0.0};
                   for (int dx = 0; dx < d1d; dx++)
                   {
-                     const real_t f = field(dx, dy, dz, vd);
+                     const real_t f = field(dx, dy, dz, c);
                      uv[0] += f * B(qx, 0, dx);
                      uv[1] += f * G(qx, 0, dx);
                   }
@@ -128,11 +170,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(dz, z, d1d)
+         MFEM_FOREACH_THREAD_DIRECT(dz, z, d1d)
          {
-            MFEM_FOREACH_THREAD(qy, y, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
-               MFEM_FOREACH_THREAD(qx, x, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
                {
                   real_t uvw[3] = {0.0, 0.0, 0.0};
                   for (int dy = 0; dy < d1d; dy++)
@@ -150,11 +192,11 @@ void map_field_to_quadrature_data_tensor_product_3d(
          }
          MFEM_SYNC_THREAD;
 
-         MFEM_FOREACH_THREAD(qz, z, q1d)
+         MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
          {
-            MFEM_FOREACH_THREAD(qy, y, q1d)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
-               MFEM_FOREACH_THREAD(qx, x, q1d)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, q1d)
                {
                   real_t uvw[3] = {0.0, 0.0, 0.0};
                   for (int dz = 0; dz < d1d; dz++)
@@ -163,19 +205,59 @@ void map_field_to_quadrature_data_tensor_product_3d(
                      uvw[1] += s3(dz, qy, qx) * B(qz, 0, dz);
                      uvw[2] += s4(dz, qy, qx) * G(qz, 0, dz);
                   }
-                  fqp(vd, 0, qx, qy, qz) = uvw[0];
-                  fqp(vd, 1, qx, qy, qz) = uvw[1];
-                  fqp(vd, 2, qx, qy, qz) = uvw[2];
+                  fqp(c, 0, qx, qy, qz) = uvw[0];
+                  fqp(c, 1, qx, qy, qz) = uvw[1];
+                  fqp(c, 2, qx, qy, qz) = uvw[2];
                }
             }
          }
          MFEM_SYNC_THREAD;
       }
+
+      /*
+      {
+         for (int c = 0; c < vdim; c++)
+         {
+            kernels::internal::LoadDofs3d(d1d, c, field, r0);
+            for (int d = 0; d < DIM; d++)
+            {
+               for (int dz = 0; dz < d1d; dz++)
+               {
+                  for (int dy = 0; dy < d1d; dy++)
+                  {
+                     for (int dx = 0; dx < d1d; dx++)
+                     {
+                        const real_t f = field(dx, dy, dz, c);
+                        assert(AlmostEq(f, r0[d][dz][dy][dx]));
+                     }
+                  }
+               }
+            }
+
+            kernels::internal::Grad3d(d1d, q1d, smem, sB, sG, r0, r1, c);
+            for (int qz = 0; qz < q1d; qz++)
+            {
+               for (int qy = 0; qy < q1d; qy++)
+               {
+                  for (int qx = 0; qx < q1d; qx++)
+                  {
+                     if (!AlmostEq(fqp(c, d, qx, qy, qz), r1[d][qz][qy][qx]))
+                     {
+                        dbg("\x1b[31m[{}:d] {} {}", c, fqp(c, d, qx, qy, qz), r1[d][qz][qy][qx]);
+                        dbg("❌❌❌"), std::exit(EXIT_FAILURE);
+                     }
+                  }
+               }
+            }
+         }
+         // dbg("✅✅✅✅✅✅✅✅✅✅✅✅✅✅✅");//, std::exit(EXIT_SUCCESS);
+      }*/
    }
    // TODO: Create separate function for clarity
    else if constexpr (
       std::is_same_v<std::decay_t<field_operator_t>, Weight>)
    {
+      // dbg("None");
       const int num_qp = integration_weights.GetShape()[0];
       // TODO: eeek
       const int q1d = (int)floor(std::pow(num_qp, 1.0/input.dim) + 0.5);
@@ -432,6 +514,9 @@ void map_fields_to_quadrature_data(
    const int &dimension,
    const bool &use_sum_factorization = false)
 {
+   // dbg();
+   assert(use_sum_factorization && "❌ use_sum_factorization required");
+
    // When the input_to_field map returns -1, this means the requested input
    // is the integration weight. Weights don't have a user defined field
    // attached to them and we create a dummy field which is not accessed
@@ -485,18 +570,19 @@ void map_field_to_quadrature_data_conditional(
    const int &dimension,
    const bool &use_sum_factorization = false)
 {
+   assert(false && "❌ condition not implemented");
    if (condition)
    {
       if (use_sum_factorization)
       {
          if (dimension == 2)
          {
-            map_field_to_quadrature_data_tensor_product_3d(
+            map_field_to_quadrature_data_tensor_product_2d(
                field_qp, dtqmap, field_e, fop, integration_weights, scratch_mem);
          }
          else if (dimension == 3)
          {
-            map_field_to_quadrature_data_tensor_product_2d(
+            map_field_to_quadrature_data_tensor_product_3d(
                field_qp, dtqmap, field_e, fop, integration_weights, scratch_mem);
          }
       }
@@ -520,6 +606,7 @@ void map_fields_to_quadrature_data_conditional(
    const std::array<bool, num_inputs> &conditions,
    const bool &use_sum_factorization = false)
 {
+   assert(false && "❌ condition not implemented");
    for_constexpr<num_inputs>([&](auto i)
    {
       map_field_to_quadrature_data_conditional(

fem/dfem/qfunction_transform.hpp

@@ -235,6 +235,8 @@ void process_qf_arg(
    }
 }
 
+// const tensor<real_t, DIM>       ∇u
+// const tensor<real_t, DIM, DIM>  D (PA_DATA)
 template <typename arg_type>
 MFEM_HOST_DEVICE inline
 void process_qf_arg(const DeviceTensor<2> &u, arg_type &arg, int qp)

fem/dfem/qp.hpp

@@ -0,0 +1,76 @@
+// 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
+
+#include "tuple.hpp"
+#include "../linalg/tensor.hpp"
+
+using namespace mfem::future;
+using mfem::future::tensor;
+
+// Helper to add dimension to tensor type
+template<typename T, int qp>
+struct AddQPDimension;
+
+// Specialization for tensor<real_t, dim>
+template<typename real_t, int dim, int qp>
+struct AddQPDimension<tensor<real_t, dim>, qp>
+{
+   using type = tensor<real_t, dim, qp>;
+};
+
+// Specialization for tensor<real_t, dim, dim>
+template<typename real_t, int dim, int qp>
+struct AddQPDimension<tensor<real_t, dim, dim>, qp>
+{
+   using type = tensor<real_t, dim, dim, qp>;
+};
+
+// Specialization for real_t (transforms to tensor<real_t, qp>)
+template<typename real_t, int qp>
+struct AddQPDimension
+{
+   using type = tensor<real_t, qp>;
+};
+
+// Helper to transform tuple
+template<typename Tuple, int qp>
+struct TransformTupleQP {};
+
+// Specialization for mfem::future::tuple
+template<int qp, typename... Types>
+struct TransformTupleQP<mfem::future::tuple<Types...>, qp>
+{
+   using type = mfem::future::tuple<typename AddQPDimension<Types, qp>::type...>;
+};
+
+template<int qp, typename... Types>
+struct TransformTupleQP<std::tuple<Types...>, qp>
+{
+   using type = std::tuple<typename AddQPDimension<Types, qp>::type...>;
+};
+
+// Function to transform tuple type with qp dimension
+template<int qp, typename qf_param_ts>
+struct add_qp_dimension
+{
+   using type = typename TransformTupleQP<qf_param_ts, qp>::type;
+};
+
+// Helper alias template for cleaner usage
+template<int qp, typename qf_param_ts>
+using add_qp_dimension_t = typename add_qp_dimension<qp, qf_param_ts>::type;
+
+// ...AddDomainIntegrator...
+// {
+//    constexpr int Q1D = 4;
+//    using qf_param_augmentd_ts = add_qp_dimension_t<Q1D, decay_tuple<qf_param_ts>>;
+// }

fem/dfem/util.hpp

@@ -34,10 +34,19 @@
 #include "parameterspace.hpp"
 #include "tuple.hpp"
 
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kSeashell
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
 namespace mfem::future
 {
 
 template<typename... Ts>
+inline MFEM_ALWAYS_INLINE
 constexpr auto to_array(const std::tuple<Ts...>& tuple)
 {
    constexpr auto get_array = [](const Ts&... x) { return std::array<typename std::common_type<Ts...>::type, sizeof...(Ts)> { x... }; };
@@ -48,6 +57,7 @@ namespace detail
 {
 
 template <typename lambda, std::size_t... i>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr(lambda&& f,
                              std::integral_constant<std::size_t, i>... Is)
 {
@@ -56,6 +66,7 @@ constexpr void for_constexpr(lambda&& f,
 
 
 template <std::size_t... n, typename lambda, typename... arg_types>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr(lambda&& f,
                              std::integer_sequence<std::size_t, n...>,
                              arg_types... args)
@@ -67,6 +78,7 @@ constexpr void for_constexpr(lambda&& f,
 }  // namespace detail
 
 template <typename lambda, std::size_t... i>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr(lambda&& f,
                              std::integer_sequence<std::size_t, i ... >)
 {
@@ -74,15 +86,18 @@ constexpr void for_constexpr(lambda&& f,
 }
 
 template <typename lambda>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr(lambda&& f, std::integer_sequence<std::size_t>) {}
 
 template <int... n, typename lambda>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr(lambda&& f)
 {
    detail::for_constexpr(f, std::make_integer_sequence<std::size_t, n> {}...);
 }
 
 template <typename lambda, typename arg_t>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg,
                                       std::integer_sequence<std::size_t>)
 {
@@ -90,6 +105,7 @@ constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg,
 }
 
 template <typename lambda, typename arg_t, std::size_t i, std::size_t... Is>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg,
                                       std::integer_sequence<std::size_t, i, Is...>)
 {
@@ -99,6 +115,7 @@ constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg,
 }
 
 template <typename lambda, typename arg_t>
+inline MFEM_ALWAYS_INLINE
 constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg)
 {
    using indices =
@@ -108,6 +125,7 @@ constexpr void for_constexpr_with_arg(lambda&& f, arg_t&& arg)
 }
 
 template <typename... input_ts, std::size_t... Is>
+inline MFEM_ALWAYS_INLINE
 auto make_dependency_map_impl(
    tuple<input_ts...> inputs,
    std::index_sequence<Is...>)
@@ -136,6 +154,7 @@ auto make_dependency_map_impl(
 // @returns an unordered_map where the keys are the field IDs and the values
 // are arrays of booleans indicating which inputs depend on each field ID.
 template <typename... input_ts>
+inline MFEM_ALWAYS_INLINE
 auto make_dependency_map(tuple<input_ts...> inputs)
 {
    return make_dependency_map_impl(inputs, std::index_sequence_for<input_ts...> {});
@@ -150,6 +169,7 @@ auto make_dependency_map(tuple<input_ts...> inputs)
 // ```
 // prints "int".
 template <typename T>
+inline MFEM_ALWAYS_INLINE
 constexpr auto get_type_name() -> std::string_view
 {
 #if defined(__clang__)
@@ -253,6 +273,7 @@ void pretty_print(const mfem::Vector& v)
 ///
 /// @param v vector of vectors to print
 template <typename T>
+inline
 void pretty_print(const mfem::Array<T>& v)
 {
    out << "[";
@@ -275,6 +296,7 @@ void pretty_print(const mfem::Array<T>& v)
 /// @tparam T type of array elements
 /// @tparam N size of array
 template<typename K, typename T, std::size_t N>
+inline
 void pretty_print(const std::unordered_map<K,std::array<T,N>>& map)
 {
    out << "{";
@@ -384,6 +406,7 @@ void pretty_print_mpi(const mfem::Vector& v)
 
 
 template <typename ... Ts>
+inline MFEM_ALWAYS_INLINE
 constexpr auto decay_types(tuple<Ts...> const &)
 -> tuple<std::remove_cv_t<std::remove_reference_t<Ts>>...>;
 
@@ -416,12 +439,14 @@ struct create_function_signature<output_t (*)(input_ts...)>
 };
 
 template <typename T>
+inline MFEM_ALWAYS_INLINE
 constexpr int GetFieldId()
 {
    return T::GetFieldId();
 }
 
 template <typename Tuple, std::size_t... Is>
+inline MFEM_ALWAYS_INLINE
 constexpr auto extract_field_ids_impl(Tuple&& t, std::index_sequence<Is...>)
 {
    return std::array<int, sizeof...(Is)>
@@ -435,6 +460,7 @@ constexpr auto extract_field_ids_impl(Tuple&& t, std::index_sequence<Is...>)
 /// @param t the tuple to extract field IDs from.
 /// @returns an array of field IDs.
 template <typename... Ts>
+inline MFEM_ALWAYS_INLINE
 constexpr auto extract_field_ids(const std::tuple<Ts...>& t)
 {
    return extract_field_ids_impl(t, std::index_sequence_for<Ts...> {});
@@ -446,6 +472,7 @@ constexpr auto extract_field_ids(const std::tuple<Ts...>& t)
 /// @param size the size of the array.
 /// @param value the value to search for.
 /// @returns true if the value is found, false otherwise.
+inline MFEM_ALWAYS_INLINE
 constexpr bool contains(const int* arr, std::size_t size, int value)
 {
    for (std::size_t i = 0; i < size; ++i)
@@ -463,6 +490,7 @@ constexpr bool contains(const int* arr, std::size_t size, int value)
 /// @param t the tuple to count unique field IDs from.
 /// @returns the number of unique field IDs.
 template <typename... Ts>
+inline MFEM_ALWAYS_INLINE
 constexpr std::size_t count_unique_field_ids(const std::tuple<Ts...>& t)
 {
    auto ids = extract_field_ids(t);
@@ -489,6 +517,7 @@ constexpr std::size_t count_unique_field_ids(const std::tuple<Ts...>& t)
 /// @param marker the marker std::array indicating which entries to get.
 /// @returns a std::vector containing the marked entries.
 template <typename T, std::size_t N>
+inline MFEM_ALWAYS_INLINE
 auto get_marked_entries(
    const std::array<T, N> &a,
    const std::array<bool, N> &marker)
@@ -509,6 +538,7 @@ auto get_marked_entries(
 /// @param t the tuple to filter fields from.
 /// @returns a tuple containing only the fields with field IDs not equal to -1.
 template <typename... Ts>
+inline MFEM_ALWAYS_INLINE
 constexpr auto filter_fields(const std::tuple<Ts...>& t)
 {
    return std::tuple_cat(
@@ -532,11 +562,13 @@ struct FieldDescriptor
    data_variant_t data;
 
    /// Default constructor
+   inline MFEM_ALWAYS_INLINE
    FieldDescriptor() :
       id(SIZE_MAX), data(data_variant_t{}) {}
 
    /// Constructor
    template <typename T>
+   inline MFEM_ALWAYS_INLINE
    FieldDescriptor(std::size_t field_id, const T* v) :
       id(field_id), data(v) {}
 };
@@ -568,11 +600,26 @@ struct ThreadBlocks
    int z = 1;
 };
 
-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+/*template <typename tag_t>
+struct ForallKernel
+{
+   template <typename func_t>
+   __global__ static void run(func_t f, int n)
+   {
+      int i = blockIdx.x;
+      extern __shared__ real_t shmem[];
+      if (i < n)
+      {
+         f(i, shmem);
+      }
+   }
+};*/
+
 template <typename func_t>
 __global__ void forall_kernel_shmem(func_t f, int n)
 {
-   int i = blockIdx.x;
+   int i = MFEM_BLOCK_ID(x);
    extern __shared__ real_t shmem[];
    if (i < n)
    {
@@ -581,27 +628,38 @@ __global__ void forall_kernel_shmem(func_t f, int n)
 }
 #endif
 
-template <typename func_t>
+template </*typename kernel_tag,*/ typename func_t>
+inline MFEM_ALWAYS_INLINE
 void forall(func_t f,
             const int &N,
             const ThreadBlocks &blocks,
             int num_shmem = 0,
             real_t *shmem = nullptr)
 {
+   // dbg("num_shmem: {}", num_shmem);
+   // dbg("blocks: {}x{}x{}", blocks.x, blocks.y, blocks.z);
+
    if (Device::Allows(Backend::CUDA_MASK) ||
        Device::Allows(Backend::HIP_MASK))
    {
 #if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-      // int gridsize = (N + Z - 1) / Z;
       int num_bytes = num_shmem * sizeof(decltype(shmem));
       dim3 block_size(blocks.x, blocks.y, blocks.z);
-      forall_kernel_shmem<<<N, block_size, num_bytes>>>(f, N);
+      // ForallKernel<kernel_tag>::run<<<N, block_size, num_bytes>>>(f, N);
+      if (num_bytes > 0)
+      {
+         forall_kernel_shmem<<<N, block_size, num_bytes>>>(f, N);
+      }
+      else
+      {
+         forall_kernel_shmem<<<N, block_size>>>(f, N);
+      }
 #if defined(MFEM_USE_CUDA)
       MFEM_GPU_CHECK(cudaGetLastError());
 #elif defined(MFEM_USE_HIP)
       MFEM_GPU_CHECK(hipGetLastError());
 #endif
-      MFEM_DEVICE_SYNC;
+      // MFEM_DEVICE_SYNC;
 #endif
    }
    else if (Device::Allows(Backend::CPU_MASK))
@@ -709,7 +767,7 @@ private:
 /// @param fields the vector of field descriptors.
 /// @returns the index of the field descriptor with the given ID,
 /// or SIZE_MAX if not found.
-inline
+inline MFEM_ALWAYS_INLINE
 std::size_t FindIdx(const std::size_t& id,
                     const std::vector<FieldDescriptor>& fields)
 {
@@ -727,7 +785,7 @@ std::size_t FindIdx(const std::size_t& id,
 ///
 /// @param f the field descriptor.
 /// @returns the vdof size of the field descriptor.
-inline
+inline MFEM_ALWAYS_INLINE
 int GetVSize(const FieldDescriptor &f)
 {
    return std::visit([](auto arg)
@@ -762,7 +820,7 @@ int GetVSize(const FieldDescriptor &f)
 /// @param f the field descriptor.
 /// @param el the element index.
 /// @param vdofs the array to store the element vdofs.
-inline
+inline MFEM_ALWAYS_INLINE
 void GetElementVDofs(const FieldDescriptor &f, int el, Array<int> &vdofs)
 {
    return std::visit([&](auto arg)
@@ -796,7 +854,7 @@ void GetElementVDofs(const FieldDescriptor &f, int el, Array<int> &vdofs)
 ///
 /// @param f the field descriptor.
 /// @returns the true dof size of the field descriptor.
-inline
+inline MFEM_ALWAYS_INLINE
 int GetTrueVSize(const FieldDescriptor &f)
 {
    return std::visit([](auto arg)
@@ -831,7 +889,7 @@ int GetTrueVSize(const FieldDescriptor &f)
 ///
 /// @param f the field descriptor.
 /// @returns the vdim of the field descriptor.
-inline
+inline MFEM_ALWAYS_INLINE
 int GetVDim(const FieldDescriptor &f)
 {
    return std::visit([](auto && arg)
@@ -863,6 +921,7 @@ int GetVDim(const FieldDescriptor &f)
 /// @tparam entity_t the entity type (see Entity).
 /// @returns the spatial dimension of the field descriptor.
 template <typename entity_t>
+inline MFEM_ALWAYS_INLINE
 int GetDimension(const FieldDescriptor &f)
 {
    return std::visit([](auto && arg)
@@ -897,7 +956,7 @@ int GetDimension(const FieldDescriptor &f)
 ///
 /// @param f the field descriptor.
 /// @returns the prolongation operator for the field descriptor.
-inline
+inline MFEM_ALWAYS_INLINE
 const Operator *get_prolongation(const FieldDescriptor &f)
 {
    return std::visit([](auto&& arg) -> const Operator*
@@ -926,7 +985,7 @@ const Operator *get_prolongation(const FieldDescriptor &f)
 /// @param o the element dof ordering.
 /// @returns the element restriction operator for the field descriptor in
 /// specified ordering.
-inline
+inline MFEM_ALWAYS_INLINE
 const Operator *get_element_restriction(const FieldDescriptor &f,
                                         ElementDofOrdering o)
 {
@@ -957,7 +1016,7 @@ const Operator *get_element_restriction(const FieldDescriptor &f,
 /// @returns the restriction operator for the field descriptor in
 /// specified ordering.
 template <typename entity_t>
-inline
+inline MFEM_ALWAYS_INLINE
 const Operator *get_restriction(const FieldDescriptor &f,
                                 const ElementDofOrdering &o)
 {
@@ -977,7 +1036,8 @@ const Operator *get_restriction(const FieldDescriptor &f,
 /// @returns a tuple containting a std::function with the transpose
 /// restriction callback and it's height.
 template <typename entity_t, typename fop_t>
-inline std::tuple<std::function<void(const Vector&, Vector&)>, int>
+inline MFEM_ALWAYS_INLINE
+std::tuple<std::function<void(const Vector&, Vector&)>, int>
 get_restriction_transpose(
    const FieldDescriptor &f,
    const ElementDofOrdering &o,
@@ -985,7 +1045,7 @@ get_restriction_transpose(
 {
    if constexpr (is_sum_fop<fop_t>::value)
    {
-      auto RT = [=](const Vector &v_e, Vector &v_l)
+      const auto RT = [=](const Vector &v_e, Vector &v_l)
       {
          v_l = v_e;
       };
@@ -994,7 +1054,8 @@ get_restriction_transpose(
    else
    {
       const Operator *R = get_restriction<entity_t>(f, o);
-      std::function<void(const Vector&, Vector&)> RT = [=](const Vector &x, Vector &y)
+      const std::function<void(const Vector&, Vector&)> RT =
+         [=](const Vector &x, Vector &y)
       {
          R->MultTranspose(x, y);
       };
@@ -1012,7 +1073,7 @@ get_restriction_transpose(
 /// @param field the field descriptor.
 /// @param x the input vector in tdofs.
 /// @param field_l the output vector in vdofs.
-inline
+inline MFEM_ALWAYS_INLINE
 void prolongation(const FieldDescriptor field, const Vector &x, Vector &field_l)
 {
    const auto P = get_prolongation(field);
@@ -1032,6 +1093,7 @@ void prolongation(const FieldDescriptor field, const Vector &x, Vector &field_l)
 /// @tparam N the number of fields.
 /// @tparam M the number of output fields.
 template <std::size_t N, std::size_t M>
+inline MFEM_ALWAYS_INLINE
 void prolongation(const std::array<FieldDescriptor, N> fields,
                   const Vector &x,
                   std::array<Vector, M> &fields_l)
@@ -1059,7 +1121,7 @@ void prolongation(const std::array<FieldDescriptor, N> fields,
 /// @param fields the array of field descriptors.
 /// @param x the input vector in tdofs.
 /// @param fields_l the array of output vectors in vdofs.
-inline
+inline MFEM_ALWAYS_INLINE
 void prolongation(const std::vector<FieldDescriptor> fields,
                   const Vector &x,
                   std::vector<Vector> &fields_l)
@@ -1086,7 +1148,7 @@ void prolongation(const std::vector<FieldDescriptor> fields,
 /// @param mpi_comm the MPI communicator.
 /// @tparam fop_t the field operator type.
 template <typename fop_t>
-inline
+inline MFEM_ALWAYS_INLINE
 std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
    const FieldDescriptor &f,
    const fop_t &fop,
@@ -1126,6 +1188,7 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
 /// @param ordering the element dof ordering.
 /// @tparam entity_t the entity type (see Entity).
 template <typename entity_t>
+inline MFEM_ALWAYS_INLINE
 void restriction(const FieldDescriptor u,
                  const Vector &u_l,
                  Vector &field_e,
@@ -1148,6 +1211,7 @@ void restriction(const FieldDescriptor u,
 /// @param offset the array index offset to start writing in fields_e.
 /// @tparam entity_t the entity type (see Entity).
 template <typename entity_t>
+inline MFEM_ALWAYS_INLINE
 void restriction(const std::vector<FieldDescriptor> u,
                  const std::vector<Vector> &u_l,
                  std::vector<Vector> &fields_e,
@@ -1167,6 +1231,7 @@ void restriction(const std::vector<FieldDescriptor> u,
 
 // TODO: keep this temporarily
 template <std::size_t N, std::size_t M>
+inline MFEM_ALWAYS_INLINE
 void element_restriction(const std::array<FieldDescriptor, N> u,
                          const std::array<Vector, N> &u_l,
                          std::array<Vector, M> &fields_e,
@@ -1190,6 +1255,7 @@ void element_restriction(const std::array<FieldDescriptor, N> u,
 /// @tparam entity_t the entity type (see Entity).
 /// @returns the number of entities of the given type.
 template <typename entity_t>
+inline MFEM_ALWAYS_INLINE
 int GetNumEntities(const mfem::Mesh &mesh)
 {
    if constexpr (std::is_same_v<entity_t, Entity::Element>)
@@ -1217,7 +1283,7 @@ int GetNumEntities(const mfem::Mesh &mesh)
 /// @param mode the mode of the DofToQuad object.
 /// @tparam entity_t the entity type (see Entity).
 template <typename entity_t>
-inline
+inline MFEM_ALWAYS_INLINE
 const DofToQuad *GetDofToQuad(const FieldDescriptor &f,
                               const IntegrationRule &ir,
                               DofToQuad::Mode mode)
@@ -1692,7 +1758,7 @@ void print_shared_memory_info(shmem_info_t &shmem_info)
 }
 
 template <std::size_t N>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 std::array<DofToQuadMap, N> load_dtq_mem(
    void *mem,
    int offset,
@@ -1755,7 +1821,7 @@ std::array<DofToQuadMap, N> load_dtq_mem(
 }
 
 template <std::size_t num_fields>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 std::array<DeviceTensor<1>, num_fields>
 load_field_mem(
    void *mem,
@@ -1789,12 +1855,12 @@ load_field_mem(
    return f;
 }
 
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 DeviceTensor<1> load_direction_mem(
    void *mem,
    int offset,
    const int &size,
-   const DeviceTensor<2> &direction,
+   const DeviceTensor<2, const real_t> &direction,
    const int &entity_idx)
 {
    int block_size = MFEM_THREAD_SIZE(x) *
@@ -1814,7 +1880,7 @@ DeviceTensor<1> load_direction_mem(
 }
 
 template <std::size_t N>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 std::array<DeviceTensor<2>, N> load_input_mem(
    void *mem,
    int offset,
@@ -1832,7 +1898,7 @@ std::array<DeviceTensor<2>, N> load_input_mem(
    return f;
 }
 
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 DeviceTensor<2> load_residual_mem(
    void *mem,
    int offset,
@@ -1844,7 +1910,7 @@ DeviceTensor<2> load_residual_mem(
 }
 
 template <std::size_t N>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 std::array<DeviceTensor<1>, 6> load_scratch_mem(
    void *mem,
    int offset,
@@ -1860,7 +1926,7 @@ std::array<DeviceTensor<1>, 6> load_scratch_mem(
 }
 
 template <typename shared_mem_info_t, std::size_t num_inputs, std::size_t num_outputs, std::size_t num_fields>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 auto unpack_shmem(
    void *shmem,
    const shared_mem_info_t &shmem_info,
@@ -1923,14 +1989,14 @@ auto unpack_shmem(
 }
 
 template <typename shared_mem_info_t, std::size_t num_inputs, std::size_t num_outputs, std::size_t num_fields>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 auto unpack_shmem(
    void *shmem,
    const shared_mem_info_t &shmem_info,
    const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
    const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
    const std::array<DeviceTensor<2>, num_fields> &wrapped_fields_e,
-   const DeviceTensor<2> &wrapped_direction_e,
+   const DeviceTensor<2, const real_t> &wrapped_direction_e,
    const int &num_qp,
    const int &e)
 {
@@ -2003,7 +2069,7 @@ auto unpack_shmem(
 }
 
 template <std::size_t... i>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 std::array<DeviceTensor<2>, sizeof...(i)> get_local_input_qp(
    const std::array<DeviceTensor<3>, sizeof...(i)> &input_qp_global, int e,
    std::index_sequence<i...>)
@@ -2018,7 +2084,7 @@ std::array<DeviceTensor<2>, sizeof...(i)> get_local_input_qp(
 }
 
 template <std::size_t N>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 void set_zero(std::array<DeviceTensor<2>, N> &v)
 {
    for (std::size_t i = 0; i < N; i++)
@@ -2033,7 +2099,7 @@ void set_zero(std::array<DeviceTensor<2>, N> &v)
 }
 
 template <std::size_t n>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 void set_zero(DeviceTensor<n> &u)
 {
    int s = 1;
@@ -2054,7 +2120,7 @@ void set_zero(DeviceTensor<n> &u)
 /// @param v destination DeviceTensor
 /// @tparam n DeviceTensor rank
 template <int n>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 void copy(DeviceTensor<n> &u, DeviceTensor<n> &v)
 {
    int s = 1;
@@ -2077,7 +2143,7 @@ void copy(DeviceTensor<n> &u, DeviceTensor<n> &v)
 /// @tparam n DeviceTensor rank
 /// @tparam m number of DeviceTensors
 template <int n, std::size_t m>
-MFEM_HOST_DEVICE inline
+MFEM_HOST_DEVICE inline MFEM_ALWAYS_INLINE
 void copy(std::array<DeviceTensor<n>, m> &u,
           std::array<DeviceTensor<n>, m> &v)
 {
@@ -2095,6 +2161,7 @@ void copy(std::array<DeviceTensor<n>, m> &u,
 /// @tparam num_fields number of fields
 /// @return array of field data wrapped in DeviceTensors
 template <std::size_t num_fields>
+inline MFEM_ALWAYS_INLINE
 std::array<DeviceTensor<2>, num_fields> wrap_fields(
    std::vector<Vector> &fields,
    std::array<int, num_fields> &field_sizes,
@@ -2131,6 +2198,7 @@ std::array<DeviceTensor<2>, num_fields> wrap_fields(
 /// size required on quadrature points using GetSizeOnQP() and adds it to the
 /// total. Non-dependent inputs contribute zero to the total size.
 template <typename input_t, std::size_t num_fields, std::size_t... i>
+inline MFEM_ALWAYS_INLINE
 int accumulate_sizes_on_qp(
    const input_t &inputs,
    std::array<bool, sizeof...(i)> &kinput_is_dependent,
@@ -2157,6 +2225,7 @@ template <
    typename field_operator_ts,
    std::size_t N = tuple_size<field_operator_ts>::value,
    std::size_t... Is>
+inline MFEM_ALWAYS_INLINE
 std::array<DofToQuadMap, N> create_dtq_maps_impl(
    field_operator_ts &fops,
    std::vector<const DofToQuad*> &dtqs,
@@ -2242,6 +2311,7 @@ template <
    typename entity_t,
    typename field_operator_ts,
    std::size_t num_fields>
+inline MFEM_ALWAYS_INLINE
 std::array<DofToQuadMap, num_fields> create_dtq_maps(
    field_operator_ts &fops,
    std::vector<const DofToQuad*> &dtqmaps,

fem/fem.hpp

@@ -51,7 +51,7 @@
 #include "hyperbolic.hpp"
 #include "bounds.hpp"
 
-#include "dfem/doperator.hpp"
+// #include "dfem/doperator.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/integ/bilininteg_diffusion_kernels.hpp

@@ -1214,20 +1214,23 @@ inline void SmemPADiffusionApply3D(const int NE,
 
 namespace
 {
-using ApplyKernelType = DiffusionIntegrator::ApplyKernelType;
-using DiagonalKernelType = DiffusionIntegrator::DiagonalKernelType;
+using DiffusionApplyKernelType =
+   DiffusionIntegrator::DiffusionApplyKernelType;
+
+using DiffusionDiagonalKernelType =
+   DiffusionIntegrator::DiffusionDiagonalKernelType;
 }
 
 template<int DIM, int T_D1D, int T_Q1D>
-ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Kernel()
+DiffusionApplyKernelType DiffusionIntegrator::DiffusionApplyPAKernel::Kernel()
 {
    if (DIM == 2) { return internal::SmemPADiffusionApply2D<T_D1D,T_Q1D>; }
    else if (DIM == 3) { return internal::SmemPADiffusionApply3D<T_D1D, T_Q1D>; }
    else { MFEM_ABORT(""); }
 }
 
-inline
-ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
+inline DiffusionApplyKernelType
+DiffusionIntegrator::DiffusionApplyPAKernel::Fallback(int DIM, int, int)
 {
    if (DIM == 2) { return internal::PADiffusionApply2D; }
    else if (DIM == 3) { return internal::PADiffusionApply3D; }
@@ -1235,15 +1238,16 @@ ApplyKernelType DiffusionIntegrator::ApplyPAKernels::Fallback(int DIM, int, int)
 }
 
 template<int DIM, int D1D, int Q1D>
-DiagonalKernelType DiffusionIntegrator::DiagonalPAKernels::Kernel()
+DiffusionDiagonalKernelType
+DiffusionIntegrator::DiffusionDiagonalPAKernel::Kernel()
 {
    if (DIM == 2) { return internal::SmemPADiffusionDiagonal2D<D1D,Q1D>; }
    else if (DIM == 3) { return internal::SmemPADiffusionDiagonal3D<D1D, Q1D>; }
    else { MFEM_ABORT(""); }
 }
 
-inline DiagonalKernelType
-DiffusionIntegrator::DiagonalPAKernels::Fallback(int DIM, int, int)
+inline DiffusionDiagonalKernelType
+DiffusionIntegrator::DiffusionDiagonalPAKernel::Fallback(int DIM, int, int)
 {
    if (DIM == 2) { return internal::PADiffusionDiagonal2D; }
    else if (DIM == 3) { return internal::PADiffusionDiagonal3D; }

fem/integ/bilininteg_diffusion_pa.cpp

@@ -16,6 +16,14 @@
 #include "../ceed/integrators/diffusion/diffusion.hpp"
 #include "bilininteg_diffusion_kernels.hpp"
 
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kNvidia
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
 namespace mfem
 {
 
@@ -31,8 +39,8 @@ void DiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
       const Array<real_t> &B = maps->B;
       const Array<real_t> &G = maps->G;
       const Vector &Dv = pa_data;
-      DiagonalPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Dv,
-                             diag, dofs1D, quad1D);
+      DiffusionDiagonalPAKernel::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Dv,
+                                     diag, dofs1D, quad1D);
    }
 }
 
@@ -67,9 +75,9 @@ void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
          MFEM_ABORT("OCCA PADiffusionApply unknown kernel!");
       }
 #endif // MFEM_USE_OCCA
-
-      ApplyPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Bt,
-                          Gt, Dv, x, y, dofs1D, quad1D);
+      db1("[DiffusionApplyPAKernel] D1D:{} Q1D:{}", dofs1D, quad1D);
+      DiffusionApplyPAKernel::Run(dim, dofs1D, quad1D, ne, symmetric, B, G, Bt,
+                                  Gt, Dv, x, y, dofs1D, quad1D);
    }
 }
 
@@ -174,9 +182,9 @@ void DiffusionIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
    abs_pa_data.Abs();
    auto abs_maps = maps->Abs();
 
-   ApplyPAKernels::Run(dim, dofs1D, quad1D, ne, symmetric,
-                       abs_maps.B, abs_maps.G, abs_maps.Bt, abs_maps.Gt,
-                       abs_pa_data, x, y, dofs1D, quad1D);
+   DiffusionApplyPAKernel::Run(dim, dofs1D, quad1D, ne, symmetric,
+                               abs_maps.B, abs_maps.G, abs_maps.Bt, abs_maps.Gt,
+                               abs_pa_data, x, y, dofs1D, quad1D);
 }
 
 void DiffusionIntegrator::AddAbsMultTransposePA(const Vector &x,

fem/kernel_dispatch.hpp

@@ -13,6 +13,7 @@
 #define MFEM_KERNEL_DISPATCH_HPP
 
 #include "../config/config.hpp"
+#include "../config/tconfig.hpp" // MFEM_ALWAYS_INLINE
 #include "kernel_reporter.hpp"
 #include <unordered_map>
 #include <tuple>
@@ -95,11 +96,13 @@ struct KernelDispatchKeyHash
 {
 private:
    template<int N>
+   inline MFEM_ALWAYS_INLINE
    size_t operator()(std::tuple<KernelParameters...> value) const { return 0; }
 
    // The hashing formula here is taken directly from the Boost library, with
    // the magic number 0x9e3779b9 chosen to minimize hashing collisions.
    template<std::size_t N, typename THead, typename... TTail>
+   inline MFEM_ALWAYS_INLINE
    size_t operator()(std::tuple<KernelParameters...> value) const
    {
       constexpr int Index = N - sizeof...(TTail) - 1;
@@ -109,6 +112,7 @@ private:
    }
 public:
    /// Returns the hash of the given @a value.
+   inline MFEM_ALWAYS_INLINE
    size_t operator()(std::tuple<KernelParameters...> value) const
    {
       return operator()<sizeof...(KernelParameters),KernelParameters...>(value);
@@ -137,7 +141,8 @@ class KernelDispatchTable<Kernels,
    /// Only valid when the function @a f is not a member function.
    template <typename F, typename... Args,
              typename std::enable_if<std::is_pointer<F>::value,bool>::type=true>
-   static void Invoke(F f, Args&&... args)
+   static inline MFEM_ALWAYS_INLINE
+   void Invoke(F f, Args&&... args)
    {
       f(std::forward<Args>(args)...);
    }
@@ -149,7 +154,8 @@ class KernelDispatchTable<Kernels,
    template <typename F, typename T, typename... Args,
              typename std::enable_if<
                 std::is_member_function_pointer<F>::value,bool>::type=true>
-   static void Invoke(F f, T&& t, Args&&... args)
+   static inline MFEM_ALWAYS_INLINE
+   void Invoke(F f, T&& t, Args&&... args)
    {
       (t.*f)(std::forward<Args>(args)...);
    }
@@ -164,7 +170,8 @@ public:
    /// If the kernel is a member function, then the first argument after @a
    /// params should be the object on which it is called.
    template<typename... Args>
-   static void Run(Params... params, Args&&... args)
+   static inline MFEM_ALWAYS_INLINE
+   void Run(Params... params, Args&&... args)
    {
       const auto &table = Kernels::Get().table;
       const std::tuple<Params...> key = std::make_tuple(params...);

fem/picojson.h

@@ -28,6 +28,8 @@
 #ifndef picojson_h
 #define picojson_h
 
+#include "../config/config.hpp"
+
 #include <algorithm>
 #include <cstdio>
 #include <cstdlib>
@@ -73,7 +75,7 @@ extern "C" {
 #endif
 
 #ifndef PICOJSON_ASSERT
-# define PICOJSON_ASSERT(e) do { if (! (e)) throw std::runtime_error(#e); } while (0)
+# define PICOJSON_ASSERT(e) do { if (! (e)) MFEM_THROW(std::runtime_error, #e); } while (0)
 #endif
 
 #ifdef _MSC_VER
@@ -202,7 +204,7 @@ namespace picojson {
         isnan(n) || isinf(n)
 #endif
         ) {
-      throw std::overflow_error("");
+      MFEM_THROW(std::runtime_error, "overflow_error");
     }
     u_.number_ = n;
   }
@@ -300,11 +302,11 @@ namespace picojson {
 #else
 #define GET(ctype, var)						\
   template <> inline const ctype& value::get<ctype>() const {	\
-    do { if (! (is<ctype>())) throw std::runtime_error("type mismatch! call is<type>() before get<type>()"); } while (0); \
+    do { if (! (is<ctype>())) MFEM_THROW(std::runtime_error, "type mismatch! call is<type>() before get<type>()"); } while (0); \
     return var;							\
   }								\
   template <> inline ctype& value::get<ctype>() {		\
-    do { if (! (is<ctype>())) throw std::runtime_error("type mismatch! call is<type>() before get<type>()"); } while (0); \
+    do { if (! (is<ctype>())) MFEM_THROW(std::runtime_error, "type mismatch! call is<type>() before get<type>()"); } while (0); \
     return var;							\
   }
 #endif

general/backends.hpp

@@ -47,6 +47,7 @@
 #define MFEM_DEVICE
 #define MFEM_HOST
 #define MFEM_LAMBDA
+#define MFEM_CONSTANT
 // #define MFEM_HOST_DEVICE // defined in config/config.hpp
 // MFEM_DEVICE_SYNC is made available for debugging purposes
 #define MFEM_DEVICE_SYNC

general/cuda.hpp

@@ -20,9 +20,10 @@
 
 #ifdef MFEM_USE_CUDA
 #define MFEM_USE_CUDA_OR_HIP
-#define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
 #define MFEM_LAMBDA __host__
+#define MFEM_DEVICE __device__
+#define MFEM_CONSTANT __constant__
 // #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(cudaDeviceSynchronize())
 #define MFEM_STREAM_SYNC MFEM_GPU_CHECK(cudaStreamSynchronize(0))

general/gecko.cpp

@@ -9,6 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
+#include "../config/config.hpp"
 #include "gecko.hpp"
 
 // This file collects the sources of the Gecko library as a single module.
@@ -392,7 +393,7 @@ Subgraph::Subgraph(Graph* g, uint n) : g(g), n(n), f(g->functional)
 {
    if (n > GECKO_WINDOW_MAX)
    {
-      throw std::out_of_range("optimization window too large");
+      MFEM_THROW(std::out_of_range, "optimization window too large");
    }
    cache = new Subnode[n << n];
 }
@@ -758,7 +759,8 @@ Graph::arc_source(Arc::Index a) const
       }
    }
    // should never get here
-   throw std::runtime_error("internal data structure corrupted");
+   MFEM_THROW(::std::runtime_error, "internal data structure corrupted");
+   return Node::null;
 }
 
 // Return reverse arc (j, i) of arc a = (i, j).

general/hip.cpp

@@ -9,6 +9,7 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
+#include <cstring>
 #include "backends.hpp"
 #include "globals.hpp"
 
@@ -175,6 +176,15 @@ void* HipMemcpyDtoHAsync(void *dst, const void *src, size_t bytes)
    return dst;
 }
 
+const void* HipMemcpyToSymbol(const void *d_sym, const void *h_src,
+                              size_t bytes)
+{
+#ifdef MFEM_USE_HIP
+   MFEM_GPU_CHECK(hipMemcpyToSymbol(d_sym, h_src, bytes));
+#endif
+   return memcpy(const_cast<void*>(d_sym), h_src, bytes);
+}
+
 void HipCheckLastError()
 {
 #ifdef MFEM_USE_HIP

general/hip.hpp

@@ -20,8 +20,9 @@
 
 #ifdef MFEM_USE_HIP
 #define MFEM_USE_CUDA_OR_HIP
-#define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
+#define MFEM_DEVICE __device__
+#define MFEM_CONSTANT __constant__
 #define MFEM_LAMBDA __host__ __device__
 // #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(hipDeviceSynchronize())
@@ -94,6 +95,10 @@ void* HipMemcpyDtoH(void *h_dst, const void *d_src, size_t bytes);
 /// Copies memory from Device to Host
 void* HipMemcpyDtoHAsync(void *h_dst, const void *d_src, size_t bytes);
 
+/// Copies data to the given symbol on the device.
+const void* HipMemcpyToSymbol(const void *d_sym, const void *h_src,
+                              size_t bytes);
+
 /// Check the error code returned by hipGetLastError(), aborting on error.
 void HipCheckLastError();
 

general/mem_manager.cpp

@@ -249,7 +249,7 @@ class Aligned32HostMemorySpace : public HostMemorySpace
 public:
    Aligned32HostMemorySpace(): HostMemorySpace() { }
    void Alloc(void **ptr, size_t bytes) override
-   { if (mfem_memalign(ptr, 32, bytes) != 0) { throw ::std::bad_alloc(); } }
+   { if (mfem_memalign(ptr, 32, bytes) != 0) { MFEM_THROW(std::bad_alloc,); } }
    void Dealloc(void *ptr) override { mfem_aligned_free(ptr); }
 };
 
@@ -259,7 +259,12 @@ class Aligned64HostMemorySpace : public HostMemorySpace
 public:
    Aligned64HostMemorySpace(): HostMemorySpace() { }
    void Alloc(void **ptr, size_t bytes) override
-   { if (mfem_memalign(ptr, 64, bytes) != 0) { throw ::std::bad_alloc(); } }
+   {
+      if (mfem_memalign(ptr, 64, bytes) != 0)
+      {
+         MFEM_THROW(::std::bad_alloc,);
+      }
+   }
    void Dealloc(void *ptr) override { mfem_aligned_free(ptr); }
 };
 
@@ -365,7 +370,7 @@ inline void MmuAlloc(void **ptr, const size_t bytes)
    const int prot = PROT_READ | PROT_WRITE;
    const int flags = MAP_ANONYMOUS | MAP_PRIVATE;
    *ptr = ::mmap(NULL, length, prot, flags, -1, 0);
-   if (*ptr == MAP_FAILED) { throw ::std::bad_alloc(); }
+   if (*ptr == MAP_FAILED) { MFEM_THROW(::std::bad_alloc,); }
 }
 
 /// MMU deallocation, through ::munmap

general/nvtx.hpp

@@ -0,0 +1 @@
+../../stash/debug/nvtx.hpp

general/zstr.hpp

@@ -75,7 +75,7 @@ inline char* check_strerror_r(char* r, char*, int)
 
 /// Overload of error-reporting function, to enable use with VS.
 /// Ref: http://stackoverflow.com/a/901316/717706
-static std::string strerror()
+[[maybe_unused]] static std::string strerror()
 {
    std::string buff(80, '\0');
 #ifdef _WIN32
@@ -148,18 +148,21 @@ struct static_method_holder
    {
       if ((mode & std::ios_base::trunc) && ! (mode & std::ios_base::out))
       {
-         throw Exception(std::string("strict_fstream: open('") + filename +
-                         "'): mode error: trunc and not out");
+         MFEM_THROW(Exception,
+                    std::string("strict_fstream: open('")
+                    + filename + "'): mode error: trunc and not out");
       }
       else if ((mode & std::ios_base::app) && ! (mode & std::ios_base::out))
       {
-         throw Exception(std::string("strict_fstream: open('") + filename +
-                         "'): mode error: app and not out");
+         MFEM_THROW(Exception,
+                    std::string("strict_fstream: open('") + filename +
+                    "'): mode error: app and not out");
       }
       else if ((mode & std::ios_base::trunc) && (mode & std::ios_base::app))
       {
-         throw Exception(std::string("strict_fstream: open('") + filename +
-                         "'): mode error: trunc and app");
+         MFEM_THROW(Exception,
+                    std::string("strict_fstream: open('") + filename +
+                    "'): mode error: trunc and app");
       }
    }
    static void check_open(std::ios * s_p, const std::string& filename,
@@ -167,26 +170,32 @@ struct static_method_holder
    {
       if (s_p->fail())
       {
-         throw Exception(std::string("strict_fstream: open('")
-                         + filename + "'," + mode_to_string(mode) + "): open failed: "
-                         + strerror());
+         MFEM_THROW(Exception,
+                    std::string("strict_fstream: open('") + filename +
+                    "'," + mode_to_string(mode) + "): open failed: " +
+                    strerror());
       }
    }
    static void check_peek(std::istream * is_p, const std::string& filename,
                           std::ios_base::openmode mode)
    {
       bool peek_failed = true;
+#ifdef MFEM_HAS_EXCEPTIONS
       try
+#endif
       {
          is_p->peek();
          peek_failed = is_p->fail();
       }
+#ifdef MFEM_HAS_EXCEPTIONS
       catch (std::ios_base::failure&) {}
+#endif
       if (peek_failed)
       {
-         throw Exception(std::string("strict_fstream: open('")
-                         + filename + "'," + mode_to_string(mode) + "): peek failed: "
-                         + strerror());
+         MFEM_THROW(Exception,
+                    std::string("strict_fstream: open('") + filename +
+                    "'," + mode_to_string(mode) + "): peek failed: " +
+                    strerror());
       }
       is_p->clear();
    }

linalg/dtensor.hpp

@@ -13,6 +13,7 @@
 #define MFEM_DTENSOR
 
 #include "../general/backends.hpp"
+#include "../config/tconfig.hpp" // MFEM_ALWAYS_INLINE
 #include <array>
 
 namespace mfem
@@ -23,8 +24,8 @@ template <int N, int Dim, typename T, typename... Args>
 class TensorInd
 {
 public:
-   MFEM_HOST_DEVICE
-   static inline int result(const int* sizes, T first, Args... args)
+   static inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   int result(const int* sizes, T first, Args... args)
    {
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<sizes[N-1],"Trying to access out of boundary.");
@@ -39,8 +40,8 @@ template <int Dim, typename T, typename... Args>
 class TensorInd<Dim, Dim, T, Args...>
 {
 public:
-   MFEM_HOST_DEVICE
-   static inline int result(const int* sizes, T first, Args... args)
+   static inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   int result(const int* sizes, T first, Args... args)
    {
 #if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
       MFEM_ASSERT(first<static_cast<T>(sizes[Dim-1]),
@@ -56,8 +57,8 @@ template <int N, int Dim, typename T, typename... Args>
 class Init
 {
 public:
-   MFEM_HOST_DEVICE
-   static inline int result(int* sizes, T first, Args... args)
+   static inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   int result(int* sizes, T first, Args... args)
    {
       sizes[N - 1] = first;
       return first * Init < N + 1, Dim, Args... >::result(sizes, args...);
@@ -69,8 +70,8 @@ template <int Dim, typename T, typename... Args>
 class Init<Dim, Dim, T, Args...>
 {
 public:
-   MFEM_HOST_DEVICE
-   static inline int result(int* sizes, T first, Args... args)
+   static inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   int result(int* sizes, T first, Args... args)
    {
       sizes[Dim - 1] = first;
       return first;
@@ -90,11 +91,12 @@ protected:
 public:
    /// Default constructor
    // DeviceTensor() = delete;
-   MFEM_HOST_DEVICE
-   DeviceTensor() {}
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   DeviceTensor() = default;
 
    /// Constructor to initialize a tensor from the Scalar array data_
-   template <typename... Args> MFEM_HOST_DEVICE
+   template <typename... Args>
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
    DeviceTensor(Scalar* data_, Args... args)
    {
       static_assert(sizeof...(args) == Dim, "Wrong number of arguments");
@@ -105,16 +107,19 @@ public:
    }
 
    /// Copy constructor (default)
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
    DeviceTensor(const DeviceTensor&) = default;
 
    /// Copy assignment (default)
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
    DeviceTensor& operator=(const DeviceTensor&) = default;
 
    /// Conversion to `Scalar *`.
-   MFEM_HOST_DEVICE inline operator Scalar *() const { return data; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE operator Scalar *() const { return data; }
 
    /// Const accessor for the data
-   template <typename... Args> MFEM_HOST_DEVICE inline
+   template <typename... Args>
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
    Scalar& operator()(Args... args) const
    {
       static_assert(sizeof...(args) == Dim, "Wrong number of arguments");
@@ -122,36 +127,37 @@ public:
    }
 
    /// Subscript operator where the tensor is viewed as a 1D array.
-   MFEM_HOST_DEVICE inline Scalar& operator[](int i) const
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE Scalar& operator[](int i) const
    {
       return data[i];
    }
 
    /// Returns the shape of the tensor.
-   MFEM_HOST_DEVICE inline auto &GetShape() const { return sizes; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto &GetShape() const { return sizes; }
 };
 
 
 /** @brief Wrap a pointer as a DeviceTensor with automatically deduced template
     parameters */
-template <typename T, typename... Dims> MFEM_HOST_DEVICE
-inline DeviceTensor<sizeof...(Dims),T> Reshape(T *ptr, Dims... dims)
+template <typename T, typename... Dims>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+DeviceTensor<sizeof...(Dims),T> Reshape(T *ptr, Dims... dims)
 {
    return DeviceTensor<sizeof...(Dims),T>(ptr, dims...);
 }
 
 
-typedef DeviceTensor<1,int> DeviceArray;
-typedef DeviceTensor<1,const int> ConstDeviceArray;
+using DeviceArray = DeviceTensor<1,int>;
+using ConstDeviceArray = DeviceTensor<1,const int>;
 
-typedef DeviceTensor<1,real_t> DeviceVector;
-typedef DeviceTensor<1,const real_t> ConstDeviceVector;
+using DeviceVector = DeviceTensor<1,real_t>;
+using ConstDeviceVector = DeviceTensor<1,const real_t>;
 
-typedef DeviceTensor<2,real_t> DeviceMatrix;
-typedef DeviceTensor<2,const real_t> ConstDeviceMatrix;
+using DeviceMatrix = DeviceTensor<2,real_t>;
+using ConstDeviceMatrix = DeviceTensor<2,const real_t>;
 
-typedef DeviceTensor<3,real_t> DeviceCube;
-typedef DeviceTensor<3,const real_t> ConstDeviceCube;
+using DeviceCube = DeviceTensor<3,real_t>;
+using ConstDeviceCube = DeviceTensor<3,const real_t>;
 
 } // mfem namespace
 

linalg/tensor.hpp

@@ -17,6 +17,7 @@
 
 #pragma once
 
+#include "../config/tconfig.hpp" // MFEM_ALWAYS_INLINE
 #include "../general/backends.hpp"
 #include "dual.hpp"
 #include <limits>
@@ -39,11 +40,14 @@ struct tensor<T>
    using type = T;
    static constexpr int ndim      = 1;
    static constexpr int first_dim = 0;
-   MFEM_HOST_DEVICE T& operator[](int /*unused*/) { return values; }
-   MFEM_HOST_DEVICE const T& operator[](int /*unused*/) const { return values; }
-   MFEM_HOST_DEVICE T& operator()(int /*unused*/) { return values; }
-   MFEM_HOST_DEVICE const T& operator()(int /*unused*/) const { return values; }
-   MFEM_HOST_DEVICE operator T() const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+   T& operator[](int /*unused*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator[](
+      int /*unused*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int /*unused*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(
+      int /*unused*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE operator T() const { return values; }
    T values;
 };
 
@@ -53,85 +57,205 @@ struct tensor<T, n0>
    using type = T;
    static constexpr int ndim      = 1;
    static constexpr int first_dim = n0;
-   MFEM_HOST_DEVICE T& operator[](int i) { return values[i]; }
-   MFEM_HOST_DEVICE const T& operator[](int i) const { return values[i]; }
-   MFEM_HOST_DEVICE T& operator()(int i) { return values[i]; }
-   MFEM_HOST_DEVICE const T& operator()(int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator[](int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator[](int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int i) const { return values[i]; }
    T values[n0];
 };
 
+template < typename T >
+struct tensor<T, 0>
+{
+   using type = T;
+   static constexpr int ndim      = 1;
+   static constexpr int first_dim = 0;
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator[](int /**/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator[](int /**/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int /**/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int /**/) const { return values; }
+   T values;
+};
+
 template < typename T, int n0, int n1 >
 struct tensor<T, n0, n1>
 {
    using type = T;
    static constexpr int ndim      = 2;
    static constexpr int first_dim = n0;
-   MFEM_HOST_DEVICE tensor< T, n1 >& operator[](int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1 >& operator[](int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n1 >& operator()(int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1 >& operator()(int i) const { return values[i]; }
-   MFEM_HOST_DEVICE T& operator()(int i, int j) { return values[i][j]; }
-   MFEM_HOST_DEVICE const T& operator()(int i, int j) const { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1 >& operator[](int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1 >& operator[](
+      int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1 >& operator()(int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1 >& operator()(
+      int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int i, int j) { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int i,
+                                                                  int j) const { return values[i][j]; }
    tensor < T, n1 > values[n0];
 };
 
+template < typename T, int n1 >
+struct tensor<T, 0, n1>
+{
+   using type = T;
+   static constexpr int ndim      = 2;
+   static constexpr int first_dim = 0;
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1 >& operator[](
+      int /**/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1 >& operator[](
+      int /**/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1 >& operator()(
+      int /**/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1 >& operator()(
+      int /**/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int /**/, int j) { return values[j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int /**/,
+                                                                  int j) const { return values[j]; }
+   tensor < T, n1 > values;
+};
+
 template < typename T, int n0, int n1, int n2 >
 struct tensor<T, n0, n1, n2>
 {
    using type = T;
    static constexpr int ndim      = 3;
    static constexpr int first_dim = n0;
-   MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator[](int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2 >& operator[](int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator()(int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2 >& operator()(int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n2 >& operator()(int i, int j) { return values[i][j]; }
-   MFEM_HOST_DEVICE const tensor< T, n2 >& operator()(int i, int j) const { return values[i][j]; }
-   MFEM_HOST_DEVICE T& operator()(int i, int j, int k) { return values[i][j][k]; }
-   MFEM_HOST_DEVICE const T& operator()(int i, int j, int k) const { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator[](
+      int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2 >&
+   operator[](int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator()(
+      int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2 >& operator()
+   (int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n2 >& operator()(int i,
+                                                                          int j) { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n2 >& operator()(
+      int i, int j) const { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int i, int j, int k) { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int i, int j,
+                                                                  int k) const { return values[i][j][k]; }
    tensor < T, n1, n2 > values[n0];
 };
 
+template < typename T, int n1, int n2 >
+struct tensor<T, 0, n1, n2>
+{
+   using type = T;
+   static constexpr int ndim      = 3;
+   static constexpr int first_dim = 0;
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator[](
+      int /*i*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2 >&
+   operator[](int /*i*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2 >& operator()(
+      int /*i*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2 >& operator()
+   (int /*i*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n2 >& operator()(
+      int /*i*/, int j) { return values[j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n2 >& operator()(
+      int i, int j) const { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int /*i*/, int j,
+                                                            int k) { return values[j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int /*i*/, int j,
+                                                                  int k) const { return values[j][k]; }
+   tensor < T, n1, n2 > values;
+};
+
 template < typename T, int n0, int n1, int n2, int n3 >
 struct tensor<T, n0, n1, n2, n3>
 {
    using type = T;
    static constexpr int ndim      = 4;
    static constexpr int first_dim = n0;
-   MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator[](int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >& operator[](int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator()(int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >& operator()(int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n2, n3 >& operator()(int i, int j) { return values[i][j]; }
-   MFEM_HOST_DEVICE const tensor< T, n2, n3 >& operator()(int i, int j) const { return values[i][j]; }
-   MFEM_HOST_DEVICE tensor< T, n3 >& operator()(int i, int j, int k) { return values[i][j][k]; }
-   MFEM_HOST_DEVICE const tensor< T, n3 >& operator()(int i, int j, int k) const { return values[i][j][k]; }
-   MFEM_HOST_DEVICE T& operator()(int i, int j, int k, int l) { return values[i][j][k][l]; }
-   MFEM_HOST_DEVICE const T&  operator()(int i, int j, int k, int l) const { return values[i][j][k][l]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator[](
+      int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >&
+   operator[](int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator()(
+      int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >&
+   operator()(int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n2, n3 >& operator()(
+      int i, int j) { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n2, n3 >& operator()
+   (int i, int j) const { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n3 >& operator()(int i,
+                                                                          int j, int k) { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n3 >& operator()(
+      int i, int j, int k) const { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int i, int j, int k,
+                                                            int l) { return values[i][j][k][l]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T&  operator()(int i, int j,
+                                                                   int k, int l) const { return values[i][j][k][l]; }
    tensor < T, n1, n2, n3 > values[n0];
 };
 
+template < typename T, int n1, int n2, int n3 >
+struct tensor<T, 0, n1, n2, n3>
+{
+   using type = T;
+   static constexpr int ndim      = 4;
+   static constexpr int first_dim = 0;
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator[](
+      int /*i*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >&
+   operator[](int /*i*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3 >& operator()(
+      int /*i*/) { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3 >&
+   operator()(int /*i*/) const { return values; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n2, n3 >& operator()(
+      int /*i*/, int j) { return values[j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n2, n3 >& operator()
+   (int /*i*/, int j) const { return values[j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n3 >& operator()(
+      int /*i*/, int j, int k) { return values[j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n3 >& operator()(
+      int /*i*/, int j,
+      int k) const { return values[j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int /*i*/, int j,
+                                                            int k, int l) { return values[j][k][l]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T&  operator()(int /*i*/,
+                                                                   int j, int k, int l) const { return values[j][k][l]; }
+   tensor < T, n1, n2, n3 > values;
+};
+
 template < typename T, int n0, int n1, int n2, int n3, int n4 >
 struct tensor<T, n0, n1, n2, n3, n4>
 {
    using type = T;
    static constexpr int ndim      = 5;
    static constexpr int first_dim = n0;
-   MFEM_HOST_DEVICE tensor< T, n1, n2, n3, n4 >& operator[](int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2, n3, n4 >& operator[](int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n1, n2, n3, n4 >& operator()(int i) { return values[i]; }
-   MFEM_HOST_DEVICE const tensor< T, n1, n2, n3, n4 >& operator()(int i) const { return values[i]; }
-   MFEM_HOST_DEVICE tensor< T, n2, n3, n4 >& operator()(int i, int j) { return values[i][j]; }
-   MFEM_HOST_DEVICE const tensor< T, n2, n3, n4 >& operator()(int i,
-                                                              int j) const { return values[i][j]; }
-   MFEM_HOST_DEVICE tensor< T, n3, n4>& operator()(int i, int j, int k) { return values[i][j][k]; }
-   MFEM_HOST_DEVICE const tensor< T, n3, n4>& operator()(int i, int j,
-                                                         int k) const { return values[i][j][k]; }
-   MFEM_HOST_DEVICE tensor< T, n4 >& operator()(int i, int j, int k, int l) { return values[i][j][k][l]; }
-   MFEM_HOST_DEVICE const tensor< T, n4 >& operator()(int i, int j, int k,
-                                                      int l) const { return values[i][j][k][l]; }
-   MFEM_HOST_DEVICE T& operator()(int i, int j, int k, int l, int m) { return values[i][j][k][l][m]; }
-   MFEM_HOST_DEVICE const T& operator()(int i, int j, int k, int l, int m) const { return values[i][j][k][l][m]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3, n4 >&
+   operator[](int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3, n4 >&
+   operator[](int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n1, n2, n3, n4 >&
+   operator()(int i) { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n1, n2, n3, n4 >&
+   operator()(int i) const { return values[i]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n2, n3, n4 >& operator()(
+      int i, int j) { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n2, n3, n4 >&
+   operator()(int i,
+              int j) const { return values[i][j]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n3, n4>& operator()(int i,
+                                                                             int j, int k) { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n3, n4>& operator()(
+      int i, int j,
+      int k) const { return values[i][j][k]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor< T, n4 >& operator()(int i,
+                                                                          int j, int k, int l) { return values[i][j][k][l]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const tensor< T, n4 >& operator()(
+      int i, int j, int k,
+      int l) const { return values[i][j][k][l]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE T& operator()(int i, int j, int k,
+                                                            int l, int m) { return values[i][j][k][l][m]; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE const T& operator()(int i, int j,
+                                                                  int k, int l, int m) const { return values[i][j][k][l][m]; }
    tensor < T, n1, n2, n3, n4 > values[n0];
 };
 
@@ -141,25 +265,25 @@ struct tensor<T, n0, n1, n2, n3, n4>
 struct zero
 {
    /** @brief `zero` is implicitly convertible to real_t with value 0.0 */
-   MFEM_HOST_DEVICE operator real_t() { return 0.0; }
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE operator real_t() { return 0.0; }
 
    /** @brief `zero` is implicitly convertible to a tensor of any shape */
    template <typename T, int... n>
-   MFEM_HOST_DEVICE operator tensor<T, n...>()
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE operator tensor<T, n...>()
    {
       return tensor<T, n...> {};
    }
 
    /** @brief `zero` can be accessed like a multidimensional array */
    template <typename... T>
-   MFEM_HOST_DEVICE zero operator()(T...)
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero operator()(T...)
    {
       return zero{};
    }
 
    /** @brief anything assigned to `zero` does not change its value and returns `zero` */
    template <typename T>
-   MFEM_HOST_DEVICE zero operator=(T)
+   inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero operator=(T)
    {
       return zero{};
    }
@@ -178,18 +302,19 @@ struct is_zero<zero> : std::true_type
 };
 
 /** @brief the sum of two `zero`s is `zero` */
-MFEM_HOST_DEVICE constexpr zero operator+(zero, zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator+(zero,
+                                                                    zero) { return zero{}; }
 
 /** @brief the sum of `zero` with something non-`zero` just returns the other value */
 template <typename T>
-MFEM_HOST_DEVICE constexpr T operator+(zero, T other)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr T operator+(zero, T other)
 {
    return other;
 }
 
 /** @brief the sum of `zero` with something non-`zero` just returns the other value */
 template <typename T>
-MFEM_HOST_DEVICE constexpr T operator+(T other, zero)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr T operator+(T other, zero)
 {
    return other;
 }
@@ -197,21 +322,22 @@ MFEM_HOST_DEVICE constexpr T operator+(T other, zero)
 /////////////////////////////////////////////////
 
 /** @brief the unary negation of `zero` is `zero` */
-MFEM_HOST_DEVICE constexpr zero operator-(zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator-(zero) { return zero{}; }
 
 /** @brief the difference of two `zero`s is `zero` */
-MFEM_HOST_DEVICE constexpr zero operator-(zero, zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator-(zero,
+                                                                    zero) { return zero{}; }
 
 /** @brief the difference of `zero` with something else is the unary negation of the other thing */
 template <typename T>
-MFEM_HOST_DEVICE constexpr T operator-(zero, T other)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr T operator-(zero, T other)
 {
    return -other;
 }
 
 /** @brief the difference of something else with `zero` is the other thing itself */
 template <typename T>
-MFEM_HOST_DEVICE constexpr T operator-(T other, zero)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr T operator-(T other, zero)
 {
    return other;
 }
@@ -219,52 +345,58 @@ MFEM_HOST_DEVICE constexpr T operator-(T other, zero)
 /////////////////////////////////////////////////
 
 /** @brief the product of two `zero`s is `zero` */
-MFEM_HOST_DEVICE constexpr zero operator*(zero, zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator*(zero,
+                                                                    zero) { return zero{}; }
 
 /** @brief the product `zero` with something else is also `zero` */
 template <typename T>
-MFEM_HOST_DEVICE constexpr zero operator*(zero, T /*other*/)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator*(zero,
+                                                                    T /*other*/)
 {
    return zero{};
 }
 
 /** @brief the product `zero` with something else is also `zero` */
 template <typename T>
-MFEM_HOST_DEVICE constexpr zero operator*(T /*other*/, zero)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator*(T /*other*/,
+                                                                    zero)
 {
    return zero{};
 }
 
 /** @brief `zero` divided by something is `zero` */
 template <typename T>
-MFEM_HOST_DEVICE constexpr zero operator/(zero, T /*other*/)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator/(zero,
+                                                                    T /*other*/)
 {
    return zero{};
 }
 
 /** @brief `zero` plus `zero` is `zero */
-MFEM_HOST_DEVICE constexpr zero operator+=(zero, zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator+=(zero,
+                                                                     zero) { return zero{}; }
 
 /** @brief `zero` minus `zero` is `zero */
-MFEM_HOST_DEVICE constexpr zero operator-=(zero, zero) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr zero operator-=(zero,
+                                                                     zero) { return zero{}; }
 
 /** @brief let `zero` be accessed like a tuple */
 template <int i>
-MFEM_HOST_DEVICE zero& get(zero& x)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero& get(zero& x)
 {
    return x;
 }
 
 /** @brief the dot product of anything with `zero` is `zero` */
 template <typename T>
-MFEM_HOST_DEVICE zero dot(const T&, zero)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero dot(const T&, zero)
 {
    return zero{};
 }
 
 /** @brief the dot product of anything with `zero` is `zero` */
 template <typename T>
-MFEM_HOST_DEVICE zero dot(zero, const T&)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero dot(zero, const T&)
 {
    return zero{};
 }
@@ -296,7 +428,8 @@ using reduced_tensor = typename std::conditional<
  *       to work around a limitation in nvcc involving __host__ __device__ lambdas with `auto` parameters.
  */
 template <typename lambda_type>
-MFEM_HOST_DEVICE constexpr auto make_tensor(lambda_type f) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE constexpr auto make_tensor(
+   lambda_type f) ->
 tensor<decltype(f())>
 {
    return {f()};
@@ -314,7 +447,7 @@ tensor<decltype(f())>
  *       to work around a limitation in nvcc involving __host__ __device__ lambdas with `auto` parameters.
  */
 template <int n1, typename lambda_type>
-MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
 tensor<decltype(f(n1)), n1>
 {
    using T = decltype(f(n1));
@@ -339,7 +472,7 @@ tensor<decltype(f(n1)), n1>
  *       to work around a limitation in nvcc involving __host__ __device__ lambdas with `auto` parameters.
  */
 template <int n1, int n2, typename lambda_type>
-MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
 tensor<decltype(f(n1, n2)), n1, n2>
 {
    using T = decltype(f(n1, n2));
@@ -368,7 +501,7 @@ tensor<decltype(f(n1, n2)), n1, n2>
  *       to work around a limitation in nvcc involving __host__ __device__ lambdas with `auto` parameters.
  */
 template <int n1, int n2, int n3, typename lambda_type>
-MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
 tensor<decltype(f(n1, n2, n3)), n1, n2, n3>
 {
    using T = decltype(f(n1, n2, n3));
@@ -401,7 +534,7 @@ tensor<decltype(f(n1, n2, n3)), n1, n2, n3>
  *       to work around a limitation in nvcc involving __host__ __device__ lambdas with `auto` parameters.
  */
 template <int n1, int n2, int n3, int n4, typename lambda_type>
-MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto make_tensor(lambda_type f) ->
 tensor<decltype(f(n1, n2, n3, n4)), n1, n2, n3, n4>
 {
    using T = decltype(f(n1, n2, n3, n4));
@@ -448,8 +581,9 @@ tensor<T, 1> get_col(tensor<T, 1, 1> A, int j)
  * @param[in] B The righthand operand
  */
 template <typename S, typename T, int... n>
-MFEM_HOST_DEVICE auto operator+(const tensor<S, n...>& A,
-                                const tensor<T, n...>& B) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator+
+(const tensor<S, n...>& A,
+ const tensor<T, n...>& B) ->
 tensor<decltype(S {} + T{}), n...>
 {
    tensor<decltype(S{} + T{}), n...> C{};
@@ -467,7 +601,8 @@ tensor<decltype(S {} + T{}), n...>
  * @param[in] A The tensor to negate
  */
 template <typename T, int... n>
-MFEM_HOST_DEVICE tensor<T, n...> operator-(const tensor<T, n...>& A)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor<T, n...> operator-
+(const tensor<T, n...>& A)
 {
    tensor<T, n...> B{};
    for (int i = 0; i < tensor<T, n...>::first_dim; i++)
@@ -486,8 +621,9 @@ MFEM_HOST_DEVICE tensor<T, n...> operator-(const tensor<T, n...>& A)
  * @param[in] B The righthand operand
  */
 template <typename S, typename T, int... n>
-MFEM_HOST_DEVICE auto operator-(const tensor<S, n...>& A,
-                                const tensor<T, n...>& B) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator-
+(const tensor<S, n...>& A,
+ const tensor<T, n...>& B) ->
 tensor<decltype(S {} + T{}), n...>
 {
    tensor<decltype(S{} + T{}), n...> C{};
@@ -509,7 +645,8 @@ tensor<decltype(S {} + T{}), n...>
 template <typename S, typename T, int... n,
           typename = typename std::enable_if<std::is_arithmetic<S>::value ||
                                              is_dual_number<S>::value>::type>
-MFEM_HOST_DEVICE auto operator*(S scale, const tensor<T, n...>& A) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator*(S scale,
+                                                          const tensor<T, n...>& A) ->
 tensor<decltype(S {} * T{}), n...>
 {
    tensor<decltype(S{} * T{}), n...> C{};
@@ -531,7 +668,8 @@ tensor<decltype(S {} * T{}), n...>
 template <typename S, typename T, int... n,
           typename = typename std::enable_if<std::is_arithmetic<S>::value ||
                                              is_dual_number<S>::value>::type>
-MFEM_HOST_DEVICE auto operator*(const tensor<T, n...>& A, S scale) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator*
+(const tensor<T, n...>& A, S scale) ->
 tensor<decltype(T {} * S{}), n...>
 {
    tensor<decltype(T{} * S{}), n...> C{};
@@ -553,7 +691,8 @@ tensor<decltype(T {} * S{}), n...>
 template <typename S, typename T, int... n,
           typename = typename std::enable_if<std::is_arithmetic<S>::value ||
                                              is_dual_number<S>::value>::type>
-MFEM_HOST_DEVICE auto operator/(S scale, const tensor<T, n...>& A) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator/(S scale,
+                                                          const tensor<T, n...>& A) ->
 tensor<decltype(S {} * T{}), n...>
 {
    tensor<decltype(S{} * T{}), n...> C{};
@@ -575,7 +714,8 @@ tensor<decltype(S {} * T{}), n...>
 template <typename S, typename T, int... n,
           typename = typename std::enable_if<std::is_arithmetic<S>::value ||
                                              is_dual_number<S>::value>::type>
-MFEM_HOST_DEVICE auto operator/(const tensor<T, n...>& A, S scale) ->
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE auto operator/
+(const tensor<T, n...>& A, S scale) ->
 tensor<decltype(T {} * S{}), n...>
 {
    tensor<decltype(T{} * S{}), n...> C{};
@@ -1314,7 +1454,8 @@ tensor<T, n, n> dev(const tensor<T, n, n>& A)
  * @return I_dim
  */
 template <int dim>
-MFEM_HOST_DEVICE tensor<real_t, dim, dim> IdentityMatrix()
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor<real_t, dim, dim>
+IdentityMatrix()
 {
    tensor<real_t, dim, dim> I{};
    for (int i = 0; i < dim; i++)
@@ -1651,13 +1792,15 @@ tensor<T, n> linear_solve(tensor<T, n, n> A, const tensor<T, n> b)
  * @note Uses a shortcut for inverting a 1x1, 2x2 and 3x3 matrix
  */
 template <typename T>
-inline MFEM_HOST_DEVICE tensor<T, 1, 1> inv(const tensor<T, 1, 1>& A)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor<T, 1, 1> inv(
+   const tensor<T, 1, 1>& A)
 {
    return tensor<T, 1, 1> {{{T{1.0} / A[0][0]}}};
 }
 
 template <typename T>
-inline MFEM_HOST_DEVICE tensor<T, 2, 2> inv(const tensor<T, 2, 2>& A)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor<T, 2, 2> inv(
+   const tensor<T, 2, 2>& A)
 {
    T inv_detA(1.0 / det(A));
 
@@ -1676,7 +1819,8 @@ inline MFEM_HOST_DEVICE tensor<T, 2, 2> inv(const tensor<T, 2, 2>& A)
  * @note Uses a shortcut for inverting a 3-by-3 matrix
  */
 template <typename T>
-inline MFEM_HOST_DEVICE tensor<T, 3, 3> inv(const tensor<T, 3, 3>& A)
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE tensor<T, 3, 3> inv(
+   const tensor<T, 3, 3>& A)
 {
    T inv_detA(1.0 / det(A));
 
@@ -1903,7 +2047,7 @@ using outer_product_t = typename detail::outer_prod<T1, T2>::type;
  * @brief Retrieves the gradient component of a real_t (which is nothing)
  * @return The sentinel, @see zero
  */
-inline MFEM_HOST_DEVICE zero get_gradient(real_t /* arg */) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero get_gradient(real_t /* arg */) { return zero{}; }
 
 /**
  * @brief get the gradient of type `tensor` (note: since its stored type is not a dual
@@ -1919,8 +2063,9 @@ MFEM_HOST_DEVICE zero get_gradient(const tensor<real_t, n...>& /* arg */)
 /**
  * @brief evaluate the change (to first order) in a function, f, given a small change in the input argument, dx.
  */
-inline MFEM_HOST_DEVICE zero chain_rule(const zero /* df_dx */,
-                                        const zero /* dx */) { return zero{}; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE zero chain_rule(
+   const zero /* df_dx */,
+   const zero /* dx */) { return zero{}; }
 
 /**
  * @overload
@@ -1948,8 +2093,8 @@ MFEM_HOST_DEVICE zero chain_rule(const T /* df_dx */,
  * @overload
  * @note for a scalar-valued function of a scalar, the chain rule is just multiplication
  */
-inline MFEM_HOST_DEVICE real_t chain_rule(const real_t df_dx,
-                                          const real_t dx) { return df_dx * dx; }
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE real_t chain_rule(const real_t df_dx,
+                                                             const real_t dx) { return df_dx * dx; }
 
 /**
  * @overload

makefile

@@ -805,7 +805,7 @@ FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
 COUT_CERR_FILES = $(foreach dir,$(DIRS),$(dir)/*.[ch]pp)
-COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp'
+COUT_CERR_EXCLUDE = '^general/error\.cpp' '^general/globals\.[ch]pp' '^general/nvtx\.hpp'
 
 DEPRECATION_WARNING := \
 "This feature is planned for removal in the next release."\

miniapps/dfem/dfem-minimal-surface.cpp

@@ -47,6 +47,8 @@
 //               visualization.
 
 #include "mfem.hpp"
+#include "fem/dfem/doperator.hpp"
+// #include <roctracer/roctx.h>
 
 using namespace mfem;
 
@@ -74,7 +76,7 @@ enum DerivativeType
 //
 // This class implements the minimal surface equation, which is a nonlinear
 // operator that provides the residual.
-template <typename dscalar_t, int dim = 2>
+template <typename dscalar_t, int dim = 3>
 class MinimalSurface : public Operator
 {
 private:
@@ -170,21 +172,28 @@ private:
          // DifferentiableOperator::AddDomainIntegrator call.
          dres_du = minsurface->res->GetDerivative(
                       SOLUTION_U, {&minsurface->u}, {mesh_nodes});
+
+         constr_op.reset(new ConstrainedOperator(dres_du.get(), minsurface->ess_tdofs));
       }
 
       void Mult(const Vector &x, Vector &y) const override
       {
-         z = x;
-         z.SetSubVector(minsurface->ess_tdofs, 0.0);
+         // MFEM_STREAM_SYNC;
+         // roctxRangePush("Operator::Mult");
+         // roctxRangePush("Essential BC PreMult");
+         // z = x;
+         // z.SetSubVector(minsurface->ess_tdofs, 0.0);
+         // roctxRangePop();
 
-         dres_du->Mult(z, y);
+         constr_op->Mult(x, y);
 
-         auto d_y = y.HostReadWrite();
-         const auto d_x = x.HostRead();
-         for (int i = 0; i < minsurface->ess_tdofs.Size(); i++)
-         {
-            d_y[minsurface->ess_tdofs[i]] = d_x[minsurface->ess_tdofs[i]];
-         }
+         // Reuse z as a temporary vector to avoid unnecessary allocations.
+         // roctxRangePush("Essential BC PostMult");
+         // x.GetSubVector(minsurface->ess_tdofs, z);
+         // y.SetSubVector(minsurface->ess_tdofs, z);
+         // roctxRangePop();
+         // MFEM_STREAM_SYNC;
+         // roctxRangePop();
       }
 
       // Pointer to the wrapped MinimalSurface operator
@@ -193,6 +202,8 @@ private:
       // Pointer to the DifferentiableOperator that computes the Jacobian
       std::shared_ptr<DerivativeOperator> dres_du;
 
+      std::unique_ptr<ConstrainedOperator> constr_op;
+
       // Temporary vector
       mutable Vector z;
    };
@@ -261,12 +272,9 @@ private:
 
          dres_du->Mult(z, y);
 
-         auto d_y = y.HostReadWrite();
-         const auto d_x = x.HostRead();
-         for (int i = 0; i < minsurface->ess_tdofs.Size(); i++)
-         {
-            d_y[minsurface->ess_tdofs[i]] = d_x[minsurface->ess_tdofs[i]];
-         }
+         // Reuse z as a temporary vector to avoid unnecessary allocations.
+         x.GetSubVector(minsurface->ess_tdofs, z);
+         y.SetSubVector(minsurface->ess_tdofs, z);
       }
 
       const MinimalSurface *minsurface = nullptr;
@@ -306,6 +314,7 @@ public:
       // Create the DifferentiableOperator on the desired mesh.
       res = std::make_shared<DifferentiableOperator>(
                solutions, parameters, *H1.GetParMesh());
+      res->UseAutomaticPA();
 
       // DifferentiableOperator::AddIntegrator consists mainly of multiple
       // components. The input and output operators and the pointwise
@@ -366,6 +375,7 @@ public:
       Array<int> ess_bdr(H1.GetParMesh()->bdr_attributes.Max());
       ess_bdr = 1;
       H1.GetEssentialTrueDofs(ess_bdr, ess_tdofs);
+      ess_tdofs.UseDevice();
    }
 
    void Mult(const Vector &x, Vector &y) const override
@@ -420,7 +430,7 @@ real_t boundary_func(const Vector &coords)
    const real_t y = coords(1);
    if (coords.Size() == 3)
    {
-      MFEM_ABORT("internal error");
+      // MFEM_ABORT("internal error");
    }
    const real_t a = 1.0e-2;
    return log(cos(a * x) / cos(a * y)) / a;
@@ -460,7 +470,7 @@ int main(int argc, char *argv[])
    if (Mpi::Root()) { device.Print(); }
 
    // 4. Create a 2D mesh on the square domain [-π/2,π/2]^2
-   Mesh mesh = Mesh::MakeCartesian2D(4, 4, Element::QUADRILATERAL);
+   Mesh mesh = Mesh::MakeCartesian3D(2, 2, 2, Element::HEXAHEDRON);
    mesh.SetCurvature(order);
 
    auto transform_mesh = [](const Vector &cold, Vector &cnew)
@@ -490,11 +500,18 @@ int main(int argc, char *argv[])
 
    // 8. Set up the integration rule
    const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(),
-                                  2 * order + 1);
+                                  2 * order + 2);
+
+   out << "Using integration rule of order: "
+       << ir->GetOrder() << "\n#qp: " << ir->GetNPoints()
+       << "\n#qp1d: " << floor(std::pow(ir->GetNPoints(), 1.0 / dim) + 0.5)
+       << std::endl;
 
    ParGridFunction u(&H1);
    Vector X(H1.GetTrueVSize());
 
+   out << "H1.GetTrueVSize():" << H1.GetTrueVSize() << std::endl;
+
    // 9. Create the nonlinear operator for the minimal surface equation
    std::unique_ptr<Operator> minsurface;
 #ifdef MFEM_USE_ENZYME
@@ -527,6 +544,57 @@ int main(int argc, char *argv[])
    krylov.SetMaxIter(500);
    krylov.SetPrintLevel(2);
 
+   auto &A = minsurface->GetGradient(X);
+   krylov.SetOperator(A);
+   Vector B(H1.GetTrueVSize());
+   B.UseDevice(true);
+   X.UseDevice(true);
+   B = 1.0;
+
+   out << "#dofs: " << H1.GetTrueVSize() << std::endl;
+
+   const int num_iterations = 1000;
+   const int num_runs = 10;
+   real_t min = std::numeric_limits<real_t>::max();
+   real_t max = 0;
+   real_t avg = 0;
+   for (int runs = 0; runs < num_runs; runs++)
+   {
+      MFEM_DEVICE_SYNC;
+      MPI_Barrier(MPI_COMM_WORLD);
+      tic();
+      for (int i = 0; i < num_iterations; i++)
+      {
+         MFEM_STREAM_SYNC;
+         // roctxRangePush("Operator::Mult");
+         A.Mult(B, X);
+         MFEM_STREAM_SYNC;
+         // roctxRangePop();
+      }
+      MFEM_DEVICE_SYNC;
+      MPI_Barrier(MPI_COMM_WORLD);
+      real_t elapsed_time = tic_toc.RealTime();
+      const real_t mdofs = H1.GetTrueVSize() * 1e-6;
+      auto cur_perf = mdofs * num_iterations / elapsed_time;
+      min = std::min(min, cur_perf);
+      max = std::max(max, cur_perf);
+      avg += cur_perf;
+      if (Mpi::Root())
+      {
+         out << "performance: " << cur_perf << " MDOFs/s" << std::endl;
+      }
+   }
+
+   avg /= num_runs;
+   if (Mpi::Root())
+   {
+      out << "min: " << min << "\n"
+          << "avg: " << avg << "\n"
+          << "max: " << max << "\n";
+   }
+
+   exit(0);
+
    // 13. Set up the nonlinear solver (Newton) for the minimal surface equation
    NewtonSolver newton(MPI_COMM_WORLD);
    newton.SetOperator(*minsurface);

tests/benchmarks/CMakeLists.txt

@@ -31,11 +31,6 @@ function(add_benchmark name)
        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
-    if (MFEM_USE_HIP)
-        set_property(SOURCE ${${NAME}_BENCH_SRCS} PROPERTY LANGUAGE 
-                     HIP_SOURCE_PROPERTY_FORMAT TRUE)
-    endif(MFEM_USE_HIP)
-
     add_executable(bench_${name} ${${NAME}_BENCH_SRCS})
     target_link_libraries(bench_${name} mfem pthread)
     add_dependencies(${MFEM_ALL_BENCHMARKS_TARGET_NAME} bench_${name})
@@ -56,6 +51,7 @@ endfunction(add_benchmark)
 #-------------------------------------------------------------------------------
 add_benchmark(assembly_levels)
 add_benchmark(ceed)
+add_benchmark(dfem)
 add_benchmark(dg_amr)
 add_benchmark(elasticity)
 add_benchmark(tmop)

tests/benchmarks/bench_ceed.cpp

@@ -119,7 +119,7 @@ struct Problem : public BakeOff<VDIM, GLL>
       cg.SetMaxIter(max_it);
       cg.SetPrintLevel(print_lvl);
       cg.iterative_mode = false;
-      MFEM_DEVICE_SYNC;
+      benchmark();
    }
 
    void benchmark() override

tests/benchmarks/bench_dfem.cpp

@@ -0,0 +1,892 @@
+// 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 "bench.hpp" // IWYU pragma: keep
+
+#ifdef MFEM_USE_BENCHMARK
+
+#include <cstdlib>
+#include <memory>
+
+#include <fem/qinterp/det.hpp>
+#include <fem/qinterp/grad.hpp> // IWYU pragma: keep
+
+#include "fem/dfem/doperator.hpp"
+#include <linalg/tensor.hpp>
+#include "tests/benchmarks/kernels_vd.hpp"
+#include "tests/benchmarks/kernels_fma.hpp"
+
+#include "fem/kernels.hpp"
+namespace ker = kernels::internal;
+
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kNvidia
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
+using namespace mfem;
+
+using mfem::future::tuple;
+using mfem::future::tensor;
+
+using future::DifferentiableOperator;
+using future::DerivativeOperator;
+using future::UniformParameterSpace;
+using future::ParameterFunction;
+using future::FieldDescriptor;
+using future::make_tensor;
+using future::Gradient;
+using future::Weight;
+using future::Identity;
+
+#ifdef MFEM_USE_ENZYME
+using dscalar_t = real_t;
+#else
+using mfem::future::dual;
+using dscalar_t = dual<real_t, real_t>;
+#endif
+
+#if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) ||       \
+     (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
+
+template<int N>
+using MQZ = std::integral_constant<int, 0>;
+
+#else
+
+template<int N>
+using MQZ = std::integral_constant<int, N>;
+
+constexpr int SetMaxOf2(int n) { return mfem::kernels::internal::NextMultipleOf<2>(n); }
+#endif
+
+/// Hall of Fame //////////////////////////////////////////////////////////////
+/// Clang 20
+/*-----------------------------------------------------------------------------------------------
+Benchmark           Time             CPU   Iterations       Dofs     MDof/s          p    version
+-------------------------------------------------------------------------------------------------
+BP3/0/6/25       17.9 ms         17.9 ms           10    15.625k  30.1555/s (!lto)   6          0
+BP3/1/6/25       12.4 ms         12.4 ms           10    15.625k  40.3044/s (lto)    6          1
+BP3/2/6/25       18.1 ms         18.1 ms           10    15.625k  27.6771/s          6          2
+BP3/3/6/25       39.7 ms         39.7 ms           10    15.625k  12.5988/s          6          3
+BP3/4/6/25       30.4 ms         30.3 ms           10    15.625k  16.4853/s          6          4
+BP3/5/6/25       14.3 ms         14.3 ms           10    15.625k  34.9011/s (lto)    6          5
+BP3/6/6/25       47.2 ms         47.1 ms           10    15.625k   10.606/s          6          6
+BP3/7/6/25       46.8 ms         46.8 ms           10    15.625k  10.6893/s          6          7
+*/
+
+/// Max number of DOFs ////////////////////////////////////////////////////////
+#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+constexpr int MAX_NDOFS = 128 * 1024;
+#ifdef MFEM_DEBUG
+constexpr int NDOFS_INC = 4;
+#else
+constexpr int NDOFS_INC = 25;
+#endif
+#else
+constexpr int MAX_NDOFS = 10 * 1024 * 1024;
+constexpr int NDOFS_INC = 25;
+#endif
+
+/// Kernels Versions //////////////////////////////////////////////////////////
+enum class kernels_t: int
+{
+   PA_STD,
+   PA_NEW,
+   PA_NEW_VDD,
+   MF_DFEM,
+   PA_DFEM,
+   PA_DFEM_NEW,
+   AUTO_PA_DFEM,
+   AUTO_PA_DFEM_NEW,
+   PA_FMA_VDD,
+   SIZE,
+} ;
+
+constexpr int PA_STD = static_cast<int>(kernels_t::PA_STD);
+constexpr int PA_NEW = static_cast<int>(kernels_t::PA_NEW);
+constexpr int PA_NEW_VDD = static_cast<int>(kernels_t::PA_NEW_VDD);
+constexpr int MF_DFEM = static_cast<int>(kernels_t::MF_DFEM);
+constexpr int PA_DFEM = static_cast<int>(kernels_t::PA_DFEM);
+constexpr int PA_DFEM_NEW = static_cast<int>(kernels_t::PA_DFEM_NEW);
+constexpr int AUTO_PA_DFEM = static_cast<int>(kernels_t::AUTO_PA_DFEM);
+constexpr int AUTO_PA_DFEM_NEW = static_cast<int>(kernels_t::AUTO_PA_DFEM_NEW);
+constexpr int PA_FMA_VDD = static_cast<int>(kernels_t::PA_FMA_VDD);
+
+constexpr auto all_kernels =
+{
+   kernels_t::PA_STD,
+   kernels_t::PA_NEW,
+   kernels_t::PA_NEW_VDD,
+   kernels_t::MF_DFEM,
+   kernels_t::PA_DFEM,
+   kernels_t::PA_DFEM_NEW,
+   kernels_t::AUTO_PA_DFEM,
+   kernels_t::AUTO_PA_DFEM_NEW,
+   // kernels_t::PA_FMA_VDD
+};
+
+static void DumpVersionInfo()
+{
+   mfem::out << "\x1b[33m";
+   mfem::out << "version " << PA_STD << ": PA std\n";
+   mfem::out << "version " << PA_NEW << ": PA new\n";
+   mfem::out << "version " << PA_NEW_VDD << ": PA new (layout by vdim)\n";
+   mfem::out << "version " << MF_DFEM << ": MF ∂fem\n";
+   mfem::out << "version " << PA_DFEM << ": PA ∂fem\n";
+   mfem::out << "version " << PA_DFEM_NEW << ": PA ∂fem new\n";
+   mfem::out << "version " << AUTO_PA_DFEM << ": Auto PA ∂fem\n";
+   mfem::out << "version " << AUTO_PA_DFEM_NEW << ": Auto PA ∂fem new\n";
+   mfem::out << "version " << PA_FMA_VDD << ": PA FMA (layout by vdim)\n";
+   mfem::out << "\x1b[m\n" << std::flush;
+}
+
+/// Benchmarks Arguments //////////////////////////////////////////////////////
+static void OrderSideVersionArgs(bmi::Benchmark *b)
+{
+   const auto est = [](int c) { return (c + 1) * (c + 1) * (c + 1); };
+   for (auto k : all_kernels)
+   {
+      for (int p = 8; p >= 1; p -= 1)
+      {
+         for (int c = NDOFS_INC; est(c) <= MAX_NDOFS; c += NDOFS_INC)
+         {
+            b->Args({ static_cast<int>(k), p, c });
+         }
+      }
+   }
+}
+
+/// Globals ///////////////////////////////////////////////////////////////////
+Device *device_ptr = nullptr;
+static int gD1D = 0, gQ1D = 0, gCGNI = 0;
+static bool use_new_kernels = false;
+static bool use_kernels_specialization = true;
+
+/// StiffnessIntegrator ///////////////////////////////////////////////////////
+template <bool layout_by_vdim, bool use_fma>
+struct StiffnessIntegrator : public BilinearFormIntegrator
+{
+   const FiniteElementSpace *fes;
+   const real_t *B, *G, *DX;
+   int ne, d1d, q1d;
+   Vector J0, dx;
+
+public:
+   StiffnessIntegrator()
+   {
+      dbg();
+      if (!use_kernels_specialization) { return; }
+      // StiffnessKernels::template Specialization<2, 3>::Add();
+      StiffnessKernels::template Specialization<3, 4>::Add();
+      // StiffnessKernels::template Specialization<4, 5>::Add();
+      StiffnessKernels::template Specialization<5, 6>::Add();
+      // StiffnessKernels::template Specialization<6, 7>::Add();
+      StiffnessKernels::template Specialization<7, 8>::Add();
+      // StiffnessKernels::template Specialization<9, 10>::Add();
+   }
+
+   void AssemblePA(const FiniteElementSpace &fespace) override
+   {
+      fes = &fespace;
+      auto *mesh = fes->GetMesh();
+      const int DIM = mesh->Dimension();
+      ne = mesh->GetNE();
+      const auto p = fes->GetFE(0)->GetOrder();
+      const auto q = 2 * p + mesh->GetElementTransformation(0)->OrderW();
+      const auto type = mesh->GetElementBaseGeometry(0);
+      const IntegrationRule &ir = IntRules.Get(type, q);
+      const int NQPT = ir.GetNPoints();
+      d1d = p + 1;
+      q1d = IntRules.Get(Geometry::SEGMENT, ir.GetOrder()).GetNPoints();
+      MFEM_VERIFY(d1d == gD1D, "D1D mismatch: " << d1d << " != " << gD1D);
+      MFEM_VERIFY(q1d == gQ1D, "Q1D mismatch: " << q1d << " != " << gQ1D);
+      MFEM_VERIFY(NQPT == q1d * q1d * q1d, "");
+      const DofToQuad *maps =
+         &fes->GetFE(0)->GetDofToQuad(ir, DofToQuad::TENSOR);
+      const GridFunction *nodes = (mesh->EnsureNodes(), mesh->GetNodes());
+      const FiniteElementSpace *nfes = nodes->FESpace();
+      const int nVDIM = nfes->GetVDim();
+      dx.SetSize(nVDIM * DIM * NQPT * ne, Device::GetDeviceMemoryType());
+      J0.SetSize(nVDIM * DIM * NQPT * ne, Device::GetDeviceMemoryType());
+      dx.UseDevice(true), J0.UseDevice(true);
+      B = maps->B.Read(), G = maps->G.Read(), DX = dx.Read();
+
+      const Operator *NR =
+         nfes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      const QuadratureInterpolator *nqi = nfes->GetQuadratureInterpolator(ir);
+      nqi->SetOutputLayout(QVectorLayout::byVDIM);
+      const int nd = nfes->GetFE(0)->GetDof();
+      Vector xe(nVDIM * nd * ne, Device::GetDeviceMemoryType());
+      NR->Mult(*nodes, (xe.UseDevice(true), xe));
+      nqi->Derivatives(xe, J0);
+
+      const int Q1D = q1d;
+      const auto w_r = ir.GetWeights().Read();
+      const auto W = Reshape(w_r, q1d, q1d, q1d);
+      const auto J = Reshape(J0.Read(), 3, 3, q1d, q1d, q1d, ne);
+      auto DX_w = Reshape(dx.Write(), 3, 3, q1d, q1d, q1d, ne);
+
+      mfem::forall_3D(ne, Q1D, Q1D, Q1D,[=] MFEM_HOST_DEVICE(int e)
+      {
+         MFEM_FOREACH_THREAD_DIRECT(qz, z, Q1D)
+         {
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  const real_t w = W(qx, qy, qz);
+                  const real_t *Jtr = &J(0, 0, qx, qy, qz, e);
+                  const real_t detJ = kernels::Det<3>(Jtr);
+                  const real_t wd = w * detJ;
+                  const real_t D[9] = { wd, 0.0, 0.0,
+                                        0.0, wd, 0.0,
+                                        0.0, 0.0, wd
+                                      };
+                  real_t Jrt[9], A[9];
+                  kernels::CalcInverse<3>(Jtr, Jrt);
+                  kernels::MultABt(3, 3, 3, D, Jrt, A);
+                  kernels::Mult(3, 3, 3, A, Jrt, &DX_w(0, 0, qx, qy, qz, e));
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+      });
+   }
+
+   ///////////////////////////////////////////////////////////////////
+   /// BP3/1/6/25 | 12.5ms | 12.5ms | 10 | 15.625k | 40.1397/s | 6 | 1
+   template <int T_D1D = 0, int T_Q1D = 0>
+   static void StiffnessMult(const int NE, const real_t *b, const real_t *g,
+                             const real_t *dx, const real_t *xe, real_t *ye,
+                             const int d1d, const int q1d)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      db1("D1D:{} Q1D:{} (NEW, no VDIM)", D1D, Q1D);
+
+      constexpr int DIM = 3, VDIM = 1;
+      const auto XE = Reshape(xe, D1D, D1D, D1D, VDIM, NE);
+      const auto DX = Reshape(dx, 3, 3, Q1D, Q1D, Q1D, NE);
+      auto YE = Reshape(ye, D1D, D1D, D1D, VDIM, NE);
+
+      mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+      {
+         // constexpr int MD1 = T_D1D > 0 ? kernels::internal::SetMaxOf(T_D1D) : 32;
+         // constexpr int MQ1 = T_Q1D > 0 ? kernels::internal::SetMaxOf(T_Q1D) : 32;
+         constexpr int MD1 = T_D1D > 0 ? T_D1D : 32;
+         constexpr int MQ1 = T_Q1D > 0 ? T_Q1D : 32;
+
+         alignas(64) MFEM_SHARED real_t smem[MQ1][MQ1];
+         alignas(64) MFEM_SHARED real_t sB[MD1][MQ1];
+         alignas(64) MFEM_SHARED real_t sG[MD1][MQ1];
+         ker::vd_regs3d_t<VDIM, DIM, MQ1> r0, r1;
+
+         ker::LoadMatrix(D1D, Q1D, b, sB);
+         ker::LoadMatrix(D1D, Q1D, g, sG);
+
+         ker::LoadDofs3d(e, D1D, XE, r0);
+         ker::Grad3d(D1D, Q1D, smem, sB, sG, r0, r1);
+
+         for (int qz = 0; qz < Q1D; qz++)
+         {
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  real_t v[3], u[3] = { r1[0][0][qz][qy][qx],
+                                        r1[0][1][qz][qy][qx],
+                                        r1[0][2][qz][qy][qx]
+                                      };
+                  const real_t *dx = &DX(0, 0, qx, qy, qz, e);
+                  kernels::Mult(3, 3, dx, u, v);
+                  r0[0][0][qz][qy][qx] = v[0];
+                  r0[0][1][qz][qy][qx] = v[1];
+                  r0[0][2][qz][qy][qx] = v[2];
+               }
+            }
+         }
+         ker::GradTranspose3d(D1D, Q1D, smem, sB, sG, r0, r1);
+         ker::WriteDofs3d(e, D1D, r1, YE);
+      });
+   }
+
+   ///////////////////////////////////////////////////////////////////
+   // BP3/1/6/25 | 49.5 ms | 49.5 ms | 10 | 15.625k | 10.0983/s | 6 | 1
+   template <int T_D1D = 0, int T_Q1D = 0>
+   static void StiffnessMultVDD(const int NE, const real_t *b, const real_t *g,
+                                const real_t *dx, const real_t *xe, real_t *ye,
+                                const int d1d, const int q1d)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      db1("D1D:{} Q1D:{} (NEW & VDIM)", D1D, Q1D);
+
+      constexpr int DIM = 3, VDIM = 1;
+      const auto XE = Reshape(xe, D1D, D1D, D1D, VDIM, NE);
+      const auto DX = Reshape(dx, 3, 3, Q1D, Q1D, Q1D, NE);
+      auto YE = Reshape(ye, D1D, D1D, D1D, VDIM, NE);
+
+      mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+      {
+         // constexpr int MD1 = T_D1D > 0 ? kernels::internal::SetMaxOf(T_D1D) : 32;
+         // constexpr int MQ1 = T_Q1D > 0 ? kernels::internal::SetMaxOf(T_Q1D) : 32;
+         constexpr int MD1 = T_D1D ? T_D1D : 32;
+         constexpr int MQ1 = T_Q1D ? T_Q1D : 32;
+
+         alignas(64) MFEM_SHARED real_t smem[MQ1][MQ1];
+         alignas(64) MFEM_SHARED real_t sB[MD1][MQ1];
+         alignas(64) MFEM_SHARED real_t sG[MD1][MQ1];
+
+         constexpr int MZ1 = MQZ<MQ1>::value;
+         alignas(64) mfem::future::tensor<real_t, MQ1, MZ1, MZ1, VDIM, DIM> v0, v1;
+
+         ker::LoadMatrix(D1D, Q1D, b, sB);
+         ker::LoadMatrix(D1D, Q1D, g, sG);
+
+         ker::vd::LoadDofs3d(e, D1D, XE, v0);
+         ker::vd::Grad3d(D1D, Q1D, smem, sB, sG, v0, v1);
+
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; qz++)
+         {
+            MFEM_UNROLL(MQ1)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               MFEM_UNROLL(MQ1)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  auto &vd_v = v0[qz][qy][qx][0];
+                  const auto &vd_u = v1[qz][qy][qx][0];
+                  const auto d = make_tensor<DIM, DIM>([&](int i, int j)
+                  {
+                     return DX(i, j, qx, qy, qz, e);
+                  });
+                  vd_v = d * vd_u;
+               }
+            }
+         }
+         ker::vd::GradTranspose3d(D1D, Q1D, smem, sB, sG, v0, v1);
+         ker::vd::WriteDofs3d(e, D1D, v1, YE);
+      });
+   }
+
+   ///////////////////////////////////////////////////////////////////
+   template <int T_D1D = 0, int T_Q1D = 0>
+   static void StiffnessMultFMA(const int NE, const real_t *b,
+                                const real_t *g,
+                                const real_t *dx, const real_t *xe, real_t *ye,
+                                const int d1d, const int q1d)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      db1("D1D:{} Q1D:{} (FMA & VDIM)", D1D, Q1D);
+
+      constexpr int DIM = 3, VDIM = 1, NBZ = 1;
+      const auto XE = Reshape(xe, D1D, D1D, D1D, VDIM, NE);
+      const auto DX = Reshape(dx, 3, 3, Q1D, Q1D, Q1D, NE);
+      auto YE = Reshape(ye, D1D, D1D, D1D, VDIM, NE);
+
+      mfem::forall_3D(NE, Q1D, Q1D, NBZ, [=] MFEM_HOST_DEVICE(int e)
+      {
+         const int tz = MFEM_THREAD_ID(z);
+         constexpr int MD1 = T_D1D ? T_D1D : 32;
+         constexpr int MQ1 = T_Q1D ? T_Q1D : 32;
+
+         alignas(64) MFEM_SHARED real_t smem[NBZ][VDIM][DIM][MQ1][MQ1][MQ1];
+         alignas(64) MFEM_SHARED real_t sB[MD1][MQ1];
+         alignas(64) MFEM_SHARED real_t sG[MD1][MQ1];
+
+         constexpr int MZ1 = MQZ<MQ1>::value;
+         alignas(64) mfem::future::tensor<real_t, MQ1, MZ1, MZ1, VDIM, DIM> r_q;
+
+         ker::LoadMatrix(D1D, Q1D, b, sB);
+         ker::LoadMatrix(D1D, Q1D, g, sG);
+
+         ker::fma::Grad3d(D1D, Q1D, smem, sB, sG, r_q, XE, e);
+
+         MFEM_UNROLL(MQ1)
+         for (int qz = 0; qz < Q1D; qz++)
+         {
+            MFEM_UNROLL(MQ1)
+            MFEM_FOREACH_THREAD_DIRECT(qy, y, Q1D)
+            {
+               MFEM_UNROLL(MQ1)
+               MFEM_FOREACH_THREAD_DIRECT(qx, x, Q1D)
+               {
+                  /*const auto d = make_tensor<3, 3>([&](int i, int j)
+                  {
+                     return DX(i, j, qx, qy, qz, e);
+                  });
+                  const auto vd_u = make_tensor<VDIM, DIM>([&](int vd, int d)
+                  {
+                     return smem[tz][vd][d][qx][qy][qz];
+                  });
+                  auto &vd_v = r_q[qz][qy][qx];
+                  // const auto &vd_u = r_q[qz][qy][qx][0];
+                  // vd_v = d * vd_u;
+
+                  // auto &vd_v = r_q[qz][qy][qx][0];
+                  // vd_v = d * v;*/
+
+                  real_t v[3], u[3] = { smem[tz][0][0][qz][qy][qx],
+                                        smem[tz][0][1][qz][qy][qx],
+                                        smem[tz][0][2][qz][qy][qx]
+                                      };
+                  const real_t *dx = &DX(0, 0, qx, qy, qz, e);
+                  kernels::Mult(3, 3, dx, u, v);
+                  r_q[qz][qy][qx][0][0] = v[0];
+                  r_q[qz][qy][qx][0][1] = v[1];
+                  r_q[qz][qy][qx][0][2] = v[2];
+               }
+            }
+         }
+         ker::fma::GradTranspose3d(D1D, Q1D, smem, sB, sG, r_q, YE, e);
+      });
+   }
+
+   using StiffnessKernelType = decltype(&StiffnessMult<>);
+   MFEM_REGISTER_KERNELS(StiffnessKernels, StiffnessKernelType, (int, int));
+
+   void AddMultPA(const Vector &x, Vector &y) const override
+   {
+      StiffnessKernels::Run(d1d, q1d,
+                            ne, B, G, DX, x.Read(), y.ReadWrite(),
+                            d1d, q1d);
+   }
+};
+
+template<bool layout_by_vdim, bool use_fma>
+template <int D1D, int Q1D>
+inline MFEM_ALWAYS_INLINE
+typename StiffnessIntegrator<layout_by_vdim, use_fma>::StiffnessKernelType
+StiffnessIntegrator<layout_by_vdim, use_fma>::StiffnessKernels::Kernel()
+{
+   if constexpr (!layout_by_vdim && !use_fma) { return StiffnessMult<D1D, Q1D>; }
+   if constexpr (layout_by_vdim && !use_fma) { return StiffnessMultVDD<D1D, Q1D>; }
+   if constexpr (layout_by_vdim &&  use_fma) { return StiffnessMultFMA<D1D, Q1D>; }
+   MFEM_ABORT("Invalid specialization for StiffnessIntegrator");
+   return nullptr; // Should never happen
+}
+
+template<bool layout_by_vdim, bool use_fma>
+inline MFEM_ALWAYS_INLINE
+typename StiffnessIntegrator<layout_by_vdim, use_fma>::StiffnessKernelType
+StiffnessIntegrator<layout_by_vdim, use_fma>::StiffnessKernels::Fallback(
+   int d1d,
+   int q1d)
+{
+   dbg("\x1b[33mFallback d1d:{} q1d:{}", d1d, q1d);
+   // if constexpr (!layout_by_vdim && !use_fma) { return StiffnessMult<>; }
+   // if constexpr (layout_by_vdim && !use_fma) { return StiffnessMultVDD<>; }
+   // if constexpr (use_fma) { return StiffnessMultFMA<>; }
+   MFEM_ABORT("Invalid specialization for StiffnessIntegrator");
+   return nullptr; // Should never happen
+}
+
+/// BakeOff ///////////////////////////////////////////////////////////////////
+template <int VDIM, bool GLL>
+struct BakeOff
+{
+   static constexpr int DIM = 3;
+   const int p, c, q, n, nx, ny, nz;
+   const bool check_x, check_y, check_z, checked;
+   Mesh smesh;
+   ParMesh pmesh;
+   H1_FECollection fec;
+   ParFiniteElementSpace pfes;
+   const Geometry::Type geom_type;
+   IntegrationRules irs;
+   const IntegrationRule *ir;
+   ConstantCoefficient one;
+   Vector uvec;
+   VectorConstantCoefficient unit_vec;
+   const int dofs;
+   ParGridFunction *nodes;
+   ParFiniteElementSpace& mfes;
+   ParGridFunction x, y;
+   ParBilinearForm a;
+   std::unique_ptr<DifferentiableOperator> dop;
+   std::shared_ptr<future::DerivativeOperator> dRdU;
+   const int elem_size, total_size, d1d, q1d;
+   UniformParameterSpace qd_ps;
+   ParameterFunction qdata;
+
+   double mdofs{};
+
+   BakeOff(int p, int side):
+      p(p), c(side), q(2 * p + (GLL ? -1 : 3)), n((assert(c >= p), c / p)),
+      nx(n + (p * (n + 1) * p * n * p * n < c * c * c ? 1 : 0)),
+      ny(n + (p * (n + 1) * p * (n + 1) * p * n < c * c * c ? 1 : 0)), nz(n),
+      check_x(p * nx * p * ny * p * nz <= c * c * c),
+      check_y(p * (nx + 1) * p * (ny + 1) * p * nz > c * c * c),
+      check_z(p * (nx + 1) * p * (ny + 1) * p * (nz + 1) > c * c * c),
+      checked((assert(check_x &&check_y &&check_z), true)),
+      smesh(Mesh::MakeCartesian3D(nx, ny, nz, Element::HEXAHEDRON)),
+      pmesh(MPI_COMM_WORLD, (smesh.EnsureNodes(), smesh)),
+      fec(p, DIM, BasisType::GaussLobatto),
+      pfes(&pmesh, &fec, VDIM),
+      geom_type(pmesh.GetTypicalElementGeometry()),
+      irs(0, GLL ? Quadrature1D::GaussLobatto : Quadrature1D::GaussLegendre),
+      ir(&irs.Get(geom_type, q)), one(1.0), uvec(DIM),
+      unit_vec((uvec = 1.0, uvec /= uvec.Norml2(), uvec)),
+      dofs(pfes.GetTrueVSize()),
+      nodes(static_cast<ParGridFunction*>(pmesh.GetNodes())),
+      mfes(*nodes->ParFESpace()),
+      x(&pfes),
+      y(&pfes),
+      a(&pfes),
+      elem_size(DIM * DIM * ir->GetNPoints()),
+      total_size(elem_size * pmesh.GetNE()),
+      d1d(p + 1),
+      q1d(IntRules.Get(Geometry::SEGMENT, ir->GetOrder()).GetNPoints()),
+      qd_ps(pmesh, *ir, DIM*DIM),
+      qdata(qd_ps)
+   {
+      smesh.Clear();
+      x = 0.0;
+
+      gD1D = d1d, gQ1D = q1d;
+      db1("D1D: {}, Q1D: {}", gD1D, gQ1D);
+      qdata.UseDevice(true);
+      MFEM_VERIFY(q1d*q1d*q1d == ir->GetNPoints(), "");
+   }
+
+   virtual void benchmark() = 0;
+
+   double SumMdofs() const { return mdofs; }
+
+   double MDofs() const { return 1e-6 * dofs; }
+};
+
+/// Diffusion /////////////////////////////////////////////////////////////////
+template <int VDIM = 1, bool GLL = false>
+struct Diffusion : public BakeOff<VDIM, GLL>
+{
+   static constexpr int DIM = 3;
+   static constexpr int U = 0, Ξ = 1, Q = 2;
+
+   const real_t rtol = 0.0;
+   const int max_it = 32, print_lvl = -1;
+
+   Array<int> ess_tdof_list, ess_bdr, all_domain_attr;
+   ParLinearForm b;
+   FieldDescriptor u_fd, Ξ_fd, q_fd;
+   std::vector<FieldDescriptor> u_sol, q_param, Ξ_q_params;
+   OperatorPtr A;
+   Operator *A_ptr;
+   Vector B, X;
+   CGSolver cg;
+
+   using BakeOff<VDIM, GLL>::a;
+   using BakeOff<VDIM, GLL>::ir;
+   using BakeOff<VDIM, GLL>::one;
+   using BakeOff<VDIM, GLL>::pmesh;
+   using BakeOff<VDIM, GLL>::pfes;
+   using BakeOff<VDIM, GLL>::mfes;
+   using BakeOff<VDIM, GLL>::x;
+   using BakeOff<VDIM, GLL>::y;
+   using BakeOff<VDIM, GLL>::mdofs;
+   using BakeOff<VDIM, GLL>::dop;
+   using BakeOff<VDIM, GLL>::dRdU;
+   using BakeOff<VDIM, GLL>::nodes;
+   using BakeOff<VDIM, GLL>::qdata;
+   using BakeOff<VDIM, GLL>::qd_ps;
+   using BakeOff<VDIM, GLL>::dofs;
+
+   Diffusion(int version, int order, int side):
+      BakeOff<VDIM, GLL>(order, side),
+      ess_bdr(pmesh.bdr_attributes.Max()),
+      all_domain_attr(pmesh.bdr_attributes.Max()),
+      b(&pfes),
+      u_fd{U, &pfes}, Ξ_fd{Ξ, &mfes}, q_fd{Q, &qd_ps},
+      u_sol{u_fd},
+      q_param {q_fd},
+      Ξ_q_params {Ξ_fd, q_fd},
+      cg(MPI_COMM_WORLD)
+   {
+      static_assert(VDIM == 1 && GLL == false);
+
+      ess_bdr = 1;
+      all_domain_attr = 1;
+      pfes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+      b.AddDomainIntegrator(new DomainLFIntegrator(this->one));
+      b.UseFastAssembly(true);
+      b.Assemble();
+
+      // PA_STD, PA_NEW, PA_NEW_VD & PA_FMA_VDD
+      if (version == PA_STD    || version == PA_NEW ||
+          version == PA_NEW_VDD || version == PA_FMA_VDD)
+      {
+         a.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+         if (version == PA_STD) { a.AddDomainIntegrator(new DiffusionIntegrator(ir)); }
+         if (version == PA_NEW) { a.AddDomainIntegrator(new StiffnessIntegrator<false, false>()); }
+         if (version == PA_NEW_VDD) { a.AddDomainIntegrator(new StiffnessIntegrator<true, false>()); }
+         if (version == PA_FMA_VDD) { a.AddDomainIntegrator(new StiffnessIntegrator<true, true>()); }
+         a.Assemble();
+         a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+         if (version == PA_STD)
+         {
+            BilinearFormIntegrator *bfi = a.GetDBFI()->operator[](0);
+            auto *di = dynamic_cast<DiffusionIntegrator*>(bfi);
+            const int d1d = di->dofs1D, q1d = di->quad1D;
+            MFEM_VERIFY(d1d == gD1D, "D1D mismatch: " << d1d << " != " << gD1D);
+            MFEM_VERIFY(q1d == gQ1D, "Q1D mismatch: " << q1d << " != " << gQ1D);
+         }
+      }
+      // ∂fem: MF
+      else if (version == MF_DFEM)
+      {
+         dbg("MF ∂fem");
+         auto solutions = std::vector{FieldDescriptor{U, &pfes}};
+         auto parameters = std::vector{FieldDescriptor{Ξ, &mfes}};
+         dop = std::make_unique<DifferentiableOperator>(solutions, parameters, pmesh);
+         dop->SetParameters({nodes});
+         const auto diffusion_mf_kernel =
+            [] MFEM_HOST_DEVICE (const tensor<real_t, DIM>& Grad_u,
+                                 const tensor<real_t, DIM, DIM>& J,
+                                 const real_t& w)
+         {
+            const auto invJ = inv(J);
+            return tuple{((Grad_u * invJ)) * transpose(invJ) * det(J) * w};
+         };
+         dop->AddDomainIntegrator(diffusion_mf_kernel,
+                                  tuple{Gradient<U>{}, Gradient<Ξ>{}, Weight{}},
+                                  tuple{Gradient<U>{}},
+                                  *ir, ess_bdr);
+         dop->FormLinearSystem(ess_tdof_list, x, b, A_ptr, X, B);
+         A.Reset(A_ptr);
+      }
+      // ∂fem: PA, PA NEW
+      else if (version == PA_DFEM || version == PA_DFEM_NEW)
+      {
+         dbg("[PA ∂fem] setup");
+         const auto pa_setup_qf =
+            [] MFEM_HOST_DEVICE(const real_t &u,
+                                const tensor<real_t, DIM, DIM> &J,
+                                const real_t &w)
+         {
+            const auto invJ = inv(J);
+            return tuple{invJ * transpose(invJ) * det(J) * w};
+         };
+         DifferentiableOperator dSetup(u_sol, Ξ_q_params, pmesh);
+         dSetup.AddDomainIntegrator(pa_setup_qf,
+                                    tuple{Identity<U> {}, Gradient<Ξ> {}, Weight{}},
+                                    tuple{Identity<Q> {}},
+                                    *ir, ess_bdr);
+         dSetup.SetParameters({nodes, &qdata});
+         X.SetSize(pfes.GetTrueVSize());
+         pfes.GetRestrictionMatrix()->Mult(x, X);
+         dSetup.Mult(X, qdata);
+
+         dbg("[PA ∂fem] apply");
+         const auto pa_apply_qf =
+            [] MFEM_HOST_DEVICE(const tensor<real_t, DIM> &Grad_u,
+                                const tensor<real_t, DIM, DIM> &Q)
+         {
+            return tuple{Grad_u * Q};
+         };
+         dop = std::make_unique<DifferentiableOperator>(u_sol, q_param, pmesh);
+         if (version == PA_DFEM_NEW) { dop->UseNewKernels(); }
+         if (use_kernels_specialization) { dop->UseKernelsSpecialization(); }
+         else { dbg("[PA ∂fem] NOT using kernels specialization"); }
+         dop->AddDomainIntegrator(pa_apply_qf,
+                                  tuple{Gradient<U> {}, Identity<Q> {}},
+                                  tuple{Gradient<U> {}},
+                                  *ir, ess_bdr);
+         dop->SetParameters({ &qdata });
+
+         dop->FormLinearSystem(ess_tdof_list, x, b, A_ptr, X, B);
+         A.Reset(A_ptr);
+      }
+      // ∂fem: Auto PA, AUTO PA NEW
+      else if (version == AUTO_PA_DFEM || version == AUTO_PA_DFEM_NEW)
+      {
+         dbg("[Auto PA ∂fem]");
+         const auto solutions = std::vector{FieldDescriptor{U, &pfes}};
+         const auto parameters = std::vector{FieldDescriptor{Ξ, &mfes}};
+         const auto derivatives = std::integer_sequence<size_t, U> {};
+
+         dop = std::make_unique<DifferentiableOperator>(solutions, parameters, pmesh);
+         if (version == AUTO_PA_DFEM_NEW) { dop->UseNewKernels(); }
+         dop->UseAutomaticPA();
+         dop->UseKernelsSpecialization();
+
+         {
+            auto stiffness_integrator = new StiffnessIntegrator<false, false>();
+            stiffness_integrator->AssemblePA(pfes);
+            dop->UsePaData(stiffness_integrator->dx.Read());
+         }
+
+         const auto diffusion_mf_kernel =
+            [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM>& Grad_u,
+                                 const tensor<real_t, DIM, DIM>& J,
+                                 const real_t& w)
+         {
+            const auto invJ = inv(J), TinJ = transpose(invJ);
+            return tuple{((Grad_u * invJ)) * TinJ * det(J) * w};
+         };
+         dop->AddDomainIntegrator(diffusion_mf_kernel,
+                                  tuple{Gradient<U>{}, Gradient<Ξ>{}, Weight{}},
+                                  tuple{Gradient<U>{}},
+                                  *ir, ess_bdr, derivatives);
+         dop->SetParameters({nodes});
+         dRdU = dop->GetDerivative(U, {&x}, {nodes});
+         dRdU->FormLinearSystem(ess_tdof_list, x, b, A_ptr, X, B);
+         A.Reset(A_ptr);
+      }
+      else { MFEM_ABORT("Invalid version"); }
+
+      cg.SetOperator(*A);
+      cg.iterative_mode = false;
+      if (dofs < 128 * 1024) // check
+      {
+#ifdef MFEM_DEBUG
+         cg.SetPrintLevel(3);
+#else
+         cg.SetPrintLevel(3/*-1*/);
+#endif
+         cg.SetMaxIter(100);//2000);
+         cg.SetRelTol(1e-8);
+         cg.SetAbsTol(0.0);
+         cg.Mult(B, X);
+         MFEM_VERIFY(cg.GetConverged(), "❌ CG did not converge.");
+         const int lCGNI = cg.GetNumIterations();
+         static const bool iniCGNumIterations = (gCGNI = lCGNI, true);
+         if (iniCGNumIterations)
+         {
+            MFEM_VERIFY(lCGNI == gCGNI, "❌ CG iters " << lCGNI << " != " << gCGNI);
+         }
+         MFEM_DEVICE_SYNC;
+      }
+      cg.SetAbsTol(0.0);
+      cg.SetRelTol(rtol);
+      cg.SetMaxIter(max_it);
+      cg.SetPrintLevel(print_lvl);
+      benchmark();
+      mdofs = 0.0;
+   }
+
+   void benchmark() override
+   {
+      cg.Mult(B, X);
+      MFEM_DEVICE_SYNC;
+      mdofs += this->MDofs() * cg.GetNumIterations();
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+#define BakeOff_Problem(i, Problem)                                  \
+   static void BP##i(bm::State &state)                               \
+   {                                                                 \
+      const auto version = static_cast<int>(state.range(0));         \
+      const auto order = static_cast<int>(state.range(1));           \
+      const auto side = static_cast<int>(state.range(2));            \
+      Problem ker(version, order, side);                             \
+      while (state.KeepRunning()) { ker.benchmark(); }               \
+      bm::Counter::Flags flags = bm::Counter::kIsRate;               \
+      state.counters["MDof/s"] = bm::Counter(ker.SumMdofs(), flags); \
+      state.counters["Dofs"] = bm::Counter(ker.dofs);                \
+      state.counters["p"] = bm::Counter(order);                      \
+      state.counters["version"] = bm::Counter(version);              \
+   }                                                                 \
+   BENCHMARK(BP##i)                                                  \
+      ->Apply(OrderSideVersionArgs)                                  \
+      ->Unit(bm::kMillisecond)
+
+BakeOff_Problem(3, Diffusion);
+
+/// Basic Kernels Specializations /////////////////////////////////////////////
+static void AddBasicKernelSpecializations()
+{
+   using Det = QuadratureInterpolator::DetKernels;
+   Det::Specialization<3, 3, 2, 2>::Add();
+   Det::Specialization<3, 3, 2, 3>::Add();
+   Det::Specialization<3, 3, 2, 5>::Add();
+   Det::Specialization<3, 3, 2, 6>::Add();
+   Det::Specialization<3, 3, 2, 7>::Add();
+
+   using Grad = QuadratureInterpolator::GradKernels;
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 3>::Add();
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 4>::Add();
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 5>::Add();
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 6>::Add();
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 7>::Add();
+   Grad::Specialization<3, QVectorLayout::byVDIM,  false, 3, 2, 8>::Add();
+   Grad::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 7>::Add();
+   Grad::Specialization<3, QVectorLayout::byNODES, false, 3, 2, 8>::Add();
+}
+
+/// main //////////////////////////////////////////////////////////////////////
+int main(int argc, char *argv[])
+{
+   dbg();
+   DumpVersionInfo();
+   AddBasicKernelSpecializations();
+   static mfem::MPI_Session mpi(argc, argv);
+
+   bm::ConsoleReporter CR;
+   bm::Initialize(&argc, argv);
+
+   // Device setup, cpu by default
+   std::string device_context = "cpu",
+               kernels_context = "std",
+               kernels_specialization = "yes";
+   const auto global_context = bmi::GetGlobalContext();
+   if (global_context != nullptr)
+   {
+      const auto device = global_context->find("device");
+      if (device != global_context->end())
+      {
+         mfem::out << device->first << " : "
+                   << device->second << std::endl;
+         device_context = device->second;
+      }
+
+      const auto kernels = global_context->find("kernels");
+      if (kernels != global_context->end())
+      {
+         mfem::out << kernels->first << " : "
+                   << kernels->second << std::endl;
+         kernels_context = kernels->second;
+         MFEM_VERIFY(kernels_context == "std" || kernels_context == "new",
+                     "Invalid kernels config: " << kernels_context);
+         use_new_kernels = (kernels_context == "new");
+      }
+
+      const auto specialization = global_context->find("specialization");
+      if (specialization != global_context->end())
+      {
+         mfem::out << specialization->first << " : "
+                   << specialization->second << std::endl;
+         kernels_specialization = specialization->second;
+         MFEM_VERIFY(kernels_specialization == "yes" || kernels_specialization == "no",
+                     "Invalid kernels specialization config: " << kernels_specialization);
+         use_kernels_specialization = (kernels_specialization == "yes");
+      }
+   }
+   dbg("device_config: {}", device_context);
+   Device device(device_context.c_str());
+   device_ptr = &device;
+   device.Print();
+
+   if (bm::ReportUnrecognizedArguments(argc, argv)) { return EXIT_FAILURE; }
+
+   bm::RunSpecifiedBenchmarks(&CR);
+
+   return EXIT_SUCCESS;
+}
+
+#endif // MFEM_USE_BENCHMARK

tests/benchmarks/kernels_fma.hpp

@@ -0,0 +1,264 @@
+// 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
+
+#include "../../config/config.hpp"
+#include "../../config/tconfig.hpp" // MFEM_ALWAYS_INLINE
+// #include "../../linalg/dtensor.hpp"
+// #include "../../linalg/tensor.hpp"
+// #include "../../general/forall.hpp" // MFEM_UNROLL
+#include "./kernels_vd.hpp"
+
+using namespace mfem::kernels::internal::vd;
+
+namespace mfem::kernels::internal::fma
+{
+
+// Grad1X
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE void
+ContractX3d(const int D1D, const int Q1D,
+            real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+            const real_t (*B)[MQ1],
+            const DeviceTensor<5, const real_t> &XE,
+            const int vd, const int d, const int e)
+{
+   const int tz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(b, y, D1D)
+   {
+      MFEM_FOREACH_THREAD(a, x, D1D)
+      {
+         // MFEM_UNROLL(Q1D)
+         for (int k = 0; k < Q1D; ++k)
+         {
+            double u = 0.0;
+            // MFEM_UNROLL(D1D)
+            for (int c = 0; c < D1D; ++c) { u += B[c][k] * XE(c,b,a, vd,e); }
+            smem[tz][vd][d][k][b][a] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+// Grad1Y
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE void
+ContractY3d(const int D1D, const int Q1D,
+            real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+            const real_t (*B)[MQ1],
+            regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+            const int vd, const int d)
+{
+   const int tz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(k, y, Q1D)
+   {
+      MFEM_FOREACH_THREAD(a, x, D1D)
+      {
+         // MFEM_UNROLL(D1D)
+         for (int b = 0; b < D1D; ++b) { r_q[k][a][b][vd][d] = smem[tz][vd][d][k][b][a]; }
+         // MFEM_UNROLL(Q1D)
+         for (int j = 0; j < Q1D; ++j)
+         {
+            double u = 0.0;
+            // MFEM_UNROLL(D1D)
+            for (int b = 0; b < D1D; ++b) { u += B[b][j] * r_q[k][a][b][vd][d]; }
+            smem[tz][vd][d][k][j][a] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+
+// Grad1Z
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void ContractZ3d(const int D1D, const int Q1D,
+                 real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+                 const real_t (*B)[MQ1],
+                 regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+                 const int vd, const int d)
+{
+   const int tz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(k,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(j,x,Q1D)
+      {
+         // MFEM_UNROLL(D1D)
+         for (int a=0; a<D1D; ++a) { r_q[k][j][a][vd][d] = smem[tz][vd][d][k][j][a]; }
+         // MFEM_UNROLL(Q1D)
+         for (int i=0; i<Q1D; ++i)
+         {
+            double u = 0.0;
+            // MFEM_UNROLL(D1D)
+            for (int a=0; a<D1D; ++a) { u += B[a][i] * r_q[k][j][a][vd][d]; }
+            smem[tz][vd][d][k][j][i] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+
+   // Flush
+   MFEM_FOREACH_THREAD(j,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(i,x,Q1D)
+      {
+         MFEM_UNROLL(Q1D)
+         for (int k = 0; k < Q1D; ++k) { r_q[k][j][i][vd][d] = 0.0; }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+/// 3D vector gradient, with component
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void Grad3d(const int D1D, const int Q1D,
+            real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+            const real_t (*B)[MQ1],
+            const real_t (*G)[MQ1],
+            regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+            const DeviceTensor<5, const real_t> &XE,
+            const int e)
+{
+   for (int vd = 0; vd < VDIM; vd++)
+   {
+      for (int d = 0; d < DIM; d++)
+      {
+         const real_t (*Bx)[MQ1] = (d == 0) ? G : B;
+         const real_t (*By)[MQ1] = (d == 1) ? G : B;
+         const real_t (*Bz)[MQ1] = (d == 2) ? G : B;
+         ContractX3d(D1D, Q1D, smem, Bx,  XE, vd, d, e);
+         ContractY3d(D1D, Q1D, smem, By, r_q, vd, d);
+         ContractZ3d(D1D, Q1D, smem, Bz, r_q, vd, d);
+      }
+   }
+}
+
+// GradZT
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void ContractZT3d(const int D1D, const int Q1D,
+                  real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+                  const real_t (*B)[MQ1],
+                  regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+                  const int vd, const int d)
+{
+   const int tz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(j,y,Q1D)
+   {
+      MFEM_FOREACH_THREAD(i,x,Q1D)
+      {
+         MFEM_UNROLL(D1D)
+         for (int c=0; c<D1D; ++c)
+         {
+            double u = 0.0;
+            MFEM_UNROLL(Q1D)
+            for (int k=0; k<Q1D; ++k) { u += B[c][k] * r_q[k][j][i][vd][d]; }
+            smem[tz][vd][d][c][j][i] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+// GradYT
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE void
+ContractYT3d(const int D1D, const int Q1D,
+             real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+             const real_t (*B)[MQ1],
+             regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+             const int vd, const int d)
+{
+   const int tz = MFEM_THREAD_ID(z);
+
+   MFEM_FOREACH_THREAD(c,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(i,x,Q1D)
+      {
+         MFEM_UNROLL(Q1D)
+         for (int j=0; j<Q1D; ++j) { r_q[c][i][j][vd][d] = smem[tz][vd][d][c][j][i]; }
+         MFEM_UNROLL(D1D)
+         for (int b=0; b<D1D; ++b)
+         {
+            double u = 0.0;
+            MFEM_UNROLL(Q1D)
+            for (int j=0; j<Q1D; ++j) { u += B[b][j] * r_q[c][i][j][vd][d]; }
+            smem[tz][vd][d][c][b][i] = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+// GradXT
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE void
+ContractXT3d(const int D1D, const int Q1D,
+             real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+             const real_t (*B)[MQ1],
+             regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+             const DeviceTensor<5, real_t> &YE,
+             const int vd, const int d, const int e)
+{
+   const int tz = MFEM_THREAD_ID(z);
+   MFEM_FOREACH_THREAD(c,y,D1D)
+   {
+      MFEM_FOREACH_THREAD(b,x,D1D)
+      {
+         MFEM_UNROLL(Q1D)
+         for (int i=0; i<Q1D; ++i) { r_q[c][b][i][vd][d] = smem[tz][vd][d][c][b][i]; }
+         MFEM_UNROLL(D1D)
+         for (int a=0; a<D1D; ++a)
+         {
+            double u = 0.0;
+            MFEM_UNROLL(Q1D)
+            for (int i=0; i<Q1D; ++i) { u += B[a][i] * r_q[c][b][i][vd][d]; }
+            // smem[tz][vd][d][c][b][a] = u;
+            YE(c,b,a,vd,e) = u;
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+/// 3D vector transposed gradient
+template <int VDIM, int DIM, int MQ1, int NBZ>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void GradTranspose3d(const int d1d, const int q1d,
+                     real_t (&smem)[NBZ][VDIM][DIM][MQ1][MQ1][MQ1],
+                     const real_t (*B)[MQ1],
+                     const real_t (*G)[MQ1],
+                     regs3d_vd_t<VDIM,DIM,MQ1> &r_q,
+                     const DeviceTensor<5, real_t> &YE,
+                     const int e)
+{
+   for (int vd = 0; vd < VDIM; vd++)
+   {
+      for (int d = 0; d < DIM; d++)
+      {
+         const real_t (*Bx)[MQ1] = (d == 0) ? G : B;
+         const real_t (*By)[MQ1] = (d == 1) ? G : B;
+         const real_t (*Bz)[MQ1] = (d == 2) ? G : B;
+         ContractZT3d(d1d, q1d, smem, Bz, r_q, vd, d);
+         ContractYT3d(d1d, q1d, smem, By, r_q, vd, d);
+         ContractXT3d(d1d, q1d, smem, Bx, r_q, YE, vd, d, e);
+      }
+   }
+}
+
+} // namespace mfem::kernels::internal::fma

tests/benchmarks/kernels_vd.hpp

@@ -0,0 +1,273 @@
+// 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
+
+#include "../../config/config.hpp"
+#include "../../config/tconfig.hpp" // MFEM_ALWAYS_INLINE
+#include "../../linalg/dtensor.hpp"
+#include "../../linalg/tensor.hpp"
+#include "../../general/forall.hpp" // MFEM_UNROLL
+// #include "../../fem/kernels.hpp"
+
+namespace mfem::kernels::internal::vd
+{
+
+#if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) ||       \
+     (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
+
+template <int VDIM, int DIM, int N>
+using regs3d_vd_t = mfem::future::tensor<real_t, N, 0, 0, VDIM, DIM>;
+
+#else
+
+template <int VDIM, int DIM, int N>
+using regs3d_vd_t = mfem::future::tensor<real_t, N, N, N, VDIM, DIM>;
+
+template <int N>
+using regs3d_t = mfem::future::tensor<real_t, N, N, N>;
+
+#endif // CUDA/HIP && DEVICE_COMPILE
+
+/// Load 3D input VDIM*DIM vector into given register tensor
+template <int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void LoadDofs3d(const int e, const int d1d,
+                const DeviceTensor<5, const real_t> &X,
+                regs3d_vd_t<VDIM, DIM, MQ1> &Y)
+{
+   for (int dz = 0; dz < d1d; ++dz)
+   {
+      MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
+      {
+         MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
+         {
+            for (int c = 0; c < VDIM; ++c)
+            {
+               for (int d = 0; d < DIM; d++)
+               {
+                  Y[dz][dy][dx][c][d] = X(dx, dy, dz, c, e);
+               }
+            }
+         }
+      }
+   }
+}
+
+/// Write 3D scalar into given device tensor, with read (i) write (j) indices
+template <int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void WriteDofs3d(const int e, const int d1d,
+                 const int i, const int j,
+                 regs3d_vd_t<VDIM, DIM, MQ1> &X,
+                 const DeviceTensor<5, real_t> &Y)
+{
+   for (int dz = 0; dz < d1d; ++dz)
+   {
+      MFEM_FOREACH_THREAD_DIRECT(dy, y, d1d)
+      {
+         MFEM_FOREACH_THREAD_DIRECT(dx, x, d1d)
+         {
+            real_t value = 0.0;
+            for (int d = 0; d < DIM; d++) { value += X[dz][dy][dx][i][d]; }
+            Y(dx, dy, dz, j, e) += value;
+         }
+      }
+   }
+}
+
+/// Write 3D VDIM*DIM vector into given device tensor
+template <int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void WriteDofs3d(const int e, const int d1d,
+                 regs3d_vd_t<VDIM, DIM, MQ1> &X,
+                 const DeviceTensor<5, real_t> &Y)
+{
+   for (int c = 0; c < VDIM; ++c) { WriteDofs3d(e, d1d, c, c, X, Y); }
+}
+
+/// 3D vector contraction, X direction
+template <bool Transpose, int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void ContractX3d(const int d1d, const int q1d,
+                 real_t (&smem)[MQ1][MQ1],
+                 const real_t (*B)[MQ1],
+                 const regs3d_vd_t<VDIM,DIM,MQ1> &X,
+                 regs3d_vd_t<VDIM,DIM,MQ1> &Y,
+                 const int c, const int d)
+{
+   for (int z = 0; z < d1d; ++z)
+   {
+      MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
+      {
+         MFEM_UNROLL(MQ1)
+         MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? q1d : d1d))
+         {
+            smem[y][x] = X[z][y][x][c][d];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD_DIRECT(y, y, d1d)
+      {
+         MFEM_UNROLL(MQ1)
+         MFEM_FOREACH_THREAD_DIRECT(x, x, (Transpose ? d1d : q1d))
+         {
+            real_t u = 0.0;
+            for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
+            {
+               u += (Transpose ? B[x][k] : B[k][x]) * smem[y][k];
+            }
+            Y[z][y][x][c][d] = u;
+         }
+      }
+      MFEM_SYNC_THREAD;
+   }
+}
+
+/// 3D vector contraction, Y direction
+template <bool Transpose, int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void ContractY3d(const int d1d, const int q1d,
+                 real_t (&smem)[MQ1][MQ1],
+                 const real_t (*B)[MQ1],
+                 const regs3d_vd_t<VDIM,DIM,MQ1> &X,
+                 regs3d_vd_t<VDIM,DIM,MQ1> &Y,
+                 const int c, const int d)
+{
+   for (int z = 0; z < d1d; ++z)
+   {
+      MFEM_UNROLL(MQ1)
+      MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? q1d : d1d))
+      {
+         MFEM_FOREACH_THREAD_DIRECT(x, x, q1d) { smem[y][x] = X[z][y][x][c][d]; }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_UNROLL(MQ1)
+      MFEM_FOREACH_THREAD_DIRECT(y, y, (Transpose ? d1d : q1d))
+      {
+         MFEM_UNROLL(MQ1)
+         MFEM_FOREACH_THREAD_DIRECT(x, x, q1d)
+         {
+            real_t u = 0.0;
+            for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
+            {
+               u += (Transpose ? B[y][k] : B[k][y]) * smem[k][x];
+            }
+            Y[z][y][x][c][d] = u;
+         }
+      }
+      MFEM_SYNC_THREAD;
+   }
+}
+
+/// 3D vector contraction, Z direction
+template <bool Transpose, int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void ContractZ3d(const int d1d, const int q1d,
+                 const real_t (*B)[MQ1],
+                 const regs3d_vd_t<VDIM,DIM,MQ1> &X,
+                 regs3d_vd_t<VDIM,DIM,MQ1> &Y,
+                 const int c, const int d)
+{
+   for (int z = 0; z < (Transpose ? d1d : q1d); ++z)
+   {
+      MFEM_UNROLL(MQ1)
+      MFEM_FOREACH_THREAD_DIRECT(y, y, q1d)
+      {
+         MFEM_UNROLL(MQ1)
+         MFEM_FOREACH_THREAD_DIRECT(x, x, q1d)
+         {
+            real_t u = 0.0;
+            for (int k = 0; k < (Transpose ? q1d : d1d); ++k)
+            {
+               u += (Transpose ? B[z][k] : B[k][z]) * X[k][y][x][c][d];
+            }
+            Y[z][y][x][c][d] = u;
+         }
+      }
+   }
+}
+
+/// 3D scalar contraction: X, Y & Z directions
+template <bool Transpose, int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void Contract3d(const int d1d, const int q1d,
+                real_t (&smem)[MQ1][MQ1],
+                const real_t (*Bx)[MQ1],
+                const real_t (*By)[MQ1],
+                const real_t (*Bz)[MQ1],
+                regs3d_vd_t<VDIM,DIM,MQ1> &X,
+                regs3d_vd_t<VDIM,DIM,MQ1> &Y,
+                const int c, const int d)
+{
+   if (!Transpose)
+   {
+      ContractX3d<false>(d1d, q1d, smem, Bx, X, Y, c, d);
+      ContractY3d<false>(d1d, q1d, smem, By, Y, X, c, d);
+      ContractZ3d<false>(d1d, q1d,       Bz, X, Y, c, d);
+   }
+   else
+   {
+      ContractZ3d<true>(d1d, q1d,       Bz, X, Y, c, d);
+      ContractY3d<true>(d1d, q1d, smem, By, Y, X, c, d);
+      ContractX3d<true>(d1d, q1d, smem, Bx, X, Y, c, d);
+   }
+}
+/// 3D vector gradient, with component
+template <int VDIM, int DIM, int MQ1, bool Transpose = false>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void Grad3d(const int d1d, const int q1d,
+            real_t (&smem)[MQ1][MQ1],
+            const real_t (*B)[MQ1],
+            const real_t (*G)[MQ1],
+            regs3d_vd_t<VDIM, DIM, MQ1> &X,
+            regs3d_vd_t<VDIM, DIM, MQ1> &Y,
+            const int c)
+{
+   for (int d = 0; d < DIM; d++)
+   {
+      const real_t (*Bx)[MQ1] = (d == 0) ? G : B;
+      const real_t (*By)[MQ1] = (d == 1) ? G : B;
+      const real_t (*Bz)[MQ1] = (d == 2) ? G : B;
+      Contract3d<Transpose>(d1d, q1d, smem, Bx, By, Bz, X, Y, c, d);
+   }
+}
+/// 3D vector gradient
+template <int VDIM, int DIM, int MQ1, bool Transpose = false>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void Grad3d(const int d1d, const int q1d,
+            real_t (&smem)[MQ1][MQ1],
+            const real_t (*B)[MQ1],
+            const real_t (*G)[MQ1],
+            regs3d_vd_t<VDIM, DIM, MQ1> &X,
+            regs3d_vd_t<VDIM, DIM, MQ1> &Y)
+{
+   for (int c = 0; c < VDIM; c++)
+   {
+      Grad3d<VDIM, DIM, MQ1, Transpose>(d1d, q1d, smem, B, G, X, Y, c);
+   }
+}
+
+/// 3D vector transposed gradient
+template <int VDIM, int DIM, int MQ1>
+inline MFEM_ALWAYS_INLINE MFEM_HOST_DEVICE
+void GradTranspose3d(const int d1d, const int q1d,
+                     real_t (&smem)[MQ1][MQ1],
+                     const real_t (*B)[MQ1],
+                     const real_t (*G)[MQ1],
+                     regs3d_vd_t<VDIM, DIM, MQ1> &X,
+                     regs3d_vd_t<VDIM, DIM, MQ1> &Y)
+{
+   Grad3d<VDIM, DIM, MQ1, true>(d1d, q1d, smem, B, G, X, Y);
+}
+
+} // namespace mfem::kernels::internal

tests/benchmarks/makefile

@@ -20,8 +20,8 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_TESTS = bench_assembly_levels bench_ceed bench_dg_amr bench_elasticity \
-            bench_tmop bench_vector bench_virtuals
+SEQ_TESTS = bench_assembly_levels bench_ceed bench_dfem bench_dg_amr \
+			bench_elasticity bench_tmop bench_vector bench_virtuals
 PAR_TESTS = 
 ifeq ($(MFEM_USE_MPI),NO)
    TESTS = $(SEQ_TESTS)

tests/unit/CMakeLists.txt

@@ -18,6 +18,7 @@ include_directories(BEFORE ${CMAKE_CURRENT_SOURCE_DIR})
 #    for d in general linalg mesh fem enzyme; do ls -1 $d/*.cpp; done
 set(UNIT_TESTS_SRCS
   dfem/test_diffusion.cpp
+  dfem/test_diffusion_q1d.cpp
   dfem/test_divergence.cpp
   dfem/test_mass.cpp
   general/test_array.cpp

tests/unit/dfem/test_diffusion.cpp

@@ -11,8 +11,17 @@
 
 #include "../unit_tests.hpp"
 #include "mfem.hpp"
+#include <fem/dfem/doperator.hpp>
 #include <utility>
 
+#if defined(__has_include) && __has_include("general/nvtx.hpp") && !defined(_WIN32)
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kMagenta
+#include "general/nvtx.hpp"
+#else
+#define dbg(...)
+#endif
+
 #ifdef MFEM_USE_MPI
 
 using namespace mfem;
@@ -48,6 +57,18 @@ template <int DIM> struct Diffusion
       }
    };
 
+   struct MFApplyNoRho
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const dvecd_t &dudxi,
+                                              const matd_t &J,
+                                              const real_t &w) const
+      {
+         const auto invJ = inv(J), TinJ = transpose(invJ);
+         return tuple{ (dudxi * invJ) * TinJ * det(J) * w };
+      }
+   };
+
+
    struct PASetup
    {
       MFEM_HOST_DEVICE inline auto operator()(const real_t u,
@@ -120,7 +141,7 @@ void DFemDiffusion(const char *filename, int p, const int r)
    blf_fa.Assemble();
    blf_fa.Finalize();
 
-   SECTION("[Partial assembly] Diffusion")
+   /*SECTION("[Partial assembly] Diffusion")
    {
       ParBilinearForm blf_pa(&pfes);
       blf_pa.AddDomainIntegrator(new DiffusionIntegrator(rho_coeff, ir));
@@ -132,7 +153,7 @@ void DFemDiffusion(const char *filename, int p, const int r)
       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);
@@ -144,7 +165,7 @@ void DFemDiffusion(const char *filename, int p, const int r)
    static constexpr int U = 0, Coords = 1, Rho = 3;
    const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
 
-   SECTION("[dFEM Matrix free] Diffusion")
+   /*SECTION("[dFEM Matrix free] Diffusion")
    {
       DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
       typename Diffusion<DIM>::MFApply mf_apply_qf;
@@ -169,9 +190,9 @@ void DFemDiffusion(const char *filename, int p, const int r)
 
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
-   }
+   }*/
 
-   SECTION("[dFEM Partial assembly] Diffusion")
+   /*SECTION("[dFEM Partial assembly] Diffusion")
    {
       static constexpr int QData = 2;
       UniformParameterSpace qd_ps(pmesh, *ir, DIM * DIM);
@@ -210,9 +231,9 @@ void DFemDiffusion(const char *filename, int p, const int r)
 
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
-   }
+   }*/
 
-   SECTION("[dFEM Linearization] Diffusion")
+   /*SECTION("[dFEM Linearization] Diffusion")
    {
       DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
       typename Diffusion<DIM>::MFApply mf_apply_qf;
@@ -232,6 +253,80 @@ void DFemDiffusion(const char *filename, int p, const int r)
       pfes.GetProlongationMatrix()->MultTranspose(y, Y);
       Y -= Z;
 
+      real_t norm_global = 0.0;
+      real_t norm_local = Y.Normlinf();
+      MPI_Allreduce(&norm_local, &norm_global, 1, MPI_DOUBLE, MPI_MAX,
+                    pmesh.GetComm());
+
+      REQUIRE(norm_global == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }*/
+
+   /*SECTION("[dFEM Linearization] Diffusion with automatic PA")
+   {
+      DOperator dop_auto_pa(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
+      dbg("Setting up automatic PA for diffusion");
+      dop_auto_pa.UseAutomaticPA();
+
+      typename Diffusion<DIM>::MFApply mf_apply_qf;
+      auto derivatives = std::integer_sequence<size_t, U> {};
+      dbg("Adding domain integrator for automatic PA");
+      dop_auto_pa.AddDomainIntegrator(mf_apply_qf,
+                                      tuple{ Gradient<U>{}, Identity<Rho>{},
+                                             Gradient<Coords>{},
+                                             Weight{} },
+                                      tuple{ Gradient<U>{} },
+                                      *ir, all_domain_attr, derivatives);
+      dop_auto_pa.SetParameters({ &rho_coeff_cv, nodes });
+      dbg("Getting derivative for automatic PA");
+      auto dRdU = dop_auto_pa.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
+
+      pfes.GetRestrictionMatrix()->Mult(x, X);
+      dbg("Mult with automatic PA operator");
+      dRdU->Mult(X, Z);
+
+      blf_fa.Mult(x, y);
+      pfes.GetProlongationMatrix()->MultTranspose(y, Y);
+      Y -= Z;
+
+      real_t norm_global = 0.0;
+      real_t norm_local = Y.Normlinf();
+      MPI_Allreduce(&norm_local, &norm_global, 1, MPI_DOUBLE, MPI_MAX,
+                    pmesh.GetComm());
+
+      REQUIRE(norm_global == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }*/
+
+   SECTION("[dFEM Linearization] Diffusion with automatic PA and new kernels")
+   {
+      DOperator dop_auto_pa(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
+      dbg("Setting up AUTO PA for diffusion");
+      dop_auto_pa.UseNewKernels();
+      dop_auto_pa.UseAutomaticPA();
+      dop_auto_pa.UseKernelsSpecialization();
+
+      typename Diffusion<DIM>::MFApply mf_apply_qf;
+      auto derivatives = std::integer_sequence<size_t, U> {};
+      dbg("Adding domain integrator for automatic PA");
+      dop_auto_pa.AddDomainIntegrator(mf_apply_qf,
+                                      tuple{ Gradient<U>{}, Identity<Rho>{},
+                                             Gradient<Coords>{},
+                                             Weight{} },
+                                      tuple{ Gradient<U>{} },
+                                      *ir, all_domain_attr, derivatives);
+      dop_auto_pa.SetParameters({ &rho_coeff_cv, nodes });
+      dbg("Getting derivative for automatic PA");
+      auto dRdU = dop_auto_pa.GetDerivative(U, {&x}, {&rho_coeff_cv, nodes});
+
+      pfes.GetRestrictionMatrix()->Mult(x, X);
+      dbg("Mult with automatic PA operator");
+      dRdU->Mult(X, Z);
+
+      blf_fa.Mult(x, y);
+      pfes.GetProlongationMatrix()->MultTranspose(y, Y);
+      Y -= Z;
+
       real_t norm_global = 0.0;
       real_t norm_local = Y.Normlinf();
       MPI_Allreduce(&norm_local, &norm_global, 1, MPI_DOUBLE, MPI_MAX,
@@ -241,7 +336,8 @@ void DFemDiffusion(const char *filename, int p, const int r)
       MPI_Barrier(MPI_COMM_WORLD);
    }
 
-   SECTION("[dFEM Matrix free] vector diffusion")
+
+   /*SECTION("[dFEM Matrix free] vector diffusion")
    {
       ParFiniteElementSpace vpfes(&pmesh, &fec, DIM);
       ParGridFunction vx(&vpfes), vy(&vpfes);
@@ -255,8 +351,8 @@ void DFemDiffusion(const char *filename, int p, const int r)
          DOperator dop_mf(vsol, {{Coords, mfes}}, pmesh);
          const auto mf_vector_diffusion_qf =
             [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM, DIM> &dudxi,
-               const tensor<real_t, DIM, DIM> &J,
-               const real_t &w)
+                                 const tensor<real_t, DIM, DIM> &J,
+                                 const real_t &w)
          {
             const auto invJ = inv(J), TinJ = transpose(invJ);
             return tuple{ (dudxi * invJ) * TinJ * det(J) * w };
@@ -281,17 +377,70 @@ void DFemDiffusion(const char *filename, int p, const int r)
                     pmesh.GetComm());
       REQUIRE(norm_global == MFEM_Approx(0.0));
       MPI_Barrier(MPI_COMM_WORLD);
-   }
+   }*/
+
+   /*SECTION("[dFEM Linearization] vector diffusion")
+   {
+      ParFiniteElementSpace vpfes(&pmesh, &fec, DIM);
+      ParGridFunction vx(&vpfes), vy(&vpfes);
+      Vector vX(vpfes.GetTrueVSize()), vY(vpfes.GetTrueVSize()),
+             vZ(vpfes.GetTrueVSize());
+
+      vX.Randomize(1), vx.SetFromTrueDofs(vX);
+
+      {
+         const auto vsol = std::vector{ FieldDescriptor{ U, &vpfes } };
+         DOperator dop_mf(vsol, {{Coords, mfes}}, pmesh);
+         dbg("Setting up automatic PA for diffusion");
+         // dop_mf.UseAutomaticPA(); // ❌ Trying to access out of boundary
+         const auto mf_vector_diffusion_qf =
+            [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM, DIM> &dudxi,
+                                 const tensor<real_t, DIM, DIM> &J,
+                                 const real_t &w)
+         {
+            const auto invJ = inv(J), TinJ = transpose(invJ);
+            return tuple{ (dudxi * invJ) * TinJ * det(J) * w };
+         };
+         auto derivatives = std::integer_sequence<size_t, U> {};
+         dbg("Adding domain integrator for automatic PA");
+         dop_mf.AddDomainIntegrator(mf_vector_diffusion_qf,
+                                    tuple{ Gradient<U>{}, Gradient<Coords>{}, Weight{} },
+                                    tuple{ Gradient<U>{} },
+                                    *ir, all_domain_attr,
+                                    derivatives);
+         dbg("Getting derivative for automatic PA");
+         auto dRdU = dop_mf.GetDerivative(U, {&vx}, {nodes});
+
+         vpfes.GetRestrictionMatrix()->Mult(vx, vX);
+         dbg("Mult with automatic PA operator");
+         dRdU->Mult(vX, vZ);
+      }
+      {
+         ConstantCoefficient one(1.0);
+         ParBilinearForm vblf_fa(&vpfes);
+         vblf_fa.AddDomainIntegrator(new VectorDiffusionIntegrator(one, ir));
+         (vblf_fa.Assemble(), vblf_fa.Finalize());
+         (vblf_fa.Mult(vx, vy), vpfes.GetProlongationMatrix()->MultTranspose(vy, vY));
+      }
+      vY -= vZ;
+      real_t norm_global = 0.0, norm_local = vY.Normlinf();
+      MPI_Allreduce(&norm_local, &norm_global, 1, MPI_DOUBLE, MPI_MAX,
+                    pmesh.GetComm());
+      REQUIRE(norm_global == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }*/
 }
 
-TEST_CASE("DFEM Diffusion", "[Parallel][DFEM]")
+TEST_CASE("DFEM Diffusion", "[Parallel][DFEM][dDiffusionLinearisation]")
 {
-   const bool all_tests = launch_all_non_regression_tests;
+   // const bool all_tests = launch_all_non_regression_tests;
 
-   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
-   const auto r = !all_tests ? 2 : GENERATE(0, 1, 2, 3);
+   // const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+   // const auto r = !all_tests ? 2 : GENERATE(0, 1, 2, 3);
+   const int p = 2, r = 1;
+   dbg("p:{} r:{}", p, r);
 
-   SECTION("2D p=" + std::to_string(p) + " r=" + std::to_string(r))
+   /*SECTION("2D p=" + std::to_string(p) + " r=" + std::to_string(r))
    {
       const auto filename =
          GENERATE("../../data/star.mesh",
@@ -300,16 +449,21 @@ TEST_CASE("DFEM Diffusion", "[Parallel][DFEM]")
                   "../../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))
    {
+#if 0
       const auto filename =
          GENERATE("../../data/fichera.mesh",
                   "../../data/fichera-q3.mesh",
                   "../../data/inline-hex.mesh",
                   "../../data/toroid-hex.mesh",
                   "../../data/periodic-cube.mesh");
+#else
+      const auto filename =
+         GENERATE("../../data/fichera.mesh");
+#endif
       DFemDiffusion<3>(filename, p, r);
    }
 }

tests/unit/dfem/test_diffusion_q1d.cpp

@@ -0,0 +1,341 @@
+// 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 <type_traits>
+#include "fem/dfem/doperator.hpp"
+#include "fem/dfem/util.hpp"
+
+#include <fem/integ/bilininteg_diffusion_kernels.hpp>
+
+#undef NVTX_COLOR
+#define NVTX_COLOR nvtx::kGold
+#include "general/nvtx.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+using mfem::future::dual;
+
+using DOperator = future::DifferentiableOperator;
+
+enum class MQ1Settings : int { kRuntime,
+                               kCompileTime,
+                               kDefault
+                             };
+
+namespace dfem_pa_kernels
+{
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T, int DIM, int T_MQ1 = 0> struct Diffusion
+{
+   using dvecd_t = tensor<T, DIM>;
+   using matd_t = tensor<real_t, DIM, DIM>;
+
+   struct MFApply
+   {
+      static constexpr int MQ1 = T_MQ1;
+      MFEM_HOST_DEVICE inline auto operator()(const dvecd_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::future::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::future::tuple{ inv(J) * transpose(inv(J)) * det(J) * w * rho };
+      }
+   };
+
+   struct PAApply
+   {
+      MFEM_HOST_DEVICE inline auto operator()(const dvecd_t &dudxi,
+                                              const matd_t &q) const
+      {
+         return mfem::future::tuple{ q * dudxi };
+      };
+   };
+};
+
+///////////////////////////////////////////////////////////////////////////////
+template <typename T, int DIM, std::size_t... MQ1s>
+struct MFDiffusionFactory
+{
+   static auto All()
+   {
+      // could also use a map instead of a tuple
+      return mfem::future::make_tuple(typename Diffusion<T, DIM, MQ1s>::MFApply{}...);
+   }
+};
+
+template <typename T, int DIM>
+using MFDiffusionFactory_1_4 = MFDiffusionFactory<T, DIM, 1, 2, 3, 4>;
+
+template <typename T, int DIM>
+class MFDiffusionQFs
+{
+   using MFApplyTuple = decltype(MFDiffusionFactory_1_4<T, DIM>::All());
+   MFApplyTuple mf_qfs;
+
+public:
+   MFDiffusionQFs(): mf_qfs(MFDiffusionFactory_1_4<T, DIM>::All()) {}
+
+   template <typename F>
+   void run(int i, F&& f)
+   {
+      MFEM_VERIFY(i >= 1, "Index must be >= 1");
+      const auto I = static_cast<size_t>(i - 1);
+      runtime_get_impl(I, std::forward<F>(f),
+                       std::make_index_sequence<mfem::future::tuple_size<MFApplyTuple>::value>());
+   }
+
+private:
+   template <typename F, size_t... I>
+   void runtime_get_impl(size_t index, F&& f, std::index_sequence<I...>)
+   {
+      using fun_ptr = std::function<void(F&&)>;
+      fun_ptr table[] = { [&](F&& f) { f(mfem::future::get<I>(mf_qfs)); } ... };
+      if (index < mfem::future::tuple_size<MFApplyTuple>::value)
+      {
+         table[index](std::forward<F>(f));
+      }
+      else
+      {
+         throw std::out_of_range("Index out of bounds");
+      }
+   }
+};
+
+///////////////////////////////////////////////////////////////////////////////
+template <int DIM>
+void DFemDiffusion(const char *filename, int p, const int r,
+                   const MQ1Settings mq1_setting)
+{
+   dbg("DIM:{}", DIM);
+   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.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.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");
+
+   ParGridFunction x(&pfes), y(&pfes), z(&pfes);
+   Vector X(pfes.GetTrueVSize()), Y(pfes.GetTrueVSize()), Z(pfes.GetTrueVSize());
+
+   X.Randomize(1);
+   x.SetFromTrueDofs(X);
+
+   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();
+
+   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, "");
+
+   UniformParameterSpace rho_ps(pmesh, *ir, 1);
+
+   static constexpr int U = 0, Coords = 1, Rho = 3;
+   const auto sol = std::vector{ FieldDescriptor{ U, &pfes } };
+
+   SECTION("DFEM Matrix free")
+   {
+      // fields = {solutions, parameters}
+      dbg("fields = {{solutions, parameters}} = {{{{U}}, {{Rho, Coords}}}}");
+      DOperator dop_mf(sol, {{Rho, &rho_ps}, {Coords, mfes}}, pmesh);
+
+      dbg("AddDomainIntegrator: {{∇U, Rho, ∇Coords, Weight}} -> {{∇U}}");
+      if (mq1_setting == MQ1Settings::kRuntime)
+      {
+         dbg("MQ1Settings::kRuntime");
+         MFEM_VERIFY(q1d == (int)floor(std::pow(ir->GetNPoints(), 1.0/DIM) + 0.5),
+                     "q1d and ir->GetNPoints() have to match");
+         auto add_domain_integrator = [&](auto &qf)
+         {
+            dbg("q1d:{} MQ1:{}", q1d, qf.MQ1);
+            MFEM_VERIFY(q1d == qf.MQ1, "q1d and qf.MQ1 have to match");
+            dop_mf.AddDomainIntegrator(qf,
+                                       tuple{ Gradient<U>{}, Identity<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                       tuple{ Gradient<U>{} }, *ir,
+                                       all_domain_attr);
+         };
+         // select the right qf from the factory
+         MFDiffusionQFs<real_t, DIM> {}.run(q1d, add_domain_integrator);
+      }
+      else if (mq1_setting == MQ1Settings::kCompileTime) // hardcoded, MQ1 = 2,3,4,5
+      {
+         dbg("MQ1Settings::kCompileTime");
+         dbg("q1d:{}", q1d);
+         if (q1d == 2)
+         {
+            typename Diffusion<real_t, DIM, 2>::MFApply mf_apply_qf;
+            MFEM_VERIFY(q1d == 2, "q1d and 2 have to match");
+            dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                       tuple{ Gradient<U>{}, Identity<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                       tuple{ Gradient<U>{} }, *ir,
+                                       all_domain_attr);
+         }
+         else if (q1d == 3)
+         {
+            typename Diffusion<real_t, DIM, 3>::MFApply mf_apply_qf;
+            MFEM_VERIFY(q1d == 3, "q1d and 3 have to match");
+            dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                       tuple{ Gradient<U>{}, Identity<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                       tuple{ Gradient<U>{} }, *ir,
+                                       all_domain_attr);
+         }
+         else if (q1d == 4)
+         {
+            typename Diffusion<real_t, DIM, 4>::MFApply mf_apply_qf;
+            MFEM_VERIFY(q1d == 4, "q1d and 4 have to match");
+            dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                       tuple{ Gradient<U>{}, Identity<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                       tuple{ Gradient<U>{} }, *ir,
+                                       all_domain_attr);
+         }
+         else if (q1d == 5)
+         {
+            typename Diffusion<real_t, DIM, 5>::MFApply mf_apply_qf;
+            MFEM_VERIFY(q1d == 5, "q1d and 5 have to match");
+            dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                       tuple{ Gradient<U>{}, Identity<Rho>{},
+                                              Gradient<Coords>{}, Weight{} },
+                                       tuple{ Gradient<U>{} }, *ir,
+                                       all_domain_attr);
+         }
+         else { MFEM_ABORT("Not supported q1d:" << q1d); }
+      }
+      else // MQ1Settings::kDefault, MQ1 = 0
+      {
+         dbg("MQ1Settings::kDefault");
+         typename Diffusion<real_t, DIM>::MFApply mf_apply_qf;
+         dop_mf.AddDomainIntegrator(mf_apply_qf,
+                                    tuple{ Gradient<U>{}, Identity<Rho>{},
+                                           Gradient<Coords>{}, Weight{} },
+                                    tuple{ Gradient<U>{} }, *ir,
+                                    all_domain_attr);
+      }
+      dop_mf.SetParameters({ &rho_coeff_cv, nodes });
+
+      pfes.GetRestrictionMatrix()->Mult(x, X);
+      dop_mf.Mult(X, Z);
+
+      blf_fa.Mult(x, y);
+      pfes.GetProlongationMatrix()->MultTranspose(y, Y);
+      Y -= Z;
+
+      real_t norm_global = 0.0;
+      real_t norm_local = Y.Normlinf();
+      MPI_Allreduce(&norm_local, &norm_global, 1, MPI_DOUBLE, MPI_MAX,
+                    pmesh.GetComm());
+
+      REQUIRE(norm_global == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+}
+
+TEST_CASE("DFEM Diffusion Q1D", "[Parallel][DFEM][MQ1]")
+{
+   // const bool all_tests = launch_all_non_regression_tests;
+
+   // const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+   // const auto r = !all_tests ? 1 : GENERATE(0, 1, 2, 3);
+   const int p = 2, r = 1;
+   dbg("p:{} r:{}", p, r);
+
+   const auto mq1_setting = MQ1Settings::kCompileTime;
+   /*const auto mq1_setting = GENERATE(MQ1Settings::kRuntime,
+                                     MQ1Settings::kCompileTime,
+                                     MQ1Settings::kDefault);*/
+
+   DiffusionIntegrator::AddSpecialization<3,3,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))
+   {
+#if 0
+      const auto filename =
+         GENERATE("../../data/fichera.mesh",
+                  "../../data/fichera-q3.mesh",
+                  "../../data/inline-hex.mesh",
+                  "../../data/toroid-hex.mesh",
+                  "../../data/periodic-cube.mesh");
+#else
+      const auto filename = "../../data/fichera.mesh";
+#endif
+      dbg("DFemDiffusion");
+      DFemDiffusion<3>(filename, p, r, mq1_setting);
+   }
+}
+
+} // namespace dfem_pa_kernels
+
+#endif

tests/unit/dfem/test_divergence.cpp

@@ -11,6 +11,7 @@
 
 #include "../unit_tests.hpp"
 #include "mfem.hpp"
+#include <fem/dfem/doperator.hpp>
 
 #ifdef MFEM_USE_MPI
 

tests/unit/dfem/test_mass.cpp

@@ -11,6 +11,7 @@
 
 #include "../unit_tests.hpp"
 #include "mfem.hpp"
+#include <fem/dfem/doperator.hpp>
 
 #ifdef MFEM_USE_MPI
 

tests/unit/fem/test_pa_kernels.cpp

@@ -844,9 +844,9 @@ TEST_CASE("Dispatch Map Specializations")
 
    DiffusionIntegrator{};
    REQUIRE_FALSE(
-      DiffusionIntegrator::ApplyPAKernels::GetDispatchTable().empty());
+      DiffusionIntegrator::DiffusionApplyPAKernel::GetDispatchTable().empty());
    REQUIRE_FALSE(
-      DiffusionIntegrator::DiagonalPAKernels::GetDispatchTable().empty());
+      DiffusionIntegrator::DiffusionDiagonalPAKernel::GetDispatchTable().empty());
 
    Mesh mesh = Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL);
    H1_FECollection fec(1, mesh.Dimension());