Fix style error

Credit to William Moses.

fem/CMakeLists.txt

@@ -213,6 +213,7 @@ set(HDRS
   dfem/qfunction_apply.hpp
   dfem/qfunction_transform.hpp
   dfem/tuple.hpp
+  dfem/univarsolvers.hpp
   dfem/util.hpp
   eltrans.hpp
   estimators.hpp

fem/dfem/univarsolvers.hpp

@@ -0,0 +1,263 @@
+// 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.
+
+/**
+ * @file univarsolvers.hpp
+ *
+ * @brief Solvers of functions of a single variable suitable for use in ∂FEM q-functions.
+ */
+
+#ifndef MFEM_UNIVARSOLVERS
+#define MFEM_UNIVARSOLVERS
+
+#include <cmath>
+#include <limits>
+
+#include "../../config/config.hpp"
+#include "../../general/enzyme.hpp"
+#include "../../general/error.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+// Currently needed to work around a bug in LLVM
+extern void __enzyme_double(void*, size_t);
+
+namespace mfem
+{
+
+namespace future
+{
+/// Representation of root search bounds
+struct Bounds
+{
+   real_t lower, upper;
+};
+
+
+/// Settings for univariate solver
+struct SolverSettings
+{
+   real_t residual_abs_tol; ///< Tolerance for convergence check on absolute value of residual
+   real_t residual_rel_tol; ///< Tolerance for convergence check on absolute value of current residual relative to absolute value of residual at initial guess
+   Bounds bounds; ///< Bounds on root
+   int max_iters = 50;
+};
+} // namespace future
+
+namespace internal
+{
+/// @cond DO_NOT_DOCUMENT
+
+using future::SolverSettings;
+
+// The noinline attribute is neccessary for Enzyme. If this function were to be
+//  inlined in the calling function, The custom derivative rules would not be
+// found (since the function they refer to would no longer exist).
+template <auto f, typename T>
+__attribute__((noinline))
+MFEM_HOST_DEVICE void SolveNewtonBisection_impl(const real_t* x0_ptr,
+                                                const T* p_ptr, const SolverSettings* settings_ptr, real_t* x_ptr)
+{
+   int max_iters = settings_ptr->max_iters;
+
+   const real_t& x0 = *x0_ptr;
+   const T& p = *p_ptr;
+   const SolverSettings& settings = *settings_ptr;
+   const real_t& left_bracket = settings.bounds.lower;
+   const real_t& right_bracket = settings.bounds.upper;
+   real_t& x = *x_ptr;
+   using std::abs;
+
+   auto fprime = [&p](real_t x)
+   {
+      real_t x_dot = 1.0;
+      return __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, x, x_dot, enzyme_const,
+                                      p);
+   };
+
+   real_t fl = f(left_bracket, p);
+   real_t fh = f(right_bracket, p);
+
+   // handle corner cases where one of the brackets is the root
+   if (abs(fl) < settings.residual_abs_tol)
+   {
+      x = left_bracket;
+      return;
+   }
+   else if (abs(fh) < settings.residual_abs_tol)
+   {
+      x = right_bracket;
+      return;
+   }
+
+   if (fl * fh > 0)
+   {
+      MFEM_WARNING("Root is not bracketed, solver may diverge.");
+   }
+
+   // clamp initial guess within root brackets
+   x = x0 > right_bracket? right_bracket : x0 < left_bracket? left_bracket : x0;
+
+   // Orient search so that f(xl) < 0
+   real_t xl = left_bracket;
+   real_t xh = right_bracket;
+   if (fl > 0.0)
+   {
+      xl = right_bracket;
+      xh = left_bracket;
+      real_t tmp = fl;
+      fl = fh;
+      fh = tmp;
+   }
+
+   real_t dx_old = abs(right_bracket - left_bracket);
+   real_t dx = dx_old;
+   x = x0;
+   real_t r = f(x, p);
+   real_t dr_dx = fprime(x);
+   real_t r0 = r;
+   for (int i = 0; i < max_iters; i++)
+   {
+      if ((((x - xh) * dr_dx - r)*((x - xl)*dr_dx - r) >= 0.0) ||
+          // Newton out of range
+          (std::abs(2.0*r) > std::abs(
+              dx_old*dr_dx)))              // Newton decreasing bracket slower than bisection
+      {
+         // Take bisection step
+         dx_old = dx;
+         dx = 0.5*(xh - xl);
+         real_t x_old = x;
+         x = xl + dx;
+         if (x == x_old) { return; }
+      }
+      else
+      {
+         // Take Newton step
+         dx_old = dx;
+         dx = -r/dr_dx;
+         real_t x_old = x;
+         x += dx;
+         if (x == x_old) { return; }
+      }
+
+      // update residual and jacobian
+      r = f(x, p);
+      dr_dx = fprime(x);
+
+      // Check convergence
+      if (abs(r) < settings.residual_rel_tol*r0 ||
+          abs(r) < settings.residual_abs_tol)
+      {
+         return;
+      }
+
+      // Update bracket
+      if (r < 0.0)
+      {
+         xl = x;
+         fl = r;
+      }
+      else
+      {
+         xh = x;
+         fh = r;
+      }
+   }
+   MFEM_ABORT("Univariate solve did not converge.");
+}
+
+template <auto f, typename T>
+void SolveNewtonBisection_impl_fwddiff(const real_t* x0,
+                                       const real_t* /* unused shadow */,
+                                       const T* p, const T* dp,
+                                       const SolverSettings* settings, const SolverSettings* /* unused shadow */,
+                                       real_t* x, real_t* dx)
+{
+   SolveNewtonBisection_impl<f>(x0, p, settings, x);
+   real_t dfdx = __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, *x, 1.0,
+                                          enzyme_const, *p);
+   real_t dfdp = __enzyme_fwddiff<real_t>((void*)+f, enzyme_const, *x, enzyme_dup,
+                                          *p, *dp);
+   *dx = -dfdp/dfdx;
+}
+
+
+template<auto f, typename T>
+void SolveNewtonBisection_impl_aug(const real_t* x0, real_t* x0_bar,
+                                   const T* p, T* p_bar,
+                                   const SolverSettings* settings, SolverSettings* settings_bar,
+                                   real_t* x, real_t* x_bar)
+{
+   SolveNewtonBisection_impl<f>(x0, p, settings, x);
+}
+
+// Change the residual function to return-by-reference so that there is a
+// slot to provide the downstream cotangent (ie the shadow for y)
+// in the reverse mode call.
+template<auto f, typename T>
+void rbr_wrapper(real_t x, T& p, real_t& y)
+{
+   y = f(x, p);
+}
+
+template<auto f, typename T>
+void SolveNewtonBisection_impl_rev(const real_t* x0, real_t* x0_bar,
+                                   const T* p, T* p_bar,
+                                   const SolverSettings* settings, SolverSettings* settings_bar,
+                                   real_t* x, real_t* x_bar)
+{
+   real_t drdx = __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, *x, 1.0,
+                                          enzyme_const, *p);
+   real_t lambda = -(*x_bar / drdx);
+   real_t r;
+   __enzyme_autodiff<void>((void*)rbr_wrapper<f, T>, enzyme_const, *x, enzyme_dup,
+                           p, p_bar, enzyme_dupnoneed, &r, &lambda);
+
+   // These are logically constants, the root has no sensitivity to these
+   *x0_bar = 0.0;
+   *settings_bar = SolverSettings{};
+}
+
+/// @endcond
+} // namespace internal
+
+
+namespace future
+{
+
+/**
+ * @brief Find the root of a univariate funtion
+ */
+template<auto f, typename T>
+MFEM_HOST_DEVICE __attribute__((always_inline)) real_t SolveNewtonBisection(
+   real_t x0, T p, SolverSettings settings)
+{
+   // We need to tell Enzyme how much memory in the SolverSettings object is
+   // used by active variables (in the sense of Enzyme activity analysis).
+   // Without this, it seems that a bug in LLVM causes this information to
+   // be lost during some optimization pass, and the Enzyme pass fails in
+   // Release builds.
+   // There are 4 real_t members in settings, which is what Enzyme will
+   // consider active.
+   // TODO: File an issue on Enzyme to remind Bill to fix this in LLVM.
+   __enzyme_double((void*)&settings, sizeof(real_t)*4);
+
+   real_t x;
+   internal::SolveNewtonBisection_impl<f>(&x0, &p, &settings, &x);
+   return x;
+}
+
+
+} // namespace future
+} // namespace mfem
+
+#endif // MFEM_USE_ENZYME
+#endif // MFEM_UNIVARSOLVERS

