Mixed derivatives for dFEM functional.

fem/dfem/doperator.cpp

@@ -40,21 +40,37 @@ using SparseAssemblyMap =
    std::map<size_t,
    std::vector<assemble_derivative_sparsematrix_callback_t>>;
 using HypreAssemblyMap =
-   std::map<size_t,
-   std::vector<assemble_derivative_hypreparmatrix_callback_t>>;
+    std::map<size_t,
+             std::vector<assemble_derivative_hypreparmatrix_callback_t>>;
 using DiagonalAssemblyMap =
-   std::map<size_t, std::vector<assemble_diagonal_callback_t>>;
+    std::map<size_t, std::vector<assemble_diagonal_callback_t>>;
+using SecondDerivativeActionMap =
+    std::map<second_derivative_key_t, std::vector<derivative_action_t>>;
+using SecondDerivativeSetupMap =
+    std::map<second_derivative_key_t, std::vector<derivative_setup_t>>;
+using SecondDerivativeFieldMap =
+    std::map<second_derivative_key_t, std::vector<FieldDescriptor>>;
+using SecondSparseAssemblyMap =
+    std::map<second_derivative_key_t,
+             std::vector<assemble_derivative_sparsematrix_callback_t>>;
+using SecondHypreAssemblyMap =
+    std::map<second_derivative_key_t,
+             std::vector<assemble_derivative_hypreparmatrix_callback_t>>;
+using SecondDiagonalAssemblyMap =
+    std::map<second_derivative_key_t, std::vector<assemble_diagonal_callback_t>>;
 
 template <typename map_t>
 const typename map_t::mapped_type &FindOrDefault(
-   const map_t &map, size_t id, const typename map_t::mapped_type &fallback)
+    const map_t &map, const typename map_t::key_type &id,
+    const typename map_t::mapped_type &fallback)
 {
    const auto it = map.find(id);
    return it == map.end() ? fallback : it->second;
 }
 
 template <typename map_t>
-typename map_t::mapped_type FindOrEmpty(const map_t &map, size_t id)
+typename map_t::mapped_type FindOrEmpty(
+   const map_t &map, const typename map_t::key_type &id)
 {
    const auto it = map.find(id);
    return it == map.end() ? typename map_t::mapped_type{} : it->second;
@@ -91,6 +107,19 @@ struct DerivativeCallbackSet
    const char *missing_action_message;
 };
 
+struct SecondDerivativeCallbackSet
+{
+   const SecondDerivativeActionMap &actions;
+   const SecondDerivativeActionMap &cached_actions;
+   const SecondDerivativeActionMap &transpose_actions;
+   const SecondDerivativeFieldMap &outfds;
+   const SecondSparseAssemblyMap &assemble_sparse;
+   const SecondHypreAssemblyMap &assemble_hypre;
+   const SecondDiagonalAssemblyMap &assemble_diagonal;
+   const SecondDerivativeSetupMap &setup;
+   const char *missing_action_message;
+};
+
 template <typename vector_t>
 std::shared_ptr<DerivativeOperator> MakeStatefulDerivativeOperator(
    size_t derivative_id,
@@ -125,6 +154,63 @@ std::shared_ptr<DerivativeOperator> MakeStatefulDerivativeOperator(
              FindOrEmpty(callbacks.assemble_diagonal, derivative_id),
              FindOrEmpty(callbacks.setup, derivative_id));
 }
