Merge branch 'master' into dfem-lvector-interface

.github/workflows/builds-and-tests.yml

@@ -132,14 +132,12 @@ jobs:
             hypre-target: int32
             precision: fp64
             enzyme: true
-            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme) LDFLAGS=-L$LLVM_PREFIX/lib/c++
+            config-opts: MFEM_USE_ENZYME=YES ENZYME_DIR=$(brew --prefix enzyme)
 
     name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}
 
     runs-on: ${{ matrix.os }}
 
-    continue-on-error: ${{ matrix.enzyme && true || false }}
-
     steps:
       # Fix 'No space left on device' errors for Ubuntu builds.
       - name: Run Actions Cleaner
@@ -294,12 +292,10 @@ jobs:
         run: |
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
-          brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
-          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
-          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
-          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
+          brew install llvm@20 enzyme
+          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
+          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV
 
       # MFEM build and test
       - name: build

fem/dfem/doperator.hpp

@@ -231,22 +231,53 @@ public:
       const std::vector<FieldDescriptor> &parameters,
       const ParMesh &mesh);
 
+   /// MultLevel enum to indicate if the T->L Operators are used in the
+   /// Mult method.
+   enum MultLevel
+   {
+      TVECTOR,
+      LVECTOR
+   };
+
+   /// @brief Set the MultLevel mode for the DifferentiableOperator.
+   /// The default is TVECTOR, which means that the Operator will use
+   /// T->L before Mult and L->T Operators after.
+   void SetMultLevel(MultLevel level)
+   {
+      mult_level = level;
+   }
+
    /// @brief Compute the action of the operator on a given vector.
    ///
-   /// @param solutions_t The solution vector in which to compute the action.
-   /// This has to be a T-dof vector.
-   /// @param result_t Result vector of the action of the operator on
-   /// solutions_t. The result is a T-dof vector.
-   void Mult(const Vector &solutions_t, Vector &result_t) const override
+   /// @param solutions_in The solution vector in which to compute the action.
+   /// This has to be a T-dof vector if MultLevel is set to TVECTOR, or L-dof
+   /// Vector if MultLevel is set to LVECTOR.
+   /// @param result_in Result vector of the action of the operator on
+   /// solutions. The result is a T-dof vector or L-dof vector depending on
+   /// the MultLevel.
+   void Mult(const Vector &solutions_in, Vector &result_in) const override
    {
       MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
-      prolongation(solutions, solutions_t, solutions_l);
-      residual_l = 0.0;
-      for (auto &action : action_callbacks)
+
+      if (mult_level == MultLevel::LVECTOR)
       {
-         action(solutions_l, parameters_l, residual_l);
+         get_lvectors(solutions, solutions_in, solutions_l);
+         result_in = 0.0;
+         for (auto &action : action_callbacks)
+         {
+            action(solutions_l, parameters_l, result_in);
+         }
+      }
+      else
+      {
+         prolongation(solutions, solutions_in, solutions_l);
+         residual_l = 0.0;
+         for (auto &action : action_callbacks)
+         {
+            action(solutions_l, parameters_l, residual_l);
+         }
+         prolongation_transpose(residual_l, result_in);
       }
-      prolongation_transpose(residual_l, result_t);
    }
 
    /// @brief Add a domain integrator to the operator.
@@ -345,6 +376,8 @@ public:
 private:
    const ParMesh &mesh;
 
+   MultLevel mult_level = TVECTOR;
+
    std::vector<action_t> action_callbacks;
    std::map<size_t,
        std::vector<derivative_action_t>> derivative_action_callbacks;
@@ -354,7 +387,6 @@ private:
        std::vector<assemble_derivative_hypreparmatrix_callback_t>>
        assemble_derivative_hypreparmatrix_callbacks;
 
-
    std::vector<FieldDescriptor> solutions;
    std::vector<FieldDescriptor> parameters;
    // solutions and parameters

fem/dfem/util.hpp

@@ -1076,6 +1076,24 @@ void prolongation(const std::vector<FieldDescriptor> fields,
    }
 }
 
+inline
+void get_lvectors(const std::vector<FieldDescriptor> fields,
+                  const Vector &x,
+                  std::vector<Vector> &fields_l)
+{
+   int data_offset = 0;
+   for (std::size_t i = 0; i < fields.size(); i++)
+   {
+      const int sz = GetVSize(fields[i]);
+      fields_l[i].SetSize(sz);
+
+      const Vector x_i(const_cast<Vector&>(x), data_offset, sz);
+      fields_l[i] = x_i;
+
+      data_offset += sz;
+   }
+}
+
 /// @brief Get a transpose prolongation callback for a field descriptor.
 ///
 /// In the special case of a one field operator, the transpose prolongation

fem/fespace.hpp

@@ -922,9 +922,6 @@ public:
    { return mesh->GetBdrElementType(i); }
 
    /// Returns ElementTransformation for the @a i-th element.
-   /// @note The returned pointer references an object owned by the associated
-   /// @a Mesh that will be modified by other calls to `GetElementTransformation`.
-   /// As such, this pointer should @b not be deleted by the caller.
    ElementTransformation *GetElementTransformation(int i) const
    { return mesh->GetElementTransformation(i); }
 

linalg/CMakeLists.txt

@@ -23,7 +23,6 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
-  functional.cpp
   symmat.cpp
   handle.cpp
   matrix.cpp