fem/fem.hpp

@@ -53,6 +53,7 @@
 #include "particleset.hpp"
 
 #include "dfem/doperator.hpp"
+#include "dfem/univarsolvers.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/prestriction.cpp

@@ -321,12 +321,17 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
    const int vd = vdim;
    const bool t = byvdim;
    const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize();
+   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
    auto d_indices1 = scatter_indices1.Read();
    auto d_indices2 = scatter_indices2.Read();
    auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   auto d_x_shared = Reshape(face_nbr_data.Read(),
-                             t?vd:nsdofs, t?nsdofs:vd);
+   const int ne_shared = nsdofs / elem_dofs;
+   const int nedof = elem_dofs;
+   // Note: the shape of face_nbr_data, as determined by
+   // ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
+   // ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
+   // element space.
+   auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
    auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
    mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
    {
@@ -346,8 +351,9 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
          }
          else if (idx2>=threshold) // shared boundary
          {
-            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
-                                              t?(idx2-threshold):c);
+            const int e_shared = (idx2 - threshold) / nedof;
+            const int i_shared = (idx2 - threshold) % nedof;
+            d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
          }
          else // true boundary
          {

fem/restriction.cpp

@@ -1398,20 +1398,17 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
    const int dim = fes.GetMesh()->Dimension();
    const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
    fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
-   Array<int> face_nbr_dofs;
-   const ParFiniteElementSpace &pfes =
-      static_cast<const ParFiniteElementSpace&>(this->fes);
-   pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);
 
    for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
    {
       const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
                                                orientation, dof1d, face_dof_elem1);
       const int volume_dof_elem2 = face_map[face_dof_elem2];
