Silenced tee.

examples/ex1p.cpp

@@ -66,6 +66,26 @@
 using namespace std;
 using namespace mfem;
 
+class VisServer
+{
+   const char *path{};
+   int port{-1};
+
+   struct Deleter { void operator()(FILE *); };
+   std::unique_ptr<FILE, Deleter> fglvis;
+
+   void Open();
+public:
+   VisServer(const char *path_, int port_ = 19916)
+      : path(path_), port(port_) { Open(); }
+
+   explicit VisServer() = default;
+
+   inline const char *GetPath() const { return path; }
+   inline int GetPort() const { return port; }
+   inline bool IsOpen() const { return (fglvis != nullptr); }
+};
+
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI and HYPRE.
@@ -82,6 +102,7 @@ int main(int argc, char *argv[])
    bool fa = false;
    const char *device_config = "cpu";
    bool visualization = true;
+   const char *visserver = "";
    bool algebraic_ceed = false;
 
    OptionsParser args(argc, argv);
@@ -106,6 +127,8 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&visserver, "-vs", "--visualization-server",
+                  "Path to GLVis binary to start own server.");
    args.Parse();
    if (!args.Good())
    {
@@ -301,6 +324,14 @@ int main(int argc, char *argv[])
    {
       char vishost[] = "localhost";
       int  visport   = 19916;
+
+      VisServer server;
+      if (strlen(visserver) > 0)
+      {
+         server = VisServer(visserver, visport);
+         visport = server.GetPort();
+      }
+
       socketstream sol_sock(vishost, visport);
       sol_sock << "parallel " << num_procs << " " << myid << "\n";
       sol_sock.precision(8);
@@ -315,3 +346,66 @@ int main(int argc, char *argv[])
 
    return 0;
 }
+
+void VisServer::Open()
+{
+   if (!path || strlen(path) <= 0) { return; }
+
+#ifdef MFEM_USE_MPI
+   if (Mpi::Root())
+#endif // MFEM_USE_MPI
+   {
+      // Open a log file
+      FILE *ftmp = popen("mktemp --tmpdir glvis-log.XXXXXX", "r");
+      constexpr size_t ssize = 256;
+      char stmp[ssize];
+      fgets(stmp, ssize, ftmp);
+      pclose(ftmp);
+      int len = strlen(stmp);
+      if (len > 0 && stmp[len-1] == '\n') { stmp[--len] = '\0'; }
+      std::cout << "Starting GLVis log in: " << stmp << std::endl;
+
+      // Start the server
+      std::stringstream ss;
+      ss << "trap '' SIGPIPE && " << path <<
+         " -no-pr -no-ex 2>&1 | tee \"" << stmp <<
+         "\" 2> /dev/null | grep -m 1 ^GLVIS_SERVER_PORT";
+      fglvis.reset(popen(ss.str().c_str(), "r"));
+
+      // Capture the port number
+      char line[ssize];
+      bool captured = false;
+      while (fgets(line, ssize, fglvis.get()))
+      {
+         if (strncmp(line, "GLVIS_SERVER_PORT=", 18) == 0)
+         {
+            sscanf(line, "GLVIS_SERVER_PORT=%d", &port);
+            captured = true;
+            break;
+         }
+      }
+      if (captured)
+      {
+         std::cout << "Started GLVis server at port: " << port << std::endl;
+      }
+      else
+      {
+         std::cerr << "GLVis server did not start normally." << std::endl;
+      }
+   }
+#ifdef MFEM_USE_MPI
+   MPI_Bcast(&port, 1, MPI_INT, 0, MPI_COMM_WORLD);
+#endif // MFEM_USE_MPI
+}
+
+void VisServer::Deleter::operator()(FILE *f)
+{
+   if (f != nullptr)
+   {
+      int ierr = pclose(f);
+      if (ierr != 0)
+      {
+         std::cerr << "GLVis server pclose() returns: " << ierr << std::endl;
+      }
+   }
+}

fem/CMakeLists.txt

@@ -213,7 +213,6 @@ 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

@@ -1,263 +0,0 @@
-// 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,7 +53,6 @@
 #include "particleset.hpp"
 
 #include "dfem/doperator.hpp"
-#include "dfem/univarsolvers.hpp"
 
 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"

fem/prestriction.cpp

@@ -321,17 +321,12 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
    const int vd = vdim;
    const bool t = byvdim;
    const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
+   const int nsdofs = pfes.GetFaceNbrVSize();
    auto d_indices1 = scatter_indices1.Read();
    auto d_indices2 = scatter_indices2.Read();
    auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs: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_x_shared = Reshape(face_nbr_data.Read(),
+                             t?vd:nsdofs, t?nsdofs:vd);
    auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
    mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
    {
@@ -351,9 +346,8 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
          }
          else if (idx2>=threshold) // shared boundary
          {
-            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);
+            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
+                                              t?(idx2-threshold):c);
          }
          else // true boundary
          {

fem/restriction.cpp

@@ -1398,17 +1398,20 @@ 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];
-      // Encode the volume DOF index and element index
-      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
+      const int global_dof_elem2 = face_nbr_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,9 +317,6 @@ 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,7 +39,6 @@ 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

@@ -1,359 +0,0 @@
-// 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,52 +117,3 @@ 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