Merge branch 'master' into dfem-transpose

fem/dfem/doperator.hpp

@@ -84,29 +84,33 @@ public:
       const std::vector<derivative_action_t> &derivative_actions,
       const FieldDescriptor &direction,
       const int &daction_l_size,
-      const std::vector<derivative_action_t> &derivative_actions_transpose,
+      const int &derivative_action_tr_l_size,
+      const std::vector<derivative_action_t> &derivative_tr_actions,
       const FieldDescriptor &transpose_direction,
       const int &daction_transpose_l_size,
       const std::vector<Vector *> &solutions_l,
       const std::vector<Vector *> &parameters_l,
       const restriction_callback_t &restriction_callback,
       const std::function<void(Vector &, Vector &)> &prolongation_transpose,
+      const std::function<void(Vector &, Vector &)> &tr_prolongation_transpose,
       const std::vector<assemble_derivative_sparsematrix_callback_t>
       &assemble_derivative_sparsematrix_callbacks,
-      const std::vector<assemble_derivative_hypreparmatrix_callback_t>
-      &assemble_derivative_hypreparmatrix_callbacks) :
+      const assemble_derivative_hypreparmatrix_callback_t
+      &assemble_derivative_hypreparmatrix_callback) :
       Operator(height, width),
       derivative_actions(derivative_actions),
       direction(direction),
       daction_l(daction_l_size),
       daction_l_size(daction_l_size),
-      derivative_actions_transpose(derivative_actions_transpose),
+      derivative_action_tr_l_size(derivative_action_tr_l_size),
+      derivative_tr_actions(derivative_tr_actions),
       transpose_direction(transpose_direction),
       prolongation_transpose(prolongation_transpose),
+      tr_prolongation_transpose(tr_prolongation_transpose),
       assemble_derivative_sparsematrix_callbacks(
          assemble_derivative_sparsematrix_callbacks),
-      assemble_derivative_hypreparmatrix_callbacks(
-         assemble_derivative_hypreparmatrix_callbacks)
+      assemble_derivative_hypreparmatrix_callback(
+         assemble_derivative_hypreparmatrix_callback)
    {
       std::vector<Vector> s_l(solutions_l.size());
       for (size_t i = 0; i < s_l.size(); i++)
@@ -156,18 +160,18 @@ public:
    /// direction_t on T-dofs.
    void MultTranspose(const Vector &direction_t, Vector &result_t) const override
    {
-      MFEM_ASSERT(!derivative_actions_transpose.empty(),
+      MFEM_ASSERT(!derivative_tr_actions.empty(),
                   "derivative can't be used to be multiplied in transpose mode");
 
-      daction_l.SetSize(width);
+      daction_l.SetSize(derivative_action_tr_l_size);
       daction_l = 0.0;
 
       prolongation(transpose_direction, direction_t, direction_l);
-      for (const auto &f : derivative_actions_transpose)
+      for (const auto &f : derivative_tr_actions)
       {
          f(fields_e, direction_l, daction_l);
       }
-      prolongation_transpose(daction_l, result_t);
+      tr_prolongation_transpose(daction_l, result_t);
    };
 
    /// @brief Assemble the derivative operator into a SparseMatrix.
@@ -183,6 +187,10 @@ public:
       {
          f(fields_e, A);
       }
+
+      // SparseMatrix A is finalized after all callbacks have contributed to
+      // it.
+      A->Finalize();
    }
 
    /// @brief Assemble the derivative operator into a HypreParMatrix.
@@ -191,13 +199,7 @@ public:
    /// be an uninitialized object.
    void Assemble(HypreParMatrix *&A)
    {
-      MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
-                  "derivative can't be assembled into a HypreParMatrix");
-
-      for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
-      {
-         f(fields_e, A);
-      }
+      assemble_derivative_hypreparmatrix_callback(fields_e, A);
    }
 
 private:
@@ -212,10 +214,12 @@ private:
 
    const int daction_l_size;
 
+   const int derivative_action_tr_l_size;
+
    /// Transpose Derivative action callbacks. Depending on the requested
    /// derivatives in DifferentiableOperator the callbacks represent certain
    /// combinations of actions of derivatives of the forward operator.
-   std::vector<derivative_action_t> derivative_actions_transpose;
+   std::vector<derivative_action_t> derivative_tr_actions;
 
    FieldDescriptor transpose_direction;
 
@@ -225,13 +229,15 @@ private:
 
    std::function<void(Vector &, Vector &)> prolongation_transpose;
 
+   std::function<void(Vector &, Vector &)> tr_prolongation_transpose;
+
    /// Callbacks that assemble derivatives into a SparseMatrix.
    std::vector<assemble_derivative_sparsematrix_callback_t>
    assemble_derivative_sparsematrix_callbacks;
 
    /// Callbacks that assemble derivatives into a HypreParMatrix.
-   std::vector<assemble_derivative_hypreparmatrix_callback_t>
-   assemble_derivative_hypreparmatrix_callbacks;
+   assemble_derivative_hypreparmatrix_callback_t
+   assemble_derivative_hypreparmatrix_callback;
 };
 
 /// Class representing a differentiable operator which acts on solution and
@@ -457,7 +463,10 @@ public:
          dir_l = s_l[derivative_idx];
       }
 
-      derivative_setup_callbacks[derivative_id][0](fields_e, dir_l);
+      for (size_t i = 0; i < derivative_setup_callbacks[derivative_id].size(); i++)
+      {
+         derivative_setup_callbacks[derivative_id][i](fields_e, dir_l);
+      }
 
       return std::make_shared<DerivativeOperator>(
                 height,
@@ -465,15 +474,17 @@ public:
                 derivative_action_callbacks[derivative_id],
                 fields[derivative_idx],
                 residual_l.Size(),
-                daction_transpose_callbacks[derivative_id],
+                derivative_action_tr_l_size[derivative_id],
+                derivative_action_tr_callbacks[derivative_id],
                 fields[test_space_field_idx],
                 GetVSize(fields[test_space_field_idx]),
                 sol_l,
                 par_l,
                 restriction_callback,
                 prolongation_transpose,
+                derivative_tr_prolongation_transpose[derivative_id],
                 assemble_derivative_sparsematrix_callbacks[derivative_id],
-                assemble_derivative_hypreparmatrix_callbacks[derivative_id]);
+                assemble_derivative_hypreparmatrix_callback[derivative_id]);
    }
 
 private:
@@ -486,13 +497,15 @@ private:
    std::map<size_t,
        std::vector<derivative_action_t>> derivative_action_callbacks;
    std::map<size_t,
-       std::vector<derivative_action_t>> daction_transpose_callbacks;
+       std::vector<derivative_action_t>> derivative_action_tr_callbacks;
+   std::map<size_t,
+       std::function<void(Vector &, Vector &)>> derivative_tr_prolongation_transpose;
+   std::map<size_t, int> derivative_action_tr_l_size;
    std::map<size_t,
        std::vector<assemble_derivative_sparsematrix_callback_t>>
        assemble_derivative_sparsematrix_callbacks;
-   std::map<size_t,
-       std::vector<assemble_derivative_hypreparmatrix_callback_t>>
-       assemble_derivative_hypreparmatrix_callbacks;
+   std::map<size_t, assemble_derivative_hypreparmatrix_callback_t>
+   assemble_derivative_hypreparmatrix_callback;
 
    std::vector<FieldDescriptor> solutions;
    std::vector<FieldDescriptor> parameters;
@@ -510,7 +523,7 @@ private:
    std::function<void(Vector &, Vector &)> output_restriction_transpose;
    restriction_callback_t restriction_callback;
 
-   std::map<size_t, Vector> derivative_qp_caches;
+   std::map<size_t, std::vector<Vector>> derivative_qp_caches;
 
    std::map<size_t, size_t> assembled_vector_sizes;
 