-      const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
+      // Encode the volume DOF index and element index
+      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
       const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
       // Trick to differentiate dof location inter/shared
-      scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
+      scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
    }
 #endif
 }

linalg/hypre.cpp

@@ -317,6 +317,9 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
 
 Vector * HypreParVector::GlobalVector() const
 {
+   MFEM_VERIFY(size > 0,
+               "GlobalVector method can only be called on vectors wherein each "
+               "process owns one or more entries");
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
    Vector *v = new Vector(hv->data, internal::to_int(hv->size));
    v->MakeDataOwner();

tests/unit/CMakeLists.txt

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

tests/unit/dfem/test_univarsolver.cpp

@@ -0,0 +1,359 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include <algorithm>
+#include <cmath>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+using mfem::real_t;
+using namespace mfem::future;
+
+MFEM_HOST_DEVICE inline real_t FlowResistance(real_t eqps, real_t sigma_y, real_t n, real_t ep_0)
+{
+    return sigma_y*(1.0 + std::pow((eqps)/ep_0, n));
+}
+
+using J2PlasticityParameters = tuple<real_t, real_t, real_t, real_t, real_t, real_t>;
+
+// Residual function that is solved in the plasticity model.
+// Made a free function to facilitate Enzyme differentiation.
+real_t J2PlasticityResidual(real_t delta_eqps, J2PlasticityParameters p)
+{
+    auto [eqps, q, G, sigma_y, n, ep_0] = p;
+    return q - 3.0*G*delta_eqps - FlowResistance(eqps + delta_eqps, sigma_y, n, ep_0);
+}
+
+struct J2Plasticity {
+  static constexpr int dim = 3;  ///< spatial dimension
+  static constexpr int N_INTERNAL_STATES = 10;
+  static constexpr real_t tol = 1e-10;  ///< relative tolerance on residual mag to judge convergence of return map
+
+  real_t E;        ///< Young's modulus
+  real_t nu;       ///< Poisson's ratio
+  real_t sigma_y;  ///< Yield strength
+  real_t n;        ///< Hardening index
+  real_t ep_0;     ///< Reference plastic strain
+  
+  /// @brief variables required to characterize the hysteresis response
+  struct InternalState {
+    tensor<real_t, dim, dim> plastic_strain;
+    real_t accumulated_plastic_strain;
+  };
+
+  /// Internal state variables in a flattened array for storing in a global field
+  using PackedInternalState = mfem::future::tensor<real_t, N_INTERNAL_STATES>;
+
+  // Unflatten internal state variables
+  MFEM_HOST_DEVICE static inline InternalState unpack_internal_state(
+      const mfem::future::tensor<real_t, N_INTERNAL_STATES>& packed_state)
+  {
+    auto plastic_strain =
+        mfem::future::make_tensor<dim, dim>([&packed_state](int i, int j) { return packed_state[dim * i + j]; });
+    real_t accumulated_plastic_strain = packed_state[N_INTERNAL_STATES - 1];
+    return {plastic_strain, accumulated_plastic_strain};
+  }
+
+  // Flatten internal state variables (for repacking into global field)
+  MFEM_HOST_DEVICE static inline PackedInternalState pack_internal_state(
+      const mfem::future::tensor<real_t, dim, dim>& plastic_strain, real_t accumulated_plastic_strain)
+  {
+    PackedInternalState packed_state{};
+    for (int i = 0, ij = 0; i < dim; i++) {
+      for (int j = 0; j < dim; j++, ij++) {
+        packed_state[ij] = plastic_strain[i][j];
+      }
+    }
+    packed_state[N_INTERNAL_STATES - 1] = accumulated_plastic_strain;
+    return packed_state;
+  }
+
+  // Compute the new stress and the internal state variables
+  MFEM_HOST_DEVICE inline tuple<tensor<real_t, dim, dim>, PackedInternalState>
+  update(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto invJ = inv(J);
+    const auto dudX = dudxi * invJ;
+    auto I = IdentityMatrix<dim>();
+    const real_t K = E / (3.0 * (1.0 - 2.0 * nu));
+    const real_t G = 0.5 * E / (1.0 + nu);
+
+    auto [plastic_strain, accumulated_plastic_strain] = unpack_internal_state(internal_state);
+    
+    auto el_strain = sym(dudX) - plastic_strain;
+    auto p = K * tr(el_strain);
+    auto s = 2.0 * G * dev(el_strain);
+    auto q = std::sqrt(1.5) * norm(s);
+    real_t denom = q > 0.0? q : 1.0;
+    auto Np = 1.5 * s / denom;
+
+    if (q > FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0)) {
+        real_t lb = 0.0;
+        real_t ub = (q - FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0))/(3*G);
+        SolverSettings settings{1e-10*sigma_y, 1e-10, {lb, ub}};
+        // Use the differentiable univariate root finder.
+        // This has custom derivatives, so it's ok to differentiate this enclosing function.
+        real_t delta_eqps = SolveNewtonBisection<J2PlasticityResidual>(
+            0.5*(lb + ub), make_tuple(accumulated_plastic_strain, q, G, sigma_y, n, ep_0), settings);
+        accumulated_plastic_strain += delta_eqps;
+        plastic_strain += delta_eqps * Np;
+        s -= 2.0 * G * delta_eqps * Np;
+    }
+    auto Q_new = pack_internal_state(plastic_strain, accumulated_plastic_strain);
+    auto stress = s + p * I;
+    const real_t dV = det(J)*w;
+    return {stress*transpose(invJ)*dV, Q_new};
+  }
+  
+  MFEM_HOST_DEVICE inline tensor<real_t, dim, dim>
+  stress(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return stress;
+  }
+
+  MFEM_HOST_DEVICE inline PackedInternalState
+  internal_state_new(tensor<real_t, dim, dim> dudxi,
+                     PackedInternalState internal_state,
+                     tensor<real_t, dim, dim> J,
+                     real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return internal_state_new;
+  }
+};
+
+
+// Register the custom derivatives for the solver.
+// This needs to be done for every residual function that the solver is applied on,
+// since the SolveNewtonBisection_impl is a function template, and we need a real
+// function with an address to specify the custom derivative.
+
+// Forward mode
+__attribute__((used))
+void *  __enzyme_register_derivative_newton_bisection_on_j2[2] = {
+  (void*) mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*) mfem::internal::SolveNewtonBisection_impl_fwddiff<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+ // Reverse mode
+__attribute__((used))
+void* __enzyme_register_gradient_SolveNewtonBisectionJ2[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+// Create free functions for Enzyme to differentiate in the tests
+// Return by value version
+tensor<real_t, 3, 3> ComputeStress(
+    J2Plasticity* material, tensor<real_t, 3, 3> dudxi, 
+    J2Plasticity::PackedInternalState Q, tensor<real_t, 3, 3> J, real_t w)
+{
+    return material->stress(dudxi, Q, J, w);
+}
+
+// Return by reference version
+void ComputeStressRef(const J2Plasticity* material, const tensor<real_t, 3, 3>& dudxi,
+                      const J2Plasticity::PackedInternalState& Q,
+                      const tensor<real_t, 3, 3>& J, real_t w,
+                      tensor<real_t, 3, 3>& sigma)
+{
+  sigma = material->stress(dudxi, Q, J, w);
+}
+
+template <int dim>
+real_t elementwise_max_norm(tensor<real_t, dim, dim> A) {
+  real_t maxval = 0;
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      maxval = std::max(std::abs(A[i][j]), maxval);
+    }
+  }
+  return maxval;
+}
+
+
+TEST_CASE("Univariate function solver in a qfunction", "[univar]")
+{
+    J2Plasticity material{70.0e3, 0.34, 240.0, 0.15, 1e-3};
+    tensor<real_t, 3, 3> H{{{0.947667  , 0.9785799 , 0.33229148},
+                            {0.46866846, 0.5698887 , 0.16550303},
+                            {0.3101946 , 0.68948054, 0.74676657}}};
+    J2Plasticity::PackedInternalState Q{};
+    const tensor<real_t, 3, 3> J = IdentityMatrix<3>();
+    const real_t w = 1.0;
+
+    SECTION("Correctness")
+    {
+      // Checks that stress after update is on the yield surface.
+      auto [stress, Q_new] = material.update(H, Q, IdentityMatrix<3>(), 1.0);
+      real_t mises = std::sqrt(1.5)*norm(dev(stress));
+      real_t eqps = Q_new[9];
+      // This test only makes sense if the displacement gradient is big enough to 
+      // cuase yielding.
+      REQUIRE(eqps > 1e-9);
+      real_t Y = FlowResistance(eqps, material.sigma_y, material.n, material.ep_0);
+      CHECK(mises == MFEM_Approx(Y, 0.0, 1e-8));
+    }
+
+    SECTION("JVP")
+    {
+      // Compare forward mode derivative to finite difference approximation
+
+      tensor<real_t, 3, 3> H_dot{{{1.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0}}};
+
+      // Enzyme directional derivative (uses custom derivative of solver)
+      auto sigma_dot = __enzyme_fwddiff<tensor<real_t, 3, 3>>((void*)ComputeStress, 
+                                                              enzyme_const, &material,
+                                                              enzyme_dup, H, H_dot,
+                                                              enzyme_const, Q,
+                                                              enzyme_const, J,
+                                                              enzyme_const, w);
+      // sigma_dot = ∂sigma / ∂H[0, 0]
+      REQUIRE(sigma_dot[0][0] > 0.0);
+
+      // Finite difference derivative approximation
+      constexpr int dim = 3;
+      real_t eps = 1e-5;
+      tensor<real_t, 3, 3> sigma = ComputeStress(&material, H, Q, J, w);
+      tensor<real_t, 3, 3> sigma_p = ComputeStress(&material, H + eps*H_dot, Q, J, w);
+      tensor<real_t, 3, 3> sigma_dot_h = (1.0/eps)*(sigma_p - sigma);
+
+      tensor<real_t, 3, 3> rel_error = sigma_dot - sigma_dot_h;
+      for (int i = 0; i < dim; i++) {
+        for (int j = 0; j < dim; j++) {
+          real_t denom = sigma[i][j] != 0? sigma[i][j] : 1.0;
+          rel_error[i][j] /= denom;
+        }
+      }
+
+      CHECK(elementwise_max_norm(rel_error) < 1e-5);
+    }
+
+    SECTION("VJP")
+    {
+      // compare reverse mode derivative to finite differences
+
+      tensor<real_t, 3, 3> sigma;
+      ComputeStressRef(&material, H, Q, J, w, sigma);
+      double epsilon = 1e-6;
+      tensor<real_t, 3, 3> dH{{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}};
+      auto H_p = H + epsilon*dH;
+      tensor<real_t, 3, 3> sigma_p;
+      ComputeStressRef(&material, H_p, Q, J, w, sigma_p);
+      auto sigma_dot_h = (sigma_p - sigma)/epsilon;
+      // Note: sigma_dot_h[i,j] = ∂sigma[i,j]/∂H[0,0]
+
+      tensor<real_t, 3, 3> sigma_bar{{{1.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0}}};
+
+      tensor<real_t, 3, 3> H_bar{};
+      J2Plasticity::PackedInternalState Q_bar{};
+      tensor<real_t, 3, 3> J_bar{};
+      __enzyme_autodiff<void>(
+        (void*)ComputeStressRef,  enzyme_const, &material, enzyme_dup, &H, &H_bar,
+        enzyme_dup, &Q, &Q_bar, enzyme_dup, &J, &J_bar, enzyme_const, w,
+        enzyme_dup, &sigma, &sigma_bar);
+
+      // H_bar[ij] = ∂sigma[0,0]/∂H[i,j]
+      // For this model, we expect the major symmetries in the tangent operator.
+      // Hence H_bar \approx sigma_dot_h
+
+      const double abs_tol = 1e-12;
+      const double rel_tol = 5e-6;
+      for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+          CHECK(H_bar[i][j] == MFEM_Approx(sigma_dot_h[i][j], abs_tol, rel_tol));
+        }
+      }
+    }
+}
+
+
+real_t nthroot_res(real_t x, tuple<real_t, real_t> p)
+{
+    auto [index, radicand] = p;
+    return std::pow(x, index) - radicand;
+}
+
+__attribute__((used))
+void* __enzyme_register_gradient_solver[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<nthroot_res, tuple<real_t, real_t>>
+};
+
+TEST_CASE("Univariate solver reverse mode derivative", "[univar]")
+{
+  auto mysqrt = [](real_t x) -> real_t
+    {
+      real_t x0 = x;
+      real_t index = 2.0;
+      real_t ub = std::max(1.0, x);
+      SolverSettings settings{1e-12, 1e-12, {0, ub}};
+      return SolveNewtonBisection<nthroot_res>(x0, make_tuple(index, x), settings);
+    };
+
+  real_t x = 2.0;
+  real_t dydx = __enzyme_autodiff<real_t>((void*)+mysqrt, enzyme_out, x);
+  CHECK(dydx == MFEM_Approx(0.5/std::sqrt(2.0)));
+}
+
+TEST_CASE("Univariate function solver robustness", "[univar]")
+{
+  SolverSettings settings{1e-12, 1e-12};
+
+  SECTION("Simple case")
+  {
+      auto Nthroot = [&settings](real_t x, real_t n) {
+          real_t x0 = std::max(x, 1.0);
+          settings.bounds = {0.0, x0};
+          return SolveNewtonBisection<nthroot_res>(x0, make_tuple(n, x), settings);
+      };
+      real_t x = 8.0;
+      real_t y = Nthroot(x, 3.0);
+      CHECK(y == MFEM_Approx(2.0));
+  }
+
+    SECTION("Stiff problem")
+    {
+        auto f = [](real_t x, real_t p) { return std::pow(x, p) - 1.0; };
+        real_t x0 = 0.1;
+        real_t p = 50;
+        settings.bounds = {0.0, 5.1};
+        real_t x = SolveNewtonBisection<+f>(x0, p, settings);
+        CHECK(x == MFEM_Approx(1.0));
+    }
+
+    SECTION("Works where Newton diverges")
+    {
+        auto f = [](double x, int) { return std::atan(x); };
+        real_t x0 = 1.5;
+        settings.bounds = {0.0, 2.0};
+        real_t x = SolveNewtonBisection<+f>(x0, int{}, settings);
+        CHECK(std::abs(x) == MFEM_Approx(0.0));
+    }
+}
+
+#endif // MFEM_USE_ENZYME