@@ -50,7 +49,6 @@ list(APPEND HDRS
   complex_operator.hpp
   constraints.hpp
   densemat.hpp
-  functional.hpp
   dinvariants.hpp
   symmat.hpp
   dtensor.hpp

linalg/functional.cpp

@@ -1,68 +0,0 @@
-#include "functional.hpp"
-
-namespace mfem
-{
-QuadraticFunctional::QuadraticFunctional(const Operator *A_,
-                                         const Vector *b_, const real_t beta_, const real_t c_)
-   : Functional(A_ ? A_->Width() : 0)
-   , A(A_), beta(beta_), b(b_), c(c_)
-   , aux(A_ ? A_->Width() : 0)
-{
-   MFEM_VERIFY(A_ != nullptr,
-               "QuadraticFunctional: A must not be nullptr. "
-               << "Use QuadraticFunctional() constructor to create an empty Quadratic functional.");
-   MFEM_VERIFY(A_->Width() == A_->Height(),
-               "QuadraticFunctional: A must be a square operator.");
-   MFEM_VERIFY(b_ == nullptr || A_->Width() == b_->Size(),
-               "QuadraticFunctional: A and b must have compatible sizes");
-}
-
-void QuadraticFunctional::SetOperator(const Operator &A_)
-{
-   MFEM_VERIFY(A_.Width() == A_.Height(),
-               "QuadraticFunctional: A must be a square operator.");
-   A = &A_;
-   width = A_.Width();
-   aux.SetSize(width);
-}
-
-void QuadraticFunctional::SetVector(const Vector &b_, const real_t beta_)
-{
-   MFEM_VERIFY(A != nullptr && A->Width() == b_.Size(),
-               "QuadraticFunctional: A and b must have compatible sizes.");
-   b = &b_;
-   beta = beta_;
-}
-void QuadraticFunctional::SetConstant(real_t c_) { c = c_; }
-
-void QuadraticFunctional::Mult(const Vector &x, Vector &y) const
-{
-   MFEM_ASSERT(A != nullptr, "QuadraticFunctional: A must not be nullptr");
-   A->Mult(x, aux); // aux = A(x)
-   if (b != nullptr) { aux.Add(2.0*beta, *b); } // aux = A(x) + 2*beta*b
-   real_t result = 0.0;
-#ifdef MFEM_USE_MPI
-   if (IsParallel())
-   {
-      result = InnerProduct(GetComm(), x, aux); // result = <A(x) + 2*beta*b, x>
-   }
-   else
-   {
-      result = InnerProduct(x, aux); // result = <A(x) + 2*beta*b, x>
-   }
-#else
-   result = InnerProduct(x, aux); // result  = <A(x) + 2*beta*b, x>
-#endif
-   y.SetSize(1);
-   y[0] = result*0.5 + c; // y = 0.5 * <A(x) + 2*beta*b, x> + c
-}
-
-void QuadraticFunctional::EvalGradient(const Vector &x,
-                                       Vector &y) const
-{
-   y.SetSize(A->Width());
-   A->Mult(x, y);
-   if (b) { y.Add(beta,*b); }
-}
-
-} // namespace mfem

linalg/functional.hpp

@@ -1,744 +0,0 @@
-#ifndef MFEM_FUNCTIONAL_HPP
-#define MFEM_FUNCTIONAL_HPP
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_MPI
-#include "../general/communication.hpp"
-#endif
-
-#include "operator.hpp"
-#include "blockvector.hpp"
-#include "solvers.hpp"
-#include <cxxabi.h>
-#include <vector>
-
-namespace mfem
-{
-/// @brief A base class for functionals F:R^n->R
-///
-/// This class provides an interface for evaluating
-/// $ F:R^n->R, \nabla F:R^n->R^n $, and $ \nabla^2 F:R^n x R^n ->R^n $.
-/// F.Mult(x, y) evaluates the functional at a point x, and stores the result in y[0]
-/// F.GetGradient() returns an operator that evaluates the gradient
-/// F.GetGradient().GetGradient(x) returns an Hessian action operator.
-///
-/// The usual Operator::GetGradient(const Vector &x) method for this method
-/// is deprecated as $ \nabla F $ only takes a single argument x.
-/// It is redundant to use F.GetGradient(x).Mult(x, y) to evaluate the gradient.
-/// Instead, use F.GetGradient().Mult(x, y) to evaluate the gradient at x.
-///
-/// The gradient and Hessian can be defined in two ways:
-/// 1. If the gradient is available as a sperate operator,
-///    then override Functional::GetGradient().
-///    In this case, Functional::HessianMult() will not be called.
-///
-/// 2. Otherwise, override Functional::EvalGradient() and Functional::HessianMult()
-///    to evaluate the gradient and Hessian action, respectively.
-///    The helper classes, GradientOperator and HessianActionOperator,
-///    will call these methods to evaluate the gradient and Hessian action.
-///    If Hessian is a seperate operator, then you can override
-///    The GradientOperator::GetGradient(x) will call Functional::GetHessian(x)
-///
-class Functional : public Operator
-{
-   Operator * riesz_map = nullptr; ///< Riesz map operator, if available
-public:
-   /// @brief Create a Functional with optional gradient and hessian
-   /// @param n number of variables
-   Functional(int n=0)
-      : Operator(1, n)
-      , grad_operator(*this)
-      , hessian_action_operator(*this)
-   { }
-
-#ifdef MFEM_USE_MPI
-   Functional(MPI_Comm comm, int n=0)
-      : Functional(n)
-   { SetComm(comm); }
-
-   void SetComm(MPI_Comm comm_)
-   {
-      parallel = comm_ != MPI_COMM_NULL;
-      comm = comm_;
-   }
-   MPI_Comm GetComm() const { return comm; }
-   bool IsParallel() const { return parallel; }
-#else
-   constexpr bool IsParallel() const { return false; }
-#endif
-
-   void SetRieszMap(Operator &op) { riesz_map = &op; }
-   /// @brief return the GradientOperator that evaluates the gradient
-   ///        input x is not used. Use GetGradient().Mult(x,y) to evaluate the gradient
-   ///        we recommend using GetGradient() instead of GetGradient(x)
-   /// Deprecated. See Functional::GetGradient()
-   MFEM_DEPRECATED
-   Operator &GetGradient(const Vector &dummy) const override final { return GetGradient(); }
-
-   /// @brief Return the GradientOperator that wraps Functional::EvalGradient() for Mult().
-   /// @note If the functional has a corresponding standalone gradient operator,
-   /// override this method to return the gradient operator.
-   virtual Operator &GetGradient() const { return grad_operator; }
-
-   /// @brief Evaluate the functional at a point x that will be called by GradientOperator::Mult()
-   /// @note This method is not meant to be called directly. See, GradientOperator
-   virtual void EvalGradient(const Vector &x, Vector &y) const
-   {
-      MFEM_ABORT("Functional::EvalGradient() not implemented");
-   }
-   /// @brief Evaluate the Hessian action at a point x and direction d
-   /// that will be called by Functional::GetGradient().GetHessian(x).Mult(d,y)
-   /// @note This method is not meant to be called directly. See, HessianActionOperator
-   virtual void HessianMult(const Vector &x, const Vector &d, Vector &y) const
-   {
-      MFEM_ABORT("Functional::HessianMult() not implemented.");
-   }
-   /// @brief Return the HessianActionOperator at evaluation point x
-   ///        that wraps Functional::HessianMult() for Mult().
-   ///        See, HessianActionOperator and Functional::HessianMult().
-   ///
-   /// @note If the Hessian is available as a seperate operator, override this method.
-   ///
-   /// @warning If GetGradient() is overridden, this method will not be used.
-   virtual Operator &GetHessian(const Vector &x) const
-   {
-      hessian_action_operator.SetX(x);
-      return hessian_action_operator;
-   }
-private:
-#ifdef MFEM_USE_MPI
-   bool parallel=false;
-#else
-   const static bool parallel=false;
-#endif
-#ifdef MFEM_USE_MPI
-   MPI_Comm comm;
-#endif
-   /// @brief A helper class to return an operator that evaluates the gradient
-   ///        using Functional::EvalGradient() method.
-   class GradientOperator : public Operator
-   {
-   private: const Functional &f; mutable Vector der;
-   public:
-      GradientOperator(const Functional &f) : Operator(f.Width()), f(f) {}
-      /// @brief Evaluate the gradient of Functional at a point x
-      void Mult(const Vector &x, Vector &y) const override final
-      {
-         if (f.riesz_map)
-         {
-            der.SetSize(f.Width());
-            f.EvalGradient(x, der);
-            f.riesz_map->Mult(der, y);
-         }
-         else
-         {
-            f.EvalGradient(x, y);
-         }
-      }
-      /// @brief Evaluate the Hessian of Functional at a point x
-      Operator &GetGradient(const Vector &x) const override final { return f.GetHessian(x); }
-   };
-   friend class GradientOperator;
-
-   /// @brief A helper class to return an operator that applies the Hessian action
-   ///        using Functional::HessianMult() method.
-   class HessianActionOperator : public Operator
-   {
-   private:
-      const Functional &f;
-      const Vector *x;
-   public:
-      HessianActionOperator(const Functional &f) : Operator(f.Width()), f(f) {}
-      void SetX(const Vector &new_x) { x = &new_x; }
-      void Mult(const Vector &d, Vector &y) const override { f.HessianMult(*x, d, y); }
-   };
-   friend class HessianActionOperator;
-
-   mutable GradientOperator grad_operator;
-   mutable HessianActionOperator hessian_action_operator;
-};
-
-/// @brief Stacked functioanl operator, [f1, ..., fk] where fi:R^n->R are functionals
-/*
-   Typical usage of this class is to provide a single operator for multiple constraints.
-   For example, consider a minimization problem with k constraints,
-   min f0(u) s.t. fi(u)=0, i=1,...,k.
-   The Lagrangian functional is
-   $ L(u, lambda) = F0(u) + sum_i lambda_i * fi(u) $
-   The first-order optimality conditions are
-   $ \nabla f0(u) + \sum_i lambda_i * grad fi(u) = 0 $
-   $ fi(u) = 0 $
-   where lambda_i are the Lagrange multipliers.
-   The StackedFunctional class can be used to represent the list of constraints fi(u).
-
-   StackedFunctional::Mult(u, y) will evaluate each functional y[i]=fi(u)
-
-   StackedFunctional::GetGradient(u) represents an operator, column-stacked gradient
-   That is, [grad f0(u), ..., grad fk(u)] in R^{n x k}
-   If you want to extract the gradient as a matrix, use
-   StackedFunctional::GetGradientMatrix(const Vector &x, DenseMatrix &grad)
-   As functionals are not assumed to return a sparse vector, the gradient is dense.
-
-   StackedFunctional::GetGradient(u).Mult(lambda, y) contract the gradients with the Lagrange multipliers
-   y = sum lambda_i * grad fi(u)
-   StackedFunctional::GetGradient(u).MultTranspose(d, y) return the directional derivative for each k
-   y[i] = <grad fi(u), d>
-
-   StackedFunctional::GetHessian(u, lambda).Mult(d, y) will return the contracted Hessian action
-   $ y = \sum_i \lambda_i * H_{fi}(u, d) $
-*/
-/// @warning Functionals should be all serial or all parallel.
-///
-class StackedFunctional : public Operator
-{
-public:
-   StackedFunctional(int n=0)
-      : Operator(0, n)
-      , funcs(0)
-      , grad_helper_op(*this)
-      , hessian_helper_op(*this)
-   {}
-   StackedFunctional(Functional &f)
-      : Operator(0, f.Width())
-      , grad_helper_op(*this)
-      , hessian_helper_op(*this)
-   { AddFunctional(f); }
-
-   StackedFunctional(const std::vector<Functional*> &funcs)
-      : Operator((int)funcs.size(), funcs[0]->Width())
-      , grad_helper_op(*this)
-      , hessian_helper_op(*this)
-   { for (auto &f : funcs) { AddFunctional(*f); } }
-
-   void AddFunctional(Functional &f)
-   {
-#ifdef MFEM_USE_MPI
-      if (funcs.empty()) { if (f.IsParallel()) { SetComm(f.GetComm()); } }
-#endif
-      MFEM_VERIFY(f.Width() == Width(),
-                  "StackedFunctional::AddFunctional: Functional width does not match with the operator.");
-      MFEM_VERIFY(parallel == f.IsParallel(),
-                  "StackedFunctional::AddFunctional: Parallelism mismatch.");
-      funcs.push_back(&f);
-      height++;
-   }
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      y.SetSize(Height());
-      Vector yview;
-      for (int i=0; i<Height(); i++)
-      {
-         yview.MakeRef(y, i, 1);
-         funcs[i]->Mult(x, yview);
-      }
-   }
-
-   Operator &GetGradient(const Vector &x) const override
-   {
-      grad_helper_op.SetX(x);
-      return grad_helper_op;
-   }
-
-   void GetGradientMatrix(const Vector &x, DenseMatrix &grads) const
-   {
-      grads.SetSize(Width(), Height());
-      Vector grad;
-      for (int i=0; i<Height(); i++)
-      {
-         grads.GetColumnReference(i, grad);
-         funcs[i]->GetGradient().Mult(x, grad);
-      }
-   }
-   Functional &GetFunctional(int i) const
-   {
-      MFEM_VERIFY(i >= 0 && i < Height(),
-                  "StackedFunctional::GetFunctional: Index out of bounds.");
-      return *funcs[i];
-   }
-
-   Operator &GetHessian(const Vector &x, const Vector &lambda) const
-   {
-      hessian_helper_op.SetX(x, lambda);
-      return hessian_helper_op;
-   }
-
-   bool parallel;
-   bool IsParallel() const { return parallel; }
-#ifdef MFEM_USE_MPI
-   void SetComm(MPI_Comm comm_)
-   {
-      parallel = comm != MPI_COMM_NULL;
-      comm = comm_;
-   }
-   MPI_Comm GetComm() const { return comm; }
-#endif
-
-protected:
-#ifdef MFEM_USE_MPI
-   MPI_Comm comm;
-#endif
-   std::vector<Functional*> funcs;
-
-   class GradientOperator : public Operator
-   {
-   public:
-      GradientOperator(const StackedFunctional &op)
-         : Operator(op.Width(), op.Height())
-         , op(op)
-         , tmp_grad(op.Width())
-      {}
-      void SetX(const Vector &x) const { x_curr = &x; }
-      Operator &GetGradient(const Vector &lambda) const override
-      {
-         op.hessian_helper_op.SetX(*x_curr, lambda);
-         return op.hessian_helper_op;
-      }
-
-      void Mult(const Vector &lambda, Vector &y) const override
-      {
-         y.SetSize(op.Width());
-         y = 0.0;
-         for (int i=0; i<op.Height(); i++)
-         {
-            op.funcs[i]->GetGradient().Mult(*x_curr, tmp_grad);
-            y.Add(lambda[i], tmp_grad);
-         }
-      }
-      void MultTranspose(const Vector &x, Vector &y) const override
-      {
-         y.SetSize(op.Height());
-         for (int i=0; i<op.Height(); i++)
-         {
-            op.funcs[i]->GetGradient().Mult(x, tmp_grad);
-            y[i] = InnerProduct(tmp_grad, *x_curr);
-         }
-#ifdef MFEM_USE_MPI
-         if (op.IsParallel())
-         {
-            MPI_Allreduce(MPI_IN_PLACE, y.GetData(), op.Height(),
-                          MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                          op.GetComm());
-         }
-#endif
-      }
-
-   private:
-      const StackedFunctional &op;
-      mutable const Vector *x_curr;
-      mutable Vector tmp_grad;
-
-   };
-   class HessianActionOperator : public Operator
-   {
-   public:
-      HessianActionOperator(const StackedFunctional &op)
-         : Operator(op.Width()), op(op)
-      {}
-      void SetX(const Vector &x, const Vector &lambda) const { x_curr = &x; lambda_curr = &lambda; }
-
-      void Mult(const Vector &d, Vector &y) const override
-      {
-         y.SetSize(op.Width());
-         y = 0.0;
-         for (int i=0; i<op.Height(); i++)
-         {
-            op.funcs[i]->GetGradient().GetGradient(*x_curr).Mult(d, tmp_hessian);
-            y.Add((*lambda_curr)[i], tmp_hessian);
-         }
-      }
-   private:
-      const StackedFunctional &op;
-      mutable Vector tmp_hessian;
-      mutable const Vector *x_curr;
-      mutable const Vector *lambda_curr;
-   };
-   friend class GradientOperator;
-   friend class HessianActionOperator;
-   mutable GradientOperator grad_helper_op;
-   mutable HessianActionOperator hessian_helper_op;
-private:
-};
-
-class ConstrainedOptimizationProblem : public Functional
-{
-public:
-   ConstrainedOptimizationProblem(Functional &objective_,
-                                  Operator *eq_constraints_=nullptr,
-                                  Operator *ineq_constraints_=nullptr)
-      : Functional(objective_.Width())
-      , objective(objective_)
-      , eq_constraints(eq_constraints_)
-      , ineq_constraints(ineq_constraints_)
-   {
-      // Check Size
-      MFEM_VERIFY((eq_constraints == nullptr ||
-                   eq_constraints->Width() == objective.Width()),
-                  "ConstrainedFunctional: Equality constraints width does not match with the objective.");
-      MFEM_VERIFY((ineq_constraints == nullptr ||
-                   ineq_constraints->Width() == objective.Width()),
-                  "ConstrainedFunctional: Inequality constraints width does not match with the objective.");
-#ifdef MFEM_USE_MPI
-      if (objective.IsParallel()) { SetComm(objective.GetComm()); }
-#endif
-   }
-
-   Functional &GetObjective() { return objective; }
-   const Functional &GetObjective() const { return objective; }
-   Operator *GetEqualityConstraints() { return eq_constraints; }
-   const Operator *GetEqualityConstraints() const { return eq_constraints; }
-   Operator *GetInequalityConstraints() { return ineq_constraints; }
-   const Operator *GetInequalityConstraints() const { return ineq_constraints; }
-protected:
-   Functional &objective;
-   Operator *eq_constraints;
-   Operator *ineq_constraints;
-};
-
-/// @brief A Lagrangian functional for
-/// min F(u)
-/// subject to C(u) = 0
-/// That is, L(u, lambda) = F(u) + <lambda, C(u)>
-///
-/// We assume that $ F:R^n -> R $ is a functional,
-///                $ C:R^n -> R^k $ is an equality constraint operator,
-/// C should return a residual. That is,
-/// C(u) = c, then C.Mult(u, y) should return y[i] = C_i(u) - c_i.
-///
-/// C.GetGradient(u):R^k -> R^n that takes lambda and returns the contracted gradient at x
-/// C.GetGradient(u).Mult(lambda, y) returns y = sum lambda_i * grad C_i(u)
-///
-/// C's gradient should support MultTranspose method
-/// That is, C.GetGradient(u).MultTranspose(d, y) returns y[i] = <grad C_i(u), d>
-///
-class LagrangianFunctional : public ConstrainedOptimizationProblem
-{
-private:
-   mutable Vector eq_residual;
-public:
-   LagrangianFunctional(Functional &objective,
-                        Operator &eq_constraints)
-      : ConstrainedOptimizationProblem(objective, &eq_constraints)
-      , eq_residual(eq_constraints.Height())
-   {
-      width = objective.Width() + eq_constraints.Height();
-
-      offsets.SetSize(3);
-      offsets[0] = 0;
-      offsets[1] = objective.Width();
-      offsets[2] = eq_constraints.Height();
-      offsets.PartialSum();
-   }
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      const BlockVector input_block(const_cast<Vector&>(x), offsets);
-      const Vector &u = input_block.GetBlock(0);
-      const Vector &lambda = input_block.GetBlock(1);
-
-      y.SetSize(1);
-      y[0] = 0.0;
-      objective.Mult(u, y);
-      eq_constraints->Mult(u, eq_residual);
-      y[0] += InnerProduct(lambda, eq_residual);
-   }
-
-   void EvalGradient(const Vector &x, Vector &y) const override
-   {
-      const BlockVector input_block(const_cast<Vector&>(x), offsets);
-      const Vector &u = input_block.GetBlock(0);
-      const Vector &lambda = input_block.GetBlock(1);
-
-      y.SetSize(Width());
-      BlockVector output_block(y, offsets);
-      Vector &opt_residual = output_block.GetBlock(0);
-      eq_residual = output_block.GetBlock(1);
-      y = 0.0;
-      // grad F(u) + \sum_i lambda_i grad C_i(u)
-      objective.GetGradient().Mult(u, opt_residual);
-      eq_constraints->GetGradient(u).AddMult(lambda, opt_residual);
-      // grad C_i(u)^T
-      eq_constraints->GetGradient(u).MultTranspose(u, eq_residual);
-   }
-
-   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
-   /// and direction d=[v, mu]
-   /// $ [H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
-   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
-   {
-      const BlockVector input_block(const_cast<Vector&>(x), offsets);
-      const Vector &u = input_block.GetBlock(0);
-      const Vector &lambda = input_block.GetBlock(1);
-
-      const BlockVector direction_block(const_cast<Vector&>(x), offsets);
-      const Vector &v = direction_block.GetBlock(0);
-      const Vector &mu = direction_block.GetBlock(1);
-
-      y.SetSize(Width());
-      BlockVector output_block(y, offsets);
-      Vector &opt_H = output_block.GetBlock(0); // Optimality Hessian
-      Vector &eq_H = output_block.GetBlock(1); // Equality Hessian
-
-      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u)
-      objective.GetGradient().GetGradient(u).Mult(v, opt_H);
-      eq_constraints->GetGradient(u).GetGradient(lambda).AddMult(v, opt_H);
-      eq_constraints->GetGradient(u).Mult(mu, eq_H);
-      // <grad C_i(u), d>
-      eq_constraints->GetGradient(u).MultTranspose(d, eq_H);
-   }
-
-protected:
-   Array<int> offsets; // offsets for [x, lambda, mu]
-};
-
-/// @brief An augmented Lagrangian functional of the form
-/// F(u) + 0.5 mu * ||C(u)||^2 + <lambda, C(u)>
-/// where F is the objective functional,
-/// C is the equality constraint operator,
-/// lambda is the Lagrange multiplier vector (initialized to zero),
-/// mu is the penalty parameter (defaults to 1.0)
-///
-/// Currently, only equality constraints are supported.
-///
-/// AugLagrangianFunctional::Update() will update the penalty and Lagrange multiplier vectors
-/// By default, lambda <- lambda + mu * C(u)
-///             mu <- mu (no update)
-class AugLagrangianFunctional : public ConstrainedOptimizationProblem
-{
-public:
-   AugLagrangianFunctional(Functional &objective_,
-                           Operator &eq_constraints_)
-      : ConstrainedOptimizationProblem(objective_, &eq_constraints_)
-      , lambda(eq_constraints_.Height())
-      , mu(1.0)
-      , eq_residual(eq_constraints_.Height())
-      , eq_dir(eq_constraints_.Height())
-   {
-      lambda = 0.0;
-   }
-
-   void SetLambda(const Vector &lambda_)
-   {
-      MFEM_VERIFY(lambda_.Size() == eq_constraints->Height(),
-                  "AugLagrangianFunctional: Lambda size does not match with the equality constraints.");
-      lambda = lambda_;
-   }
-
-   void SetPenalty(real_t mu_)
-   {
-      MFEM_VERIFY(mu_ >= 0.0,
-                  "AugLagrangianFunctional: Penalty parameter mu must be non-negative.");
-      mu = mu_;
-   }
-
-   virtual void Update(const Vector &x)
-   {
-      // Update the Lagrange multipliers
-      eq_constraints->AddMult(x, lambda, mu);
-      // Update the penalty parameter
-      // Do nothing
-   }
-   const Vector &GetLambda() const { return lambda; }
-   real_t GetPenalty() const { return mu; }
-
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      y.SetSize(1);
-      objective.Mult(x, y);
-      eq_constraints->Mult(x, eq_residual);
-      y[0] += lambda*eq_residual;
-      y[0] += 0.5 * mu * (eq_residual*eq_residual);
-   }
-
-   void EvalGradient(const Vector &x, Vector &y) const override
-   {
-      y.SetSize(Width());
-      // grad F(x) + \sum_i (lambda_i + mu * C_i(x)) grad C_i(x)
-      Vector curr_lambda = lambda; // store lambda + mu * C(x)
-      objective.GetGradient().Mult(x, y);
-      eq_constraints->Mult(x, eq_residual);
-      curr_lambda.Add(mu, eq_residual);
-      eq_constraints->GetGradient(x).AddMult(curr_lambda, y);
-   }
-
-   /// @brief Evaluate the Hessian action at a point x=[u, lambda]
-   /// and direction d=[v, mu]
-   /// $ H_F(u,d) + \sum_i \lambda_i H_{C_i}(u, d) + <\nabla C_i(u), mu> $
-   void HessianMult(const Vector &x, const Vector &d, Vector &y) const override
-   {
-      // H_F(u,d) + \sum_i lambda_i H_{C_i}(u, d) + mu_i grad C_i(u) <grad C_i(u), d>
-      objective.GetGradient().GetGradient(x).Mult(d, y);
-
-      Vector curr_lambda = lambda;
-      eq_constraints->Mult(x, eq_residual);
-      curr_lambda.Add(mu, eq_residual);
-
-      eq_constraints->GetGradient(x).GetGradient(curr_lambda).AddMult(d, y);
-      // eq_dir = <grad C_i(u), d>
-      eq_constraints->GetGradient(x).MultTranspose(d, eq_dir);
-      // mu_i <grad C_i(u), eq_dir>
-      eq_constraints->GetGradient(x).AddMult(eq_dir, y, mu);
-   }
-
-protected:
-   Vector lambda;
-   real_t mu;
-   mutable Vector eq_residual; // residual of the equality constraints, R^k
-   // directional derivative of the equality constraints, R^k
-   mutable Vector eq_dir;
-};
-
-/// @brief Quadratic functional of the form
-///        f(u) = 0.5 * <A u, u> + beta<b, u> + c
-///        where A is a square (possibly nonlinear) operator,
-///        beta is a scalar (defaults to 1.0, not used when b is nullptr),
-///        b is a vector (independent of u, optional),
-///        c is a constant (independent of u, optional).
-///        GetHessian() returns the operator A.
-///
-class QuadraticFunctional : public Functional
-{
-public:
-   QuadraticFunctional()
-      : Functional(0)
-      , A(nullptr), b(nullptr), c(0.0)
-   {}
-   QuadraticFunctional(const Operator *A_, const Vector *b_=nullptr,
-                       const real_t beta_=1.0,
-                       const real_t c_=0.0);
-#ifdef MFEM_USE_MPI
-   QuadraticFunctional(MPI_Comm comm_)
-      : QuadraticFunctional()
-   { SetComm(comm_); }
-
-   QuadraticFunctional(MPI_Comm comm_, const Operator *A_,
-                       const Vector *b_=nullptr,
-                       const real_t beta_=1.0, const real_t c_=0.0)
-      : QuadraticFunctional(A_, b_, beta_, c_)
-   { SetComm(comm_); }
-#endif
-
-   void SetOperator(const Operator &A_);
-   void SetVector(const Vector &b_, const real_t beta_=1.0);
-   void SetConstant(real_t c_);
-
-   void Mult(const Vector &x, Vector &y) const override;
-   void EvalGradient(const Vector &x, Vector &y) const override;
-protected:
-   const Operator *A;
-   real_t beta;
-   const Vector *b;
-   real_t c;
-   mutable Vector aux;
-protected:
-   /// @brief return the underlying Operator A
-   /// @warning Modifying the returned operator leads to undefined behavior.
-   Operator& GetHessian(const Vector &dummy) const override
-   {
-      return const_cast<Operator&>(*A);
-   }
-};
-
-class Optimizer : public IterativeSolver
-{
-public:
-   Optimizer() : IterativeSolver(), f(nullptr) { }
-#ifdef MFEM_USE_MPI
-   Optimizer(MPI_Comm comm) : IterativeSolver(comm), f(nullptr) { }
-#endif
-   // @brief Set the subproblem functional operator
-   // @param op the functional operator
-   // @note The functional will be stored in subproblem, and oper will be set to the gradient of the functional.
-   void SetOperator(const Functional &f_)
-   {
-      f = &f_;
-      IterativeSolver::SetOperator(f_.GetGradient());
-   }
-   virtual void SetLinearSolver(Solver &prec) { IterativeSolver::SetPreconditioner(prec); }
-
-   /// @brief This will abort. Should be called only with a Functional operator.
-   void SetOperator(const Operator &op) override
-   {
-      MFEM_ABORT("OptSolver::SetOperator() should not be called directly. Use SetFunctional() instead.");
-   }
-protected:
-   const Functional * f;
-};
-
-class NewtonOptimizer : public Optimizer
-{
-private:
-   real_t step_size = 1.0; // default step size
-public:
-   NewtonOptimizer() : Optimizer() { }
-#ifdef MFEM_USE_MPI
-   NewtonOptimizer(MPI_Comm comm) : Optimizer(comm) { }
-#endif
-   void SetStepSize(real_t step_size_) { step_size = step_size_; }
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      dx.SetSize(x.Size());
-      y.SetSize(x.Size());
-      y = x;
-      MFEM_ASSERT(f != nullptr,
-                  "NewtonOptimizer::Mult() called without a functional operator.");
-      MFEM_ASSERT(prec != nullptr,
-                  "NewtonOptimizer::Mult() called without a linear solver.");
-      for (int i=0; i<max_iter; i++)
-      {
-         oper->Mult(y, grad);
-         Operator &hess = oper->GetGradient(y);
-         prec->SetOperator(hess);
-         prec->Mult(grad, dx);
-         y.Add(-step_size, dx);
-         if (Dot(dx, dx) < abs_tol*abs_tol)
-         {
-            break;
-         }
-      }
-   }
-private:
-   mutable Vector grad;
-   mutable Vector dx;
-};
-
-class GradientDescentOptimizer : public Optimizer
-{
-private:
-   real_t step_size = 1.0; // default step size
-public:
-   GradientDescentOptimizer() : Optimizer() { }
-#ifdef MFEM_USE_MPI
-   GradientDescentOptimizer(MPI_Comm comm) : Optimizer(comm) { }
-#endif
-   void SetStepSize(real_t step_size_) { step_size = step_size_; }
-
-   void Mult(const Vector &x, Vector &y) const override
-   {
-      grad.SetSize(x.Size());
-      y.SetSize(x.Size());
-      y = x;
-      MFEM_ASSERT(f != nullptr,
-                  "NewtonOptimizer::Mult() called without a functional operator.");
-      MFEM_ASSERT(prec != nullptr,
-                  "NewtonOptimizer::Mult() called without a linear solver.");
-      for (int i=0; i<max_iter; i++)
-      {
-         oper->Mult(y, grad);
-         y.Add(-step_size, grad);
-         if (Dot(grad, grad) < abs_tol*abs_tol)
-         {
-            break;
-         }
-      }
-   }
-private:
-   mutable Vector grad;
-};
-
-
-} // namespace mfem
-#endif // MFEM_FUNCTIONAL_HPP