@@ -769,9 +782,10 @@ void DifferentiableOperator::AddIntegrator(
    auto input_size_on_qp =
       get_input_size_on_qp(inputs, std::make_index_sequence<num_inputs> {});
 
+   // printf("calculate shmem action info\n");
    auto action_shmem_info =
       get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>
-      (input_dtq_maps, output_dtq_maps, fields, num_entities, inputs, num_qp,
+      (input_dtq_maps, output_dtq_maps, fields, num_entities, num_qp,
        input_size_on_qp, residual_size_on_qp, element_dof_ordering);
 
    Vector shmem_cache(action_shmem_info.total_size);
@@ -894,9 +908,10 @@ void DifferentiableOperator::AddIntegrator(
          const int da_size_on_qp =
             GetSizeOnQP<entity_t>(output_fop, fields[test_space_field_idx]);
 
+         // printf("calculate shmem derivative action info\n");
          auto shmem_info =
             get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>(
-               input_dtq_maps, output_dtq_maps, fields, num_entities, inputs,
+               input_dtq_maps, output_dtq_maps, fields, num_entities,
                num_qp, input_size_on_qp, residual_size_on_qp,
                element_dof_ordering, d_field_idx);
 
@@ -968,9 +983,12 @@ void DifferentiableOperator::AddIntegrator(
          // Quadrature point local derivative cache for each element, with data
          // layout:
          // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         derivative_qp_caches[derivative_id].push_back(
+            Vector(test_vdim * test_op_dim * trial_vdim * total_trial_op_dim * num_qp *
+                   num_entities));
+
+         const int cache_index = this->derivative_qp_caches[derivative_id].size() - 1;
+
          // Create local references for MSVC lambda capture compatibility
          auto& fields_ref = this->fields;
          auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
@@ -1012,6 +1030,8 @@ void DifferentiableOperator::AddIntegrator(
                trial_vdim,
                inputs_trial_op_dim,
 
+               qpdc_idx = cache_index,
+
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref
             ](std::vector<Vector> &f_e, const Vector &dir_l) mutable
@@ -1024,7 +1044,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.ReadWrite(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].ReadWrite(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1095,6 +1115,7 @@ void DifferentiableOperator::AddIntegrator(
                inputs_trial_op_dim,
                total_trial_op_dim,
                trial_vdim,
+               qpdc_idx = cache_index,
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref,
                &or_transpose
@@ -1112,7 +1133,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1157,6 +1178,154 @@ void DifferentiableOperator::AddIntegrator(
             or_transpose(derivative_action_e, der_action_l);
          });
 
+         // This prevents Sum/Identity. These are invalid
+         // as input FieldOperators anyways.
+         constexpr auto dummy_fop = Value<0> {};
+
+         auto [input_rt,
+               input_e_sz] = get_restriction_transpose<entity_t>
+                             (fields[d_field_idx],
+                              element_dof_ordering, dummy_fop);
+
+         const auto input_restriction_transpose = input_rt;
+
+         derivative_tr_prolongation_transpose[derivative_id] =
+            get_prolongation_transpose(
+               fields[d_field_idx], dummy_fop, mesh.GetComm());
+
+         const auto d_tr_field_idx = test_space_field_idx;
+         const auto direction_tr = fields[d_tr_field_idx];
+
+         auto output_size_on_qp =
+            get_input_size_on_qp(outputs, std::make_index_sequence<num_outputs> {});
+
+         const int residual_tr_size_on_qp = trial_vdim * total_trial_op_dim;
+
+         auto shmem_tr_info =
+            get_shmem_info<entity_t, num_fields, num_outputs, num_inputs>(
+               output_dtq_maps, input_dtq_maps, fields, num_entities,
+               num_qp, output_size_on_qp, residual_tr_size_on_qp,
+               element_dof_ordering, test_space_field_idx);
+
+         // print_shared_memory_info(shmem_tr_info);
+
+         // TODO: this is a hack to extend the shared memory with a known
+         // offset for a temp variable
+         Vector shmem_tr_cache(shmem_tr_info.total_size + residual_tr_size_on_qp *
+                               num_qp);
+
+         Vector direction_tr_e(get_restriction<entity_t>(
+                                  fields[test_space_field_idx],
+                                  element_dof_ordering)->Height());
+
+         derivative_action_tr_l_size[derivative_id] =
+            get_restriction<entity_t>(fields[d_field_idx],
+                                      element_dof_ordering)->Width();
+
+         Vector derivative_action_tr_e(input_e_sz);
+         derivative_action_tr_e = 0.0;
+
+         derivative_action_tr_callbacks[derivative_id].push_back(
+            [
+               // capture by copy:
+               dimension,             // int
+               num_entities,          // int
+               num_trial_dof,         // int
+               num_qp,                // int
+               q1d,                   // int
+               test_vdim,             // int (= output_fop.vdim)
+               test_op_dim,           // int (derived from output_fop)
+               inputs,                // mfem::future::tuple
+               outputs,               // mfem::future::tuple
+               attributes,            // Array<int>
+               ir_weights,            // DeviceTensor
+               use_sum_factorization, // bool
+               input_dtq_maps,        // std::array<DofToQuadMap, num_fields>
+               output_dtq_maps,       // std::array<DofToQuadMap, num_fields>
+               // output_fop,            // class derived from FieldOperator
+               thread_blocks,         // ThreadBlocks
+               shmem_tr_cache,           // Vector (local)
+               shmem_tr_info,            // SharedMemoryInfo
+               // TODO: make this Array<int> a member of the DifferentiableOperator
+               //       and capture it by ref.
+               elem_attributes,       // Array<int>
+
+               input_is_dependent,
+               direction_tr,             // FieldDescriptor
+               direction_tr_e,           // Vector
+               derivative_action_tr_e,   // Vector
+               element_dof_ordering,  // ElementDofOrdering
+               inputs_trial_op_dim,
+               total_trial_op_dim,
+               trial_vdim,
+               input_restriction_transpose,
+               qpdc_idx = cache_index,
+               // capture by ref:
+               &qpdc_mem = derivative_qp_caches_ref
+            ](
+               std::vector<Vector> &f_e, const Vector &dir_tr_l,
+               Vector &derivative_action_tr_l) mutable
+         {
+            restriction<entity_t>(direction_tr, dir_tr_l, direction_tr_e,
+                                  element_dof_ordering);
+            auto ye = Reshape(derivative_action_tr_e.ReadWrite(), num_trial_dof,
+                              trial_vdim, num_entities);
+            auto wrapped_fields_e = wrap_fields(f_e, shmem_tr_info.field_sizes,
+                                                num_entities);
+            auto wrapped_direction_e = Reshape(direction_tr_e.ReadWrite(),
+                                               shmem_tr_info.direction_size,
+                                               num_entities);
+
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
+                                trial_vdim, total_trial_op_dim, num_qp, num_entities);
+
+            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
+
+            const bool has_attr = attributes.Size() > 0;
+            const auto d_attr = attributes.Read();
+            const auto d_elem_attr = elem_attributes->Read();
+
+            derivative_action_tr_e = 0.0;
+            forall([=] MFEM_HOST_DEVICE (int e, real_t *shmem)
+            {
+               if (has_attr && !d_attr[d_elem_attr[e] - 1]) { return; }
+
+               auto [output_dtq_shmem, input_dtq_shmem, fields_shmem,
+                                       direction_shmem, input_shmem,
+                                       shadow_shmem_, residual_shmem,
+                                       scratch_shmem] =
+                        unpack_shmem(shmem, shmem_tr_info, output_dtq_maps, input_dtq_maps,
+                                     wrapped_fields_e, wrapped_direction_e, num_qp, e);
+               auto &shadow_shmem = shadow_shmem_;
+
+               std::array<bool, num_outputs> all_true{true};
+               map_direction_to_quadrature_data_conditional(
+                  shadow_shmem, direction_shmem, output_dtq_shmem, outputs,
+                  ir_weights, scratch_shmem, all_true, dimension,
+                  use_sum_factorization);
+
+               auto fhat = Reshape(&residual_shmem(0, 0), trial_vdim,
+                                   total_trial_op_dim, num_qp);
+
+               auto qpdce = Reshape(&qpdc(0, 0, 0, 0, 0, e), test_vdim, test_op_dim,
+                                    trial_vdim, total_trial_op_dim, num_qp);
+
+               constexpr bool transpose = true;
+               apply_qpdc(fhat, shadow_shmem, qpdce, itod, q1d, dimension,
+                          use_sum_factorization, transpose);
+
+               auto y = Reshape(&ye(0, 0, e), num_trial_dof, trial_vdim);
+               auto fi_shmem = Reshape(shmem + shmem_tr_info.total_size, trial_vdim,
+                                       total_trial_op_dim, num_qp);
+
+               map_quadrature_data_to_fields_conditional(
+                  y, fhat, inputs, itod, input_dtq_shmem, scratch_shmem, fi_shmem,
+                  input_is_dependent, dimension, use_sum_factorization);
+            }, num_entities, thread_blocks, shmem_tr_info.total_size,
+            shmem_tr_cache.ReadWrite());
+            input_restriction_transpose(derivative_action_tr_e, derivative_action_tr_l);
+         });
+
          assemble_derivative_sparsematrix_callbacks[derivative_id].push_back(
             [
                // capture by copy:
@@ -1190,7 +1359,7 @@ void DifferentiableOperator::AddIntegrator(
                inputs_trial_op_dim,
                Ae_mem,
                output_to_field,
-
+               qpdc_idx = cache_index,
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref,
                &fields = fields_ref
@@ -1202,7 +1371,7 @@ void DifferentiableOperator::AddIntegrator(
                                                shmem_info.direction_size,
                                                num_entities);
 
-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                 trial_vdim, total_trial_op_dim, num_qp, num_entities);
 
             auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1250,7 +1419,10 @@ void DifferentiableOperator::AddIntegrator(
             auto test_fes = *std::get_if<const ParFiniteElementSpace *>
                             (&fields[output_to_field[0]].data);
 
-            A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            if (A == nullptr)
+            {
+               A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            }
 
             auto tmp = Reshape(Ae_mem.HostReadWrite(), num_test_dof * test_vdim,
                                num_trial_dof * trial_vdim, num_entities);
@@ -1321,14 +1493,15 @@ void DifferentiableOperator::AddIntegrator(
                   A->AddSubMatrix(test_vdofs, trial_vdofs, Aee, 1);
                }
             }
-            A->Finalize();
+            // Don't finalize here since multiple callbacks might contribute to the same matrix
+            // A->Finalize() will be called after all callbacks have contributed
          });
 
          // Create local references for MSVC lambda capture compatibility
          auto& assemble_derivative_sparsematrix_callbacks_ref =
             this->assemble_derivative_sparsematrix_callbacks[derivative_id];
 
-         assemble_derivative_hypreparmatrix_callbacks[derivative_id].push_back(
+         assemble_derivative_hypreparmatrix_callback[derivative_id] =
             [
                input_is_dependent,
                input_to_field,
@@ -1342,6 +1515,7 @@ void DifferentiableOperator::AddIntegrator(
             {
                f(f_e, spmat);
             }
+            spmat->Finalize();
 
             if (spmat == nullptr)
             {
@@ -1395,7 +1569,7 @@ void DifferentiableOperator::AddIntegrator(
                        trial_fes->Dof_TrueDof_Matrix());
             }
             delete spmat;
-         });
+         };
       }, derivative_ids);
    }
 }

