Added dfem-reaction-diffusion (single field, multiple outputs example)

.gitignore

@@ -214,7 +214,7 @@ miniapps/adjoint/adjoint_advection_diffusion
 
 miniapps/dfem/dfem-minimal-surface
 miniapps/dfem/dfem-minimal-surface-output
-miniapps/dfem/dfem-hyperelasticity_energy
+miniapps/dfem/dfem-hyperelasticity
 
 miniapps/electromagnetics/volta
 miniapps/electromagnetics/tesla

fem/dfem/backends/local_qf/derivative_assemble.hpp

@@ -16,8 +16,10 @@
 #include "kernels.hpp"
 #include "util.hpp"
 
+#include <algorithm>
 #include <array>
 #include <type_traits>
+#include <vector>
 
 namespace ker = mfem::kernels::internal;
 
@@ -592,10 +594,10 @@ class DerivativeAssemble
    const std::array<DofToQuadMap, n_outputs> output_dtq_maps;
    const std::array<bool, n_inputs> input_is_dependent;
    const size_t trial_field_uf;
-   const size_t test_field_uf;
-   const ParFiniteElementSpace *test_fes;
+   /// Column space of every row block. GetDerivative differentiates w.r.t. one
+   /// field, so there is exactly one trial space. Null when that field is not
+   /// an FE space, in which case nothing can be assembled.
    const ParFiniteElementSpace *trial_fes;
-   const int test_vdim;
    /// Per-output row geometry of the quadrature point cache. DerivativeSetup
    /// lays that cache out over every output FieldOperator, so reading it needs
    /// all of them.
@@ -603,7 +605,19 @@ class DerivativeAssemble
    const std::array<int, n_outputs> out_op_dim;
    const std::array<int, n_outputs> out_offsets;
    const int output_size_on_qp;
-   const int num_test_dof;
+   /// Output field ids, in order of first appearance among @a outputs. Each one
+   /// is a row block of the derivative: see @ref compute_group_field_ids.
+   const std::vector<int> group_field_ids;
+   /// Output FieldOperator -> row block it contributes to.
+   const std::array<int, n_outputs> out_group;
+   /// Row block -> position in ctx.outfds (-1 if the field is not an output of
+   /// the operator, which should not happen).
+   const std::vector<int> group_outfd_idx;
+   /// Row block -> test space, null for spaces without a basis.
+   const std::vector<const ParFiniteElementSpace *> group_fes;
+   const std::vector<int> group_test_vdim;
+   const std::vector<int> group_num_test_dof;
+   const std::vector<bool> group_assemblable;
    const int trial_vdim;
    const int trial_op_dim;
    const int num_trial_dof;
@@ -611,7 +625,36 @@ class DerivativeAssemble
    const int num_trial_dof_1d;
    const int total_trial_op_dim;
    mutable Vector inputs_trial_op_dim;
-   mutable Vector Ae_mem;
+   /// One element matrix bank per row block. Outputs on the same test field
+   /// share a bank and accumulate into it; outputs on different test fields are
+   /// different row blocks and need banks of their own, because their element
+   /// matrices have different shapes and are filled through different
+   /// ElementRestrictions.
+   mutable std::vector<Vector> group_Ae_mem;
+
+   /// @brief Distinct output field ids, in order of first appearance.
+   ///
+   /// Outputs sharing a field id form one row block: they are summed into a
+   /// single element matrix and produce a single assembled matrix. This is the
+   /// multi-output case, e.g. Outputs<Value<U>, Gradient<U>> giving mass plus
+   /// diffusion. Outputs on different field ids are separate row blocks of
+   ///
+   ///     dR/dU = [ dR_U/dU ; dR_Y/dU ]
+   ///
+   /// and are assembled into one matrix each.
+   static std::vector<int> compute_group_field_ids(const outputs_t &outs)
+   {
+      std::vector<int> ids;
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         const int fid = get<o>(outs).GetFieldId();
+         if (std::find(ids.begin(), ids.end(), fid) == ids.end())
+         {
+            ids.push_back(fid);
+         }
+      });
+      return ids;
+   }
 
 public:
    DerivativeAssemble() = delete;
@@ -658,29 +701,18 @@ public:
                       ctx_in.ir)),
    input_is_dependent(compute_input_is_dependent(inputs, derivative_id)),
    trial_field_uf(find_union_field_index(ctx_in, derivative_id)),
-   test_field_uf(
-      find_union_field_index(ctx_in, get<0>(outputs).GetFieldId())),
-   test_fes(
-      [&]
-   {
-      const auto *fes = std::get_if<const ParFiniteElementSpace *>(
-         &ctx_in.unionfds[test_field_uf].data);
-      MFEM_ASSERT(fes != nullptr && *fes != nullptr,
-                  "LocalQFBackend: test space is not a ParFiniteElementSpace");
-      return *fes;
-   }()),
    trial_fes(
-      [&]
+      [&]() -> const ParFiniteElementSpace *
    {
+      // Not every field is an FE space: a derivative w.r.t. a ParameterSpace or
+      // a QuadratureSpace has no basis to assemble columns against. That is
+      // only an error if assembly is actually requested, so keep it null here
+      // and report it in operator() instead of aborting registration.
+      if (trial_field_uf >= ctx_in.unionfds.size()) { return nullptr; }
       const auto *fes = std::get_if<const ParFiniteElementSpace *>(
          &ctx_in.unionfds[trial_field_uf].data);
-      MFEM_ASSERT(fes != nullptr && *fes != nullptr,
-                  "LocalQFBackend: trial space is not a ParFiniteElementSpace");
-      return *fes;
+      return fes ? *fes : nullptr;
    }()),
-   // All outputs are attached to the same test field, so vdim is common to
-   // them; only the operator dimension differs, and that lives in out_op_dim.
-   test_vdim(get<0>(outputs).vdim),
    out_vdim(get_vdim(outputs)),
    out_op_dim(compute_out_op_dim(outputs)),
    out_offsets(compute_out_offsets(out_vdim, out_op_dim)),
@@ -691,7 +723,96 @@ public:
       for_constexpr<n_outputs>([&](auto o) { s += get<o>(outputs).size_on_qp; });
       return s;
    }()),
-   num_test_dof(test_fes->GetFE(0)->GetDof()),
+   group_field_ids(compute_group_field_ids(outputs_in)),
+   out_group(
+      [&]
+   {
+      std::array<int, n_outputs> g {};
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         const int fid = get<o>(outputs_in).GetFieldId();
+         const auto &ids = group_field_ids;
+         g[o] = static_cast<int>(std::find(ids.begin(), ids.end(), fid)
+                                 - ids.begin());
+      });
+      return g;
+   }()),
+   group_outfd_idx(
+      [&]
+   {
+      std::vector<int> idx(group_field_ids.size(), -1);
+      for (size_t g = 0; g < group_field_ids.size(); g++)
+      {
+         for (size_t f = 0; f < ctx_in.outfds.size(); f++)
+         {
+            if (static_cast<int>(ctx_in.outfds[f].id) == group_field_ids[g])
+            {
+               idx[g] = static_cast<int>(f);
+               break;
+            }
+         }
+      }
+      return idx;
+   }()),
+   group_fes(
+      [&]
+   {
+      std::vector<const ParFiniteElementSpace *> v(group_field_ids.size(),
+                                                   nullptr);
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         if (group_outfd_idx[g] < 0) { continue; }
+         const auto *fes = std::get_if<const ParFiniteElementSpace *>(
+                              &ctx_in.outfds[group_outfd_idx[g]].data);
+         v[g] = fes ? *fes : nullptr;
+      }
+      return v;
+   }()),
+   group_test_vdim(
+      [&]
+   {
+      // Outputs in a group share a field, hence a vdim; only the operator
+      // dimension differs between them, and that lives in out_op_dim.
+      std::vector<int> v(group_field_ids.size(), 0);
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         v[out_group[o]] = get<o>(outputs_in).vdim;
+      });
+      return v;
+   }()),
+   group_num_test_dof(
+      [&]
+   {
+      std::vector<int> v(group_field_ids.size(), 0);
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         if (group_fes[g] == nullptr) { continue; }
+         v[g] = group_fes[g]->GetFE(0)->GetDof();
+      }
+      return v;
+   }()),
+   group_assemblable(
+      [&]
+   {
+      // A row block can only be assembled when both its trial and test fields are
+      // ParFiniteElementSpaces. Quadrature and parameter spaces have no element
+      // basis or ElementRestriction through which to form a SparseMatrix.
+      //
+      // Identity outputs are also excluded for now as it has no supported mapping from
+      // its pointwise quadrature rows into that FE row space.
+      std::vector<bool> v(group_field_ids.size(), false);
+      if (trial_fes == nullptr) { return v; }
+      for (size_t g = 0; g < v.size(); g++) { v[g] = group_fes[g] != nullptr; }
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         using output_fop_t = std::decay_t<decltype(get<o>(outputs_in))>;
+         if constexpr (is_identity_fop_v<output_fop_t>)
+         {
+            v[out_group[o]] = false;
+         }
+      });
+      return v;
+   }()),
    trial_vdim(compute_trial_vdim(inputs, derivative_id)), trial_op_dim(
       [&]
    {
@@ -705,7 +826,7 @@ public:
       });
       return top;
    }()),
-   num_trial_dof(trial_fes->GetFE(0)->GetDof()),
+   num_trial_dof(trial_fes ? trial_fes->GetFE(0)->GetDof() : 0),
    dim(ctx_in.mesh.Dimension()), ne(ctx_in.nentities),
    nq(ctx_in.ir.GetNPoints()), q1d(tensor_1d_size(nq, dim)),
    num_trial_dof_1d(tensor_1d_size(num_trial_dof, dim)), total_trial_op_dim(
@@ -716,14 +837,12 @@ public:
       return compute_total_trial_op_dim(
          inputs, input_is_dependent, in_qp_sizes);
    }()),
-   inputs_trial_op_dim(), Ae_mem()
+   inputs_trial_op_dim(), group_Ae_mem()
    {
       MFEM_ASSERT(ctx.unionfds.size() == nfields,
                   "LocalQFBackend: unionfds size mismatch");
       MFEM_ASSERT(trial_field_uf != SIZE_MAX,
                   "DerivativeAssemble: trial field not found in unionfds");
-      MFEM_ASSERT(test_field_uf != SIZE_MAX,
-                  "DerivativeAssemble: test field not found in unionfds");
 
       MFEM_ASSERT(trial_vdim > 0,
                   "LocalQFBackend: could not determine trial vdim");
@@ -741,79 +860,114 @@ public:
             : 0;
       });
 
-      const int elem_mat_size =
-         num_test_dof * test_vdim * num_trial_dof * trial_vdim;
-      Ae_mem.SetSize(elem_mat_size * ne, Device::GetDeviceMemoryType());
-      Ae_mem.UseDevice(true);
-      Ae_mem = 0.0;
+      group_Ae_mem.resize(group_field_ids.size());
+      for (size_t g = 0; g < group_Ae_mem.size(); g++)
+      {
+         if (!group_assemblable[g]) { continue; }
+         const int elem_mat_size = group_num_test_dof[g] * group_test_vdim[g] *
+                                   num_trial_dof * trial_vdim;
+         group_Ae_mem[g].SetSize(elem_mat_size * ne,
+                                 Device::GetDeviceMemoryType());
+         group_Ae_mem[g].UseDevice(true);
+         group_Ae_mem[g] = 0.0;
+      }
    }
 
-   void operator()(SparseMatrix *&A) const
+   /// @brief Assemble one SparseMatrix per assemblable output field.
+   ///
+   /// @a A is indexed by position in ctx.outfds, so A[f] receives the row block
+   /// belonging to output field f. Slots this integrator writes nothing
+   /// assemblable to -- quadrature or parameter spaces, or fields it does not
+   /// touch at all -- are left alone, which lets several integrators fill
+   /// different row blocks of the same vector. Ownership of the matrices passes
+   /// to the caller.
+   void operator()(std::vector<SparseMatrix *> &A) const
    {
-      // Every output is contracted into one element matrix Ae, sized from the
-      // test space of get<0>(outputs), and filled through a single test
-      // ElementRestriction.
-      //
-      // WIP:
-      // This takes care of single-field, multiple-outputs case.
-      // For a multiple fields case, outputs on a second field would need a second
-      // matrix -- the derivative then eould be a block column with one row block per
-      // test space.
-      //
-      // For now we just add a check that all outputs are attached to the same test field, and abort if not.
-      for_constexpr<n_outputs>([&](auto o)
-      {
-         MFEM_VERIFY(get<o>(outputs).GetFieldId() == get<0>(outputs).GetFieldId(),
-                     "DerivativeAssemble: every output FieldOperator has to be "
-                     "attached to the same test field; assembling outputs that "
-                     "span several fields is not supported");
-      });
-
       if (ctx.attr.Size() == 0) { return; }
 
+      const bool any_assemblable =
+         std::find(group_assemblable.begin(), group_assemblable.end(), true) !=
+         group_assemblable.end();
+      if (!any_assemblable) { return; }
+
       if (!(use_sum_factorization && (dim == 2 || dim == 3)))
       {
          MFEM_ABORT("DerivativeAssemble optimized path is implemented "
                     "for tensor-product 2D/3D elements only");
       }
 
-      DerivativeAssembleHO::Run(dim,
-                                q1d,
-                                ctx,
-                                qp_cache,
-                                Ae_mem,
-                                inputs,
-                                outputs,
-                                input_dtq_maps,
-                                output_dtq_maps,
-                                out_vdim,
-                                out_op_dim,
-                                out_offsets,
-                                output_size_on_qp,
-                                inputs_trial_op_dim,
-                                test_vdim,
-                                num_test_dof,
-                                num_trial_dof,
-                                num_trial_dof_1d,
-                                trial_vdim,
-                                total_trial_op_dim,
-                                nq,
-                                ne,
-                                q1d,
-                                dim);
+      MFEM_VERIFY(trial_fes != nullptr,
+                  "DerivativeAssemble: the differentiated field is not a "
+                  "ParFiniteElementSpace, so the columns of the derivative "
+                  "have no basis to be assembled against");
 
-      A = new SparseMatrix;
-      A->OverrideSize(test_fes->GetVSize(), trial_fes->GetVSize());
-
-      const auto *test_restr = dynamic_cast<const ElementRestriction *>(
-                                  test_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC));
       const auto *trial_restr = dynamic_cast<const ElementRestriction *>(
-                                   trial_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC));
-      MFEM_VERIFY(test_restr != nullptr && trial_restr != nullptr,
+                                   trial_fes->GetElementRestriction(
+                                      ElementDofOrdering::LEXICOGRAPHIC));
+      MFEM_VERIFY(trial_restr != nullptr,
                   "DerivativeAssemble SparseMatrix assembly requires "
                   "H1/conforming ElementRestriction spaces");
 
-      test_restr->FillSparseMatrix(Ae_mem, *A, *trial_restr);
+      if (A.size() < ctx.outfds.size())
+      {
+         A.resize(ctx.outfds.size(), nullptr);
+      }
+
+      for (size_t g = 0; g < group_Ae_mem.size(); g++)
+      {
+         if (!group_assemblable[g]) { continue; }
+
+         // The kernel accumulates into Ae, so the bank has to start clean on
+         // every call: otherwise assembling twice (the next Newton step, say)
+         // would double every entry.
+         group_Ae_mem[g] = 0.0;
+
+         DerivativeAssembleHO::Run(dim,
+                                   q1d,
+                                   ctx,
+                                   qp_cache,
+                                   group_Ae_mem[g],
+                                   inputs,
+                                   outputs,
+                                   input_dtq_maps,
+                                   output_dtq_maps,
+                                   out_vdim,
+                                   out_op_dim,
+                                   out_offsets,
+                                   out_group,
+                                   static_cast<int>(g),
+                                   output_size_on_qp,
+                                   inputs_trial_op_dim,
+                                   group_test_vdim[g],
+                                   group_num_test_dof[g],
+                                   num_trial_dof,
+                                   num_trial_dof_1d,
+                                   trial_vdim,
+                                   total_trial_op_dim,
+                                   nq,
+                                   ne,
+                                   q1d,
+                                   dim);
+
+         const ParFiniteElementSpace *test_fes = group_fes[g];
+         const auto *test_restr = dynamic_cast<const ElementRestriction *>(
+                                     test_fes->GetElementRestriction(
+                                        ElementDofOrdering::LEXICOGRAPHIC));
+         MFEM_VERIFY(test_restr != nullptr,
+                     "DerivativeAssemble SparseMatrix assembly requires "
+                     "H1/conforming ElementRestriction spaces");
+
+         const int f = group_outfd_idx[g];
+         MFEM_VERIFY(A[f] == nullptr,
+                     "DerivativeAssemble: output field already carries an "
+                     "assembled matrix; two integrators contributing to the "
+                     "same row block cannot be assembled into one matrix");
+
+         auto *M = new SparseMatrix;
+         M->OverrideSize(test_fes->GetVSize(), trial_fes->GetVSize());
+         test_restr->FillSparseMatrix(group_Ae_mem[g], *M, *trial_restr);
+         A[f] = M;
+      }
    }
 
    template<typename backend_t = LocalQFHOBackend<3>, int T_Q1D = 0>
@@ -828,6 +982,8 @@ public:
       const std::array<int, n_outputs> &out_vdim,
       const std::array<int, n_outputs> &out_op_dim,
       const std::array<int, n_outputs> &out_offsets,
+      const std::array<int, n_outputs> &out_group,
+      const int group,
       const int output_size_on_qp,
       const Vector &inputs_trial_op_dim,
       const int test_vdim,
@@ -853,6 +1009,7 @@ public:
                   "DerivativeAssemble: nq exceeds backend quadrature capacity");
       for_constexpr<n_outputs>([&](auto o)
       {
+         if (out_group[o] != group) { return; }
          MFEM_VERIFY(out_op_dim[o] <= DIM,
                      "DerivativeAssemble: test_op_dim exceeds spatial DIM");
       });
@@ -896,14 +1053,25 @@ public:
 
          // Each output contributes its own rows of the cache, contracted
          // against its own test basis operation; map_quadrature_data_to_fields
-         // accumulates, so the element matrix is the sum over outputs for the
-         // same field.
+         // accumulates, so the element matrix is the sum over the outputs that
+         // belong to this row block. Outputs on other test fields are skipped
+         // here and picked up by their own launch, which has its own Ae.
+         //
+         // A launch only ever happens for an assemblable row block, and such a
+         // block contains no Identity output, so the is_identity_fop_v test
+         // below never rejects anything at run time. It is there to keep
+         // assemble_element_mat_sumfact from being instantiated for a fop it
+         // cannot handle.
          for_constexpr<n_outputs>([&](auto o)
          {
             using output_fop_t = std::decay_t<decltype(get<o>(outputs))>;
 
             if constexpr (!is_identity_fop_v<output_fop_t>)
             {
+               // Uniform across the thread block, so returning early here
+               // cannot desynchronise the barriers below.
+               if (out_group[o] != group) { return; }
+
                // The outputs share fhat_storage, so one has to be done with it
                // before the next zeroes it.
                MFEM_SYNC_THREAD;

fem/dfem/backends/local_qf/derivative_assemble_diagonal.hpp

@@ -15,12 +15,18 @@
 #include "kernels.hpp"
 #include "util.hpp"
 
+#include <algorithm>
 #include <array>
+#include <vector>
 
 namespace mfem::future::LocalQFImpl
 {
 
-// Assemble diagonal of cached Jacobian (square trial == test, tensor 2D/3D)
+// Assemble the diagonal of one row block of a cached Jacobian (tensor 2D/3D).
+//
+// The derivative is a block column, one row block per output field. Only a
+// block whose test space is the trial space is square, and only a square block
+// has a diagonal at all, so the row block is chosen per call and checked.
 
 template<int derivative_id,
          typename qfunc_t,
@@ -47,22 +53,47 @@ class DerivativeAssembleDiagonal
    const std::array<DofToQuadMap, n_outputs> output_dtq_maps;
    const std::array<bool, n_inputs> input_is_dependent;
    const size_t trial_field_uf;
-   const size_t test_field_uf;
-   const bool is_square;
-   const int test_vdim;
+   /// Column space of every row block; null when the differentiated field is
+   /// not an FE space, in which case no block has a diagonal.
+   const ParFiniteElementSpace *trial_fes;
    const std::array<int, n_outputs> out_vdim;
    const std::array<int, n_outputs> out_op_dim;
    const std::array<int, n_outputs> out_offsets;
    const int output_size_on_qp;
-   const int num_test_dof;
-   const int num_test_dof_1d;
+   /// Row blocks: the distinct output field ids, in order of first appearance.
+   /// Outputs sharing a field id are summed into one diagonal, which is how
+   /// Value<U> + Gradient<U> becomes mass plus diffusion.
+   const std::vector<int> group_field_ids;
+   const std::array<int, n_outputs> out_group;
+   const std::vector<const ParFiniteElementSpace *> group_fes;
+   const std::vector<int> group_test_vdim;
+   const std::vector<int> group_num_test_dof;
+   const std::vector<int> group_num_test_dof_1d;
+   /// Whether a row block has a diagonal: its test space has to be an FE space
+   /// and has to *be* the trial space, and no output on it may be an Identity.
+   const std::vector<bool> group_has_diagonal;
    const int trial_vdim;
    const int total_trial_op_dim;
    const int num_trial_dof_1d;
    const int residual_size_on_qp;
    const int dim, ne, nq, q1d;
    const std::array<int, n_inputs> inputs_trial_op_dim;
-   mutable Vector Ye_mem;
+   mutable std::vector<Vector> group_Ye_mem;
+
+   /// Distinct output field ids, in order of first appearance.
+   static std::vector<int> compute_group_field_ids(const outputs_t &outs)
+   {
+      std::vector<int> ids;
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         const int fid = get<o>(outs).GetFieldId();
+         if (std::find(ids.begin(), ids.end(), fid) == ids.end())
+         {
+            ids.push_back(fid);
+         }
+      });
+      return ids;
+   }
 
 public:
    DerivativeAssembleDiagonal() = delete;
@@ -109,19 +140,14 @@ public:
                       ctx_in.ir)),
    input_is_dependent(compute_input_is_dependent(inputs, derivative_id)),
    trial_field_uf(find_union_field_index(ctx_in, derivative_id)),