+
+const std::vector<derivative_action_t> &SelectSecondDerivativeActionCallbacks(
+    const std::vector<derivative_action_t> &direct_actions,
+    const SecondDerivativeActionMap &cached_actions,
+    second_derivative_key_t derivative_key,
+    bool use_cached_setup)
+{
+   if (use_cached_setup)
+   {
+      const auto it_apply = cached_actions.find(derivative_key);
+      if (it_apply != cached_actions.end() && !it_apply->second.empty())
+      {
+         return it_apply->second;
+      }
+   }
+
+   return direct_actions;
+}
+
+template <typename vector_t>
+std::shared_ptr<DerivativeOperator> MakeStatefulSecondDerivativeOperator(
+    size_t gradient_id,
+    size_t direction_id,
+    const vector_t &x,
+    const std::vector<FieldDescriptor> &infds,
+    const std::vector<FieldDescriptor> &default_outfds,
+    const SecondDerivativeCallbackSet &callbacks,
+    bool use_cached_setup)
+{
+   const second_derivative_key_t derivative_key{gradient_id, direction_id};
+   const auto it_action = callbacks.actions.find(derivative_key);
+   MFEM_ASSERT(it_action != callbacks.actions.end(),
+               callbacks.missing_action_message << "(" << gradient_id << ", "
+                                                << direction_id << ")");
+
+   const size_t dfidx = FindIdx(direction_id, infds);
+   const auto &doutfds =
+       FindOrDefault(callbacks.outfds, derivative_key, default_outfds);
+   const auto &mult_callbacks =
+       SelectSecondDerivativeActionCallbacks(
+           it_action->second, callbacks.cached_actions, derivative_key,
+           use_cached_setup);
+
+   return std::make_shared<DerivativeOperator>(
+       GetTotalTrueVSize(doutfds),
+       GetTrueVSize(infds[dfidx]),
+       mult_callbacks,
+       FindOrEmpty(callbacks.transpose_actions, derivative_key),
+       infds[dfidx],
+       x,
+       infds,
+       doutfds,
+       FindOrEmpty(callbacks.assemble_sparse, derivative_key),
+       FindOrEmpty(callbacks.assemble_hypre, derivative_key),
+       FindOrEmpty(callbacks.assemble_diagonal, derivative_key),
+       FindOrEmpty(callbacks.setup, derivative_key));
+}
 }
 
 DifferentiableOperator::DifferentiableOperator(
@@ -241,12 +327,18 @@ std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetDerivative(
 
 std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetSecondDerivative(
    size_t derivative_id, const Vector &x)
+{
+   return GetSecondDerivative(derivative_id, derivative_id, x);
+}
+
+std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetSecondDerivative(
+   size_t gradient_id, size_t direction_id, const Vector &x)
 {
    MFEM_ASSERT(has_functional_integrator,
                "second derivatives are available only for functionals");
 
-   return MakeStatefulDerivativeOperator(
-             derivative_id, x, infds, outfds,
+   return MakeStatefulSecondDerivativeOperator(
+             gradient_id, direction_id, x, infds, outfds,
    {
       second_derivative_action_callbacks,
       second_derivative_apply_callbacks,
@@ -263,12 +355,19 @@ std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetSecondDerivative(
 
 std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetSecondDerivative(
    size_t derivative_id, const MultiVector &x, const bool use_cached_setup)
+{
+   return GetSecondDerivative(derivative_id, derivative_id, x, use_cached_setup);
+}
+
+std::shared_ptr<DerivativeOperator> DifferentiableOperator::GetSecondDerivative(
+   size_t gradient_id, size_t direction_id, const MultiVector &x,
+   const bool use_cached_setup)
 {
    MFEM_ASSERT(has_functional_integrator,
                "second derivatives are available only for functionals");
 
-   return MakeStatefulDerivativeOperator(
-             derivative_id, x, infds, outfds,
+   return MakeStatefulSecondDerivativeOperator(
+             gradient_id, direction_id, x, infds, outfds,
    {
       second_derivative_action_callbacks,
       second_derivative_apply_callbacks,

fem/dfem/doperator.hpp

@@ -14,6 +14,7 @@
 
 #ifdef MFEM_USE_MPI
 #include <memory>
+#include <utility>
 
 #include "../../linalg/multivector.hpp"
 
@@ -126,6 +127,8 @@ using assemble_diagonal_callback_t = std::function<void(Vector &)>;
 using restriction_callback_t =
    std::function<void(std::vector<Vector> &, std::vector<Vector> &)>;
 
+using second_derivative_key_t = std::pair<size_t, size_t>;
+
 namespace detail
 {
 template <typename T, typename = void>
@@ -821,6 +824,24 @@ public:
 
    /// @brief Create a second derivative operator for a functional.
    ///
+   /// Returns the derivative of grad_{gradient_id} f in the direction of
+   /// @a direction_id. For example, GetSecondDerivative(U, Rho, x) applies
+   /// d/dRho(grad_U f) = d^2 f / dU dRho to the current state @a x.
+   ///
+   /// This overload accepts the state as a T-vector BlockVector and uses direct
+   /// derivative-action callbacks.
+   ///
+   /// @param gradient_id The derivative ID of the gradient to be differentiated.
+   /// @param direction_id The derivative ID of the direction in which to differentiate the gradient.
+   /// @param x Current state as a BlockVector stored through the Vector
+   /// interface.
+   /// @return A shared pointer to the configured DerivativeOperator.
+   std::shared_ptr<DerivativeOperator> GetSecondDerivative(
+      size_t gradient_id, size_t direction_id, const Vector &x);
+
+   /// @brief Create a second derivative operator for a functional.
+   ///        Shorthand of the above, expands to GetSecondDerivative(derivative_id, derivative_id, x).
+   ///
    /// Returns a DerivativeOperator representing the second derivative of a
    /// functional with respect to the given derivative ID. This is available for
    /// integrators that are declared as a functional.
@@ -851,6 +872,21 @@ public:
       size_t derivative_id, const MultiVector &x,
       const bool use_cached_setup = false);
 
+   /// @brief Create a mixed second derivative operator for a functional.
+   ///
+   /// This overload accepts the state as a MultiVector. 
+   ///
+   /// @param gradient_id The derivative ID of the gradient to be differentiated.
+   /// @param direction_id The derivative ID of the direction in which to differentiate the gradient.
+   /// @param x Current state as a BlockVector stored through the Vector
+   /// interface.
+   /// @param use_cached_setup Whether to prefer cached derivative-apply
+   /// callbacks over direct derivative actions.
+   /// @return A shared pointer to the configured DerivativeOperator.
+   std::shared_ptr<DerivativeOperator> GetSecondDerivative(
+      size_t gradient_id, size_t direction_id, const MultiVector &x,
+      const bool use_cached_setup = false);
+
    template <typename qfunc_t>
    qfunc_t *GetDerivativeActionQFunction(size_t derivative_id,
                                          size_t integrator = 0)
@@ -913,33 +949,39 @@ private:
        std::vector<derivative_action_t>> daction_transpose_callbacks;
    std::map<size_t, std::vector<FieldDescriptor>> derivative_outfds;
    std::map<size_t, std::vector<FieldDescriptor>> derivative_unionfds;
-   std::map<size_t,
+   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,
-       std::vector<assemble_diagonal_callback_t>>
+   std::map<size_t, std::vector<assemble_diagonal_callback_t>>
        assemble_diagonal_callbacks;
-   std::map<size_t, std::vector<derivative_setup_t>>
-                                                  second_derivative_setup_callbacks;
-   std::map<size_t,
-       std::vector<derivative_action_t>> second_derivative_action_callbacks;
-   std::map<size_t,
-       std::vector<derivative_action_t>> second_derivative_apply_callbacks;
-   std::map<size_t,
-       std::vector<derivative_action_t>> second_daction_transpose_callbacks;
-   std::map<size_t, std::vector<FieldDescriptor>> second_derivative_outfds;
-   std::map<size_t, std::vector<FieldDescriptor>> second_derivative_unionfds;
-   std::map<size_t,
-       std::vector<assemble_derivative_sparsematrix_callback_t>>
+   std::map<second_derivative_key_t, std::vector<derivative_setup_t>>
+       second_derivative_setup_callbacks;
+   std::map<second_derivative_key_t,
+            std::vector<derivative_action_t>>
+       second_derivative_action_callbacks;
+   std::map<second_derivative_key_t,
+            std::vector<derivative_action_t>>
+       second_derivative_apply_callbacks;
+   std::map<second_derivative_key_t,
+            std::vector<derivative_action_t>>
+       second_daction_transpose_callbacks;
+   std::map<second_derivative_key_t,
+            std::vector<FieldDescriptor>>
+       second_derivative_outfds;
+   std::map<second_derivative_key_t,
+            std::vector<FieldDescriptor>>
+       second_derivative_unionfds;
+   std::map<second_derivative_key_t,
+            std::vector<assemble_derivative_sparsematrix_callback_t>>
        assemble_second_derivative_sparsematrix_callbacks;
-   std::map<size_t,
-       std::vector<assemble_derivative_hypreparmatrix_callback_t>>
+   std::map<second_derivative_key_t,
+            std::vector<assemble_derivative_hypreparmatrix_callback_t>>
        assemble_second_derivative_hypreparmatrix_callbacks;
-   std::map<size_t,
-       std::vector<assemble_diagonal_callback_t>>
+   std::map<second_derivative_key_t,
+            std::vector<assemble_diagonal_callback_t>>
        assemble_second_derivative_diagonal_callbacks;
 
    std::vector<FieldDescriptor> infds;
@@ -1143,12 +1185,35 @@ void DifferentiableOperator::AddIntegrator(
    };
 
    [[maybe_unused]] auto set_second_derivative_fds =
+      [&](second_derivative_key_t derivative_key,
+          const std::vector<FieldDescriptor> &out,
+          const std::vector<FieldDescriptor> &all)
+      -> IntegratorContext
+   {
+      auto &stored_out = second_derivative_outfds[derivative_key];
+      stored_out = stored_out.empty() ? out : stored_out;
+      MFEM_VERIFY(stored_out == out,
+                     "inconsistent second derivative output FieldDescriptors");
+      auto &stored_union = second_derivative_unionfds[derivative_key];
+      stored_union = stored_union.empty() ? all : stored_union;
+      MFEM_VERIFY(stored_union == all,
+                  "inconsistent second derivative union FieldDescriptors");
+
+      return IntegratorContext
+      {
+         mesh, elem_attributes, attributes, num_entities,
+         infds, stored_out, stored_union, integration_rule,
+         in_qlayouts, out_qlayouts
+      };
+   };
+
+   [[maybe_unused]] auto set_functional_derivative_fds =
       [&](size_t derivative_id,
           const std::vector<FieldDescriptor> &out,
           const std::vector<FieldDescriptor> &all)
       -> IntegratorContext
    {
-      auto &stored_out = second_derivative_outfds[derivative_id];
+      auto &stored_out = derivative_outfds[derivative_id];
       if (stored_out.empty())
       {
          stored_out = out;
@@ -1156,10 +1221,10 @@ void DifferentiableOperator::AddIntegrator(
       else
       {
          MFEM_VERIFY(stored_out == out,
-                     "inconsistent second derivative output FieldDescriptors");
+                     "inconsistent first derivative output FieldDescriptors");
       }
 
-      auto &stored_union = second_derivative_unionfds[derivative_id];
+      auto &stored_union = derivative_unionfds[derivative_id];
       if (stored_union.empty())
       {
          stored_union = all;
@@ -1167,7 +1232,7 @@ void DifferentiableOperator::AddIntegrator(
       else
       {
          MFEM_VERIFY(stored_union == all,
-                     "inconsistent second derivative union FieldDescriptors");
+                     "inconsistent first derivative union FieldDescriptors");
       }
 
       return IntegratorContext
@@ -1179,19 +1244,20 @@ void DifferentiableOperator::AddIntegrator(
    };
 
    for_constexpr([&](auto i)
-   {
-      integrator_qp_caches.emplace_back(std::make_unique<Vector>());
-      Vector &qp_cache = *integrator_qp_caches.back();
-
-      auto create_callbacks = [&](auto &setup_callbacks,
+                 {
+      auto create_callbacks = [&](auto derivative_id,
+                      auto callback_key,
+                      auto &setup_callbacks,
                                   auto &apply_callbacks,
                                   auto &transpose_callbacks,
                                   auto &assemble_sparsematrix_callbacks,
                                   auto &assemble_diagonal_cbs,
                                   auto &action_cbs,
+                                  Vector &callback_qp_cache,
                                   const IntegratorContext &callback_ctx,
                                   auto &qf, auto outputs)
       {
+         constexpr size_t derivative_idx = decltype(derivative_id)::value;
          using callback_outputs_t = std::decay_t<decltype(outputs)>;
 
          bool disable_assemble = false;
@@ -1205,41 +1271,41 @@ void DifferentiableOperator::AddIntegrator(
          }, std::make_index_sequence<tuple_size<callback_outputs_t>::value> {});
 
          // Setup the qp cache for the derivative
-         setup_callbacks[i].push_back(
+         setup_callbacks[callback_key].push_back(
             MakeDerivativeSetupCallback(
-               backend_t::template MakeDerivativeSetup<i>(
-                  callback_ctx, qf, inputs, outputs, qp_cache)));
+               backend_t::template MakeDerivativeSetup<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs, callback_qp_cache)));
 
          // Apply the derivative to the qp cache
-         apply_callbacks[i].push_back(
+         apply_callbacks[callback_key].push_back(
             derivative_action_t(
-               backend_t::template MakeDerivativeApply<i>(
-                  callback_ctx, qf, inputs, outputs, qp_cache)));
+               backend_t::template MakeDerivativeApply<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs, callback_qp_cache)));
 
          // Apply the transpose of the derivative to the qp cache
-         transpose_callbacks[i].push_back(
+         transpose_callbacks[callback_key].push_back(
             derivative_action_t(
-               backend_t::template MakeDerivativeApplyTranspose<i>(
-                  callback_ctx, qf, inputs, outputs, qp_cache)));
+               backend_t::template MakeDerivativeApplyTranspose<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs, callback_qp_cache)));
 
          if (!disable_assemble)
          {
             // Assemble the derivative into a SparseMatrix
-            assemble_sparsematrix_callbacks[i].push_back(
-               backend_t::template MakeDerivativeAssemble<i>(
-                  callback_ctx, qf, inputs, outputs, qp_cache));
+            assemble_sparsematrix_callbacks[callback_key].push_back(
+               backend_t::template MakeDerivativeAssemble<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs, callback_qp_cache));
 
             // Assemble the diagonal of the derivative into an L-vector
-            assemble_diagonal_cbs[i].push_back(
-               backend_t::template MakeDerivativeAssembleDiagonal<i>(
-                  callback_ctx, qf, inputs, outputs, qp_cache));
+            assemble_diagonal_cbs[callback_key].push_back(
+               backend_t::template MakeDerivativeAssembleDiagonal<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs, callback_qp_cache));
          }
 
          // Apply the derivative
-         action_cbs[i].push_back(
+         action_cbs[callback_key].push_back(
             MakeDerivativeActionCallback(
-               backend_t::template MakeDerivativeAction<i>(callback_ctx, qf,
-                                                           inputs, outputs)));
+               backend_t::template MakeDerivativeAction<derivative_idx>(
+                  callback_ctx, qf, inputs, outputs)));
       };
 
 #ifdef MFEM_USE_ENZYME
@@ -1276,17 +1342,29 @@ void DifferentiableOperator::AddIntegrator(
 
          const auto derivative_all_fds =
             make_union_fds(infds, derivative_outputs_fds);
-         const auto derivative_ctx =
-            set_second_derivative_fds(idx, derivative_outputs_fds,
-                                      derivative_all_fds);
+         const auto first_derivative_ctx =
+            set_functional_derivative_fds(idx, derivative_outputs_fds,
+                                          derivative_all_fds);
+         for_constexpr([&](auto j)
+         {
+            constexpr size_t direction_idx = decltype(j)::value;
+            integrator_qp_caches.emplace_back(std::make_unique<Vector>());
+            Vector &second_qp_cache = *integrator_qp_caches.back();
+            const second_derivative_key_t derivative_key{idx, direction_idx};
+            const auto derivative_ctx =
+               set_second_derivative_fds(derivative_key, derivative_outputs_fds,
+                                         derivative_all_fds);
 
-         create_callbacks(second_derivative_setup_callbacks,
-                          second_derivative_apply_callbacks,
-                          second_daction_transpose_callbacks,
-                          assemble_second_derivative_sparsematrix_callbacks,
-                          assemble_second_derivative_diagonal_callbacks,
-                          second_derivative_action_callbacks,
-                          derivative_ctx, dqfunc, first_derivative_outputs);
+            create_callbacks(j, derivative_key,
+                             second_derivative_setup_callbacks,
+                             second_derivative_apply_callbacks,
+                             second_daction_transpose_callbacks,
+                             assemble_second_derivative_sparsematrix_callbacks,
+                             assemble_second_derivative_diagonal_callbacks,
+                             second_derivative_action_callbacks,
+                             second_qp_cache,
+                             derivative_ctx, dqfunc, first_derivative_outputs);
+         }, derivative_ids);
 
          // The first derivative (gradient) of the functional is the plain
          // action of the reverse-mode-differentiated energy dqfunc. Register
@@ -1295,7 +1373,7 @@ void DifferentiableOperator::AddIntegrator(
          // gradient is a function of the captured state only, so the direction
          // is ignored here.
          auto grad_action =
-            backend_t::MakeAction(derivative_ctx, dqfunc, inputs,
+            backend_t::MakeAction(first_derivative_ctx, dqfunc, inputs,
                                   first_derivative_outputs);
          derivative_action_callbacks[idx].push_back(
             [grad_action](const std::vector<Vector *> &xe,
@@ -1305,31 +1383,24 @@ void DifferentiableOperator::AddIntegrator(
             grad_action(xe, ye);
          });
 
-         auto &stored_first_out = derivative_outfds[idx];
-         if (stored_first_out.empty())
-         {
-            stored_first_out = derivative_outputs_fds;
-         }
-         else
-         {
-            MFEM_VERIFY(stored_first_out == derivative_outputs_fds,
-                        "inconsistent first derivative output FieldDescriptors");
-         }
 #else
          MFEM_ABORT("functional integrators require Enzyme support to compute derivatives");
 #endif
       }
       else
       {
-         create_callbacks(derivative_setup_callbacks,
+         integrator_qp_caches.emplace_back(std::make_unique<Vector>());
+         Vector &qp_cache = *integrator_qp_caches.back();
+         create_callbacks(i, idx,
+                          derivative_setup_callbacks,
                           derivative_apply_callbacks,
                           daction_transpose_callbacks,
                           assemble_derivative_sparsematrix_callbacks,
                           assemble_diagonal_callbacks,
                           derivative_action_callbacks,
+                          qp_cache,
                           ctx, qfunc, outputs);
-      }
-   }, derivative_ids);
+      } }, derivative_ids);
 }
 
 } // namespace mfem::future

tests/unit/dfem/test_second_derivative.cpp

@@ -111,6 +111,34 @@ struct MinimalSurfaceHessianAction
    }
 };
 
+template <typename dscalar_t, int dim>
+struct DensityWeightedQuadraticFunctional
+{
+   MFEM_HOST_DEVICE inline __attribute__((always_inline))
+   auto operator()(const dscalar_t &u,
+                   const dscalar_t &rho,
+                   const tensor<real_t, dim, dim> &J,
+                   const real_t &w,
+                   real_t &f) const
+   {
+      f = rho * u * u * det(J) * w;
+   }
+};
+
+template <typename dscalar_t, int dim>
+struct DensityWeightedQuadraticMixedAction
+{
+   MFEM_HOST_DEVICE inline __attribute__((always_inline))
+   auto operator()(const real_t &drho,
+                   const dscalar_t &u,
+                   const tensor<real_t, dim, dim> &J,
+                   const real_t &w,
+                   real_t &v) const
+   {
+      v = 2.0_r * u * drho * det(J) * w;
+   }
+};
+
 template <int dim>
 class MyFunctional
 {
@@ -385,6 +413,118 @@ void second_derivative(const char *filename, int p)
    // pretty_print(Hv);
 }
 
+template <int DIM>
+void mixed_second_derivative(const char *filename, int p)
+{
+   static constexpr int U = 0, Rho = 1, Coords = 2, Q = 3, DRho = 4;
+   CAPTURE(filename, DIM, p);
+
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+
+   pmesh.EnsureNodes();
+   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
+   ParFiniteElementSpace *mfes = nodes->ParFESpace();
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace fes(&pmesh, &fec);
+   const int tvsize = fes.GetTrueVSize();
+
+   const IntegrationRule &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;
+   }
+
+   // Use smooth, non-constant fields so the mixed action is nontrivial while
+   // keeping the exact expression simple: d/dRho(grad_U J)[drho] = 2 u drho.
+   // The current rho state is intentionally non-constant, even though this
+   // particular mixed block does not depend on rho itself.
+   ParGridFunction u_gf(&fes), rho_gf(&fes), drho_gf(&fes);
+   FunctionCoefficient u_coeff(
+      [](const auto &x)
+   {
+      return 1.0_r + x[0] + 0.25_r * x[1];
+   });
+   FunctionCoefficient rho_coeff(
+      [](const auto &x)
+   {
+      return 0.5_r + 0.2_r * x[0] * x[0] + 0.1_r * x[1];
+   });
+   FunctionCoefficient drho_coeff(
+      [](const auto &x)
+   {
+      return sin(M_PI * x[0]) + 0.5_r * x[1];
+   });
+   u_gf.ProjectCoefficient(u_coeff);
+   rho_gf.ProjectCoefficient(rho_coeff);
+   drho_gf.ProjectCoefficient(drho_coeff);
+
+   Vector u(tvsize), rho(tvsize), drho(tvsize), coords;
+   u_gf.GetTrueDofs(u);
+   rho_gf.GetTrueDofs(rho);
+   drho_gf.GetTrueDofs(drho);
+   pmesh.GetNodes()->GetTrueDofs(coords);
+
+   const auto functional_in = std::vector
+   {
+      FieldDescriptor{U, &fes},
+      FieldDescriptor{Rho, &fes},
+      FieldDescriptor{Coords, mfes}
+   };
+   QuadratureSpace qspace(pmesh, ir);
+   VectorQuadratureSpace qspace_vec(qspace, 1);
+   const auto functional_out = std::vector
+   {
+      FieldDescriptor{Q, &qspace_vec}
+   };
+
+   DifferentiableOperator functional_dop(functional_in, functional_out, pmesh);
+   DensityWeightedQuadraticFunctional<real_t, DIM> functional;
+   functional_dop.AddDomainIntegrator<LocalQFBackend>(
+      functional,
+      Inputs<Value<U>, Value<Rho>, Gradient<Coords>, Weight> {},
+      Outputs<FunctionalValue<Q>> {},
+      ir, all_domain_attr,
+      Derivatives<U, Rho> {});
+
+   const auto exact_in = std::vector
+   {
+      FieldDescriptor{DRho, &fes},
+      FieldDescriptor{U, &fes},
+      FieldDescriptor{Coords, mfes}
+   };
+   const auto exact_out = std::vector
+   {
+      FieldDescriptor{U, &fes}
+   };
+   DifferentiableOperator exact_dop(exact_in, exact_out, pmesh);
+   DensityWeightedQuadraticMixedAction<real_t, DIM> exact_action;
+   exact_dop.AddDomainIntegrator<LocalQFBackend>(
+      exact_action,
+      Inputs<Value<DRho>, Value<U>, Gradient<Coords>, Weight> {},
+      Outputs<Value<U>> {},
+      ir, all_domain_attr);
+
+   MultiVector X{u, rho, coords};
+   Vector mixed(tvsize);
+   MultiVector Mixed{mixed};
+   functional_dop.GetSecondDerivative(U, Rho, X)->Mult(drho, Mixed);
+
+   MultiVector XExact{drho, u, coords};
+   Vector exact(tvsize);
+   MultiVector Exact{exact};
+   exact_dop.Mult(XExact, Exact);
+
+   Vector diff(mixed);
+   diff -= exact;
+   REQUIRE(MFEM_Approx(diff.Norml2()) == 0.0);
+}
+
 } // namespace second_derivative_test
 
 TEST_CASE("dFEM functional second derivative action matches mfem",
@@ -419,6 +559,19 @@ TEST_CASE("dFEM functional second derivative action matches mfem",
    // }
 }
 
+TEST_CASE("dFEM functional mixed second derivative action matches exact action",
+          "[Parallel][dFEM][second-derivative]")
+{
+   const bool all_tests = launch_all_non_regression_tests;
+   const auto p = !all_tests ? 1 : GENERATE(1, 2, 3);
+
+   SECTION("2d")
+   {
+      const auto f = GENERATE("../../data/inline-quad.mesh");
+      second_derivative_test::mixed_second_derivative<2>(f, p);
+   }
+}
+
 #endif // MFEM_USE_MPI
 
 #endif // MFEM_USE_ENZYME