fem/dfem/integrate.hpp

@@ -530,4 +530,65 @@ void map_quadrature_data_to_fields(
    }
 }
 
+template <size_t N, typename field_operator_ts>
+MFEM_HOST_DEVICE
+void map_quadrature_data_to_fields_conditional(
+   DeviceTensor<2, real_t> &y,
+   const DeviceTensor<3, real_t> &f,
+   const field_operator_ts &fops,
+   const DeviceTensor<1, const real_t> &op_dims,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   std::array<DeviceTensor<1>, 6> &scratch_mem,
+   const DeviceTensor<3> &fi_shmem,
+   const std::array<bool, N> &conditions,
+   const int &dimension,
+   const bool &use_sum_factorization)
+{
+   int offset = 0;
+   for_constexpr<N>([&](auto i)
+   {
+      if (conditions[i])
+      {
+         [[maybe_unused]] const auto [K, unused, M] = f.GetShape();
+         const int L = static_cast<int>(op_dims(static_cast<size_t>(i)));
+         auto fi = Reshape(&fi_shmem(0, 0, 0), K, L, M);
+         for (int k = 0; k < K; k++)
+         {
+            for (int l = 0; l < L; l++)
+            {
+               for (int m = 0; m < M; m++)
+               {
+                  fi(k, l, m) = f(k, l + offset, m);
+               }
+            }
+         }
+
+         if (use_sum_factorization)
+         {
+            if (dimension == 1)
+            {
+               map_quadrature_data_to_fields_tensor_impl_1d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 2)
+            {
+               map_quadrature_data_to_fields_tensor_impl_2d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 3)
+            {
+               map_quadrature_data_to_fields_tensor_impl_3d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else { MFEM_ABORT_KERNEL("dimension not supported"); }
+         }
+         else
+         {
+            map_quadrature_data_to_fields_impl(y, fi, get<i>(fops), dtqmaps[i]);
+         }
+         offset += L;
+      }
+   });
+}
+
 } // namespace mfem::future

fem/dfem/interpolate.hpp

@@ -505,13 +505,13 @@ void map_field_to_quadrature_data(
    }
 }
 
-template <typename field_operator_ts, size_t num_inputs, size_t num_fields>
+template <typename field_operator_ts, size_t N, size_t M>
 MFEM_HOST_DEVICE inline
 void map_fields_to_quadrature_data(
-   std::array<DeviceTensor<2>, num_inputs> &fields_qp,
-   const std::array<DeviceTensor<1>, num_fields> &fields_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
-   const std::array<size_t, num_inputs> &input_to_field,
+   std::array<DeviceTensor<2>, N> &fields_qp,
+   const std::array<DeviceTensor<1>, M> &fields_e,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   const std::array<size_t, N> &input_to_field,
    const field_operator_ts &fops,
    const DeviceTensor<1, const real_t> &integration_weights,
    const std::array<DeviceTensor<1>, 6> &scratch_mem,
@@ -523,7 +523,7 @@ void map_fields_to_quadrature_data(
    // attached to them and we create a dummy field which is not accessed
    // inside the functions it is passed to.
    const auto dummy_field_weight = DeviceTensor<1>(nullptr, 0);
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
    {
       const DeviceTensor<1> &field_e =
          (input_to_field[i] == SIZE_MAX) ? dummy_field_weight :
@@ -549,12 +549,7 @@ void map_fields_to_quadrature_data(
                fields_qp[i], dtqmaps[i], field_e, get<i>(fops),
                integration_weights, scratch_mem);
          }
-         else
-         {
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-            MFEM_ABORT("unsupported dimension");
-#endif
-         }
+         else { MFEM_ABORT_KERNEL("unsupported dimension"); }
       }
       else
       {
@@ -627,20 +622,20 @@ void map_fields_to_quadrature_data_conditional(
    });
 }
 