-   test_field_uf(
-      find_union_field_index(ctx_in, get<0>(outputs).GetFieldId())),
-   is_square(
-      [&]
+   trial_fes(
+      [&]() -> const ParFiniteElementSpace *
    {
-      const auto *test_fes = std::get_if<const ParFiniteElementSpace *>(
-         &ctx_in.unionfds[test_field_uf].data);
-      const auto *trial_fes = std::get_if<const ParFiniteElementSpace *>(
+      if (trial_field_uf >= ctx_in.unionfds.size()) { return nullptr; }
+      const auto *fes = std::get_if<const ParFiniteElementSpace *>(
          &ctx_in.unionfds[trial_field_uf].data);
-      return test_fes && trial_fes && *test_fes && *trial_fes &&
-      (*test_fes == *trial_fes);
+      return fes ? *fes : nullptr;
    }()),
-   test_vdim(get<0>(outputs).vdim),
    out_vdim(get_vdim(outputs_in)),
    out_op_dim(compute_out_op_dim(outputs_in)),
    out_offsets(compute_out_offsets(out_vdim, out_op_dim)),
@@ -132,16 +158,96 @@ public:
       for_constexpr<n_outputs>([&](auto o)
       { s += get<o>(outputs_in).size_on_qp; });
       return s;
-   }()), num_test_dof(
+   }()),
+   group_field_ids(compute_group_field_ids(outputs_in)),
+   out_group(
       [&]
    {
-      const auto *test_fes = std::get_if<const ParFiniteElementSpace *>(
-         &ctx_in.unionfds[test_field_uf].data);
-      MFEM_ASSERT(test_fes != nullptr && *test_fes != nullptr,
-                  "LocalQFBackend: test space is not a ParFiniteElementSpace");
-      return (*test_fes)->GetFE(0)->GetDof();
+      std::array<int, n_outputs> g {};
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         const int fid = get<o>(outputs_in).GetFieldId();
+         const auto &ids = group_field_ids;
+         g[o] = static_cast<int>(std::find(ids.begin(), ids.end(), fid)
+                                 - ids.begin());
+      });
+      return g;
+   }()),
+   group_fes(
+      [&]
+   {
+      std::vector<const ParFiniteElementSpace *> v(group_field_ids.size(),
+                                                   nullptr);
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         const size_t uf = find_union_field_index(ctx_in, group_field_ids[g]);
+         if (uf >= ctx_in.unionfds.size()) { continue; }
+         const auto *fes = std::get_if<const ParFiniteElementSpace *>(
+                              &ctx_in.unionfds[uf].data);
+         v[g] = fes ? *fes : nullptr;
+      }
+      return v;
+   }()),
+   group_test_vdim(
+      [&]
+   {
+      std::vector<int> v(group_field_ids.size(), 0);
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         v[out_group[o]] = get<o>(outputs_in).vdim;
+      });
+      return v;
+   }()),
+   group_num_test_dof(
+      [&]
+   {
+      std::vector<int> v(group_field_ids.size(), 0);
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         if (group_fes[g] == nullptr) { continue; }
+         v[g] = group_fes[g]->GetFE(0)->GetDof();
+      }
+      return v;
+   }()),
+   group_num_test_dof_1d(
+      [&]
+   {
+      std::vector<int> v(group_field_ids.size(), 0);
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         if (group_num_test_dof[g] > 0)
+         {
+            v[g] = tensor_1d_size(group_num_test_dof[g],
+                                  ctx_in.mesh.Dimension());
+         }
+      }
+      return v;
+   }()),
+   group_has_diagonal(
+      [&]
+   {
+      // A diagonal needs row space == column space, so only a row block on the
+      // trial space qualifies. Squareness alone cannot pick a block when
+      // several output fields share that space, which is why the caller names
+      // the row. Identity outputs are quadrature point data and are excluded
+      // for the same reason as in DerivativeAssemble: they cannot be
+      // contracted, and every output on a field lands in the same block.
+      std::vector<bool> v(group_field_ids.size(), false);
+      if (trial_fes == nullptr) { return v; }
+      for (size_t g = 0; g < v.size(); g++)
+      {
+         v[g] = (group_fes[g] != nullptr) && (group_fes[g] == trial_fes);
+      }
+      for_constexpr<n_outputs>([&](auto o)
+      {
+         using output_fop_t = std::decay_t<decltype(get<o>(outputs_in))>;
+         if constexpr (is_identity_fop_v<output_fop_t>)
+         {
+            v[out_group[o]] = false;
+         }
+      });
+      return v;
    }()),
-   num_test_dof_1d(tensor_1d_size(num_test_dof, ctx_in.mesh.Dimension())),
    trial_vdim(compute_trial_vdim(inputs, derivative_id)), total_trial_op_dim(
       [&]
    {
@@ -151,15 +257,9 @@ public:
          inputs, input_is_dependent, input_size_on_qp);
    }()),
    num_trial_dof_1d(
-      [&]
-   {
-      const auto *trial_fes = std::get_if<const ParFiniteElementSpace *>(
-         &ctx_in.unionfds[trial_field_uf].data);
-      MFEM_ASSERT(trial_fes != nullptr && *trial_fes != nullptr,
-                  "LocalQFBackend: trial space is not a ParFiniteElementSpace");
-      const int num_trial_dof = (*trial_fes)->GetFE(0)->GetDof();
-      return tensor_1d_size(num_trial_dof, ctx_in.mesh.Dimension());
-   }()),
+      trial_fes ? tensor_1d_size(trial_fes->GetFE(0)->GetDof(),
+                                 ctx_in.mesh.Dimension())
+      : 0),
    residual_size_on_qp(output_size_on_qp * trial_vdim * total_trial_op_dim),
    dim(ctx_in.mesh.Dimension()), ne(ctx_in.nentities),
    nq(ctx_in.ir.GetNPoints()), q1d(tensor_1d_size(nq, dim)),
@@ -175,53 +275,62 @@ public:
       });
       return itod;
    }()),
-   Ye_mem()
+   group_Ye_mem()
    {
       MFEM_ASSERT(ctx.unionfds.size() == nfields,
                   "LocalQFBackend: unionfds size mismatch");
       MFEM_ASSERT(
          trial_field_uf != SIZE_MAX,
          "DerivativeAssembleDiagonal: trial field not found in unionfds");
-      MFEM_ASSERT(
-         test_field_uf != SIZE_MAX,
-         "DerivativeAssembleDiagonal: test field not found in unionfds");
       MFEM_ASSERT(trial_vdim > 0,
                   "LocalQFBackend: could not determine trial vdim");
       MFEM_ASSERT(total_trial_op_dim > 0,
                   "LocalQFBackend: no dependent inputs found");
       for_constexpr<n_outputs>([&](auto o)
       {
-         MFEM_CONTRACT_VAR(o);
-         MFEM_ASSERT(out_vdim[o] == test_vdim,
-                     "DerivativeAssembleDiagonal: all outputs must share the "
-                     "test field vdim");
+         MFEM_ASSERT(out_vdim[o] == group_test_vdim[out_group[o]],
+                     "DerivativeAssembleDiagonal: outputs on one field must "
+                     "share its vdim");
       });
 
-      if (is_square)
+      group_Ye_mem.resize(group_field_ids.size());
+      for (size_t g = 0; g < group_Ye_mem.size(); g++)
       {
-         Ye_mem.SetSize(num_test_dof * test_vdim * ne);
-         Ye_mem.UseDevice(true);
+         if (!group_has_diagonal[g]) { continue; }
+         group_Ye_mem[g].SetSize(group_num_test_dof[g] * group_test_vdim[g] *
+                                 ne);
+         group_Ye_mem[g].UseDevice(true);
       }
    }
 
+   /// Index of the row block for output field @a field_id, or -1.
+   int FindGroup(int field_id) const
+   {
+      const auto &ids = group_field_ids;
+      const auto it = std::find(ids.begin(), ids.end(), field_id);
+      return (it == ids.end()) ? -1 : static_cast<int>(it - ids.begin());
+   }
+
    template<typename Backend>
-   void run_kernels() const
+   void run_kernels(const int g) const
    {
       Backend::Run(dim,
                    q1d,
                    ctx,
                    qp_cache,
-                   Ye_mem,
+                   group_Ye_mem[g],
                    inputs,
                    outputs,
                    output_dtq_maps,
                    input_dtq_maps,
-                   test_vdim,
+                   out_group,
+                   g,
+                   group_test_vdim[g],
                    out_op_dim,
                    out_offsets,
                    output_size_on_qp,
-                   num_test_dof,
-                   num_test_dof_1d,
+                   group_num_test_dof[g],
+                   group_num_test_dof_1d[g],
                    trial_vdim,
                    total_trial_op_dim,
                    residual_size_on_qp,
@@ -232,9 +341,14 @@ public:
                    dim);
    }
 
-   void operator()(Vector &diag_e) const
+   /// Add this integrator's contribution to the diagonal of the row block of
+   /// output field @a out_field_id. Adds nothing if the integrator writes no
+   /// square, basis-backed block for that field; the caller is responsible for
+   /// rejecting a row that no integrator can serve.
+   void operator()(const int out_field_id, Vector &diag_e) const
    {
-      if (!is_square) { return; }
+      const int g = FindGroup(out_field_id);
+      if (g < 0 || !group_has_diagonal[g]) { return; }
       if (ctx.attr.Size() == 0) { return; }
 
       if (!(use_sum_factorization && (dim == 2 || dim == 3)))
@@ -242,27 +356,28 @@ public:
          MFEM_ABORT("DerivativeAssembleDiagonal optimized path is implemented "
                     "for tensor-product 2D/3D elements only");
       }
-      MFEM_VERIFY(num_test_dof_1d == num_trial_dof_1d,
+      MFEM_VERIFY(group_num_test_dof_1d[g] == num_trial_dof_1d,
                   "DerivativeAssembleDiagonal requires matching tensor dofs");
-      MFEM_VERIFY(num_test_dof_1d <= DeviceDofQuadLimits::Get().MAX_D1D, "");
-      MFEM_VERIFY(q1d <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+      const auto &limits = DeviceDofQuadLimits::Get();
+      MFEM_VERIFY(group_num_test_dof_1d[g] <= limits.MAX_D1D, "");
+      MFEM_VERIFY(q1d <= limits.MAX_Q1D, "");
 
-      Ye_mem = 0.0;
+      group_Ye_mem[g] = 0.0;
 
       if (q1d <= LocalQFLOBackendMQ1())
       {
-         run_kernels<DerivativeAssembleDiagonalLO>();
+         run_kernels<DerivativeAssembleDiagonalLO>(g);
       }
       else if (q1d <= LocalQFHOBackendMQ1())
       {
-         run_kernels<DerivativeAssembleDiagonalHO>();
+         run_kernels<DerivativeAssembleDiagonalHO>(g);
       }
       else
       {
          MFEM_ABORT("Unsupported quadrature order for LocalQF backend");
       }
 
-      diag_e += Ye_mem;
+      diag_e += group_Ye_mem[g];
    }
 
    template<typename backend_t = LocalQFLOBackend<3>, int T_Q1D = 0>
@@ -274,6 +389,8 @@ public:
       const outputs_t &outputs,
       const std::array<DofToQuadMap, n_outputs> &output_dtq_maps,
       const std::array<DofToQuadMap, n_inputs> &input_dtq_maps,
+      const std::array<int, n_outputs> &out_group,
+      const int row_group,
       const int test_vdim,
       const std::array<int, n_outputs> &out_op_dim,
       const std::array<int, n_outputs> &out_offsets,
@@ -339,85 +456,92 @@ public:
             }
             MFEM_SYNC_THREAD;
 
-            // Accumulate every (output o, test op k, dependent input s,
-            // trial op m) block of the cached Jacobian into the diagonal via
-            // the backend driver.
+            // Accumulate every output belonging to the requested row block.
+            // This sums multiple contributions, such as Value<U> +
+            // Gradient<U>, while skipping outputs on the other row blocks.
+            // The row is a run time choice, so unlike the field id it cannot
+            // gate the instantiation; is_identity_fop_v still does, since
+            // eval_test has no meaning for quadrature point data.
             for_constexpr<n_outputs>([&](auto o)
             {
                using test_fop_t = std::decay_t<decltype(get<o>(outputs))>;
-               const auto &out_dtq = output_dtq_maps[o];
-               const int test_op_dim = out_op_dim[static_cast<int>(o)];
-
-               // Test-basis factor along a spatial axis
-               const auto eval_test =
-                  [&](const int k, const int axis, const int q, const int d)
+               if constexpr (!is_identity_fop_v<test_fop_t>)
                {
-                  const auto &B = out_dtq.B;
-                  const auto &G = out_dtq.G;
-                  if constexpr (is_value_fop<test_fop_t>::value)
-                  {
-                     return (k == 0) ? B(q, 0, d) : 0.0;
-                  }
-                  else if constexpr (is_gradient_fop<test_fop_t>::value)
-                  {
-                     return (k == axis) ? G(q, 0, d) : B(q, 0, d);
-                  }
-                  else
-                  {
-                     return 0.0;
-                  }
-               };
+                  if (out_group[static_cast<int>(o)] != row_group) { return; }
+                  const auto &out_dtq = output_dtq_maps[o];
+                  const int test_op_dim = out_op_dim[static_cast<int>(o)];
 
-               for (int k = 0; k < test_op_dim; k++)
-               {
-                  const int row =
-                     out_offsets[static_cast<int>(o)] + vd * test_op_dim + k;
-                  int m_offset = 0;
-                  for_constexpr<n_inputs>([&](auto s)
+                  // Test-basis factor along a spatial axis
+                  const auto eval_test =
+                     [&](const int k, const int axis, const int q, const int d)
                   {
-                     using fop_t = std::decay_t<decltype(get<s>(inputs))>;
-                     const int trial_op_dim =
-                        inputs_trial_op_dim[static_cast<int>(s)];
-                     if (trial_op_dim == 0) { return; }
-
-                     const auto &in_dtq = input_dtq_maps[s];
-                     const auto eval_input =
-                        [&](const int m, const int axis, const int q,
-                            const int d)
+                     const auto &B = out_dtq.B;
+                     const auto &G = out_dtq.G;
+                     if constexpr (is_value_fop<test_fop_t>::value)
                      {
-                        if constexpr (is_value_fop<fop_t>::value)
-                        {
-                           return (m == 0) ? in_dtq.B(q, 0, d) : 0.0;
-                        }
-                        else if constexpr (is_gradient_fop<fop_t>::value)
-                        {
-                           return (m == axis) ? in_dtq.G(q, 0, d)
-                                  : in_dtq.B(q, 0, d);
-                        }
-                        else
-                        {
-                           return 0.0;
-                        }
-                     };
-
-                     for (int m = 0; m < trial_op_dim; m++)
-                     {
-                        const int col = m_offset + m;
-                        backend_t::DiagContract(
-                           s_diag,
-                           num_test_dof_1d,
-                           q1d,
-                           nz_dof,
-                           [&](int axis, int q, int d)
-                        { return eval_test(k, axis, q, d); },
-                        [&](int axis, int q, int d)
-                        { return eval_input(m, axis, q, d); },
-                        [&](int q) { return qpdc(q, col, vd, row); },
-                        [&](int dx, int dy, int dz, real_t u)
-                        { Y(dx, dy, dz) += u; });
+                        return (k == 0) ? B(q, 0, d) : 0.0;
                      }
-                     m_offset += trial_op_dim;
-                  });
+                     else if constexpr (is_gradient_fop<test_fop_t>::value)
+                     {
+                        return (k == axis) ? G(q, 0, d) : B(q, 0, d);
+                     }
+                     else
+                     {
+                        return 0.0;
+                     }
+                  };
+
+                  for (int k = 0; k < test_op_dim; k++)
+                  {
+                     const int row = out_offsets[static_cast<int>(o)] +
+                                     vd * test_op_dim + k;
+                     int m_offset = 0;
+                     for_constexpr<n_inputs>([&](auto s)
+                     {
+                        using fop_t = std::decay_t<decltype(get<s>(inputs))>;
+                        const int trial_op_dim =
+                           inputs_trial_op_dim[static_cast<int>(s)];
+                        if (trial_op_dim == 0) { return; }
+
+                        const auto &in_dtq = input_dtq_maps[s];
+                        const auto eval_input =
+                           [&](const int m, const int axis, const int q,
+                               const int d)
+                        {
+                           if constexpr (is_value_fop<fop_t>::value)
+                           {
+                              return (m == 0) ? in_dtq.B(q, 0, d) : 0.0;
+                           }
+                           else if constexpr (is_gradient_fop<fop_t>::value)
+                           {
+                              return (m == axis) ? in_dtq.G(q, 0, d)
+                                     : in_dtq.B(q, 0, d);
+                           }
+                           else
+                           {
+                              return 0.0;
+                           }
+                        };
+
+                        for (int m = 0; m < trial_op_dim; m++)
+                        {
+                           const int col = m_offset + m;
+                           backend_t::DiagContract(
+                              s_diag,
+                              num_test_dof_1d,
+                              q1d,
+                              nz_dof,
+                              [&](int axis, int q, int d)
+                           { return eval_test(k, axis, q, d); },
+                           [&](int axis, int q, int d)
+                           { return eval_input(m, axis, q, d); },
+                           [&](int q) { return qpdc(q, col, vd, row); },
+                           [&](int dx, int dy, int dz, real_t u)
+                           { Y(dx, dy, dz) += u; });
+                        }
+                        m_offset += trial_op_dim;
+                     });
+                  }
                }
             });
          }

fem/dfem/backends/local_qf/derivative_setup.hpp

@@ -384,6 +384,8 @@ public:
                            for_constexpr<n_outputs>([&](auto oc)
                            {
                               constexpr size_t o = oc.value, ao = n_inputs + o;
+                              using out_fop_t =
+                                 std::decay_t<tuple_element_t<o, outputs_t>>;
                               const auto &tangent = get<ao>(shadow_args);
                               const int tv = out_vdim[o], to = out_op_dim[o];
                               for (int i = 0; i < tv; i++)
@@ -394,8 +396,26 @@ public:
                                     const int cache_idx =
                                        row * trial_vdim * total_trial_op_dim +
                                        j * total_trial_op_dim + col_m;
-                                    cache_tensor(q, cache_idx, e) =
-                                       qf_value_at(tangent, i, k);
+                                    // An Identity output is flat quadrature
+                                    // point data: its FieldOperator vdim counts
+                                    // components, not rows of the q-function
+                                    // argument's shape. The two index form
+                                    // assumes vdim == extents[0] and would run
+                                    // off the end of, say, a
+                                    // tensor<real_t, DIM, DIM> bound to a vdim
+                                    // DIM*DIM space, so read it flat with the
+                                    // same column major packing that
+                                    // identity_qp_write_value writes.
+                                    if constexpr (is_identity_fop_v<out_fop_t>)
+                                    {
+                                       cache_tensor(q, cache_idx, e) =
+                                          qf_flat_value(tangent, i * to + k);
+                                    }
+                                    else
+                                    {
+                                       cache_tensor(q, cache_idx, e) =
+                                          qf_value_at(tangent, i, k);
+                                    }
                                  }
                               }
                            });