tests/unit/CMakeLists.txt

@@ -20,6 +20,7 @@ set(UNIT_TESTS_SRCS
   dfem/test_diffusion.cpp
   dfem/test_divergence.cpp
   dfem/test_mass.cpp
+  dfem/test_lvector_interface.cpp
   general/test_array.cpp
   general/test_reduction.cpp
   general/test_arrays_by_name.cpp
@@ -144,11 +145,11 @@ set(UNIT_TESTS_SRCS
 # SERIAL CPU TESTS: unit_tests
 #-----------------------------------------------------------
 if (MFEM_USE_CUDA)
-   set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                 PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-  set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
+    set_property(SOURCE unit_test_main.cpp ${UNIT_TESTS_SRCS}
                PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 
@@ -175,7 +176,7 @@ COMMAND ${CMAKE_COMMAND} -E copy_directory
 #   make unit_tests
 #   ctest -R unit_tests [-V]
 if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-  add_test(NAME unit_tests COMMAND unit_tests)
+    add_test(NAME unit_tests COMMAND unit_tests)
 endif()
 
 #-----------------------------------------------------------
@@ -183,16 +184,16 @@ endif()
 #-----------------------------------------------------------
 # Create CUDA executable and test
 if (MFEM_USE_CUDA)
-   # gpu_unit_tests
-   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-   set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-   target_link_libraries(gpu_unit_tests mfem)
-   add_dependencies(gpu_unit_tests copy_data)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-   endif()
+    # gpu_unit_tests
+    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+    set_property(SOURCE ${GPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+    target_link_libraries(gpu_unit_tests mfem)
+    add_dependencies(gpu_unit_tests copy_data)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -200,15 +201,15 @@ endif()
 #-----------------------------------------------------------
 # Create HIP 'gpu_unit_tests' executable and test
 if (MFEM_USE_HIP)
-   # gpu_unit_tests
-   set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
-   mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-   target_link_libraries(gpu_unit_tests mfem)
-   add_dependencies(gpu_unit_tests copy_data)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
-   if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
-   endif()
+    # gpu_unit_tests
+    set(GPU_UNIT_TESTS_SRCS gpu_unit_test_main.cpp)
+    mfem_add_executable(gpu_unit_tests ${GPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+    target_link_libraries(gpu_unit_tests mfem)
+    add_dependencies(gpu_unit_tests copy_data)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} gpu_unit_tests)
+    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+        add_test(NAME gpu_unit_tests COMMAND gpu_unit_tests)
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -224,7 +225,7 @@ function(add_serial_miniapp_test name test_uvm)
 
     set(${NAME}_TESTS_SRCS miniapps/test_${name}.cpp)
     if (MFEM_USE_CUDA)
-       set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+        set_property(SOURCE ${${NAME}_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
     endif(MFEM_USE_CUDA)
 
     mfem_add_executable(${name}_tests_cpu ${${NAME}_TESTS_SRCS})
@@ -244,25 +245,25 @@ function(add_serial_miniapp_test name test_uvm)
     endif()
 
     if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-       mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
-       target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
-       target_link_libraries(${name}_tests_gpu mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
-       if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-           add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
-       endif()
+        mfem_add_executable(${name}_tests_gpu ${${NAME}_TESTS_SRCS})
+        target_compile_definitions(${name}_tests_gpu PUBLIC MFEM_${NAME}_DEVICE="gpu")
+        target_link_libraries(${name}_tests_gpu mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ${name}_tests_gpu)
+        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+            add_test(NAME ${name}_tests_gpu COMMAND ${name}_tests_gpu)
+        endif()
 
-       if (test_uvm)
-           mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
-           target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
+        if (test_uvm)
+            mfem_add_executable(${name}_tests_gpu_uvm ${${NAME}_TESTS_SRCS})
+            target_compile_definitions(${name}_tests_gpu_uvm PUBLIC
                MFEM_${NAME}_DEVICE="gpu:uvm")
-           target_link_libraries(${name}_tests_gpu_uvm mfem)
-           add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
+            target_link_libraries(${name}_tests_gpu_uvm mfem)
+            add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME}
                ${name}_tests_gpu_uvm)
-           if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-               add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
-           endif()
-       endif()
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME ${name}_tests_gpu_uvm COMMAND ${name}_tests_gpu_uvm)
+            endif()
+        endif()
     endif()
 endfunction(add_serial_miniapp_test)
 
@@ -277,25 +278,25 @@ add_dependencies(tmop_pa_tests_cpu copy_miniapps_meshing_data)
 #-----------------------------------------------------------
 # Add 'ceed_tests' executable and test; add extra tests 'ceed_test_*'
 if (MFEM_USE_CEED)
-   set(CEED_TESTS_SRCS
+    set(CEED_TESTS_SRCS
       ceed/test_ceed.cpp
       ceed/test_ceed_main.cpp)
-   if (MFEM_USE_CUDA)
-      set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-   endif(MFEM_USE_CUDA)
-   mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
-   target_link_libraries(ceed_tests mfem)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
-   # Add CEED tests
-   add_test(NAME ceed_tests COMMAND ceed_tests)
-   if (MFEM_USE_CUDA)
-      add_test(NAME ceed_tests_cuda_ref
+    if (MFEM_USE_CUDA)
+        set_property(SOURCE ${CEED_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+    endif(MFEM_USE_CUDA)
+    mfem_add_executable(ceed_tests ${CEED_TESTS_SRCS})
+    target_link_libraries(ceed_tests mfem)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} ceed_tests)
+    # Add CEED tests
+    add_test(NAME ceed_tests COMMAND ceed_tests)
+    if (MFEM_USE_CUDA)
+        add_test(NAME ceed_tests_cuda_ref
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/ref)
-      add_test(NAME ceed_tests_cuda_shared
+        add_test(NAME ceed_tests_cuda_shared
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/shared)
-      add_test(NAME ceed_tests_cuda_gen
+        add_test(NAME ceed_tests_cuda_gen
                COMMAND ceed_tests --device ceed-cuda:/gpu/cuda/gen)
-   endif()
+    endif()
 endif()
 
 #-----------------------------------------------------------
@@ -303,54 +304,54 @@ endif()
 #-----------------------------------------------------------
 # Define executables and tests
 if (MFEM_USE_MPI)
-   # punit_tests
-   if (MFEM_USE_CUDA)
-      set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
-   endif()
-   mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
-   target_link_libraries(punit_tests mfem)
-   add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
-   foreach(np 1 ${MFEM_MPI_NP})
-      if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-         add_test(NAME punit_tests_np=${np}
+    # punit_tests
+    if (MFEM_USE_CUDA)
+        set_property(SOURCE punit_test_main.cpp PROPERTY LANGUAGE CUDA)
+    endif()
+    mfem_add_executable(punit_tests punit_test_main.cpp ${UNIT_TESTS_SRCS})
+    target_link_libraries(punit_tests mfem)
+    add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} punit_tests)
+    foreach(np 1 ${MFEM_MPI_NP})
+        if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+            add_test(NAME punit_tests_np=${np}
                   COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                   ${MPIEXEC_PREFLAGS} $<TARGET_FILE:punit_tests>
                   ${MPIEXEC_POSTFLAGS})
-      endif()
-   endforeach()
-   if (MFEM_USE_CUDA)
-       # pgpu_unit_tests
-       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-       set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
-       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-       add_dependencies(pgpu_unit_tests copy_data)
-       target_link_libraries(pgpu_unit_tests mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-       foreach(np 1 ${MFEM_MPI_NP})
-          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-             add_test(NAME pgpu_unit_tests_np=${np}
+        endif()
+    endforeach()
+    if (MFEM_USE_CUDA)
+        # pgpu_unit_tests
+        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+        set_property(SOURCE ${PGPU_UNIT_TESTS_SRCS} PROPERTY LANGUAGE CUDA)
+        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+        add_dependencies(pgpu_unit_tests copy_data)
+        target_link_libraries(pgpu_unit_tests mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+        foreach(np 1 ${MFEM_MPI_NP})
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-          endif()
-       endforeach()
+            endif()
+        endforeach()
     endif()
     if (MFEM_USE_HIP)
-       # pgpu_unit_tests
-       set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
-       mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
-       add_dependencies(pgpu_unit_tests copy_data)
-       target_link_libraries(pgpu_unit_tests mfem)
-       add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
-       foreach(np 1 ${MFEM_MPI_NP})
-          if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-             add_test(NAME pgpu_unit_tests_np=${np}
+        # pgpu_unit_tests
+        set(PGPU_UNIT_TESTS_SRCS pgpu_unit_test_main.cpp)
+        mfem_add_executable(pgpu_unit_tests ${PGPU_UNIT_TESTS_SRCS} ${UNIT_TESTS_SRCS})
+        add_dependencies(pgpu_unit_tests copy_data)
+        target_link_libraries(pgpu_unit_tests mfem)
+        add_dependencies(${MFEM_ALL_TESTS_TARGET_NAME} pgpu_unit_tests)
+        foreach(np 1 ${MFEM_MPI_NP})
+            if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
+                add_test(NAME pgpu_unit_tests_np=${np}
                       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${np}
                       ${MPIEXEC_PREFLAGS} $<TARGET_FILE:pgpu_unit_tests>
                       ${MPIEXEC_POSTFLAGS})
-          endif()
-       endforeach()
-   endif()
+            endif()
+        endforeach()
+    endif()
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -424,8 +425,8 @@ endfunction(add_parallel_miniapp_test)
 
 # Additional MPI unit tests
 if (MFEM_USE_MPI)
-   add_parallel_miniapp_test(sedov TRUE)
-   add_parallel_miniapp_test(tmop_pa FALSE)
+    add_parallel_miniapp_test(sedov TRUE)
+    add_parallel_miniapp_test(tmop_pa FALSE)
 endif(MFEM_USE_MPI)
 
 #-----------------------------------------------------------
@@ -434,10 +435,10 @@ endif(MFEM_USE_MPI)
 #-----------------------------------------------------------
 set(DEBUG_DEVICE_SRCS miniapps/test_debug_device.cpp)
 if (MFEM_USE_CUDA)
-   set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
+    set_property(SOURCE ${DEBUG_DEVICE_SRCS} PROPERTY LANGUAGE CUDA)
 endif()
 if (MFEM_USE_HIP)
-   set_property(SOURCE ${DEBUG_DEVICE_SRCS}
+    set_property(SOURCE ${DEBUG_DEVICE_SRCS}
                 PROPERTY HIP_SOURCE_PROPERTY_FORMAT TRUE)
 endif()
 mfem_add_executable(debug_device_tests ${DEBUG_DEVICE_SRCS})

tests/unit/dfem/test_diffusion.cpp

@@ -255,8 +255,8 @@ void DFemDiffusion(const char *filename, int p, const int r)
          DOperator dop_mf(vsol, {{Coords, mfes}}, pmesh);
          const auto mf_vector_diffusion_qf =
             [] MFEM_HOST_DEVICE (const tensor<dscalar_t, DIM, DIM> &dudxi,
-               const tensor<real_t, DIM, DIM> &J,
-               const real_t &w)
+                                 const tensor<real_t, DIM, DIM> &J,
+                                 const real_t &w)
          {
             const auto invJ = inv(J), TinJ = transpose(invJ);
             return tuple{ (dudxi * invJ) * TinJ * det(J) * w };

tests/unit/dfem/test_lvector_interface.cpp

@@ -0,0 +1,98 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../unit_tests.hpp"
+#include "mfem.hpp"
+#include <utility>
+
+#ifdef MFEM_USE_MPI
+
+using namespace mfem;
+using namespace mfem::future;
+using mfem::future::tensor;
+
+constexpr int DIM = 3;
+
+namespace kernels
+{
+struct MFApply
+{
+   MFEM_HOST_DEVICE inline auto operator()(const tensor<real_t, DIM> &dudxi,
+                                           const tensor<real_t, DIM, DIM> &J,
+                                           const real_t &w) const
+   {
+      const auto invJ = inv(J);
+      return tuple{ (dudxi * invJ) * transpose(invJ) * det(J) * w };
+   }
+};
+}
+
+TEST_CASE("DFEM L-Vector interface", "[Parallel][DFEM]")
+{
+   constexpr int p = 2; // Polynomial order
+   constexpr int r = 1;
+   constexpr int q = 2 * p + r;
+
+   const auto filename = GENERATE("../../data/fichera.mesh");
+   Mesh smesh(filename);
+   ParMesh pmesh(MPI_COMM_WORLD, smesh);
+   MFEM_VERIFY(pmesh.Dimension() == DIM, "Mesh dimension mismatch");
+
+   pmesh.EnsureNodes();
+   auto *nodes = static_cast<ParGridFunction *>(pmesh.GetNodes());
+   smesh.Clear();
+
+   Array<int> all_domain_attr;
+   if (pmesh.attributes.Size() > 0)
+   {
+      all_domain_attr.SetSize(pmesh.attributes.Max());
+      all_domain_attr = 1;
+   }
+
+   H1_FECollection fec(p, DIM);
+   ParFiniteElementSpace pfes(&pmesh, &fec);
+   ParFiniteElementSpace *mfes = nodes->ParFESpace();
+
+   const auto *ir = &IntRules.Get(pmesh.GetTypicalElementGeometry(), q);
+
+   ParGridFunction x(&pfes), y(&pfes), z(&pfes);
+   Vector X(pfes.GetTrueVSize()), Y(pfes.GetTrueVSize()), Z(pfes.GetTrueVSize());
+
+   X.Randomize(1);
+   x.SetFromTrueDofs(X);
+
+   ParBilinearForm blf_fa(&pfes);
+   blf_fa.AddDomainIntegrator(new DiffusionIntegrator(ir));
+   blf_fa.Assemble();
+   blf_fa.Finalize();
+
+   static constexpr int U = 0, Coords = 1;
+
+   const auto solution = std::vector{FieldDescriptor{U, &pfes}};
+   DifferentiableOperator dop(solution, {{Coords, mfes}}, pmesh);
+
+   kernels::MFApply mf_apply_qf;
+   dop.AddDomainIntegrator(mf_apply_qf,
+                           tuple{Gradient<U>{}, Gradient<Coords>{}, Weight{}},
+                           tuple{Gradient<U>{}}, *ir, all_domain_attr);
+
+   // Use the L-vector interface to multiply
+   dop.SetMultLevel(DifferentiableOperator::MultLevel::LVECTOR);
+   dop.SetParameters({nodes});
+   dop.Mult(x, z);
+
+   blf_fa.Mult(x, y);
+
+   z -= y;
+   REQUIRE(z.Normlinf() == MFEM_Approx(0.0));
+}
+
+#endif