-template <size_t num_inputs, typename field_operator_ts>
+template <size_t N, typename field_operator_ts>
 MFEM_HOST_DEVICE
 void map_direction_to_quadrature_data_conditional(
-   std::array<DeviceTensor<2>, num_inputs> &directions_qp,
+   std::array<DeviceTensor<2>, N> &directions_qp,
    const DeviceTensor<1> &direction_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
+   const std::array<DofToQuadMap, N> &dtqmaps,
    field_operator_ts fops,
    const DeviceTensor<1, const real_t> &integration_weights,
    const std::array<DeviceTensor<1>, 6> &scratch_mem,
-   const std::array<bool, num_inputs> &conditions,
+   const std::array<bool, N> &conditions,
    const int &dimension,
    const bool &use_sum_factorization)
 {
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
    {
       if (conditions[i])
       {

fem/dfem/qfunction_apply.hpp

@@ -379,54 +379,102 @@ namespace detail
 /// @param shadow_shmem the shadow shared memory.
 /// @param qpdc the quadrature point data cache holding the resulting
 /// Jacobians on each quadrature point.
-/// @param itod inputs trial operator dimension.
-/// If input is dependent the value corresponds to the spatial dimension, otherwise
-/// a zero indicates non-dependence on the variable.
+/// @param op_dims operator dimensions.
+/// If an operator is dependent, the value corresponds to the spatial dimension.
+/// Otherwise a zero indicates indepence on the variable.
 /// @param q the current quadrature point index.
-template <size_t num_fields>
+/// @param transpose switch to use transpose action.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
    DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
    const DeviceTensor<5, const real_t> &qpdc,
-   const DeviceTensor<1, const real_t> &itod,
-   const int &q)
+   const DeviceTensor<1, const real_t> &op_dims,
+   const int &q,
+   bool transpose)
 {
+   const size_t num_ops = op_dims.GetShape()[0];
+
    const int test_vdim = qpdc.GetShape()[0];
    const int test_op_dim = qpdc.GetShape()[1];
    const int trial_vdim = qpdc.GetShape()[2];
-   const int num_qp = qpdc.GetShape()[4];
-   const size_t num_inputs = itod.GetShape()[0];
+   const int total_trial_op_dim = qpdc.GetShape()[3];
 
-   for (int i = 0; i < test_vdim; i++)
+   const int num_qp = qpdc.GetShape()[4];
+
+   if (transpose)
    {
-      for (int k = 0; k < test_op_dim; k++)
+      for (int j = 0; j < trial_vdim; j++)
       {
-         real_t sum = 0.0;
-         int m_offset = 0;
-         for (size_t s = 0; s < num_inputs; s++)
+         for (int m = 0; m < total_trial_op_dim; m++)
          {
-            const int trial_op_dim = static_cast<int>(itod(s));
-            if (trial_op_dim == 0)
-            {
-               continue;
-            }
-            const auto d_qp =
-               Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
-            for (int j = 0; j < trial_vdim; j++)
-            {
-               for (int m = 0; m < trial_op_dim; m++)
-               {
-                  sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
-               }
-            }
-            m_offset += trial_op_dim;
+            fhat(j, m, q) = 0.0;
+         }
+      }
+
+      // Since we don't support more than output space right now
+      // shadow_shmem will always be of size 1.
+      constexpr int shadow_idx_tr = 0;
+      auto d_qp = Reshape(&(shadow_shmem[shadow_idx_tr])[0], test_vdim, test_op_dim,
+                          num_qp);
+
+      int m_offset = 0;
+      for (size_t s = 0; s < num_ops; s++)
+      {
+         const int trial_op_dim = static_cast<int>(op_dims(s));
+         if (trial_op_dim == 0) { continue; }
+
+         for (int j = 0; j < trial_vdim; j++)
+         {
+            for (int m = 0; m < trial_op_dim; m++)
+            {
+               real_t sum = 0.0;
+               for (int i = 0; i < test_vdim; i++)
+               {
+                  for (int k = 0; k < test_op_dim; k++)
+                  {
+                     const real_t contrib = qpdc(i, k, j, m + m_offset, q) * d_qp(i, k, q);
+                     sum += contrib;
+                  }
+               }
+               fhat(j, m + m_offset, q) += sum;
+            }
+         }
+         m_offset += trial_op_dim;
+      }
+   }
+   else
+   {
+      for (int i = 0; i < test_vdim; i++)
+      {
+         for (int k = 0; k < test_op_dim; k++)
+         {
+            real_t sum = 0.0;
+            int m_offset = 0;
+            for (size_t s = 0; s < num_ops; s++)
+            {
+               const int trial_op_dim = static_cast<int>(op_dims(s));
+               if (trial_op_dim == 0) { continue; }
+
+               const auto d_qp =
+                  Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
+               for (int j = 0; j < trial_vdim; j++)
+               {
+                  for (int m = 0; m < trial_op_dim; m++)
+                  {
+                     sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
+                  }
+               }
+               m_offset += trial_op_dim;
+            }
+            fhat(i, k, q) = sum;
          }
-         fhat(i, k, q) = sum;
       }
    }
 }
-}
+
+} // namespace detail
 
 /// @brief Apply the quadrature point data cache (qpdc) to a vector
 /// (usually a direction).
@@ -445,16 +493,18 @@ void apply_qpdc(
 /// @param q1d number of quadrature points in 1D.
 /// @param dimension spatial dimension.
 /// @param use_sum_factorization whether to use sum factorization.
-template <size_t num_fields>
+/// @param T switch to use transpose application.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
    DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
    const DeviceTensor<5, const real_t> &qpdc,
    const DeviceTensor<1, const real_t> &itod,
    const int &q1d,
    const int &dimension,
-   const bool &use_sum_factorization)
+   const bool &use_sum_factorization,
+   const bool T = false)
 {
    if (use_sum_factorization)
    {
@@ -462,7 +512,7 @@ void apply_qpdc(
       {
          MFEM_FOREACH_THREAD_DIRECT(q, x, q1d)
          {
-            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
          }
       }
       else if (dimension == 2)
@@ -472,7 +522,7 @@ void apply_qpdc(
             MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
             {
                const int q = qx + q1d * qy;
-               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
             }
          }
       }
@@ -485,7 +535,7 @@ void apply_qpdc(
                MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
                {
                   const int q = qx + q1d * (qy + q1d * qz);
-                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
                }
             }
          }
@@ -500,7 +550,7 @@ void apply_qpdc(
       const int num_qp = qpdc.GetShape()[4];
       MFEM_FOREACH_THREAD_DIRECT(q, x, num_qp)
       {
-         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
       }
    }
 }

fem/dfem/qfunction_transform.hpp

@@ -243,6 +243,35 @@ void process_qf_arg(
    }
 }
 
+template <typename T, int n>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n> &arg)
+{
+   for (int i = 0; i < n; i++)
+   {
+      arg(i) = u(i);
+   }
+}
+
+template <typename T, int n, int m>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n, m> &arg)
+{
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         arg(j, i) = u((i * n) + j);
+      }
+   }
+}
+
 template <typename arg_type>
 MFEM_HOST_DEVICE inline
 void process_qf_arg(const DeviceTensor<2> &u, arg_type &arg, int qp)
@@ -327,20 +356,4 @@ void process_qf_result(
    }
 }
 
-template <typename T, int n, int m>
-MFEM_HOST_DEVICE inline
-void process_qf_arg(
-   const DeviceTensor<1, T> &u,
-   const DeviceTensor<1, T> &v,
-   tensor<T, n, m> &arg)
-{
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < n; j++)
-      {
-         arg(j, i) = u((i * n) + j);
-      }
-   }
-}
-
 } // namespace mfem::future

fem/dfem/util.hpp