fem/dfem/doperator.hpp

@@ -109,19 +109,29 @@ struct derivative_action_t
    std::function<void *()> qfunc_shadow;
 };
 
-/// @brief Type alias for a function that assembles the SparseMatrix of a
-/// derivative operator
+/// @brief Type alias for a function that assembles the SparseMatrix row blocks
+/// of a derivative operator
+///
+/// The derivative of a residual with several output (test) fields is a block
+/// column, one row block per output field, so the callback fills a vector
+/// indexed by position in the operator's output FieldDescriptors rather than a
+/// single matrix.
 using assemble_derivative_sparsematrix_callback_t =
-   std::function<void(SparseMatrix *&)>;
+   std::function<void(std::vector<SparseMatrix *> &)>;
 
-/// @brief Type alias for a function that assembles the HypreParMatrix of a
-/// derivative operator
+/// @brief Type alias for a function that assembles the HypreParMatrix row
+/// blocks of a derivative operator
+///
+/// @see assemble_derivative_sparsematrix_callback_t for the indexing.
 using assemble_derivative_hypreparmatrix_callback_t =
-   std::function<void(HypreParMatrix *&)>;
+   std::function<void(std::vector<HypreParMatrix *> &)>;
 
-/// @brief Type alias for a function that assembles the diagonal of a derivative
-/// operator into an E-vector
-using assemble_diagonal_callback_t = std::function<void(Vector &)>;
+/// @brief Type alias for a function that assembles the diagonal of one row
+/// block of a derivative operator into an E-vector
+///
+/// The first argument names the output (test) field whose row block is wanted;
+/// only a square block has a diagonal at all.
+using assemble_diagonal_callback_t = std::function<void(int, Vector &)>;
 
 /// @brief Type alias for a function that applies the appropriate restriction to
 /// the solution and parameters
@@ -248,60 +258,73 @@ MakeDerivativeHypreParMatrixAssemble(
    std::vector<assemble_derivative_sparsematrix_callback_t> &sparse_callbacks,
    const IntegratorContext &ctx)
 {
-   return [derivative_idx, &sparse_callbacks, ctx](HypreParMatrix *&A)
+   using blocks_t = std::vector<HypreParMatrix *>;
+   return [derivative_idx, &sparse_callbacks, ctx](blocks_t &A)
    {
-      MFEM_VERIFY(ctx.outfds.size() == 1,
-                  "HypreParMatrix assembly requires a single output field");
-
       const size_t trial_field_idx = FindIdx(derivative_idx, ctx.unionfds);
       MFEM_VERIFY(trial_field_idx != SIZE_MAX,
                   "derivative field not found for HypreParMatrix assembly");
 
-      const auto *test_fes_ptr =
-         std::get_if<const ParFiniteElementSpace *>(&ctx.outfds[0].data);
       const auto *trial_fes_ptr =
          std::get_if<const ParFiniteElementSpace *>(
             &ctx.unionfds[trial_field_idx].data);
-      MFEM_VERIFY(test_fes_ptr && *test_fes_ptr,
-                  "HypreParMatrix assembly requires a ParFiniteElementSpace "
-                  "output field");
       MFEM_VERIFY(trial_fes_ptr && *trial_fes_ptr,
                   "HypreParMatrix assembly requires a ParFiniteElementSpace "
                   "derivative field");
-
-      const ParFiniteElementSpace *test_fes = *test_fes_ptr;
       const ParFiniteElementSpace *trial_fes = *trial_fes_ptr;
-      MFEM_VERIFY(test_fes->GetComm() == trial_fes->GetComm(),
-                  "test and trial spaces must use the same MPI communicator");
 
-      SparseMatrix *spmat = nullptr;
+      // Every integrator fills the row blocks it contributes to; the local
+      // blocks are then RAP'd one by one, each with its own test space.
+      std::vector<SparseMatrix *> spmat(ctx.outfds.size(), nullptr);
       for (const auto &f : sparse_callbacks)
       {
          f(spmat);
       }
 
-      MFEM_VERIFY(spmat != nullptr,
+      bool any = false;
+      for (const auto *m : spmat) { any = any || (m != nullptr); }
+      MFEM_VERIFY(any,
                   "internal error: sparse derivative assembly returned NULL");
-      MFEM_VERIFY(spmat->Finalized(),
-                  "local derivative matrix must be finalized");
 
-      if (test_fes == trial_fes)
+      if (A.size() < ctx.outfds.size())
       {
-         HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
-                           test_fes->GetDofOffsets(), spmat);
-         A = RAP(&dA, test_fes->Dof_TrueDof_Matrix());
-      }
-      else
-      {
-         HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
-                           trial_fes->GlobalVSize(),
-                           test_fes->GetDofOffsets(),
-                           trial_fes->GetDofOffsets(), spmat);
-         A = RAP(test_fes->Dof_TrueDof_Matrix(), &dA,
-                 trial_fes->Dof_TrueDof_Matrix());
+         A.resize(ctx.outfds.size(), nullptr);
       }
 
-      delete spmat;
+      for (size_t i = 0; i < spmat.size(); i++)
+      {
+         if (spmat[i] == nullptr) { continue; }
+         MFEM_VERIFY(spmat[i]->Finalized(),
+                     "local derivative matrix must be finalized");
+
+         const auto *test_fes_ptr =
+            std::get_if<const ParFiniteElementSpace *>(&ctx.outfds[i].data);
+         MFEM_VERIFY(test_fes_ptr && *test_fes_ptr,
+                     "HypreParMatrix assembly requires a ParFiniteElementSpace "
+                     "output field");
+         const ParFiniteElementSpace *test_fes = *test_fes_ptr;
+         MFEM_VERIFY(test_fes->GetComm() == trial_fes->GetComm(),
+                     "test and trial spaces must use the same "
+                     "MPI communicator");
+
+         if (test_fes == trial_fes)
+         {
+            HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
+                              test_fes->GetDofOffsets(), spmat[i]);
+            A[i] = RAP(&dA, test_fes->Dof_TrueDof_Matrix());
+         }
+         else
+         {
+            HypreParMatrix dA(test_fes->GetComm(), test_fes->GlobalVSize(),
+                              trial_fes->GlobalVSize(),
+                              test_fes->GetDofOffsets(),
+                              trial_fes->GetDofOffsets(), spmat[i]);
+            A[i] = RAP(test_fes->Dof_TrueDof_Matrix(), &dA,
+                       trial_fes->Dof_TrueDof_Matrix());
+         }
+
+         delete spmat[i];
+      }
    };
 }
 
@@ -563,17 +586,64 @@ public:
       }
    }
 
+   /// @brief Assemble the row blocks of the derivative operator into
+   /// SparseMatrices.
+   ///
+   /// A q-function may write to several output (test) fields. Differentiating
+   /// w.r.t. one field then gives a block column with one row block per output
+   /// field,
+   ///
+   ///     dR/dU = [ dR_U/dU ; dR_Y/dU ],
+   ///
+   /// all sharing the trial space of the differentiated field. Each row block
+   /// is an ordinary test x trial matrix and is assembled separately; the
+   /// stacked matrix, if wanted, is one HypreParMatrixFromBlocks call away.
+   ///
+   /// @param A Resized to the number of output fields and indexed the same way,
+   /// so A[f] is the row block of output field f. Fields whose rows cannot be
+   /// materialised -- quadrature and parameter spaces, which have no basis to
+   /// contract against -- are left as nullptr. Ownership passes to the caller.
+   void Assemble(std::vector<SparseMatrix *> &A)
+   {
+      MFEM_VERIFY(!assemble_derivative_sparsematrix_callbacks.empty(),
+                  "derivative can't be assembled into a SparseMatrix");
+      EnsureQpCache();
+
+      A.assign(outfds.size(), nullptr);
+      for (const auto &f : assemble_derivative_sparsematrix_callbacks)
+      {
+         f(A);
+      }
+   }
+
    /// @brief Assemble the derivative operator into a SparseMatrix.
    ///
+   /// Convenience overload for the common case of a derivative with a single
+   /// assemblable row block. Use Assemble(std::vector<SparseMatrix *> &) when
+   /// the q-function writes to more than one test field.
+   ///
    /// @param A The SparseMatrix to assemble the derivative operator into. Can
    /// be an uninitialized object.
    void Assemble(SparseMatrix *&A)
    {
-      MFEM_ASSERT(!assemble_derivative_sparsematrix_callbacks.empty(),
-                  "derivative can't be assembled into a SparseMatrix");
+      std::vector<SparseMatrix *> blocks;
+      Assemble(blocks);
+      A = SingleBlock(blocks);
+   }
+
+   /// @brief Assemble the row blocks of the derivative operator into
+   /// HypreParMatrices.
+   ///
+   /// @see Assemble(std::vector<SparseMatrix *> &) for the indexing and the
+   /// block structure.
+   void Assemble(std::vector<HypreParMatrix *> &A)
+   {
+      MFEM_VERIFY(!assemble_derivative_hypreparmatrix_callbacks.empty(),
+                  "derivative can't be assembled into a HypreParMatrix");
       EnsureQpCache();
 
-      for (const auto &f : assemble_derivative_sparsematrix_callbacks)
+      A.assign(outfds.size(), nullptr);
+      for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
       {
          f(A);
       }
@@ -581,18 +651,17 @@ public:
 
    /// @brief Assemble the derivative operator into a HypreParMatrix.
    ///
+   /// Convenience overload for the common case of a derivative with a single
+   /// assemblable row block. Use Assemble(std::vector<HypreParMatrix *> &) when
+   /// the q-function writes to more than one test field.
+   ///
    /// @param A The HypreParMatrix to assemble the derivative operator into. Can
    /// be an uninitialized object.
    void Assemble(HypreParMatrix *&A)
    {
-      MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
-                  "derivative can't be assembled into a HypreParMatrix");
-      EnsureQpCache();
-
-      for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
-      {
-         f(A);
-      }
+      std::vector<HypreParMatrix *> blocks;
+      Assemble(blocks);
+      A = SingleBlock(blocks);
    }
 
    /// @brief Assemble the derivative of a functional into a Vector.
@@ -635,35 +704,89 @@ public:
       }
    }
 
-   /// @brief Assemble the diagonal of the derivative operator into a T-vector.
+   /// @brief Assemble the diagonal of one row block of the derivative
+   /// operator into a T-vector.
    ///