tests/unit/fem/test_face_restriction.cpp

@@ -117,3 +117,52 @@ TEST_CASE("Vector FE Face Restriction", "[FaceRestriction]")
    gf2 -= gf;
    REQUIRE(gf2.Normlinf() == MFEM_Approx(0.0));
 }
+
+#ifdef MFEM_USE_MPI
+
+TEST_CASE("L2 Face Restriction", "[FaceRestriction][Parallel]")
+{
+   const int dim = GENERATE(2, 3);
+   constexpr int nx = 3;
+   constexpr int order = 2;
+   constexpr int vdim = 2;
+   const Ordering::Type ordering = GENERATE(Ordering::byNODES, Ordering::byVDIM);
+
+   Mesh serial_mesh = MakeCartesianMesh(nx, dim);
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+
+   L2_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&mesh, &fec, vdim, ordering);
+
+   auto *R = fes.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC,
+                                    FaceType::Interior);
+
+   Vector vals({1.0, 2.0});
+   VectorConstantCoefficient coeff(vals);
+
+   ParGridFunction gf(&fes);
+   gf.ProjectCoefficient(coeff);
+
+   Vector face_vec(R->Height());
+   R->Mult(gf, face_vec);
+
+   const int nf = mesh.GetNFbyType(FaceType::Interior);
+   const int face_dofs = fes.GetTypicalTraceElement()->GetDof();
+   auto h_face_vec = Reshape(face_vec.HostRead(), face_dofs, vdim, 2, nf);
+
+   for (int f = 0; f < nf; ++f)
+   {
+      for (int m = 0; m < 2; ++m)
+      {
+         for (int c = 0; c < vdim; ++c)
+         {
+            for (int i = 0; i < face_dofs; ++i)
+            {
+               REQUIRE(h_face_vec(i, c, m, f) == vals[c]);
+            }
+         }
+      }
+   }
+}
+
+#endif