@@ -1202,7 +1202,14 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
    const Operator *P = get_prolongation(f);
    auto PT = [=](const Vector &r_local, Vector &y)
    {
-      P->MultTranspose(r_local, y);
+      if (P)
+      {
+         P->MultTranspose(r_local, y);
+      }
+      else
+      {
+         y = r_local;
+      }
    };
    return PT;
 }
@@ -1580,14 +1587,13 @@ struct SharedMemoryInfo
    std::array<int, 6> temp_sizes;
 };
 
-template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs, typename input_t>
+template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs>
 SharedMemoryInfo<num_fields, num_inputs, num_outputs>
 get_shmem_info(
    const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
    const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
    const std::vector<FieldDescriptor> &fields,
    const int &num_entities,
-   const input_t &inputs,
    const int &num_qp,
    const std::vector<int> &input_size_on_qp,
    const int &residual_size_on_qp,

tests/unit/CMakeLists.txt

@@ -39,6 +39,8 @@ set(UNIT_TESTS_SRCS
   dfem/test_divergence.cpp
   dfem/test_lvector_interface.cpp
   dfem/test_mass.cpp
+  dfem/test_transpose.cpp
+  dfem/test_multiple_integrators.cpp
   general/test_array.cpp
   general/test_scan.cpp
   general/test_arrays_by_name.cpp

tests/unit/dfem/test_mass.cpp

@@ -218,7 +218,7 @@ template <int DIM> void mass_mat_mixed(const char* filename, int p)
 
    SECTION("spmat")
    {
-      SparseMatrix *A;
+      SparseMatrix *A = nullptr;
       ddopdu->Assemble(A);
       TestSameMatrices(*A, blf.SpMat());
       delete A;

tests/unit/dfem/test_multiple_integrators.cpp

@@ -0,0 +1,202 @@
+// 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 "../linalg/test_same_matrices.hpp"
+#include "mfem.hpp"
+#include "fem/dfem/doperator.hpp"
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+#ifdef MFEM_USE_ENZYME
+using dscalar_t = real_t;
+#else
+using mfem::future::dual;
+using dscalar_t = dual<real_t, real_t>;
+#endif
+
+template <int DIM>
+void mult_integ(const char *filename, int p)
+{
+   CAPTURE(filename, DIM, p);
+
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   pmesh.EnsureNodes();
+   auto* nodes = static_cast<ParGridFunction*>(pmesh.GetNodes());
+   p = std::max(p, pmesh.GetNodalFESpace()->GetMaxElementOrder());
+   smesh.Clear();
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace fes(&pmesh, &fec);
+
+   ParGridFunction x(&fes), y(&fes), z(&fes);
+   Vector X(fes.GetTrueVSize()), Y(fes.GetTrueVSize()), Z(fes.GetTrueVSize());
+
+   X.Randomize(1);
+   x.SetFromTrueDofs(X);
+
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * p);
+
+   Array<int> all_domain_attr;
+   if (pmesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   ParBilinearForm blf(&fes);
+   blf.AddDomainIntegrator(new DiffusionIntegrator(ir));
+   blf.AddDomainIntegrator(new MassIntegrator(ir));
+   blf.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   blf.Assemble();
+   blf.Mult(x, y);
+   fes.GetProlongationMatrix()->MultTranspose(y, Y);
+
+   const auto mass_qf =
+      [] MFEM_HOST_DEVICE(
+         const dscalar_t &u,
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{u * w * det(J)};
+   };
+
+   const auto diffusion_qf =
+      [] MFEM_HOST_DEVICE(
+         const tensor<dscalar_t, DIM> &dudxi,
+         const tensor<real_t, DIM, DIM> &J,
+         const real_t &w)
+   {
+      return tuple{(dudxi * inv(J)) * transpose(inv(J)) * w * det(J)};
+   };
+
+   static constexpr int U = 0, Coords = 1;
+   const auto sol = std::vector{ FieldDescriptor{ U, &fes } };
+   DifferentiableOperator dop(sol, {{Coords, nodes->ParFESpace()}}, pmesh);
+
+   auto derivatives = std::integer_sequence<size_t, U> {};
+
+   dop.AddDomainIntegrator(diffusion_qf,
+                           tuple{ Gradient<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Gradient<U>{} },
+                           *ir, all_domain_attr, derivatives);
+
+   dop.AddDomainIntegrator(mass_qf,
+                           tuple{ Value<U>{}, Gradient<Coords>{}, Weight{} },
+                           tuple{ Value<U>{} },
+                           *ir, all_domain_attr, derivatives);
+
+   SECTION("action")
+   {
+      dop.SetParameters({ nodes });
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      dop.Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   SECTION("linearized action")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      ddopdu->Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+   }
+
+   SECTION("linearized assembled SparseMatrix")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      SparseMatrix *A = nullptr;
+      ddopdu->Assemble(A);
+
+      A->Mult(x, z);
+      fes.GetProlongationMatrix()->MultTranspose(z, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+      delete A;
+   }
+
+   SECTION("linearized assembled HypreParMatrix")
+   {
+      auto ddopdu = dop.GetDerivative(U, {&x}, {nodes});
+
+      HypreParMatrix *A = nullptr;
+      ddopdu->Assemble(A);
+
+      fes.GetRestrictionMatrix()->Mult(x, X);
+      A->Mult(X, Z);
+
+      Y -= Z;
+      real_t norm_g, norm_l = Y.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+      MPI_Barrier(MPI_COMM_WORLD);
+      delete A;
+   }
+
+}
+
+// no GPU tag to avoid failing 'hypre parallel mat' section
+TEST_CASE("dFEM Multiple Integrators", "[Parallel][dFEM][XXX]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+
+   SECTION("2d")
+   {
+      const auto filename2d =
+         GENERATE(
+            "../../data/star.mesh",
+            "../../data/star-q3.mesh",
+            "../../data/rt-2d-q3.mesh",
+            "../../data/inline-quad.mesh",
+            "../../data/periodic-square.mesh"
+         );
+      mult_integ<2>(filename2d, p);
+   }
+
+   SECTION("3d")
+   {
+      const auto filename3d =
+         GENERATE(
+            "../../data/fichera.mesh",
+            "../../data/fichera-q3.mesh",
+            "../../data/inline-hex.mesh",
+            "../../data/toroid-hex.mesh",
+            "../../data/periodic-cube.mesh"
+         );
+      mult_integ<3>(filename3d, p);
+   }
+}
+
+#endif // MFEM_USE_MPI

tests/unit/dfem/test_transpose.cpp

@@ -0,0 +1,409 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../unit_tests.hpp"
+#include "mfem.hpp"
+
+#include "fem/dfem/doperator.hpp" // TODO: remove before merge
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+#ifdef MFEM_USE_ENZYME
+using dscalar_t = real_t;
+#else
+using mfem::future::dual;
+using dscalar_t = dual<real_t, real_t>;
+#endif
+
+template <int DIM>
+void transpose(const char *filename, int p)
+{
+   CAPTURE(filename, DIM, p);
+
+   Mesh serial_mesh(filename);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+   mesh.EnsureNodes();
+   auto* nodes = static_cast<ParGridFunction*>(mesh.GetNodes());
+   p = std::max(p, mesh.GetNodalFESpace()->GetMaxElementOrder());
+
+   Array<int> all_domain_attr;
+   if (mesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(mesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace scalar_fes(&mesh, &fec);
+   ParFiniteElementSpace vector_fes(&mesh, &fec, DIM);
+
+   ParGridFunction sgf(&scalar_fes);
+
+   auto f0 = [](const Vector &x)
+   {
+      if constexpr (DIM == 3)
+      {
+         return M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]) * sin(M_PI*x[2]);
+      }
+      return M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]);
+   };
+
+   FunctionCoefficient f0_coeff(f0);
+   sgf.ProjectCoefficient(f0_coeff);
+
+   ParGridFunction vgf(&vector_fes);
+
+   auto gradf1 = [](const Vector &x, Vector &u)
+   {
+      if constexpr (DIM == 3)
+      {
+         u(0) = M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]) * sin(M_PI*x[2]);
+         u(1) = M_PI*sin(M_PI*x[0]) * cos(M_PI*x[1]) * sin(M_PI*x[2]);
+         u(2) = M_PI*sin(M_PI*x[0]) * sin(M_PI*x[1]) * cos(M_PI*x[2]);
+         return;
+      }
+      u(0) = M_PI*cos(M_PI*x[0]) * sin(M_PI*x[1]);
+      u(1) = M_PI*sin(M_PI*x[0]) * cos(M_PI*x[1]);
+   };
+
+   VectorFunctionCoefficient gradf1_coeff(DIM, gradf1);
+   vgf.ProjectCoefficient(gradf1_coeff);
+
+   const auto* ir = &IntRules.Get(mesh.GetTypicalElementGeometry(), 2 * p);
+
+   SECTION("Mass Transpose Action")
+   {
+      ParBilinearForm Mblf(&scalar_fes);
+      auto mass_integ = new MassIntegrator;
+      mass_integ->SetIntegrationRule(*ir);
+      Mblf.AddDomainIntegrator(mass_integ);
+      Mblf.Assemble();
+      Mblf.Finalize();
+      auto Mmat = Mblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector{FieldDescriptor{COORDINATES, nodes->ParFESpace()}};
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const dscalar_t &u,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         return tuple{u * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Value<SCALAR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<SCALAR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector S, T, U;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      U.SetSize(scalar_fes.GetTrueVSize());
+
+      sgf.GetTrueDofs(S);
+
+      Mmat->MultTranspose(S, T);
+
+      auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {nodes});
+      ddop->MultTranspose(S, U);
+
+      T -= U;
+      real_t norm_g, norm_l = T.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Mmat;
+   }
+
+   SECTION("Vector Mass Transpose Action")
+   {
+      ParBilinearForm Mvblf(&vector_fes);
+      auto mass_integ = new VectorMassIntegrator;
+      mass_integ->SetIntegrationRule(*ir);
+      Mvblf.AddDomainIntegrator(mass_integ);
+      Mvblf.Assemble();
+      Mvblf.Finalize();
+      auto Mvmat = Mvblf.ParallelAssemble();
+
+      static constexpr int VECTOR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{VECTOR, &vector_fes}};
+      const auto par = std::vector{FieldDescriptor{COORDINATES, nodes->ParFESpace()}};
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const tensor<dscalar_t, DIM> &u,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         return tuple{u * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, VECTOR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Value<VECTOR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VECTOR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector V, W, Z;
+      V.SetSize(vector_fes.GetTrueVSize());
+      W.SetSize(vector_fes.GetTrueVSize());
+      Z.SetSize(vector_fes.GetTrueVSize());
+
+      vgf.GetTrueDofs(V);
+
+      Mvmat->MultTranspose(V, W);
+
+      auto ddop = dop.GetDerivative(VECTOR, {&vgf}, {nodes});
+      ddop->MultTranspose(V, Z);
+
+      W -= Z;
+      real_t norm_g, norm_l = W.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Mvmat;
+   }
+
+   SECTION("Discrete Gradient Transpose Action")
+   {
+      ParMixedBilinearForm Gblf(&scalar_fes, &vector_fes);
+      auto grad_integ = new GradientIntegrator;
+      grad_integ->SetIntegrationRule(*ir);
+      Gblf.AddDomainIntegrator(grad_integ);
+      Gblf.Assemble();
+      Gblf.Finalize();
+      auto Gmat = Gblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, VECTOR = 2, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{VECTOR, &vector_fes},
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+      const auto gradient_qf = [] MFEM_HOST_DEVICE(
+                                  const tensor<dscalar_t, DIM> &dudxi,
+                                  const tensor<real_t, DIM, DIM> &J,
+                                  const real_t &w)
+      {
+         const auto dudx = dudxi * inv(J);
+         return tuple{dudx * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(gradient_qf,
+                              tuple{Gradient<SCALAR>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VECTOR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({&vgf, nodes});
+
+      Vector S, T, V;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      vgf.GetTrueDofs(V);
+
+      Gmat->MultTranspose(V, S);
+
+      auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {&vgf, nodes});
+      ddop->MultTranspose(V, T);
+
+      S -= T;
+      real_t norm_g, norm_l = S.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+
+      delete Gmat;
+   }
+
+   SECTION("Scalar Convection Transpose Action")
+   {
+      auto b_func = [](const Vector &x, Vector &b)
+      {
+         b(0) = cos(x[0] * 2.0 * M_PI);
+         b(1) = 1.0 + cos(x[1] * 2.0 * M_PI);
+         if constexpr (DIM == 3)
+         {
+            b(2) = 2.0 + cos(x[2] * 2.0 * M_PI);
+         }
+      };
+      VectorFunctionCoefficient b_coeff(DIM, b_func);
+
+      ParBilinearForm Gblf(&scalar_fes);
+      auto conv_integ = new ConvectionIntegrator(b_coeff);
+      conv_integ->SetIntegrationRule(*ir);
+      Gblf.AddDomainIntegrator(conv_integ);
+      Gblf.Assemble();
+      Gblf.Finalize();
+      auto Gmat = Gblf.ParallelAssemble();
+
+      static constexpr int SCALAR = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{SCALAR, &scalar_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+
+      const auto convection_qf =
+         [] MFEM_HOST_DEVICE(
+            const tensor<dscalar_t, DIM> &dudxi,
+            const tensor<real_t, DIM> &x,
+            const tensor<real_t, DIM, DIM> &J,
+            const real_t &w)
+      {
+         const auto dudx = dudxi * inv(J);
+         tensor<dscalar_t, DIM> b{};
+         b(0) = cos(x[0] * 2.0 * M_PI);
+         b(1) = 1.0 + cos(x[1] * 2.0 * M_PI);
+         if constexpr (DIM == 3)
+         {
+            b(2) = 2.0 + cos(x[2] * 2.0 * M_PI);
+         }
+         return tuple{dot(b, dudx) * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, SCALAR> {};
+      dop.AddDomainIntegrator(convection_qf,
+                              tuple{Gradient<SCALAR>{}, Value<COORDINATES>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<SCALAR>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      Vector S, T, U;
+      S.SetSize(scalar_fes.GetTrueVSize());
+      T.SetSize(scalar_fes.GetTrueVSize());
+      U.SetSize(scalar_fes.GetTrueVSize());
+      U.Randomize(1);
+
+      {
+         Gmat->MultTranspose(U, S);
+
+         auto ddop = dop.GetDerivative(SCALAR, {&sgf}, {nodes});
+         ddop->MultTranspose(U, T);
+
+         S -= T;
+         real_t norm_g, norm_l = S.Normlinf();
+         MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+         REQUIRE(norm_g == MFEM_Approx(0.0));
+      }
+
+      delete Gmat;
+   }
+
+   SECTION("Nonlinear VectorConvection Transpose Action")
+   {
+      auto b_func = [](const Vector &x, Vector &b)
+      {
+         b(0) = cos(x[0]) * sin(x[0]) * x[1];
+         b(1) = cos(x[1]) * sin(x[1]) * x[0];
+         if constexpr (DIM == 3)
+         {
+            b(2) = cos(x[2]) * sin(x[2]) * x[0];
+         }
+      };
+      VectorFunctionCoefficient b_coeff(DIM, b_func);
+
+      ParGridFunction ugf(&vector_fes);
+      ugf.ProjectCoefficient(b_coeff);
+
+      Vector U(vector_fes.GetTrueVSize());
+      ugf.GetTrueDofs(U);
+
+      ParNonlinearForm nlf(&vector_fes);
+      const auto vcinteg = new VectorConvectionNLFIntegrator();
+      vcinteg->SetIntegrationRule(*ir);
+      nlf.AddDomainIntegrator(vcinteg);
+      HypreParMatrix &Nmat = dynamic_cast<HypreParMatrix&>(nlf.GetGradient(U));
+
+      static constexpr int VELOCITY = 0, COORDINATES = 1;
+      const auto sol = std::vector{FieldDescriptor{VELOCITY, &vector_fes}};
+      const auto par = std::vector
+      {
+         FieldDescriptor{COORDINATES, nodes->ParFESpace()}
+      };
+      DifferentiableOperator dop(sol, par, mesh);
+
+      const auto nlconvection_qf =
+         [] MFEM_HOST_DEVICE(
+            const tensor<dscalar_t, DIM> &u,
+            const tensor<dscalar_t, DIM, DIM> &dudxi,
+            const tensor<real_t, DIM, DIM> &J,
+            const real_t &w)
+      {
+         const auto invJ = inv(J);
+         const auto dudx = dudxi * invJ;
+         return tuple{dot(dudx, u) * w * det(J)};
+      };
+
+      auto derivatives = std::integer_sequence<size_t, VELOCITY> {};
+      dop.AddDomainIntegrator(nlconvection_qf,
+                              tuple{Value<VELOCITY>{}, Gradient<VELOCITY>{}, Gradient<COORDINATES>{}, Weight{}},
+                              tuple{Value<VELOCITY>{}},
+                              *ir, all_domain_attr, derivatives);
+      dop.SetParameters({nodes});
+
+      auto ddop = dop.GetDerivative(VELOCITY, {&ugf}, {nodes});
+
+      Vector S(U.Size()), T(U.Size()), Se(vector_fes.GetVSize());
+
+      Nmat.MultTranspose(U, S);
+      ddop->MultTranspose(U, T);
+
+      S -= T;
+      real_t norm_g, norm_l = S.Normlinf();
+      MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX, mesh.GetComm());
+      REQUIRE(norm_g == MFEM_Approx(0.0));
+   }
+}
+
+TEST_CASE("dFEM Transpose", "[Parallel][dFEM][XXX]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+
+   const auto p = !all_tests ? 2 : GENERATE(1, 2, 3);
+
+   SECTION("2d")
+   {
+      const auto filename2d =
+         GENERATE(
+            "../../data/star.mesh",
+            "../../data/star-q3.mesh",
+            "../../data/rt-2d-q3.mesh",
+            "../../data/inline-quad.mesh"
+         );
+      transpose<2>(filename2d, p);
+   }
+
+   SECTION("3d")
+   {
+      const auto filename3d =
+         GENERATE(
+            "../../data/fichera.mesh",
+            "../../data/fichera-q3.mesh",
+            "../../data/inline-hex.mesh",
+            "../../data/toroid-hex.mesh"
+         );
+      transpose<3>(filename3d, p);
+   }
+}
+
+#endif