-   /// @param diag The vector to receive the diagonal (must be T-dof sized).
-   void AssembleDiagonal(Vector &diag) const override
+   /// The derivative is a block column, one row block per output field, all
+   /// sharing the trial space of the differentiated field. Only a block whose
+   /// test space *is* that trial space is square, and only a square block has a
+   /// diagonal. Squareness alone cannot pick the block when several output
+   /// fields live on the trial space, so the row is named rather than inferred.
+   /// Output field ids may differ from the differentiated field's id while
+   /// sharing its ParFiniteElementSpace, in which case several row blocks are
+   /// square at once. Naming the output field id resolves the ambiguity.
+   ///
+   /// @param out_field_id The output field whose row block is wanted.
+   /// @param diag The vector to receive the diagonal, resized to that field's
+   /// T-dof size.
+   void AssembleDiagonal(size_t out_field_id, Vector &diag) const
    {
-      MFEM_ASSERT(!assemble_diagonal_callbacks.empty(),
+      MFEM_VERIFY(!assemble_diagonal_callbacks.empty(),
                   "derivative can't assemble diagonal");
       EnsureQpCache();
-      MFEM_ASSERT(outfds.size() == 1,
-                  "AssembleDiagonal currently requires a single output field");
+
+      const size_t diagonal_idx = FindIdx(out_field_id, outfds);
+      MFEM_VERIFY(diagonal_idx != SIZE_MAX,
+                  "AssembleDiagonal: field " << out_field_id << " is not an "
+                  "output field of this operator");
 
       const auto *test_pf =
-         std::get_if<const ParFiniteElementSpace *>(&outfds[0].data);
+         std::get_if<const ParFiniteElementSpace *>(&outfds[diagonal_idx].data);
       MFEM_VERIFY(test_pf && *test_pf,
                   "AssembleDiagonal: test field must be a ParFiniteElementSpace");
+      const auto *trial_pf =
+         std::get_if<const ParFiniteElementSpace *>(&direction.data);
+      MFEM_VERIFY(trial_pf && *trial_pf,
+                  "AssembleDiagonal: the differentiated field must be a "
+                  "ParFiniteElementSpace");
+      MFEM_VERIFY(*test_pf == *trial_pf,
+                  "AssembleDiagonal: the requested row block is not square and "
+                  "so has no diagonal; its test space is not the trial space "
+                  "of the differentiated field");
 
       prepare_residual(outfds, out_rcache, daction_e);
       for (auto *v : daction_e) { *v = 0.0; }
 
       for (const auto &f : assemble_diagonal_callbacks)
       {
-         f(*daction_e[0]);
+         f(static_cast<int>(out_field_id), *daction_e[diagonal_idx]);
       }
 
       restriction_transpose(outfds, out_rcache, daction_e, daction_l);
-      prolongation_transpose(outfds[0], *daction_l[0], diag);
+      prolongation_transpose(outfds[diagonal_idx],
+                             *daction_l[diagonal_idx], diag);
+   }
+
+   /// @brief Assemble the diagonal of the derivative operator into a T-vector.
+   ///
+   /// Uses the row block of the differentiated field, i.e. the usual diagonal
+   /// block dR_U/dU for GetDerivative(U). Call
+   /// AssembleDiagonal(size_t, Vector &) for any other row block.
+   ///
+   /// @param diag The vector to receive the diagonal.
+   void AssembleDiagonal(Vector &diag) const override
+   {
+      AssembleDiagonal(direction.id, diag);
    }
 
 private:
+   /// The single non-null row block of @a blocks, for the scalar Assemble
+   /// overloads.
+   template <typename mat_t>
+   static mat_t *SingleBlock(const std::vector<mat_t *> &blocks)
+   {
+      mat_t *single = nullptr;
+      int count = 0;
+      for (auto *m : blocks)
+      {
+         if (m != nullptr) { single = m; count++; }
+      }
+      MFEM_VERIFY(count == 1,
+                  "the derivative has " << count << " assemblable row blocks, "
+                  "not one; assemble it into a std::vector instead");
+      return single;
+   }
+
    /// Derivative action callbacks. Depending on the requested derivatives in
    /// DifferentiableOperator the callbacks represent certain combinations of
    /// actions of derivatives of the forward operator.
@@ -1479,16 +1602,23 @@ void DifferentiableOperator::AddIntegrator(
          constexpr size_t derivative_idx = decltype(derivative_id)::value;
          using callback_outputs_t = std::decay_t<decltype(outputs)>;
 
-         bool disable_assemble = false;
+
+         // NOTE: before disable_assemble was looking for any identity outputs,
+         // But this would silently disable assembly for other valid outputs across
+         // different fields.
+         bool any_assemblable_output = false;
          for_constexpr([&](auto j)
          {
-            using output_fop_t = tuple_element_t<j, callback_outputs_t>;
-            if constexpr (is_identity_fop_v<std::decay_t<output_fop_t>>)
+            using output_fop_t =
+               std::decay_t<tuple_element_t<j, callback_outputs_t>>;
+            if constexpr (!is_identity_fop_v<output_fop_t>)
             {
-               disable_assemble = true;
+               any_assemblable_output = true;
             }
          }, std::make_index_sequence<tuple_size<callback_outputs_t>::value> {});
 
+         const bool disable_assemble = !any_assemblable_output;
+
          // Setup the qp cache for the derivative
          setup_callbacks[callback_key].push_back(
             MakeDerivativeSetupCallback(
@@ -1514,13 +1644,21 @@ void DifferentiableOperator::AddIntegrator(
                backend_t::template MakeDerivativeAssemble<derivative_idx>(
                   callback_ctx, qf, inputs, outputs, callback_qp_cache));
 
-            // Assemble the derivative into a HypreParMatrix
-            assemble_hypreparmatrix_callbacks[callback_key].push_back(
-               MakeDerivativeHypreParMatrixAssemble(
-                  derivative_idx,
-                  assemble_sparsematrix_callbacks[callback_key],
-                  callback_ctx));
+            // Assemble the derivative into a HypreParMatrix. This one runs
+            // every sparse callback registered under the key, so it is
+            // registered once and not once per integrator.
+            if (assemble_hypreparmatrix_callbacks[callback_key].empty())
+            {
+               assemble_hypreparmatrix_callbacks[callback_key].push_back(
+                  MakeDerivativeHypreParMatrixAssemble(
+                     derivative_idx,
+                     assemble_sparsematrix_callbacks[callback_key],
+                     callback_ctx));
+            }
+         }
 
+         if (!disable_assemble)
+         {
             // Assemble the diagonal of the derivative into an L-vector
             assemble_diagonal_cbs[callback_key].push_back(
                backend_t::template MakeDerivativeAssembleDiagonal<derivative_idx>(

miniapps/dfem/CMakeLists.txt

@@ -14,8 +14,20 @@ if (MFEM_USE_MPI)
         MAIN dfem-minimal-surface.cpp
         LIBRARIES mfem)
 
-    add_mfem_miniapp(dfem-hyperelasticity-energy
-        MAIN dfem-hyperelasticity_energy.cpp
+    add_mfem_miniapp(dfem-hyperelasticity
+        MAIN dfem-hyperelasticity.cpp
+        LIBRARIES mfem)
+
+    add_mfem_miniapp(dfem-nonlinear-poisson
+        MAIN dfem-nonlinear-poisson.cpp
+        LIBRARIES mfem)
+
+    add_mfem_miniapp(dfem-hyperbolic-heat
+        MAIN dfem-hyperbolic-heat.cpp
+        LIBRARIES mfem)
+
+    add_mfem_miniapp(dfem-reaction-diffusion
+        MAIN dfem-reaction-diffusion.cpp
         LIBRARIES mfem)
 
     # Add parallel tests.
@@ -25,9 +37,35 @@ if (MFEM_USE_MPI)
         ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-minimal-surface> -der 0 -o 1 -r 2 -no-vis
         ${MPIEXEC_POSTFLAGS})
 
-        add_test(NAME dfem-hyperelasticity-energy=${MFEM_MPI_NP}
+        add_test(NAME dfem-hyperelasticity=${MFEM_MPI_NP}
         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperelasticity-energy> -o 1 -rs 0 -no-vis
+        ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperelasticity> -o 1 -rs 0 -no-vis
         ${MPIEXEC_POSTFLAGS})
+
+        add_test(NAME dfem-nonlinear-poisson=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-nonlinear-poisson>
+        -k grad -f energy -o 1 -r 2 -no-vis
+        ${MPIEXEC_POSTFLAGS})
+
+        add_test(NAME dfem-nonlinear-poisson-residual=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-nonlinear-poisson>
+        -k sol -f residual -o 1 -r 2 -no-vis
+        ${MPIEXEC_POSTFLAGS})
+
+        add_test(NAME dfem-reaction-diffusion=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-reaction-diffusion>
+        -b 1 -o 1 -r 2 -pc 2 -no-vis
+        ${MPIEXEC_POSTFLAGS})
+
+        foreach(AL partial legacy full)
+            add_test(NAME dfem-hyperbolic-heat-${AL}=${MFEM_MPI_NP}
+            COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+            ${MPIEXEC_PREFLAGS} $<TARGET_FILE:dfem-hyperbolic-heat>
+            -al ${AL} -o 1 -r 0 -tf 0.1 -dt 0.05 -no-vis
+            ${MPIEXEC_POSTFLAGS})
+        endforeach()
     endif()
 endif()

miniapps/dfem/dfem-hyperbolic-heat.cpp

@@ -0,0 +1,839 @@
+// 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.
+//
+//                --------------------------------------------
+//                Nonlinear Hyperbolic Heat Equation with dFEM
+//                --------------------------------------------
+//
+// Compile with: make dfem-hyperbolic-heat
+//
+// Sample runs:  mpirun -np 4 dfem-hyperbolic-heat
+//               mpirun -np 4 dfem-hyperbolic-heat -b 0
+//               mpirun -np 4 dfem-hyperbolic-heat -tau 0.02
+//               mpirun -np 4 dfem-hyperbolic-heat -q 5 -rc 4
+//               mpirun -np 4 dfem-hyperbolic-heat -r 1 -tf 0.6
+//               mpirun -np 4 dfem-hyperbolic-heat -dt 0.2 -tf 2
+//               mpirun -np 4 dfem-hyperbolic-heat -al legacy
+//
+// Description:  This miniapp solves the hyperbolic (Maxwell-Cattaneo) heat
+//               equation with a nonlinear conductivity,
+//
+//                 tau d^2T/dt^2 + rho c dT/dt = div( k(T) grad T ) + Q,
+//
+//               on the unit cube (set the compile-time constant dim to 2 for
+//               the unit square), starting from a Gaussian hot spot at rest.
+//               The relaxation term tau d^2T/dt^2 turns the parabolic heat
+//               equation into a damped wave equation to comply with the
+//               finite speed sqrt(k/tau) of heat propagation.
+//               Small tau recovers diffusive behavior (Fourier's law),
+//               large tau a wave-like behavior.
+//
+//               The conductivity k(T) = k0 (1 + beta T^2) is the only
+//               nonlinearity. It is written as a dFEM q-function, so both the
+//               diffusion residual
+//
+//                 F(T) = int k(T) grad T . grad w dx
+//
+//               and its tangent J_F(T) come from the same pointwise
+//               description; nothing here differentiates k(T) by hand.
+//
+//               We recommend viewing examples 16 and 23, and dfem miniapps
+//               before viewing this example.
+
+#include "mfem.hpp"
+#include "../../fem/dfem/doperator.hpp"
+#include "../../fem/dfem/backends/local_qf/prelude.hpp"
+
+#include <memory>
+
+using namespace mfem;
+
+// This example code demonstrates the use of new features in MFEM that are in
+// development but exposed through the mfem::future namespace. All features
+// under this namespace might change their interface or behavior in upcoming
+// releases until they have stabilized.
+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
+
+// Space dimension. The q-function and the diffusion operator are templated on
+// it, so it is fixed at compile time; set it to 2 to run the same problem on
+// the unit square.
+constexpr int dim = 3;
+
+// Field IDs used by the dFEM integrator.
+static constexpr int Temperature = 0;
+static constexpr int Coords = 1;
+
+enum class PreconditionerType { None, Diagonal };
+
+/// Map an AssemblyLevel name onto the MFEM assembly level used for the mass operators.
+/// Only partial assembly keeps the whole solve matrix free; the others are here
+/// so the cost of that choice can be measured against the alternatives.
+static AssemblyLevel ParseAssemblyLevel(const char *name)
+{
+   const std::string s(name);
+   if (s == "partial") { return AssemblyLevel::PARTIAL; }
+   if (s == "legacy")  { return AssemblyLevel::LEGACY; }
+   if (s == "full")    { return AssemblyLevel::FULL; }
+   MFEM_ABORT("Unknown assembly level '" << name << "'. Available levels: "
+              "partial, legacy, full.");
+   return AssemblyLevel::LEGACY;
+}
+
+// ----------------------------------------------------------------------------
+// Pointwise description of the nonlinear diffusion term.
+// ----------------------------------------------------------------------------
+
+/// Quadrature point kernel for F(T) = int k(T) grad T . grad w dx with the
+/// temperature dependent conductivity k(T) = k0 (1 + beta T^2).
+template <int DIM>
+struct NonlinearDiffusionQFunction
+{
+   real_t k0 = 1.0;
+   real_t beta = 1.0;
+
+   MFEM_HOST_DEVICE inline
+   auto operator()(const dscalar_t &T,
+                   const tensor<dscalar_t, DIM> &dTdxi,
+                   const tensor<real_t, DIM, DIM> &J,
+                   const real_t &w,
+                   tensor<dscalar_t, DIM> &dwdx) const
+   {
+      const auto invJ = inv(J);
+      const auto dTdx = dTdxi * invJ;
+      const dscalar_t k = k0 * (1.0_r + beta * T * T);
+      dwdx = k * dTdx * transpose(invJ) * det(J) * w;
+   }
+};
+
+/// dFEM wrapper for the nonlinear diffusion residual F(T) and its tangent
+/// J_F(T).
+///
+/// Neither F nor J_F applies any essential-dof treatment. That belongs to the
+/// operator that assembles the full residual, since the boundary condition
+/// constrains the acceleration, not the diffusion term by itself.
+template <int DIM>
+class NonlinearDiffusionOperator
+{
+public:
+   NonlinearDiffusionOperator(ParFiniteElementSpace &fes,
+                              const IntegrationRule &ir,
+                              real_t k0, real_t beta)
+   {
+      auto &mesh_nodes =
+         *static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
+      mesh_nodes_fes = mesh_nodes.ParFESpace();
+      mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
+
+      const std::vector<FieldDescriptor> inputs =
+      {
+         {Temperature, &fes}, {Coords, mesh_nodes_fes}
+      };
+      // The output operator is the gradient of the test function basis, which
+      // completes the weak form B^T D(B T, B x, w).
+      const std::vector<FieldDescriptor> outputs = {{Temperature, &fes}};
+
+      dop = std::make_shared<DifferentiableOperator>(
+               inputs, outputs, *fes.GetParMesh());
+
+      Array<int> all_domain_attr;
+      if (fes.GetMesh()->attributes.Size() > 0)
+      {
+         all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
+         all_domain_attr = 1;
+      }
+
+      NonlinearDiffusionQFunction<DIM> qf;
+      qf.k0 = k0;
+      qf.beta = beta;
+
+      // Requesting the derivative with respect to Temperature is what makes
+      // J_F available later through GetDerivative.
+      auto derivatives = std::integer_sequence<size_t, Temperature> {};
+      dop->AddDomainIntegrator<LocalQFBackend>(
+         qf,
+         Inputs < Value<Temperature>, Gradient<Temperature>,
+         Gradient<Coords>, Weight > {},
+         Outputs<Gradient<Temperature>> {},
+         ir, all_domain_attr, derivatives);
+   }
+
+   /// F(T), the nonlinear diffusion residual.
+   void Mult(const Vector &T, Vector &F) const
+   {
+      MultiVector X{T, mesh_nodes_tdofs};
+      MultiVector Y{F};
+      dop->Mult(X, Y);
+   }
+
+   /// J_F(T), the tangent of F linearized at @a T. The returned operator
+   /// captures @a T, so it has to be rebuilt whenever the linearization point
+   /// moves.
+   std::shared_ptr<DerivativeOperator> GetGradient(const Vector &T) const
+   {
+      MultiVector X{T, mesh_nodes_tdofs};
+      return dop->GetDerivative(Temperature, X);
+   }
+
+private:
+   ParFiniteElementSpace *mesh_nodes_fes = nullptr;
+   Vector mesh_nodes_tdofs;
+   std::shared_ptr<DifferentiableOperator> dop;
+};
+
+// ----------------------------------------------------------------------------
+// Newton residual for one implicit stage.
+// ----------------------------------------------------------------------------
+///
+/// Residual of the semi-discrete system, written in the acceleration unknown
+/// a = d^2T/dt^2 that MFEM's second order ODE solvers solve for. With the stage
+/// quantities requested by the integrator,
+///
+///    T_stage = T + fac0 a,      v_stage = dT/dt + fac1 a,
+///
+/// the weak form M_tau a + M_rhoc v + F(T) - Q = 0 becomes
+///
+///    R(a) = M_tau a + M_rhoc v_stage + F(T_stage) - Q.
+///
+class ResidualOperator : public Operator
+{
+public:
+   /// Newton tangent
+   ///
+   ///    J_R(a) = M_tau + fac1 M_rhoc + fac0 J_F(T + fac0 a).
+   ///
+   /// None of the three terms is a matrix: the mass operators are partially
+   /// assembled and J_F is a dFEM derivative action, so this class only ever
+   /// composes actions and diagonals.
+   class JacobianOperator : public Operator
+   {
+   public:
+      JacobianOperator(const ResidualOperator &oper, const Vector &T_stage) :
+         Operator(oper.Height()), oper(oper), z(oper.Height()), w(oper.Height())
+      {
+         diffusion_tangent = oper.diffusion.GetGradient(T_stage);
+      }
+
+      void Mult(const Vector &x, Vector &y) const override
+      {
+         // Essential directions are removed before the apply and restored as
+         // identity rows afterwards.
+         z = x;
+         z.SetSubVector(oper.ess_tdofs, 0.0);
+
+         oper.M_tau.Mult(z, y);
+         oper.M_rhoc.AddMult(z, y, oper.fac1);
+
+         MultiVector W{w};
+         diffusion_tangent->Mult(z, W);
+         y.Add(oper.fac0, w);
+
+         auto d_y = y.ReadWrite();
+         const auto d_x = x.Read();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_y[d_dofs[i]] = d_x[d_dofs[i]];
+         });
+      }
+
+      void AssembleDiagonal(Vector &diag) const override
+      {
+         oper.M_tau.AssembleDiagonal(diag);
+
+         oper.M_rhoc.AssembleDiagonal(w);
+         diag.Add(oper.fac1, w);
+
+         diffusion_tangent->AssembleDiagonal(w);
+         diag.Add(oper.fac0, w);
+
+         auto d_diag = diag.ReadWrite();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_diag[d_dofs[i]] = 1.0;
+         });
+      }
+
+   private:
+      const ResidualOperator &oper;
+      mutable Vector z, w;
+      std::shared_ptr<DerivativeOperator> diffusion_tangent;
+   };
+
+   /// @a M_tau and @a M_rhoc are the mass operators as returned by
+   /// ParBilinearForm::FormSystemMatrix, i.e. constrained true-dof operators.
+   ResidualOperator(const Operator &M_tau, const Operator &M_rhoc,
+                    const NonlinearDiffusionOperator<dim> &diffusion,
+                    const Array<int> &ess_tdofs) :
+      Operator(M_tau.Height()),
+      M_tau(M_tau), M_rhoc(M_rhoc), diffusion(diffusion), ess_tdofs(ess_tdofs),
+      T_stage(height), v_stage(height), z(height) { }
+
+   /// Describe the stage the time integrator is asking about. @a T, @a dTdt and
+   /// @a Q are referenced, not copied, so they have to outlive the Newton solve.
+   void SetParameters(real_t fac0_, real_t fac1_, const Vector *T_,
+                      const Vector *dTdt_, const Vector *Q_)
+   {
+      fac0 = fac0_;
+      fac1 = fac1_;
+      T = T_;
+      dTdt = dTdt_;
+      Q = Q_;
+   }
+
+   void Mult(const Vector &a, Vector &R) const override
+   {
+      MFEM_ASSERT(T && dTdt && Q, "call SetParameters() first");
+
+      add(*T, fac0, a, T_stage);
+      add(*dTdt, fac1, a, v_stage);
+
+      M_tau.Mult(a, R);
+      M_rhoc.AddMult(v_stage, R);
+
+      diffusion.Mult(T_stage, z);
+      R += z;
+
+      R -= *Q;
+
+      R.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   Operator &GetGradient(const Vector &a) const override
+   {
+      // The tangent has to be linearized at the stage temperature, not at the
+      // temperature of the previous step.
+      add(*T, fac0, a, T_stage);
+      jacobian = std::make_shared<JacobianOperator>(*this, T_stage);
+      return *jacobian;
+   }
+
+private:
+   const Operator &M_tau;
+   const Operator &M_rhoc;
+   const NonlinearDiffusionOperator<dim> &diffusion;
+   const Array<int> &ess_tdofs;
+
+   real_t fac0 = 0.0, fac1 = 0.0;
+   const Vector *T = nullptr;
+   const Vector *dTdt = nullptr;
+   const Vector *Q = nullptr;
+
+   mutable Vector T_stage, v_stage, z;
+   mutable std::shared_ptr<JacobianOperator> jacobian;
+};
+
+// ----------------------------------------------------------------------------
+// The second order time dependent operator.
+// ----------------------------------------------------------------------------
+
+/// Semi-discrete hyperbolic heat equation
+///
+///     M_tau a + M_rhoc v + F(T) - Q = 0,
+///
+class HyperbolicHeatOperator : public SecondOrderTimeDependentOperator
+{
+public:
+   HyperbolicHeatOperator(ParFiniteElementSpace &fes,
+                          const IntegrationRule &ir,
+                          const Array<int> &ess_bdr,
+                          real_t tau, real_t rho_c, real_t k0, real_t beta,
+                          real_t source, PreconditionerType pc_type,
+                          AssemblyLevel assembly);
+
+   /// Compute a = M_tau^{-1} ( Q - M_rhoc v - F(T) ).
+   using SecondOrderTimeDependentOperator::Mult;
+   void Mult(const Vector &T, const Vector &dTdt,
+             Vector &d2Tdt2) const override;
+
+   /// Solve R(a) = M_tau a + M_rhoc (dTdt + fac1 a) + F(T + fac0 a) - Q = 0 for
+   /// the unknown acceleration a = d2Tdt2.
+   using SecondOrderTimeDependentOperator::ImplicitSolve;
+   void ImplicitSolve(const real_t fac0, const real_t fac1,
+                      const Vector &T, const Vector &dTdt,
+                      Vector &d2Tdt2) override;
+
+   const Array<int> &GetEssentialTrueDofs() const { return ess_tdof_list; }
+
+   /// Newton iterations taken by the last ImplicitSolve.
+   int GetNewtonIterations() const { return newton_solver.GetNumIterations(); }
+
+   /// Compute the physical L2 norm directly from a true-dof vector.
+   real_t ComputeL2Norm(const Vector &x) const;
+
+   ~HyperbolicHeatOperator() override;
+
+private:
+   /// Set up one mass operator weighted by the constant coefficient @a w at the
+   /// requested @a assembly level, and hand back its constrained true-dof form
+   /// in @a M.
+   void SetupMass(ParBilinearForm &m, ConstantCoefficient &w,
+                  const IntegrationRule &ir, AssemblyLevel assembly,
+                  OperatorHandle &M);
+
+   ParFiniteElementSpace &fes;
+   Array<int> ess_tdof_list; // empty for pure Neumann boundary conditions
+
+   // Mass operators, weighted by tau and by rho*c. The forms have to outlive
+
+   // the handles: under partial assembly the operator in the handle refers to
+   // the quadrature point data held by the form.
+   ConstantCoefficient tau_coeff, rhoc_coeff;
+   real_t tau;
+   ParBilinearForm m_tau, m_rhoc;
+   OperatorHandle M_tau, M_rhoc;
+
+   std::unique_ptr<NonlinearDiffusionOperator<dim>> diffusion;
+   Vector Q;
+
+   CGSolver M_solver;
+   OperatorJacobiSmoother M_prec;
+
+   std::unique_ptr<ResidualOperator> residual_op;
+   NewtonSolver newton_solver;
+   std::unique_ptr<IterativeSolver> krylov;
+   std::unique_ptr<Solver> krylov_prec;
+
+   mutable Vector cached_accel;
+   mutable Vector rhs, z;
+};
+
+static std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
+{
+   switch (type)
+   {
+      case PreconditionerType::None:     return nullptr;
+      case PreconditionerType::Diagonal:
+         return std::make_unique<OperatorJacobiSmoother>();
+      default:
+         MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
+   }
+   return nullptr;
+}
+
+void HyperbolicHeatOperator::SetupMass(ParBilinearForm &m,
+                                       ConstantCoefficient &w,
+                                       const IntegrationRule &ir,
+                                       AssemblyLevel assembly,
+                                       OperatorHandle &M)
+{
+   m.SetAssemblyLevel(assembly);
+   m.AddDomainIntegrator(new MassIntegrator(w, &ir));
+   m.Assemble();
+   m.FormSystemMatrix(ess_tdof_list, M);
+}
+
+HyperbolicHeatOperator::HyperbolicHeatOperator(ParFiniteElementSpace &fes,
+                                               const IntegrationRule &ir,
+                                               const Array<int> &ess_bdr,
+                                               real_t tau, real_t rho_c,
+                                               real_t k0, real_t beta,
+                                               real_t source,
+                                               PreconditionerType pc_type,
+                                               AssemblyLevel assembly)
+   : SecondOrderTimeDependentOperator(fes.GetTrueVSize(), (real_t) 0.0),
+     fes(fes),
+     tau_coeff(tau), rhoc_coeff(rho_c),
+     tau(tau),
+     m_tau(&fes), m_rhoc(&fes),
+     M_solver(fes.GetComm()),
+     newton_solver(fes.GetComm()),
+     cached_accel(height), rhs(height), z(height)
+{
+   fes.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+
+   /// Setup the mass operators for the acceleration and velocity terms.
+   SetupMass(m_tau, tau_coeff, ir, assembly, M_tau);
+   SetupMass(m_rhoc, rhoc_coeff, ir, assembly, M_rhoc);
+
+   /// Create the nonlinear diffusion operator as a dFEM operator
+   diffusion = std::make_unique<NonlinearDiffusionOperator<dim>>(fes, ir, k0,
+                                                                 beta);
+
+   /// The source term is a constant coefficient, so we can assemble it once and for all.
+   ConstantCoefficient source_coeff(source);
+   ParLinearForm b(&fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(source_coeff, &ir));
+   b.Assemble();
+   Q.SetSize(height);
+   b.ParallelAssemble(Q);
+
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-10);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(200);
+   M_solver.SetPrintLevel(0);
+   M_solver.SetPreconditioner(M_prec);
+   M_solver.SetOperator(*M_tau);
+
+   /// Set up the nonlinear residual operator used by the Newton solver.
+   residual_op = std::make_unique<ResidualOperator>(*M_tau, *M_rhoc, *diffusion,
+                                                    ess_tdof_list);
+
+   /// Set up the Newton solver and Krylov solver for linearized systems.
+   if (beta == 0.0)
+   {
+      krylov = std::make_unique<CGSolver>(fes.GetComm());
+   }
+   else
+   {
+      krylov = std::make_unique<GMRESSolver>(fes.GetComm());
+   }
+   krylov->SetRelTol(1e-10);
+   krylov->SetAbsTol(0.0);
+   krylov->SetMaxIter(500);
+   krylov->SetPrintLevel(0);
+
+   krylov_prec = MakePreconditioner(pc_type);
+   if (krylov_prec) { krylov->SetPreconditioner(*krylov_prec); }
+
+   newton_solver.SetOperator(*residual_op);
+   newton_solver.SetSolver(*krylov);
+#ifdef MFEM_USE_SINGLE
+   newton_solver.SetRelTol(1e-5);
+   newton_solver.SetAbsTol(1e-8);
+#else
+   newton_solver.SetRelTol(1e-8);
+   newton_solver.SetAbsTol(1e-12);
+#endif
+   newton_solver.SetMaxIter(25);
+   newton_solver.SetPrintLevel(0);
+
+   cached_accel = 0.0;
+}
+
+real_t HyperbolicHeatOperator::ComputeL2Norm(const Vector &x) const
+{
+   M_tau->Mult(x, z);
+   return sqrt(InnerProduct(fes.GetComm(), x, z) / tau);
+}
+
+void HyperbolicHeatOperator::Mult(const Vector &T, const Vector &dTdt,
+                                  Vector &d2Tdt2) const
+{
+   diffusion->Mult(T, rhs);
+   M_rhoc->Mult(dTdt, z);
+   rhs += z;
+   rhs.Neg();
+   rhs += Q;
+
+   rhs.SetSubVector(ess_tdof_list, 0.0);
+   M_solver.Mult(rhs, d2Tdt2);
+   d2Tdt2.SetSubVector(ess_tdof_list, 0.0);
+
+   cached_accel = d2Tdt2;
+}
+
+void HyperbolicHeatOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
+                                           const Vector &T, const Vector &dTdt,
+                                           Vector &d2Tdt2)
+{
+   residual_op->SetParameters(fac0, fac1, &T, &dTdt, &Q);
+
+   d2Tdt2 = cached_accel;
+   d2Tdt2.SetSubVector(ess_tdof_list, 0.0);
+
+   Vector zero;
+   newton_solver.Mult(zero, d2Tdt2);
+   MFEM_VERIFY(newton_solver.GetConverged(), "Nonlinear solve failed.");
+
+   cached_accel = d2Tdt2;
+}
+
+HyperbolicHeatOperator::~HyperbolicHeatOperator() = default;
+
+static real_t InitialTemperature(const Vector &x)
+{
+   real_t r2 = 0.0;
+   for (int d = 0; d < x.Size(); d++)
+   {
+      const real_t s = x(d) - 0.5;
+      r2 += s * s;
+   }
+   return exp(-30.0 * r2);
+}
+
+static real_t InitialRate(const Vector &) { return 0.0; }
+
+
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // 2. Parse command-line options
+   int order = 2;
+   int refinements = 0;
+   const char *device_config = "cpu";
+   int ode_solver_type = 11;
+   real_t t_final = 0.4;
+   real_t dt = 0.02;
+   real_t tau = 1.0;
+   real_t rho_c = 1.0;
+   real_t k0 = 1.0;
+   real_t beta = 1.0;
+   real_t source = 0.0;
+   int prec_type = static_cast<int>(PreconditionerType::Diagonal);
+   const char *assembly_name = "partial";
+   bool dirichlet = true;
+   bool visualization = false;
+   bool paraview = false;
+   int vis_steps = 5;
+   int visport = 19916;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   // FEM
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&refinements, "-r", "--refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&prec_type, "-pc", "--preconditioner",
+                  "Preconditioner for the Newton tangent: 0 = none, "
+                  "1 = diagonal/Jacobi.");
+   args.AddOption(&assembly_name, "-al", "--assembly-level",
+                  "Assembly level of the two mass operators: partial, legacy, or full.");
+   // Time integration
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  SecondOrderODESolver::Types.c_str());
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step", "Time step.");
+   // Parameters/Physics
+   args.AddOption(&tau, "-tau", "--tau",
+                  "Thermal relaxation coefficient of d^2T/dt^2. The wave speed "
+                  "is sqrt(k/tau); tau -> 0 recovers the parabolic limit.");
+   args.AddOption(&rho_c, "-rc", "--rho-c",
+                  "Volumetric heat capacity rho*c, the damping of the wave.");
+   args.AddOption(&k0, "-k", "--conductivity",
+                  "Conductivity k0 in k(T) = k0 (1 + beta T^2).");
+   args.AddOption(&beta, "-b", "--beta",
+                  "Nonlinearity beta in k(T) = k0 (1 + beta T^2). "
+                  "beta = 0 gives a linear, symmetric problem.");
+   args.AddOption(&source, "-q", "--source", "Constant volumetric source Q.");
+   args.AddOption(&dirichlet, "-dir", "--dirichlet", "-neu", "--neumann",
+                  "BC switch: T = 0 or zero flux on the whole boundary.");
+   // Visualization
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or disable ParaView DataCollection output.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
+   args.ParseCheck();
+
+   MFEM_VERIFY(tau > 0.0, "-tau must be positive: the hyperbolic term is what "
+               "makes this a second order problem in time.");
+   MFEM_VERIFY(rho_c >= 0.0, "-rc must be non-negative.");
+   MFEM_VERIFY(k0 > 0.0, "-k must be positive.");
+   MFEM_VERIFY(beta >= 0.0, "-b must be non-negative to keep k(T) positive.");
+
+   const AssemblyLevel assembly = ParseAssemblyLevel(assembly_name);
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   if (Mpi::Root())
+   {
+      mfem::out << "k(T) = " << k0 << " (1 + " << beta << " T^2)"
+                << ", wave speed sqrt(k0/tau) = " << sqrt(k0 / tau)
+                << ", relaxation time tau/(rho c) = "
+                << (rho_c > 0.0 ? tau / rho_c
+                    : std::numeric_limits<real_t>::infinity())
+                << ", tangent: "
+                << (beta == 0.0 ? "symmetric (CG)" : "non-symmetric (GMRES)")
+                << ", mass assembly: " << assembly_name << std::endl;
+   }
+
+   // 4. Define the second order ODE solver used for time integration
+   std::unique_ptr<SecondOrderODESolver> ode_solver(
+      SecondOrderODESolver::Select(ode_solver_type));
+
+   // 5. Create the unit square/cube mesh and refine it
+   Mesh mesh = (dim == 2)
+               ? Mesh::MakeCartesian2D(8, 8, Element::QUADRILATERAL)
+               : Mesh::MakeCartesian3D(8, 8, 8, Element::HEXAHEDRON);
+   for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
+   mesh.SetCurvature(order);
+
+   // 6. Define a parallel mesh. dFEM reads the mesh coordinates as a field, so
+   //    the nodes have to exist.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   pmesh.EnsureNodes();
+
+   // 7. Define a scalar H1 space for the temperature. GlobalTrueVSize is
+   //    collective, so every rank has to reach it.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fes(&pmesh, &fec, 1);
+   const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      mfem::out << "Number of temperature unknowns: " << global_dofs
+                << std::endl;
+   }
+
+   // 8. Select quadrature for the nonlinear weak term
+   //    k(T) grad(T).grad(w). For k(T) = k0 (1 + beta T^2), its
+   //    polynomial degree on affine elements is up to 4*p - 2.
+   const int integration_order =
+      (beta == 0.0) ? 2 * order : 4 * order - 2;
+   const IntegrationRule &ir =
+      IntRules.Get(pmesh.GetTypicalElementGeometry(), integration_order);
+
+   // 9. Set the initial conditions. We assume a gaussian hot spot at rest.
+   ParGridFunction T_gf(&fes), dTdt_gf(&fes);
+
+   FunctionCoefficient T_0(InitialTemperature);
+   T_gf.ProjectCoefficient(T_0);
+   Vector T;
+   T_gf.GetTrueDofs(T);
+
+   FunctionCoefficient dTdt_0(InitialRate);
+   dTdt_gf.ProjectCoefficient(dTdt_0);
+   Vector dTdt;
+   dTdt_gf.GetTrueDofs(dTdt);
+
+   // 10. Initialize the hyperbolic heat operator
+   Array<int> ess_bdr;
+   if (pmesh.bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmesh.bdr_attributes.Max());
+      ess_bdr = dirichlet ? 1 : 0;
+   }
+
+   HyperbolicHeatOperator oper(fes, ir, ess_bdr, tau, rho_c, k0, beta, source,
+                               static_cast<PreconditionerType>(prec_type), assembly);
+
+   // The projected initial data does not necessarily satisfy T = 0 exactly on
+   // the boundary, and the operator keeps the acceleration at zero there, so the
+   // constraint is imposed on the initial state instead.
+   T.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+   dTdt.SetSubVector(oper.GetEssentialTrueDofs(), 0.0);
+   T_gf.SetFromTrueDofs(T);
+   dTdt_gf.SetFromTrueDofs(dTdt);
+
+   // 11. Set up visualization
+   socketstream sout;
+   if (visualization)
+   {
+      sout.open("localhost", visport);
+      if (!sout)
+      {
+         if (Mpi::Root())
+         {
+            mfem::out << "Unable to connect to GLVis server at localhost:"
+                      << visport << "\nGLVis visualization disabled.\n";
+         }
+         visualization = false;
+      }
+      else
+      {
+         sout.precision(8);
+         sout << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
+              << "\n";
+         sout << "solution\n" << pmesh << T_gf;
+         sout << "window_title 'Temperature'\n" << std::flush;
+      }
+   }
+
+   // Create a ParaView data collection
+   ParaViewDataCollection pd("dfem-hyperbolic-heat-output", &pmesh);
+   if (paraview)
+   {
+      pd.RegisterField("temperature", &T_gf);
+      pd.RegisterField("rate", &dTdt_gf);
+      pd.SetDataFormat(VTKFormat::BINARY);
+      if (order > 1)
+      {
+         pd.SetHighOrderOutput(true);
+         pd.SetLevelsOfDetail(order);
+      }
+      pd.SetCycle(0);
+      pd.SetTime(0.0);
+      pd.Save();
+   }
+
+   {
+      const real_t T_norm = oper.ComputeL2Norm(T);
+      const real_t dTdt_norm = oper.ComputeL2Norm(dTdt);
+      if (Mpi::Root())
+      {
+         mfem::out << "step        t      ||T||_L2   ||dT/dt||_L2   newton\n";
+         mfem::out << "   0  " << std::setw(8) << 0.0
+                   << "  " << std::setw(11) << T_norm
+                   << "  " << std::setw(13) << dTdt_norm
+                   << "  " << std::setw(7) << "-" << std::endl;
+      }
+   }
+
+   // 12. Perform time integration
+   ode_solver->Init(oper);
+   real_t t = 0.0;
+
+   bool last_step = false;
+   for (int ti = 1; !last_step; ti++)
+   {
+      if (t + dt >= t_final - dt / 2) { last_step = true; }
+
+      ode_solver->Step(T, dTdt, t, dt);
+
+      const real_t T_norm = oper.ComputeL2Norm(T);
+      const real_t dTdt_norm = oper.ComputeL2Norm(dTdt);
+      if (Mpi::Root())
+      {
+         mfem::out << std::setw(4) << ti << "  " << std::setw(8) << t
+                   << "  " << std::setw(11) << T_norm
+                   << "  " << std::setw(13) << dTdt_norm
+                   << "  " << std::setw(7) << oper.GetNewtonIterations()
+                   << std::endl;
+      }
+
+      if ((last_step || (ti % vis_steps) == 0) && (visualization || paraview))
+      {
+         T_gf.SetFromTrueDofs(T);
+         dTdt_gf.SetFromTrueDofs(dTdt);
+
+         if (visualization)
+         {
+            sout << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
+                 << "\n";
+            sout << "solution\n" << pmesh << T_gf << std::flush;
+         }
+
+         if (paraview)
+         {
+            pd.SetCycle(ti);
+            pd.SetTime(t);
+            pd.Save();
+         }
+      }
+   }
+
+   return 0;
+}
+

miniapps/dfem/dfem-hyperelasticity.cpp

@@ -18,7 +18,7 @@
 // Sample runs:  mpirun -np 4 dfem-hyperelasticity -o 1 -rs 0 -no-vis
 //               mpirun -np 4 dfem-hyperelasticity -mat linear-elastic -no-vis
 //               mpirun -np 4 dfem-hyperelasticity -mat mooney-rivlin -no-vis
-//               mpirun -np 4 dfem-hyperelasticity -mat holzapfel -no-vis
+//               mpirun -np 4 dfem-hyperelasticity -o 2 -rs 1 -pc 2 -no-vis
 //
 // Description:  This miniapp solves a quasistatic solid mechanics problem on
 //               the 3D beam used by the Hooke miniapp. The material response is
@@ -26,6 +26,9 @@
 //               residual is obtained with DifferentiableOperator::GetDerivative
 //               and Newton's Hessian-vector products are obtained with the new
 //               DifferentiableOperator::GetSecondDerivative functionality.
+//               The same second derivative can also be assembled into a
+//               HypreParMatrix, which -pc 2 uses to build a BoomerAMG
+//               preconditioner for the matrix-free Newton tangent.
 
 #include "mfem.hpp"
 #include "../../fem/dfem/doperator.hpp"
@@ -42,9 +45,6 @@ using mfem::future::dual;
 using dscalar_t = dual<real_t, real_t>;
 #endif
 
-#ifdef MFEM_USE_ENZYME
-
-
 constexpr int dim = 3;
 constexpr int Displacement = 0;
 constexpr int Coords = 1;
@@ -55,13 +55,13 @@ enum class MaterialType
    NeoHookean,
    LinearElastic,
    MooneyRivlin,
-   Holzapfel
 };
 
 enum class PreconditionerType
 {
    None,
-   Diagonal
+   Diagonal,
+   AMG
 };
 
 MaterialType ParseMaterial(const char *material)
@@ -79,13 +79,8 @@ MaterialType ParseMaterial(const char *material)
    {
       return MaterialType::MooneyRivlin;
    }
-   if (name == "holzapfel" || name == "fiber" || name == "fiber-reinforced")
-   {
-      return MaterialType::Holzapfel;
-   }
    MFEM_ABORT("Unknown material '" << name
-              << "'. Available materials: neo-hookean, linear-elastic, "
-              << "mooney-rivlin, holzapfel.");
+              << "'. Available materials: neo-hookean, linear-elastic, mooney-rivlin.");
    return MaterialType::NeoHookean;
 }
 
@@ -187,61 +182,11 @@ struct MooneyRivlinEnergy :
    }
 };
 
-template <typename dscalar_t>
-struct HolzapfelEnergy :
-   HyperelasticEnergyQFunction<HolzapfelEnergy<dscalar_t>, dscalar_t>
-{
-   real_t c = 50.0;
-   real_t kappa = 100.0;
-   real_t k1 = 10.0;
-   real_t k2 = 20.0;
-   //tensor<real_t, dim> a0 = {1.0, 0.0, 0.0};
-   tensor<real_t, dim> a0 = {0.7071067811865475, 0.7071067811865475, 0.0}; // Normalized fiber direction in the x-y plane at 45 degrees a = (1, 1, 0) / sqrt(2)
-   tensor<real_t, dim, dim> A = make_tensor<dim, dim>(
-   [&](int i, int j) { return a0(i) * a0(j); });
-
-   // Holzapfel-type transversely reinforced model with one fiber family a0:
-   // Ψ(F) = c/2 (Ī₁ - dim) + κ/2 log(J)^2
-   //      + k1/(2 k2) (exp(k2 (Ī₄ - 1)^2) - 1),
-   // where Ī₁ = J^(-2/3) tr(C), I₄ = A : C = a0.C.a0
-   // A = a0 ⊗ a0 is the fiber direction tensor.
-   // We use the full I4 invariant, to avoid auxetic behavior
-   //
-   // In this case we assume that the fiber direction is oriented 45 degrees in the x-y plane, i.e. a0 = (1, 1, 0)/sqrt(2).
-   // For a more general case one could start from an external fiber "field", and provide it as an input
-   // to the q-function, and then compute A = a0 ⊗ a0 at each quadrature point.
-   // This would require a different q-function signature including the fiber direction as well.
-   MFEM_HOST_DEVICE inline
-   dscalar_t psi(const tensor<dscalar_t, dim, dim> &F,
-                 const tensor<dscalar_t, dim, dim> & /* dudx */) const
-   {
-      // Kinematic quantities
-      const auto C = transpose(F) * F;
-      const auto J = det(F);
-      const auto Jm23 = pow(J, -2.0_r / 3.0_r);
-
-      // Strain invariants
-      const auto I1_bar = Jm23 * tr(C);
-      const auto I4 = ddot(C, A);
-      const auto fiber_strain = I4 - 1.0_r;
-      const auto log_J = log(J);
-
-      // Strain energy density components
-      const auto psi_vol = 0.5_r * kappa * log_J * log_J;
-      const auto psi_iso = 0.5_r * c * (I1_bar - real_t(dim));
-      const auto psi_aniso = (k1 / (2.0_r * k2)) * (exp(k2 * fiber_strain *
-                                                        fiber_strain) - 1.0_r);
-      // NOTE: in practice the anisotropic term should contribute only in tension, i.e. when I4 > 1, or fiber_strain > 0.
-      // mathematically this would introduce a non-smoothness in the energy functional that needs to taken care of.
-
-      return psi_vol + psi_iso + psi_aniso;
-   }
-};
-
 
 
 class HyperelasticOperator : public Operator
 {
+public:
    // Matrix-free Hessian-vector product used by Newton's method. This wraps the
    // functional second-derivative interface and applies the same essential-dof
    // treatment as the Hooke elasticity Jacobian operator.
@@ -293,6 +238,17 @@ class HyperelasticOperator : public Operator
          });
       }
 
+      // Matrix counterpart of AssembleDiagonal, for preconditioners that need
+      // a real matrix. @a A can be uninitialized; it is allocated by dFEM and
+      // ownership is passed to the caller. Eliminating the essential rows and
+      // columns puts 1.0 on their diagonal, matching Mult above.
+      void AssembleHessian(HypreParMatrix *&A) const
+      {
+         hessian->Assemble(A);
+         auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
+         delete Ae;
+      }
+
    private:
       const HyperelasticOperator &oper;
       Vector state;
@@ -310,9 +266,11 @@ public:
                         MaterialType material) :
       Operator(fes.GetTrueVSize()),
       fes(fes),
+      ir(ir),
       qspace(*fes.GetParMesh(), ir),
       qspace_vec(qspace, 1),
-      q(qspace_vec)
+      q(qspace_vec),
+      material(material)
    {
       auto &mesh_nodes =
          *static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
@@ -380,19 +338,6 @@ public:
                ir, all_domain_attr, derivatives, second_derivatives);
             break;
          }
-         case MaterialType::Holzapfel:
-         {
-            // Fiber-reinforced Holzapfel-type material with fibers aligned to
-            // the beam axis. The functional registration lets dFEM derive both
-            // the residual and the Hessian-vector product from the energy.
-            HolzapfelEnergy<dscalar_t> energy;
-            internal_energy_dop->AddDomainIntegrator<LocalQFBackend>(
-               energy,
-               Inputs<Gradient<Displacement>, Gradient<Coords>, Weight> {},
-               Outputs<FunctionalValue<Energy>> {},
-               ir, all_domain_attr, derivatives, second_derivatives);
-            break;
-         }
       }
 
       // The first variation of a functional is exposed as a stateless
@@ -439,9 +384,11 @@ public:
 private:
    ParFiniteElementSpace &fes;
    ParFiniteElementSpace *mesh_nodes_fes = nullptr;
+   const IntegrationRule &ir;
    QuadratureSpace qspace;
    VectorQuadratureSpace qspace_vec;
    QuadratureFunction q;
+   MaterialType material;
    Vector mesh_nodes_tdofs;
    Array<int> ess_tdofs;
    Array<int> prescribed_tdofs;
@@ -455,21 +402,78 @@ private:
    mutable std::shared_ptr<HessianOperator> hessian;
 };
 
-#endif // MFEM_USE_ENZYME
+// BoomerAMG preconditioner for the matrix-free Hessian. HypreBoomerAMG needs a
+// real HypreParMatrix, so on every SetOperator we let dFEM assemble the second
+// derivative and setup AMG.
+class HessianAMG : public Solver
+{
+public:
+   HessianAMG(ParFiniteElementSpace *fes_) : fes(fes_)
+   {
+      amg.SetPrintLevel(0);
+   }
+
+   void SetOperator(const Operator &op) override
+   {
+      const auto *H =
+         dynamic_cast<const HyperelasticOperator::HessianOperator *>(&op);
+      MFEM_VERIFY(H, "HessianAMG requires a HessianOperator");
+      height = width = op.Height();
+
+      delete A;
+      A = nullptr;
+      H->AssembleHessian(A);
+      amg.SetOperator(*A);
+      // Tell BoomerAMG this is a dim-component displacement system rather than
+      // a scalar one. order_bynodes MUST be true: the state space is built as
+      // ParFiniteElementSpace(&pmesh, &fec, dim, Ordering::byNODES), while
+      // hypre's default assumes byVDIM. Getting this flag wrong hands hypre a
+      // bogus dof -> function map and converges *worse* than passing no systems
+      // options at all. Must follow SetOperator: the dof map is sized from
+      // height, which SetOperator establishes.
+      amg.SetSystemsOptions(fes->GetVDim(), /*order_bynodes=*/true);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
+
+   ~HessianAMG() { delete A; }
+
+private:
+   ParFiniteElementSpace *fes;
+   HypreParMatrix *A = nullptr;
+   HypreBoomerAMG amg;
+};
+
+// Build the preconditioner selected by -pc. Newton hands the current tangent to
+// the Krylov solver on every iteration, and IterativeSolver::SetOperator
+// forwards it to the preconditioner, so the AMG hierarchy is rebuilt from the
+// freshly assembled Hessian at each Newton step.
+std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type,
+                                           ParFiniteElementSpace *fes)
+{
+   switch (type)
+   {
+      case PreconditionerType::None:
+         return nullptr;
+      case PreconditionerType::Diagonal:
+         return std::make_unique<OperatorJacobiSmoother>();
+      case PreconditionerType::AMG:
+         return std::make_unique<HessianAMG>(fes);
+      default:
+         MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
+   }
+   return nullptr;
+}
+
 
 int main(int argc, char *argv[])
 {
    Mpi::Init(argc, argv);
+   const int num_procs = Mpi::WorldSize();
+   const int myid = Mpi::WorldRank();
    Hypre::Init();
 
-#ifndef MFEM_USE_ENZYME
-   if (Mpi::Root())
-   {
-      mfem::out << "This miniapp requires MFEM_USE_ENZYME=YES because it uses "
-                << "dFEM functional second derivatives.\n";
-   }
-   return 0;
-#else
+
    int order = 1;
    const char *device_config = "cpu";
    int serial_refinement_levels = 0;
@@ -479,7 +483,6 @@ int main(int argc, char *argv[])
    bool visualization = true;
    bool paraview = false;
    int visport = 19916;
-   const char *outfolder = "./Output";
 
    OptionsParser args(argc, argv);
    args.AddOption(&order, "-o", "--order",
@@ -489,9 +492,9 @@ int main(int argc, char *argv[])
    args.AddOption(&serial_refinement_levels, "-rs", "--ref-serial",
                   "Number of uniform refinements on the serial mesh.");
    args.AddOption(&material_name, "-mat", "--material",
-                  "Material: neo-hookean, linear-elastic, mooney-rivlin, or holzapfel.");
+                  "Material: neo-hookean, linear-elastic, mooney-rivlin.");
    args.AddOption(&prec_type, "-pc", "--preconditioner",
-                  "Preconditioner: 0=none, 1=diagonal.");
+                  "Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
    args.AddOption(&cg_tol, "-tol", "--cg-tol",
                   "Relative tolerance for the CG solver.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -501,8 +504,6 @@ int main(int argc, char *argv[])
                   "--no-paraview",
                   "Enable or disable ParaView DataCollection output.");
    args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
-   args.AddOption(&outfolder, "-of", "--output-folder",
-                  "Output folder for ParaView DataCollection files.");
    args.ParseCheck();
 
    const MaterialType material = ParseMaterial(material_name);
@@ -569,18 +570,9 @@ int main(int argc, char *argv[])
    cg.SetMaxIter(10000);
    cg.SetPrintLevel(2);
 
-   std::unique_ptr<Solver> pc;
-   switch (static_cast<PreconditionerType>(prec_type))
-   {
-      case PreconditionerType::None:
-         break;
-      case PreconditionerType::Diagonal:
-         pc = std::make_unique<OperatorJacobiSmoother>();
-         cg.SetPreconditioner(*pc);
-         break;
-      default:
-         MFEM_ABORT("Unknown preconditioner type: " << prec_type);
-   }
+   std::unique_ptr<Solver> pc =
+      MakePreconditioner(static_cast<PreconditionerType>(prec_type), &fes);
+   if (pc) { cg.SetPreconditioner(*pc); }
 
    NewtonSolver newton(MPI_COMM_WORLD);
    newton.SetSolver(cg);
@@ -604,18 +596,15 @@ int main(int argc, char *argv[])
    {
       char vishost[] = "localhost";
       socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
-               << "\n";
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
       sol_sock.precision(8);
       sol_sock << "solution\n" << pmesh << U_gf << std::flush;
    }
 
    if (paraview)
    {
-      // Create a ParaView data collection. Save() creates the output directory
-      // tree under the prefix path itself, with the ranks synchronized.
-      ParaViewDataCollection pd("dfem-hyperelasticity", &pmesh);
-      pd.SetPrefixPath(outfolder);
+      // Create a ParaView data collection
+      ParaViewDataCollection pd("dfem-hyperelasticity-output", &pmesh);
       pd.RegisterField("displacement", &U_gf);
       pd.SetDataFormat(VTKFormat::BINARY);
       if (order > 1)
@@ -629,5 +618,4 @@ int main(int argc, char *argv[])
    }
 
    return 0;
-#endif // MFEM_USE_ENZYME
 }

miniapps/dfem/dfem-nonlinear-poisson.cpp

@@ -0,0 +1,622 @@
+// 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.
+//
+//               ----------------------------------------------
+//               Nonlinear Poisson (Diffusion) Problem with dFEM
+//               ----------------------------------------------
+//
+// Compile with: make dfem-nonlinear-poisson
+//
+// Sample runs:  mpirun -np 4 dfem-nonlinear-poisson -k grad -f energy -no-vis
+//               mpirun -np 4 dfem-nonlinear-poisson -k grad -f residual -no-vis
+//               mpirun -np 4 dfem-nonlinear-poisson -k sol -f residual -no-vis
+//               mpirun -np 4 dfem-nonlinear-poisson -o 2 -r 4 -pc 2 -no-vis
+//
+// Description:  This miniapp solves the nonlinear diffusion problem
+//
+//                   -div( kappa(u, grad u) grad u ) = f   in Omega
+//                                                 u = 0   on dOmega
+//
+//               on the unit square (set the compile-time constant dim to 3 for
+//               the unit cube). It shows the two ways a nonlinear
+//               problem can be given to dFEM, and the rule that decides which
+//               one is available: the coefficient (-k) determines whether an
+//               energy exists, and the formulation (-f) picks how it is
+//               written.
+//
+//               -k grad : kappa = 1 + |grad u|^2 depends only on the gradient,
+//                         so the operator is the first variation of
+//                         Pi(u) = int psi(|grad u|^2) dx with 2 psi' = kappa.
+//                         The tangent is symmetric positive definite, and both
+//                         formulations apply.
+//
+//               -k sol  : kappa = 1 + u^2 depends on the solution, so the
+//                         operator is NOT the gradient of any functional: the
+//                         tangent gains int kappa'(u) du (grad u . grad v) dx,
+//                         which is not symmetric. Only -f residual applies, and
+//                         the linear solves need GMRES rather than CG.
+//
+//               -f energy   : the q-function returns the energy density psi as
+//                             a FunctionalValue. The residual is GetDerivative,
+//                             the Newton tangent is GetSecondDerivative.
+//
+//               -f residual : the q-function returns the pointwise residual, as
+//                             in the minimal surface miniapp, and the tangent
+//                             is GetDerivative.
+//
+//               With -k grad both formulations describe the same problem and
+//               produce the same answer. Convergence of both against a
+//               manufactured solution is verified in
+//               tests/unit/dfem/test_nonlinear_poisson.cpp.
+
+#include "mfem.hpp"
+#include "../../fem/dfem/doperator.hpp"
+#include "../../fem/dfem/backends/local_qf/prelude.hpp"
+
+using namespace mfem;
+
+// This example code demonstrates the use of new features in MFEM that are in
+// development but exposed through the mfem::future namespace. All features
+// under this namespace might change their interface or behavior in upcoming
+// releases until they have stabilized.
+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
+
+// Space dimension. The q-functions and the operator are templated on it, so it
+// is fixed at compile time; set it to 3 to run the same problem on the unit
+// cube.
+constexpr int dim = 2;
+
+// Field IDs shared by every integrator in this miniapp.
+static constexpr int Solution = 0;
+static constexpr int Coords = 1;
+static constexpr int Energy = 2;
+
+// The coefficient determines whether an energy formulation exists at all.
+enum class KappaType
+{
+   GradientDependent, // kappa = 1 + |grad u|^2, symmetric, energy available
+   SolutionDependent  // kappa = 1 + u^2, non-symmetric, residual only
+};
+
+enum class FormType
+{
+   Energy,  // q-function returns psi, tangent = GetSecondDerivative
+   Residual // q-function returns the residual, tangent = GetDerivative
+};
+
+enum class PreconditionerType { None, Diagonal, AMG };
+
+// ----------------------------------------------------------------------------
+// Pointwise description of the problem: the coefficients and one q-function per
+// formulation.
+// ----------------------------------------------------------------------------
+template <int DIM>
+struct NonlinearPoisson
+{
+   // kappa = 1 + |grad u|^2, as a function of s = |grad u|^2.
+   template <typename T>
+   MFEM_HOST_DEVICE static inline T Kappa_grad(const T &s)
+   {
+      return 1.0_r + s;
+   }
+
+   // kappa = 1 + u^2.
+   template <typename T>
+   MFEM_HOST_DEVICE static inline T Kappa_sol(const T &u)
+   {
+      return 1.0_r + u * u;
+   }
+
+   // Energy density with 2 Psi_grad'(s) = Kappa_grad(s), so that the first
+   // variation of int psi dx is int kappa grad u . grad v dx.
+   template <typename T>
+   MFEM_HOST_DEVICE static inline T Psi_grad(const T &s)
+   {
+      return 0.5_r * s + 0.25_r * s * s;
+   }
+
+   // The q-function supplies only the energy density; dFEM forms both the first
+   // variation (the residual) and the second variation (the Newton tangent)
+   // from it. Only Kappa_grad has a potential, so only it can be written here.
+   struct EnergyBased
+   {
+      MFEM_HOST_DEVICE inline
+      auto operator()(const tensor<dscalar_t, DIM> &dudxi,
+                      const tensor<real_t, DIM, DIM> &J,
+                      const real_t &w,
+                      dscalar_t &energy) const
+      {
+         const auto dudx = dudxi * inv(J);
+         // det(J) * w is the quadrature form of dx in Pi(u) = int psi dx.
+         energy = Psi_grad(sqnorm(dudx)) * det(J) * w;
+      }
+   };
+
+   // KAPPA selects the coefficient family. Both variants take the same inputs,
+   // so they register identically; the solution value is interpolated either
+   // way even though only Kappa_sol consumes it.
+   template <KappaType KAPPA>
+   struct ResidualBased
+   {
+      MFEM_HOST_DEVICE inline
+      auto operator()(const dscalar_t &u,
+                      const tensor<dscalar_t, DIM> &dudxi,
+                      const tensor<real_t, DIM, DIM> &J,
+                      const real_t &w,
+                      tensor<dscalar_t, DIM> &dvdx) const
+      {
+         const auto invJ = inv(J);
+         const auto dudx = dudxi * invJ;
+         dscalar_t kappa;
+         if constexpr (KAPPA == KappaType::GradientDependent)
+         {
+            kappa = Kappa_grad(sqnorm(dudx));
+         }
+         else
+         {
+            kappa = Kappa_sol(u);
+         }
+         dvdx = kappa * dudx * transpose(invJ) * det(J) * w;
+      }
+   };
+};
+
+// ----------------------------------------------------------------------------
+// Depending on @a form, the DifferentiableOperator below holds either the
+// energy Pi(u), whose first derivative is the residual, or the residual R(u)
+// directly. Both expose a tangent through the same DerivativeOperator
+// interface, so the Newton plumbing is shared.
+// ----------------------------------------------------------------------------
+template <int DIM>
+class NonlinearPoissonOperator : public Operator
+{
+public:
+   // Matrix-free Newton tangent. Essential directions are removed before the
+   // apply and restored as identity rows afterwards.
+   class JacobianOperator : public Operator
+   {
+   public:
+      JacobianOperator(const NonlinearPoissonOperator &oper,
+                       const Vector &state) :
+         Operator(oper.Height()), oper(oper), z(oper.Height())
+      {
+         MultiVector X{state, oper.mesh_nodes_tdofs};
+         // In the energy form the residual is already the first derivative of
+         // the energy, so its tangent is the second derivative. In the residual
+         // form the residual is differentiated directly.
+         tangent = (oper.form == FormType::Energy)
+                   ? oper.dop->GetSecondDerivative(Solution, X)
+                   : oper.dop->GetDerivative(Solution, X);
+      }
+
+      void Mult(const Vector &x, Vector &y) const override
+      {
+         z = x;
+         z.SetSubVector(oper.ess_tdofs, 0.0);
+
+         MultiVector Y{y};
+         tangent->Mult(z, Y);
+
+         auto d_y = y.ReadWrite();
+         const auto d_x = x.Read();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_y[d_dofs[i]] = d_x[d_dofs[i]];
+         });
+      }
+
+      void AssembleDiagonal(Vector &diag) const override
+      {
+         tangent->AssembleDiagonal(diag);
+         auto d_diag = diag.ReadWrite();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_diag[d_dofs[i]] = 1.0;
+         });
+      }
+
+      // Matrix counterpart of AssembleDiagonal, for preconditioners that need a
+      // real matrix. @a A can be uninitialized; it is allocated by dFEM and
+      // ownership is passed to the caller. Eliminating the essential rows and
+      // columns puts 1.0 on their diagonal, matching Mult above.
+      void AssembleJacobian(HypreParMatrix *&A) const
+      {
+         tangent->Assemble(A);
+         auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
+         delete Ae;
+      }
+
+   private:
+      const NonlinearPoissonOperator &oper;
+      mutable Vector z;
+      std::shared_ptr<DerivativeOperator> tangent;
+   };
+
+   NonlinearPoissonOperator(ParFiniteElementSpace &fes,
+                            const IntegrationRule &ir,
+                            FormType form, KappaType kappa) :
+      Operator(fes.GetTrueVSize()),
+      fes(fes),
+      qspace(*fes.GetParMesh(), ir),
+      qspace_vec(qspace, 1),
+      form(form)
+   {
+      MFEM_VERIFY(form == FormType::Residual ||
+                  kappa == KappaType::GradientDependent,
+                  "the energy formulation requires -k grad: a solution "
+                  "dependent coefficient is not the gradient of a functional");
+
+      auto &mesh_nodes =
+         *static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
+      mesh_nodes_fes = mesh_nodes.ParFESpace();
+      mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
+
+      const std::vector<FieldDescriptor> inputs =
+      {
+         {Solution, &fes}, {Coords, mesh_nodes_fes}
+      };
+
+      Array<int> all_domain_attr;
+      if (fes.GetMesh()->attributes.Size() > 0)
+      {
+         all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
+         all_domain_attr = 1;
+      }
+
+      auto derivatives = std::integer_sequence<size_t, Solution> {};
+
+      if (form == FormType::Energy)
+      {
+         // A single scalar per quadrature point, summed by dFEM into the value
+         // of the functional.
+         const std::vector<FieldDescriptor> outputs = {{Energy, &qspace_vec}};
+         dop = std::make_shared<DifferentiableOperator>(
+                  inputs, outputs, *fes.GetParMesh());
+
+         // Requesting the second derivative here is what makes the Newton
+         // tangent available later through GetSecondDerivative.
+         auto second_derivatives = SecondDerivatives<Pairs::All> {};
+         typename NonlinearPoisson<DIM>::EnergyBased qf;
+         dop->AddDomainIntegrator<LocalQFBackend>(
+            qf,
+            Inputs<Gradient<Solution>, Gradient<Coords>, Weight> {},
+            Outputs<FunctionalValue<Energy>> {},
+            ir, all_domain_attr, derivatives, second_derivatives);
+
+         // The first variation of a functional is stateless: the current
+         // solution is passed to gradient->Mult(X, Y) on every residual
+         // evaluation, so this wrapper is built once and reused.
+         gradient = dop->GetDerivative(Solution);
+      }
+      else
+      {
+         // The output operator is the gradient of the test function basis,
+         // completing the diffusion-like weak form B^T D(B u, B x, w).
+         const std::vector<FieldDescriptor> outputs = {{Solution, &fes}};
+         dop = std::make_shared<DifferentiableOperator>(
+                  inputs, outputs, *fes.GetParMesh());
+
+         const auto in = Inputs < Value<Solution>, Gradient<Solution>,
+                    Gradient<Coords>, Weight > {};
+         const auto out = Outputs<Gradient<Solution>> {};
+
+         if (kappa == KappaType::GradientDependent)
+         {
+            typename NonlinearPoisson<DIM>::template
+            ResidualBased<KappaType::GradientDependent> qf;
+            dop->AddDomainIntegrator<LocalQFBackend>(
+               qf, in, out, ir, all_domain_attr, derivatives);
+         }
+         else
+         {
+            typename NonlinearPoisson<DIM>::template
+            ResidualBased<KappaType::SolutionDependent> qf;
+            dop->AddDomainIntegrator<LocalQFBackend>(
+               qf, in, out, ir, all_domain_attr, derivatives);
+         }
+      }
+   }
+
+   void SetEssentialAttributes(const Array<int> &ess_bdr)
+   {
+      fes.GetEssentialTrueDofs(ess_bdr, ess_tdofs);
+   }
+
+   // The load stays outside the dFEM functional: the full potential is
+   // Pi(u) - (b, u), whose first variation is the residual below.
+   void SetLoad(const Vector &b) { load = b; }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      MultiVector X{x, mesh_nodes_tdofs};
+      MultiVector Y{y};
+      // For a functional, GetDerivative returns the gradient action directly,
+      // including the pointwise reverse seed for the summed energy.
+      if (form == FormType::Energy) { gradient->Mult(X, Y); }
+      else { dop->Mult(X, Y); }
+
+      y -= load;
+      y.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   Operator& GetGradient(const Vector &x) const override
+   {
+      jacobian = std::make_shared<JacobianOperator>(*this, x);
+      return *jacobian;
+   }
+
+private:
+   ParFiniteElementSpace &fes;
+   ParFiniteElementSpace *mesh_nodes_fes = nullptr;
+   QuadratureSpace qspace;
+   VectorQuadratureSpace qspace_vec;
+   FormType form;
+   Vector mesh_nodes_tdofs, load;
+   Array<int> ess_tdofs;
+
+   std::shared_ptr<DifferentiableOperator> dop;
+   std::shared_ptr<DerivativeOperator> gradient;
+   mutable std::shared_ptr<JacobianOperator> jacobian;
+};
+
+// BoomerAMG for the matrix-free tangent. HypreBoomerAMG needs a real
+// HypreParMatrix, so on every SetOperator we let dFEM assemble the tangent and
+// set up AMG. This works for both formulations, since GetDerivative and
+// GetSecondDerivative both return a DerivativeOperator that can assemble itself.
+template <int DIM>
+class JacobianAMG : public Solver
+{
+   using JacobianOperator =
+      typename NonlinearPoissonOperator<DIM>::JacobianOperator;
+
+public:
+   JacobianAMG() { amg.SetPrintLevel(0); }
+
+   void SetOperator(const Operator &op) override
+   {
+      const auto *Jop = dynamic_cast<const JacobianOperator *>(&op);
+      MFEM_VERIFY(Jop, "JacobianAMG requires a JacobianOperator");
+      height = width = op.Height();
+
+      delete A;
+      A = nullptr;
+      Jop->AssembleJacobian(A);
+      amg.SetOperator(*A);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
+
+   ~JacobianAMG() { delete A; }
+
+private:
+   HypreParMatrix *A = nullptr;
+   HypreBoomerAMG amg;
+};
+
+template <int DIM>
+std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
+{
+   switch (type)
+   {
+      case PreconditionerType::None:     return nullptr;
+      case PreconditionerType::Diagonal: return
+            std::make_unique<OperatorJacobiSmoother>();
+      case PreconditionerType::AMG:      return std::make_unique<JacobianAMG<DIM>>();
+      default:
+         MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
+   }
+   return nullptr;
+}
+
+// ----------------------------------------------------------------------------
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // 2. Parse command-line options
+   int order = 1;
+   int refinements = 3;
+   const char *device_config = "cpu";
+   const char *kappa_name = "grad";
+   const char *form_name = "energy";
+   int prec_type = static_cast<int>(PreconditionerType::AMG);
+   real_t krylov_tol = 1e-8;
+   real_t source = 1.0;
+   bool visualization = true;
+   bool paraview = false;
+   int visport = 19916;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&refinements, "-r", "--refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&kappa_name, "-k", "--kappa",
+                  "Diffusion coefficient: grad = 1 + |grad u|^2, "
+                  "sol = 1 + u^2.");
+   args.AddOption(&form_name, "-f", "--form",
+                  "dFEM formulation: energy (requires -k grad), residual.");
+   args.AddOption(&prec_type, "-pc", "--preconditioner",
+                  "Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
+   args.AddOption(&krylov_tol, "-tol", "--krylov-tol",
+                  "Relative tolerance for the linear solver.");
+   args.AddOption(&source, "-s", "--source", "Constant source term.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or disable ParaView DataCollection output.");
+   args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
+   args.ParseCheck();
+
+   const std::string kappa_str(kappa_name), form_str(form_name);
+   MFEM_VERIFY(kappa_str == "grad" || kappa_str == "sol",
+               "-k must be grad or sol");
+   MFEM_VERIFY(form_str == "energy" || form_str == "residual",
+               "-f must be energy or residual");
+
+   const KappaType kappa = (kappa_str == "grad")
+                           ? KappaType::GradientDependent
+                           : KappaType::SolutionDependent;
+   const FormType form = (form_str == "energy") ? FormType::Energy
+                         : FormType::Residual;
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   if (Mpi::Root())
+   {
+      mfem::out << "coefficient: "
+                << (kappa == KappaType::GradientDependent ? "1 + |grad u|^2"
+                    : "1 + u^2")
+                << ", formulation: " << form_str
+                << ", tangent: "
+                << (kappa == KappaType::GradientDependent
+                    ? "symmetric (CG)" : "non-symmetric (GMRES)")
+                << std::endl;
+   }
+
+   // 4. Create the unit square/cube mesh and refine it
+   Mesh mesh = (dim == 2)
+               ? Mesh::MakeCartesian2D(4, 4, Element::QUADRILATERAL)
+               : Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
+   for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
+   mesh.SetCurvature(order);
+
+   // 5. Define a parallel mesh
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   pmesh.EnsureNodes();
+
+   // 6. Define a scalar H1 space on the parallel mesh. GlobalTrueVSize is
+   //    collective, so every rank has to reach it.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fes(&pmesh, &fec, 1);
+   const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
+   if (Mpi::Root()) { mfem::out << "#dofs: " << global_dofs << std::endl; }
+
+   // 7. Set up the integration rule. It has to resolve the nonlinear
+   //    coefficient, not just the bilinear part, hence the extra order
+   //    compared to a linear problem.
+   const IntegrationRule &ir =
+      IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * order + 2);
+
+   // 8. Create the nonlinear operator for the chosen coefficient and
+   //    formulation
+   NonlinearPoissonOperator<dim> poisson_op(fes, ir, form, kappa);
+
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 1;
+   poisson_op.SetEssentialAttributes(ess_bdr);
+
+   // 9. Assemble the constant source term. The boundary data is homogeneous.
+   ConstantCoefficient source_coeff(source);
+   ParLinearForm b(&fes);
+   b.AddDomainIntegrator(new DomainLFIntegrator(source_coeff, &ir));
+   b.Assemble();
+
+   Vector B(fes.GetTrueVSize());
+   b.ParallelAssemble(B);
+   poisson_op.SetLoad(B);
+
+   ParGridFunction u_gf(&fes);
+   u_gf = 0.0;
+
+   Vector U(fes.GetTrueVSize());
+   u_gf.GetTrueDofs(U);
+
+   // 10. Set up the linear solver used within Newton's method. The
+   //     solution-dependent coefficient gives a non-symmetric tangent, so CG
+   //     does not apply there.
+   std::unique_ptr<IterativeSolver> krylov;
+   if (kappa == KappaType::GradientDependent)
+   {
+      krylov = std::make_unique<CGSolver>(MPI_COMM_WORLD);
+   }
+   else
+   {
+      krylov = std::make_unique<GMRESSolver>(MPI_COMM_WORLD);
+   }
+   krylov->SetAbsTol(0.0);
+   krylov->SetRelTol(krylov_tol);
+   krylov->SetMaxIter(2000);
+   krylov->SetPrintLevel(0);
+
+   std::unique_ptr<Solver> pc =
+      MakePreconditioner<dim>(static_cast<PreconditionerType>(prec_type));
+   if (pc) { krylov->SetPreconditioner(*pc); }
+
+   // 11. Set up the nonlinear solver (Newton)
+   NewtonSolver newton(MPI_COMM_WORLD);
+   newton.SetOperator(poisson_op);
+   newton.SetSolver(*krylov);
+   newton.SetAbsTol(0.0);
+#ifdef MFEM_USE_SINGLE
+   newton.SetRelTol(1e-6);
+#else
+   newton.SetRelTol(1e-10);
+#endif
+   newton.SetMaxIter(30);
+   newton.SetPrintLevel(1);
+
+   // 12. Solve the nonlinear system using Newton's method
+   Vector zero;
+   newton.Mult(zero, U);
+   MFEM_VERIFY(newton.GetConverged(), "Newton did not converge");
+
+   u_gf.Distribute(U);
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
+               << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << u_gf << std::flush;
+   }
+
+   // 14. Save the solution in parallel using ParaView format
+   if (paraview)
+   {
+      ParaViewDataCollection pd("dfem-nonlinear-poisson-output", &pmesh);
+      pd.RegisterField("solution", &u_gf);
+      pd.SetDataFormat(VTKFormat::BINARY);
+      if (order > 1)
+      {
+         pd.SetHighOrderOutput(true);
+         pd.SetLevelsOfDetail(order);
+      }
+      pd.SetCycle(0);
+      pd.SetTime(0.0);
+      pd.Save();
+   }
+
+   return 0;
+}

miniapps/dfem/dfem-reaction-diffusion.cpp

@@ -0,0 +1,592 @@
+// 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.
+//
+//              ------------------------------------------------
+//              Nonlinear Reaction-Diffusion with dFEM (two
+//              outputs from a single quadrature point function)
+//              ------------------------------------------------
+//
+// Compile with: make dfem-reaction-diffusion
+//
+// Sample runs:  mpirun -np 4 dfem-reaction-diffusion -b 1 -no-vis
+//               mpirun -np 4 dfem-reaction-diffusion -b 1 -o 2 -r 4 -pc 2 -no-vis
+//               mpirun -np 4 dfem-reaction-diffusion -b 10 -a 0 -no-vis
+//               mpirun -np 4 dfem-reaction-diffusion -b 0 -check -no-vis
+//
+// Description:  This miniapp solves the nonlinear reaction-diffusion problem
+//
+//                   -div( kappa grad u ) + alpha u + beta u^3 = f   in Omega
+//                                                          u = 0   on dOmega
+//
+//               on the unit square (set the compile-time constant dim to 3 for
+//               the unit cube). 
+//
+//               This miniapp demonstrates the use of dFEM in multiple-output mode. 
+//               The q-function writes two outputs onto the same test field, e.g.:
+//
+//                   Outputs<Gradient<Solution>, Value<Solution>>
+//
+//               dFEM integrates each output against its own test basis
+//               operation and accumulates all of them into a single output
+//               vector,
+//
+//                   R_i = int (kappa grad u) . grad phi_i dx
+//                       + int (alpha u + beta u^3 - f) phi_i dx
+//
+//               so the diffusion and the reaction term share one kernel and one
+//               restriction of u. 
+//
+
+#include "mfem.hpp"
+#include "../../fem/dfem/doperator.hpp"
+#include "../../fem/dfem/backends/local_qf/prelude.hpp"
+
+#include <memory>
+
+using namespace mfem;
+
+// This example code demonstrates the use of new features in MFEM that are in
+// development but exposed through the mfem::future namespace. All features
+// under this namespace might change their interface or behavior in upcoming
+// releases until they have stabilized.
+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
+
+// Space dimension. The q-function and the operator are templated on it, so it
+// is fixed at compile time; set it to 3 to run the same problem on the unit
+// cube.
+constexpr int dim = 2;
+
+// Field IDs used by the dFEM integrator.
+static constexpr int Solution = 0;
+static constexpr int Coords = 1;
+
+enum class PreconditionerType { None, Diagonal, AMG };
+
+// ----------------------------------------------------------------------------
+// Pointwise description of the problem.
+// ----------------------------------------------------------------------------
+
+/// Quadrature point kernel for
+///
+///   R(u; v) = int kappa grad u . grad v dx
+///           + int (alpha u + beta u^3 - f) v dx
+///
+/// The two integrands leave through two different outputs: @a dvdx is paired
+/// with Gradient<Solution> and @a v with Value<Solution>.
+/// However, both are attached to the same fieldID, so dFEM sums their contributions
+/// into one residual vector.
+template <int DIM>
+struct ReactionDiffusion
+{
+   real_t kappa = 1.0;
+   real_t alpha = 1.0;
+   real_t beta = 1.0;
+   real_t source = 1.0;
+
+   MFEM_HOST_DEVICE inline
+   auto operator()(const dscalar_t &u,
+                   const tensor<dscalar_t, DIM> &dudxi,
+                   const tensor<real_t, DIM, DIM> &J,
+                   const real_t &w,
+                   tensor<dscalar_t, DIM> &dvdx,
+                   dscalar_t &v) const
+   {
+      const auto invJ = inv(J);
+      const auto dudx = dudxi * invJ;
+      const real_t dxw = det(J) * w;
+
+      dvdx = kappa * dudx * transpose(invJ) * dxw;
+      v = (alpha * u + beta * u * u * u - source) * dxw;
+   }
+};
+
+// ----------------------------------------------------------------------------
+// The residual operator. A single DifferentiableOperator holds both terms, so
+// the Newton tangent below is the derivative of the whole q-function.
+// ----------------------------------------------------------------------------
+template <int DIM>
+class ReactionDiffusionOperator : public Operator
+{
+public:
+   /// Matrix-free Newton tangent K + (alpha + 3 beta u^2) M. Essential
+   /// directions are removed before the apply and restored as identity rows
+   /// afterwards.
+   class JacobianOperator : public Operator
+   {
+   public:
+      JacobianOperator(const ReactionDiffusionOperator &oper,
+                       const Vector &state) :
+         Operator(oper.Height()), oper(oper), z(oper.Height())
+      {
+         MultiVector X{state, oper.mesh_nodes_tdofs};
+         // Differentiating the q-function linearizes both outputs at once,
+         // so no term has to be added by hand here.
+         tangent = oper.dop->GetDerivative(Solution, X);
+      }
+
+      void Mult(const Vector &x, Vector &y) const override
+      {
+         z = x;
+         z.SetSubVector(oper.ess_tdofs, 0.0);
+
+         MultiVector Y{y};
+         tangent->Mult(z, Y);
+
+         auto d_y = y.ReadWrite();
+         const auto d_x = x.Read();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_y[d_dofs[i]] = d_x[d_dofs[i]];
+         });
+      }
+
+      void AssembleDiagonal(Vector &diag) const override
+      {
+         tangent->AssembleDiagonal(diag);
+         auto d_diag = diag.ReadWrite();
+         const auto d_dofs = oper.ess_tdofs.Read();
+         mfem::forall(oper.ess_tdofs.Size(), [=] MFEM_HOST_DEVICE (int i)
+         {
+            d_diag[d_dofs[i]] = 1.0;
+         });
+      }
+
+      // Matrix counterpart of AssembleDiagonal, for preconditioners that need a
+      // real matrix. @a A can be uninitialized; it is allocated by dFEM and
+      // ownership is passed to the caller. Eliminating the essential rows and
+      // columns puts 1.0 on their diagonal, matching Mult above.
+      void AssembleJacobian(HypreParMatrix *&A) const
+      {
+         tangent->Assemble(A);
+         auto Ae = A->EliminateRowsCols(oper.ess_tdofs);
+         delete Ae;
+      }
+
+   private:
+      const ReactionDiffusionOperator &oper;
+      mutable Vector z;
+      std::shared_ptr<DerivativeOperator> tangent;
+   };
+
+   ReactionDiffusionOperator(ParFiniteElementSpace &fes,
+                             const IntegrationRule &ir,
+                             real_t kappa, real_t alpha, real_t beta,
+                             real_t source) :
+      Operator(fes.GetTrueVSize()),
+      fes(fes)
+   {
+      MFEM_VERIFY(kappa > 0.0, "the diffusion coefficient has to be positive");
+      MFEM_VERIFY(alpha >= 0.0 && beta >= 0.0,
+                  "negative reaction coefficients make the tangent indefinite; "
+                  "this miniapp assumes -a and -b are non-negative so that the "
+                  "problem stays convex and CG applies");
+
+      auto &mesh_nodes =
+         *static_cast<ParGridFunction *>(fes.GetParMesh()->GetNodes());
+      mesh_nodes_fes = mesh_nodes.ParFESpace();
+      mesh_nodes.GetTrueDofs(mesh_nodes_tdofs);
+
+      const std::vector<FieldDescriptor> inputs =
+      {
+         {Solution, &fes}, {Coords, mesh_nodes_fes}
+      };
+      // One output field, even though the q-function writes two outputs onto
+      // it: the number of FieldOperators does not set the number of blocks.
+      const std::vector<FieldDescriptor> outputs = {{Solution, &fes}};
+
+      Array<int> all_domain_attr;
+      if (fes.GetMesh()->attributes.Size() > 0)
+      {
+         all_domain_attr.SetSize(fes.GetMesh()->attributes.Max());
+         all_domain_attr = 1;
+      }
+
+      dop = std::make_shared<DifferentiableOperator>(
+               inputs, outputs, *fes.GetParMesh());
+
+      ReactionDiffusion<DIM> qf{kappa, alpha, beta, source};
+      dop->AddDomainIntegrator<LocalQFBackend>(
+         qf,
+         Inputs < Value<Solution>, Gradient<Solution>,
+         Gradient<Coords>, Weight > {},
+         Outputs<Gradient<Solution>, Value<Solution>> {},
+         ir, all_domain_attr, Derivatives<Solution> {});
+   }
+
+   void SetEssentialAttributes(const Array<int> &ess_bdr)
+   {
+      fes.GetEssentialTrueDofs(ess_bdr, ess_tdofs);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      MultiVector X{x, mesh_nodes_tdofs};
+      MultiVector Y{y};
+      // Both terms and the load are inside the q-function, so there is nothing
+      // to add to the residual here.
+      dop->Mult(X, Y);
+      y.SetSubVector(ess_tdofs, 0.0);
+   }
+
+   Operator& GetGradient(const Vector &x) const override
+   {
+      jacobian = std::make_shared<JacobianOperator>(*this, x);
+      return *jacobian;
+   }
+
+private:
+   ParFiniteElementSpace &fes;
+   ParFiniteElementSpace *mesh_nodes_fes = nullptr;
+   Vector mesh_nodes_tdofs;
+   Array<int> ess_tdofs;
+
+   std::shared_ptr<DifferentiableOperator> dop;
+   mutable std::shared_ptr<JacobianOperator> jacobian;
+};
+
+// BoomerAMG for the matrix-free tangent. HypreBoomerAMG needs a real
+// HypreParMatrix, so on every SetOperator we let dFEM assemble the tangent and
+// set up AMG. The assembled matrix carries both outputs, so the reaction term
+// is present in the preconditioner as well.
+template <int DIM>
+class JacobianAMG : public Solver
+{
+   using JacobianOperator =
+      typename ReactionDiffusionOperator<DIM>::JacobianOperator;
+
+public:
+   JacobianAMG() { amg.SetPrintLevel(0); }
+
+   void SetOperator(const Operator &op) override
+   {
+      const auto *Jop = dynamic_cast<const JacobianOperator *>(&op);
+      MFEM_VERIFY(Jop, "JacobianAMG requires a JacobianOperator");
+      height = width = op.Height();
+
+      delete A;
+      A = nullptr;
+      Jop->AssembleJacobian(A);
+      amg.SetOperator(*A);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override { amg.Mult(x, y); }
+
+   ~JacobianAMG() { delete A; }
+
+private:
+   HypreParMatrix *A = nullptr;
+   HypreBoomerAMG amg;
+};
+
+template <int DIM>
+std::unique_ptr<Solver> MakePreconditioner(PreconditionerType type)
+{
+   switch (type)
+   {
+      case PreconditionerType::None:     return nullptr;
+      case PreconditionerType::Diagonal: return
+            std::make_unique<OperatorJacobiSmoother>();
+      case PreconditionerType::AMG:      return std::make_unique<JacobianAMG<DIM>>();
+      default:
+         MFEM_ABORT("Unknown preconditioner type: " << static_cast<int>(type));
+   }
+   return nullptr;
+}
+
+// ----------------------------------------------------------------------------
+// Verification of the two-output path against a classical assembly.
+// ----------------------------------------------------------------------------
+
+/// Largest absolute entry of @a v across all ranks.
+static real_t GlobalNormlinf(MPI_Comm comm, const Vector &v)
+{
+   real_t local = v.Normlinf(), global = 0.0;
+   MPI_Allreduce(&local, &global, 1, MPITypeMap<real_t>::mpi_type, MPI_MAX,
+                 comm);
+   return global;
+}
+
+/// With beta = 0 and f = 0 the operator is exactly
+/// DiffusionIntegrator(kappa) + MassIntegrator(alpha), which a single
+/// ParBilinearForm can build. Comparing against it checks that the Gradient and
+/// the Value output really do accumulate into the same residual vector: drop
+/// either FieldOperator from the Outputs tuple and these comparisons fail.
+///
+/// The operator used here is built without essential attributes, so no
+/// elimination happens on either side and the raw action and matrix are
+/// compared.
+template <int DIM>
+static void CheckLinearLimit(ParFiniteElementSpace &fes,
+                             const IntegrationRule &ir,
+                             real_t kappa, real_t alpha)
+{
+   MPI_Comm comm = fes.GetComm();
+   const int tvsize = fes.GetTrueVSize();
+
+   ReactionDiffusionOperator<DIM> lin_op(fes, ir, kappa, alpha, 0.0, 0.0);
+
+   ConstantCoefficient kappa_coeff(kappa), alpha_coeff(alpha);
+   ParBilinearForm a(&fes);
+   a.AddDomainIntegrator(new DiffusionIntegrator(kappa_coeff, &ir));
+   a.AddDomainIntegrator(new MassIntegrator(alpha_coeff, &ir));
+   a.Assemble();
+   a.Finalize();
+
+   Vector U(tvsize);
+   U.Randomize(0x9e3779b9);
+   ParGridFunction u_gf(&fes), y_gf(&fes);
+   u_gf.SetFromTrueDofs(U);
+
+   // 1. Action of the two-output q-function against the two integrators.
+   Vector Y(tvsize), Yref(tvsize);
+   lin_op.Mult(U, Y);
+   a.Mult(u_gf, y_gf);
+   fes.GetProlongationMatrix()->MultTranspose(y_gf, Yref);
+
+   Vector diff(Y);
+   diff -= Yref;
+   const real_t action_err = GlobalNormlinf(comm, diff);
+
+   // 2. Assembled tangent and its diagonal. Both assembly paths walk the
+   //    quadrature point cache written by DerivativeSetup, which is laid out
+   //    over every output FieldOperator, so both have to pick up the reaction
+   //    term as well as the diffusion term. The matrix comparison is done
+   //    through matrix-vector products, so that it does not depend on the
+   //    sparsity pattern or the ordering of the two assemblies.
+   const auto &J = lin_op.GetGradient(U);
+   const auto *Jop =
+      dynamic_cast<const typename ReactionDiffusionOperator<DIM>::JacobianOperator *>
+      (&J);
+   MFEM_VERIFY(Jop, "expected a JacobianOperator");
+
+   HypreParMatrix *A = nullptr;
+   Jop->AssembleJacobian(A);
+   HypreParMatrix *Aref = a.ParallelAssemble();
+
+   real_t matrix_err = 0.0;
+   Vector v(tvsize), w1(tvsize), w2(tvsize);
+   for (int k = 0; k < 3; k++)
+   {
+      v.Randomize(0x01000193 + k);
+      A->Mult(v, w1);
+      Aref->Mult(v, w2);
+      w1 -= w2;
+      matrix_err = std::max(matrix_err, GlobalNormlinf(comm, w1));
+   }
+
+   Vector d1(tvsize), d2(tvsize);
+   Jop->AssembleDiagonal(d1);
+   Aref->GetDiag(d2);
+   d1 -= d2;
+   const real_t diag_err = GlobalNormlinf(comm, d1);
+
+   delete A;
+   delete Aref;
+
+   if (Mpi::Root())
+   {
+      mfem::out << "linear limit (beta = 0, f = 0) vs "
+                << "DiffusionIntegrator + MassIntegrator:\n"
+                << "  action   max |dFEM - reference| = " << action_err << '\n'
+                << "  tangent  max |dFEM - reference| = " << matrix_err << '\n'
+                << "  diagonal max |dFEM - reference| = " << diag_err
+                << std::endl;
+   }
+
+   // Loose enough for the accumulated round-off of two different assembly
+   // paths, tight enough that a missing output term cannot slip through.
+   const real_t tol = 1e-10 * std::max(1.0_r, kappa + alpha);
+   MFEM_VERIFY(action_err < tol, "two-output action does not match the "
+               "DiffusionIntegrator + MassIntegrator reference");
+   MFEM_VERIFY(matrix_err < tol, "assembled two-output tangent does not match "
+               "the DiffusionIntegrator + MassIntegrator reference");
+   MFEM_VERIFY(diag_err < tol, "diagonal of the two-output tangent does not "
+               "match the DiffusionIntegrator + MassIntegrator reference");
+}
+
+// ----------------------------------------------------------------------------
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // 2. Parse command-line options
+   int order = 1;
+   int refinements = 3;
+   const char *device_config = "cpu";
+   real_t kappa = 1.0;
+   real_t alpha = 1.0;
+   real_t beta = 1.0;
+   real_t source = 1.0;
+   int prec_type = static_cast<int>(PreconditionerType::AMG);
+   real_t krylov_tol = 1e-8;
+   bool check = false;
+   bool visualization = true;
+   bool paraview = false;
+   int visport = 19916;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&refinements, "-r", "--refinements",
+                  "Number of uniform refinements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&kappa, "-k", "--kappa", "Diffusion coefficient, positive.");
+   args.AddOption(&alpha, "-a", "--alpha",
+                  "Linear reaction coefficient, non-negative.");
+   args.AddOption(&beta, "-b", "--beta",
+                  "Cubic reaction coefficient, non-negative. Zero gives the "
+                  "linear screened Poisson problem.");
+   args.AddOption(&source, "-s", "--source", "Constant source term.");
+   args.AddOption(&prec_type, "-pc", "--preconditioner",
+                  "Preconditioner: 0 = none, 1 = diagonal/Jacobi, 2 = AMG.");
+   args.AddOption(&krylov_tol, "-tol", "--krylov-tol",
+                  "Relative tolerance for the linear solver.");
+   args.AddOption(&check, "-check", "--check-linear", "-no-check",
+                  "--no-check-linear",
+                  "Compare the beta = 0 limit against a ParBilinearForm "
+                  "carrying DiffusionIntegrator and MassIntegrator.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&paraview, "-pv", "--paraview", "-no-pv", "--no-paraview",
+                  "Enable or disable ParaView DataCollection output.");
+   args.AddOption(&visport, "-p", "--send-port", "Socket for GLVis.");
+   args.ParseCheck();
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   if (Mpi::Root())
+   {
+      mfem::out << "operator: -div(" << kappa << " grad u) + " << alpha
+                << " u + " << beta << " u^3 = " << source
+                << ", outputs: Gradient<Solution>, Value<Solution>"
+                << std::endl;
+   }
+
+   // 4. Create the unit square/cube mesh and refine it
+   Mesh mesh = (dim == 2)
+               ? Mesh::MakeCartesian2D(4, 4, Element::QUADRILATERAL)
+               : Mesh::MakeCartesian3D(4, 4, 4, Element::HEXAHEDRON);
+   for (int l = 0; l < refinements; l++) { mesh.UniformRefinement(); }
+   mesh.SetCurvature(order);
+
+   // 5. Define a parallel mesh
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   pmesh.EnsureNodes();
+
+   // 6. Define a scalar H1 space on the parallel mesh. GlobalTrueVSize is
+   //    collective, so every rank has to reach it.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fes(&pmesh, &fec, 1);
+   const HYPRE_BigInt global_dofs = fes.GlobalTrueVSize();
+   if (Mpi::Root()) { mfem::out << "#dofs: " << global_dofs << std::endl; }
+
+   // 7. Set up the integration rule. It has to resolve the cubic reaction
+   //    term, not just the bilinear part.
+   const IntegrationRule &ir =
+      IntRules.Get(pmesh.GetTypicalElementGeometry(), 2 * order + 2);
+
+   // 8. Optionally verify the two-output path against a classical assembly.
+   //    This uses its own operator, so it is independent of -b and -s.
+   if (check) { CheckLinearLimit<dim>(fes, ir, kappa, alpha); }
+
+   // 9. Create the nonlinear operator
+   ReactionDiffusionOperator<dim> rd_op(fes, ir, kappa, alpha, beta, source);
+
+   Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+   ess_bdr = 1;
+   rd_op.SetEssentialAttributes(ess_bdr);
+
+   ParGridFunction u_gf(&fes);
+   u_gf = 0.0;
+
+   Vector U(fes.GetTrueVSize());
+   u_gf.GetTrueDofs(U);
+
+   // 10. Set up the linear solver used within Newton's method. For
+   //     non-negative alpha and beta the tangent is symmetric positive
+   //     definite, so CG works.
+   CGSolver krylov(MPI_COMM_WORLD);
+   krylov.SetAbsTol(0.0);
+   krylov.SetRelTol(krylov_tol);
+   krylov.SetMaxIter(2000);
+   krylov.SetPrintLevel(0);
+
+   std::unique_ptr<Solver> pc =
+      MakePreconditioner<dim>(static_cast<PreconditionerType>(prec_type));
+   if (pc) { krylov.SetPreconditioner(*pc); }
+
+   // 11. Set up the nonlinear solver (Newton)
+   NewtonSolver newton(MPI_COMM_WORLD);
+   newton.SetOperator(rd_op);
+   newton.SetSolver(krylov);
+   newton.SetAbsTol(0.0);
+#ifdef MFEM_USE_SINGLE
+   newton.SetRelTol(1e-6);
+#else
+   newton.SetRelTol(1e-10);
+#endif
+   newton.SetMaxIter(30);
+   newton.SetPrintLevel(1);
+
+   // 12. Solve the nonlinear system using Newton's method
+   Vector zero;
+   newton.Mult(zero, U);
+   MFEM_VERIFY(newton.GetConverged(), "Newton did not converge");
+
+   u_gf.Distribute(U);
+
+   // 13. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank()
+               << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << u_gf << std::flush;
+   }
+
+   // 14. Save the solution in parallel using ParaView format
+   if (paraview)
+   {
+      ParaViewDataCollection pd("dfem-reaction-diffusion-output", &pmesh);
+      pd.RegisterField("solution", &u_gf);
+      pd.SetDataFormat(VTKFormat::BINARY);
+      if (order > 1)
+      {
+         pd.SetHighOrderOutput(true);
+         pd.SetLevelsOfDetail(order);
+      }
+      pd.SetCycle(0);
+      pd.SetTime(0.0);
+      pd.Save();
+   }
+
+   return 0;
+}

miniapps/dfem/makefile

@@ -21,7 +21,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_EXAMPLES =
-PAR_EXAMPLES = dfem-minimal-surface dfem-hyperelasticity_energy
+PAR_EXAMPLES = dfem-minimal-surface dfem-hyperelasticity dfem-nonlinear-poisson \
+   dfem-hyperbolic-heat dfem-reaction-diffusion
 ifeq ($(MFEM_USE_MPI),NO)
    EXAMPLES = $(SEQ_EXAMPLES)
 else
@@ -64,4 +65,5 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf dfem-minimal-surface-output
+	@rm -rf dfem-minimal-surface-output dfem-nonlinear-poisson-output \
+	   dfem-hyperbolic-heat-output dfem-reaction-diffusion-output Output

tests/unit/dfem/test_multiple_outputs.cpp

@@ -174,6 +174,56 @@ struct mass_diffusion_local_qf
    }
 };
 
+// Three outputs across two test fields: V carries mass + diffusion (two
+// outputs on one field), P carries a diffusion scaled by kappa. The
+// two blocks are different, so we should spot if the row blocks
+// get mismatched/corrupted inadvertedly
+constexpr real_t kappa = 2.0;
+
+struct two_field_local_qf
+{
+   inline MFEM_HOST_DEVICE
+   void operator()(
+      const dscalar_t &u,
+      const tensor<dscalar_t, DIM> &dudxi,
+      const tensor<real_t, DIM, DIM> &J,
+      const real_t &w,
+      dscalar_t &out_v,
+      tensor<dscalar_t, DIM> &out_dv,
+      tensor<dscalar_t, DIM> &out_dp) const
+   {
+      const auto invJ = inv(J);
+      const auto detJ = det(J);
+      const auto grad = (dudxi * invJ) * transpose(invJ);
+      out_v = u * detJ * w;
+      out_dv = grad * (detJ * w);
+      out_dp = grad * (kappa * detJ * w);
+   }
+};
+
+// Mass + diffusion on an FE test field, plus quadrature point data on a
+// VectorQuadratureSpace. The second row block has no basis to contract against,
+// so it is not assemblable while the first one still is.
+struct mass_diffusion_qdata_local_qf
+{
+   inline MFEM_HOST_DEVICE
+   void operator()(
+      const dscalar_t &u,
+      const tensor<dscalar_t, DIM> &dudxi,
+      const tensor<real_t, DIM, DIM> &J,
+      const real_t &w,
+      dscalar_t &out1,
+      tensor<dscalar_t, DIM> &out2,
+      tensor<real_t, DIM, DIM> &out3) const
+   {
+      const auto invJ = inv(J);
+      const auto detJ = det(J);
+      out1 = u * detJ * w;
+      out2 = (dudxi * invJ) * transpose(invJ) * (detJ * w);
+      out3 = J;
+   }
+};
+
 TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
 {
    const bool all_tests = launch_all_non_regression_tests;
@@ -496,7 +546,7 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
          ParBilinearForm blf_fa(&fes);
          blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
          blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
-         blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
+         blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACY);
          blf_fa.Assemble();
          blf_fa.Finalize();
 
@@ -506,6 +556,9 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
             {COORDINATES, nodes->ParFESpace()},
          };
 
+         // The test field is named V even though it is the same space as the
+         // trial field U, which is the usual "V is the test function" idiom.
+         // The diagonal of dR_V/dU is then named explicitly.
          const std::vector<FieldDescriptor> out_fds
          {
             {V, &fes},
@@ -542,7 +595,7 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
          SECTION("Multiple Outputs Assemble Diagonal")
          {
             Vector diag(fes.GetTrueVSize()), diag_ref_l(fes.GetVSize());
-            dRdU->AssembleDiagonal(diag);
+            dRdU->AssembleDiagonal(V, diag);
             blf_fa.SpMat().GetDiag(diag_ref_l);
 
             Vector diag_ref(fes.GetTrueVSize());
@@ -556,6 +609,312 @@ TEST_CASE("dFEM Multiple Outputs", "[Parallel][dFEM][GPU]")
             MPI_Barrier(MPI_COMM_WORLD);
          }
       }
+
+      // Outputs spanning two test fields. dR/dU is a block column, one row
+      // block per output field, all sharing the trial space of U, and each row
+      // block assembles into its own matrix.
+      // This would look smth like:
+      //
+      // dR/dU = [ dR_P/dU; dR_U/dU ]
+      //
+      {
+         static constexpr int P = 5;
+
+         ParBilinearForm blf_u(&fes);
+         blf_u.AddDomainIntegrator(new MassIntegrator(ir));
+         blf_u.AddDomainIntegrator(new DiffusionIntegrator(ir));
+         blf_u.SetAssemblyLevel(AssemblyLevel::LEGACY);
+         blf_u.Assemble();
+         blf_u.Finalize();
+
+         ConstantCoefficient kappa_coeff(kappa);
+         ParBilinearForm blf_p(&fes);
+         blf_p.AddDomainIntegrator(new DiffusionIntegrator(kappa_coeff, ir));
+         blf_p.SetAssemblyLevel(AssemblyLevel::LEGACY);
+         blf_p.Assemble();
+         blf_p.Finalize();
+
+         const std::vector<FieldDescriptor> in_fds
+         {
+            {U, &fes},
+            {COORDINATES, nodes->ParFESpace()},
+         };
+
+         // U is deliberately *not* the first output field: A[f] and the
+         // blocks of Mult are indexed by position in out_fds, and
+         // AssembleDiagonal has to pick the block whose field is the
+         // differentiated one rather than simply the first.
+         const std::vector<FieldDescriptor> out_fds
+         {
+            {P, &fes},
+            {U, &fes},
+         };
+
+         DifferentiableOperator dop(in_fds, out_fds, pmesh);
+
+         auto qf = two_field_local_qf{};
+         dop.AddDomainIntegrator<LocalQFBackend>(
+            qf,
+            tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
+            tuple{Value<U>{}, Gradient<U>{}, Gradient<P>{}},
+            *ir, all_domain_attr, Derivatives<U> {});
+
+         Vector nodestv;
+         nodes->GetTrueDofs(nodestv);
+         fes.GetRestrictionMatrix()->Mult(x, xtvec);
+         MultiVector X{xtvec, nodestv};
+         auto dRdU = dop.GetDerivative(U, X);
+
+         REQUIRE(dRdU->Height() == 2 * fes.GetTrueVSize());
+         REQUIRE(dRdU->Width() == fes.GetTrueVSize());
+
+         SECTION("Multiple Fields SparseMatrix")
+         {
+            std::vector<SparseMatrix *> A;
+            dRdU->Assemble(A);
+
+            REQUIRE(A.size() == 2);
+            REQUIRE(A[0] != nullptr);
+            REQUIRE(A[1] != nullptr);
+
+            // TestSameMatrices only goes thru the first matrix' sparsity pattern,
+            // so if we compare both ways we can catch entries missing from either side.
+            TestSameMatrices(*A[0], blf_p.SpMat());
+            TestSameMatrices(blf_p.SpMat(), *A[0]);
+            TestSameMatrices(*A[1], blf_u.SpMat());
+            TestSameMatrices(blf_u.SpMat(), *A[1]);
+
+            delete A[0];
+            delete A[1];
+
+            // The element matrix banks are reused, so assembling a second time
+            // (what a Newton loop does every iteration) has to give the same
+            // matrices instead of accumulating on top of the first ones.
+            dRdU->Assemble(A);
+            TestSameMatrices(*A[0], blf_p.SpMat());
+            TestSameMatrices(*A[1], blf_u.SpMat());
+
+            delete A[0];
+            delete A[1];
+            MPI_Barrier(MPI_COMM_WORLD);
+         }
+
+         SECTION("Multiple Fields HypreParMatrix")
+         {
+            std::vector<HypreParMatrix *> A;
+            dRdU->Assemble(A);
+
+            REQUIRE(A.size() == 2);
+            REQUIRE(A[0] != nullptr);
+            REQUIRE(A[1] != nullptr);
+
+            std::unique_ptr<HypreParMatrix> ref_u(blf_u.ParallelAssemble());
+            std::unique_ptr<HypreParMatrix> ref_p(blf_p.ParallelAssemble());
+
+            // The assembled blocks have to agree with the matrix free action
+            // block by block, which is also what pins A[f] to output field f.
+            Vector dir(fes.GetTrueVSize());
+            dir.Randomize(7);
+
+            Vector yp(fes.GetTrueVSize()), yu(fes.GetTrueVSize());
+            MultiVector Y{yp, yu};
+            dRdU->Mult(dir, Y);
+
+            Vector a(fes.GetTrueVSize());
+            A[0]->Mult(dir, a);
+            a -= Y[0];
+            real_t norm_l = a.Normlinf(), norm_g = norm_l;
+            MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
+                          pmesh.GetComm());
+            REQUIRE(norm_g == MFEM_Approx(0.0));
+            Vector r(fes.GetTrueVSize());
+            ref_p->Mult(dir, r);
+            a = Y[0];
+            a -= r;
+            norm_l = a.Normlinf();
+            MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
+                          pmesh.GetComm());
+            REQUIRE(norm_g == MFEM_Approx(0.0));
+
+            A[1]->Mult(dir, a);
+            a -= Y[1];
+            norm_l = a.Normlinf();
+            MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
+                          pmesh.GetComm());
+            REQUIRE(norm_g == MFEM_Approx(0.0));
+            ref_u->Mult(dir, r);
+            a = Y[1];
+            a -= r;
+            norm_l = a.Normlinf();
+            MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
+                          pmesh.GetComm());
+            REQUIRE(norm_g == MFEM_Approx(0.0));
+
+            delete A[0];
+            delete A[1];
+            MPI_Barrier(MPI_COMM_WORLD);
+         }
+
+         SECTION("Multiple Fields Assemble Diagonal")
+         {
+            // Both row blocks are square here, since P and U share a space, so
+            // squareness alone cannot pick one. They carry different operators
+            // though, so reading the wrong row block cannot pass unnoticed.
+            const auto check_diag = [&](const Vector &diag,
+                                        ParBilinearForm &ref)
+            {
+               REQUIRE(diag.Size() == fes.GetTrueVSize());
+
+               Vector diag_ref_l(fes.GetVSize());
+               ref.SpMat().GetDiag(diag_ref_l);
+               Vector diag_ref(fes.GetTrueVSize());
+               fes.GetProlongationMatrix()->MultTranspose(diag_ref_l, diag_ref);
+
+               Vector d(diag);
+               d -= diag_ref;
+               real_t norm_l = d.Normlinf(), norm_g = norm_l;
+               MPI_Allreduce(&norm_l, &norm_g, 1, MPI_DOUBLE, MPI_MAX,
+                             pmesh.GetComm());
+               REQUIRE(norm_g == MFEM_Approx(0.0));
+            };
+
+            // Default: the row block of the differentiated field, dR_U/dU.
+            Vector diag_u;
+            dRdU->AssembleDiagonal(diag_u);
+            check_diag(diag_u, blf_u);
+
+            // Named: the other row block, dR_P/dU.
+            Vector diag_p;
+            dRdU->AssembleDiagonal(P, diag_p);
+            check_diag(diag_p, blf_p);
+
+            MPI_Barrier(MPI_COMM_WORLD);
+         }
+      }
+
+      // The same two field structure, but now the two output fields live on
+      // different spaces. That is what makes AssembleDiagonal's choice of row
+      // block observable at all: with both fields on one space, reading the
+      // wrong block still lands on an identically sized vector holding the
+      // same numbers, so the check above cannot tell the two apart.
+      {
+         static constexpr int P = 5;
+
+         H1_FECollection fec_p(p + 1, DIM);
+         ParFiniteElementSpace fes_p(&pmesh, &fec_p);
+
+         ParBilinearForm blf_u(&fes);
+         blf_u.AddDomainIntegrator(new MassIntegrator(ir));
+         blf_u.AddDomainIntegrator(new DiffusionIntegrator(ir));
+         blf_u.SetAssemblyLevel(AssemblyLevel::LEGACY);
+         blf_u.Assemble();
+         blf_u.Finalize();
+
+         const std::vector<FieldDescriptor> in_fds
+         {
+            {U, &fes},
+            {COORDINATES, nodes->ParFESpace()},
+         };
+
+         // U is the second output field on purpose, and fes_p is a different
+         // space, so picking out_fds[0] would give a differently sized vector.
+         const std::vector<FieldDescriptor> out_fds
+         {
+            {P, &fes_p},
+            {U, &fes},
+         };
+
+         DifferentiableOperator dop(in_fds, out_fds, pmesh);
+
+         auto qf = two_field_local_qf{};
+         dop.AddDomainIntegrator<LocalQFBackend>(
+            qf,
+            tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
+            tuple{Value<U>{}, Gradient<U>{}, Gradient<P>{}},
+            *ir, all_domain_attr, Derivatives<U> {});
+
+         Vector nodestv;
+         nodes->GetTrueDofs(nodestv);
+         fes.GetRestrictionMatrix()->Mult(x, xtvec);
+         MultiVector X{xtvec, nodestv};
+         auto dRdU = dop.GetDerivative(U, X);
+
+         SECTION("Assemble Diagonal Picks The Square Block")
+         {
+            Vector diag;
+            dRdU->AssembleDiagonal(diag);
+
+            // Sized by U's space, not by the first output field's.
+            REQUIRE(diag.Size() == fes.GetTrueVSize());
+
+            Vector diag_ref_l(fes.GetVSize());
+            blf_u.SpMat().GetDiag(diag_ref_l);
+            Vector diag_ref(fes.GetTrueVSize());
+            fes.GetProlongationMatrix()->MultTranspose(diag_ref_l, diag_ref);
+
+            diag -= diag_ref;
+            real_t norm_l = diag.Normlinf(), norm_g = norm_l;
+            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);
+         }
+      }
+
+      // One assemblable row block and one that is not: S lives on quadrature
+      // points, so it has no basis to contract against and stays null, while
+      // the mass + diffusion block on V assembles as usual.
+      {
+         ParBilinearForm blf_fa(&fes);
+         blf_fa.AddDomainIntegrator(new MassIntegrator(ir));
+         blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
+         blf_fa.SetAssemblyLevel(AssemblyLevel::LEGACYFULL);
+         blf_fa.Assemble();
+         blf_fa.Finalize();
+
+         const std::vector<FieldDescriptor> in_fds
+         {
+            {U, &fes},
+            {COORDINATES, nodes->ParFESpace()},
+         };
+
+         const std::vector<FieldDescriptor> out_fds
+         {
+            {V, &fes},
+            {S, &vqs},
+         };
+
+         DifferentiableOperator dop(in_fds, out_fds, pmesh);
+
+         auto qf = mass_diffusion_qdata_local_qf{};
+         dop.AddDomainIntegrator<LocalQFBackend>(
+            qf,
+            tuple{Value<U>{}, Gradient<U>{}, Gradient<COORDINATES>{}, Weight{}},
+            tuple{Value<V>{}, Gradient<V>{}, Identity<S>{}},
+            *ir, all_domain_attr, Derivatives<U> {});
+
+         Vector nodestv;
+         nodes->GetTrueDofs(nodestv);
+         fes.GetRestrictionMatrix()->Mult(x, xtvec);
+         MultiVector X{xtvec, nodestv};
+         auto dRdU = dop.GetDerivative(U, X);
+
+         SECTION("Partly Assemblable Outputs SparseMatrix")
+         {
+            std::vector<SparseMatrix *> A;
+            dRdU->Assemble(A);
+
+            REQUIRE(A.size() == 2);
+            REQUIRE(A[0] != nullptr);
+            REQUIRE(A[1] == nullptr);
+
+            TestSameMatrices(*A[0], blf_fa.SpMat());
+            TestSameMatrices(blf_fa.SpMat(), *A[0]);
+
+            delete A[0];
+            MPI_Barrier(MPI_COMM_WORLD);
+         }
+      }
    }
 }