Updated explanation from PR

Action of HypreSmoother::MultTranspose

CHANGELOG

@@ -49,6 +49,12 @@ Version 4.2.1 (development)
   deprecated EvalSymmetric in MatrixCoefficient. Added DiagonalMatrixCoefficient
   for clarity, which is a typedef of VectorCoefficient.
 
+- Implemented an adaptive linear solver tolerance option for NewtonSolver based
+  on the algorithm of Eisenstat and Walker.
+
+- Extending support for L2 basis functions using MapTypes VALUE and INTEGRAL in
+  linear interpolators and GridFunction "GetValue" methods.
+
 
 Version 4.2, released on October 30, 2020
 =========================================

config/defaults.mk

@@ -161,7 +161,7 @@ endif
 ZLIB_DIR =
 ZLIB_OPT = $(if $(ZLIB_DIR),-I$(ZLIB_DIR)/include)
 ZLIB_LIB = $(if $(ZLIB_DIR),$(ZLIB_RPATH) -L$(ZLIB_DIR)/lib ,)-lz
-ZLIB_RPATH = -Wl,-rpath,$(ZLIB_DIR)/lib
+ZLIB_RPATH = $(XLINKER)-rpath,$(ZLIB_DIR)/lib
 
 LIBUNWIND_OPT = -g
 LIBUNWIND_LIB = $(if $(NOTMAC),-lunwind -ldl,)
@@ -231,19 +231,21 @@ MESQUITE_LIB = -L$(MESQUITE_DIR)/lib -lmesquite
 LIB_RT = $(if $(NOTMAC),-lrt,)
 SUITESPARSE_DIR = @MFEM_DIR@/../SuiteSparse
 SUITESPARSE_OPT = -I$(SUITESPARSE_DIR)/include
-SUITESPARSE_LIB = -Wl,-rpath,$(SUITESPARSE_DIR)/lib -L$(SUITESPARSE_DIR)/lib\
- -lklu -lbtf -lumfpack -lcholmod -lcolamd -lamd -lcamd -lccolamd\
- -lsuitesparseconfig $(LIB_RT) $(METIS_LIB) $(LAPACK_LIB)
+SUITESPARSE_LIB = $(XLINKER)-rpath,$(SUITESPARSE_DIR)/lib\
+ -L$(SUITESPARSE_DIR)/lib -lklu -lbtf -lumfpack -lcholmod -lcolamd -lamd -lcamd\
+ -lccolamd -lsuitesparseconfig $(LIB_RT) $(METIS_LIB) $(LAPACK_LIB)
 
 # SuperLU library configuration
 ifeq ($(MFEM_USE_SUPERLU5),YES)
    SUPERLU_DIR = @MFEM_DIR@/../SuperLU_DIST_5.1.0
    SUPERLU_OPT = -I$(SUPERLU_DIR)/include
-   SUPERLU_LIB = -Wl,-rpath,$(SUPERLU_DIR)/lib -L$(SUPERLU_DIR)/lib -lsuperlu_dist_5.1.0
+   SUPERLU_LIB = $(XLINKER)-rpath,$(SUPERLU_DIR)/lib -L$(SUPERLU_DIR)/lib\
+      -lsuperlu_dist_5.1.0
 else
    SUPERLU_DIR = @MFEM_DIR@/../SuperLU_DIST_6.3.1
    SUPERLU_OPT = -I$(SUPERLU_DIR)/include
-   SUPERLU_LIB = -Wl,-rpath,$(SUPERLU_DIR)/lib64 -L$(SUPERLU_DIR)/lib64 -lsuperlu_dist -lblas
+   SUPERLU_LIB = $(XLINKER)-rpath,$(SUPERLU_DIR)/lib64 -L$(SUPERLU_DIR)/lib64\
+      -lsuperlu_dist -lblas
 endif
 
 # SCOTCH library configuration (required by STRUMPACK <= v2.1.0, optional in
@@ -269,7 +271,7 @@ MPI_FORTRAN_LIB = -lmpifort
 # MUMPS library configuration
 MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.2.0
 MUMPS_OPT = -I$(MUMPS_DIR)/include
-MUMPS_LIB = -Wl,-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\
+MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\
  -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
 
 # STRUMPACK library configuration
@@ -299,8 +301,8 @@ GNUTLS_LIB = -lgnutls
 NETCDF_DIR = $(HOME)/local
 HDF5_DIR   = $(HOME)/local
 NETCDF_OPT = -I$(NETCDF_DIR)/include -I$(HDF5_DIR)/include $(ZLIB_OPT)
-NETCDF_LIB = -Wl,-rpath,$(NETCDF_DIR)/lib -L$(NETCDF_DIR)/lib\
- -Wl,-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib\
+NETCDF_LIB = $(XLINKER)-rpath,$(NETCDF_DIR)/lib -L$(NETCDF_DIR)/lib\
+ $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib\
  -lnetcdf -lhdf5_hl -lhdf5 $(ZLIB_LIB)
 
 # PETSc library configuration (version greater or equal to 3.8 or the dev branch)
@@ -312,9 +314,10 @@ PETSC_INC_VAR = PETSC_CC_INCLUDES
 PETSC_LIB_VAR = PETSC_EXTERNAL_LIB_BASIC
 ifeq ($(PETSC_FOUND),YES)
    PETSC_OPT := $(shell sed -n "s/$(PETSC_INC_VAR) = *//p" $(PETSC_VARS))
-   PETSC_LIB := $(shell sed -n "s/$(PETSC_LIB_VAR) = *//p" $(PETSC_VARS))
-   PETSC_LIB := -Wl,-rpath,$(abspath $(PETSC_DIR))/lib\
-      -L$(abspath $(PETSC_DIR))/lib -lpetsc $(PETSC_LIB)
+   PETSC_DEP := $(shell sed -n "s/$(PETSC_LIB_VAR) = *//p" $(PETSC_VARS))
+   PETSC_LIB = $(XLINKER)-rpath,$(abspath $(PETSC_DIR))/lib\
+      -L$(abspath $(PETSC_DIR))/lib -lpetsc\
+      $(subst $(CXX_XLINKER),$(XLINKER),$(PETSC_DEP))
 endif
 
 SLEPC_DIR := $(MFEM_DIR)/../slepc
@@ -326,9 +329,10 @@ ifeq ($(SLEPC_FOUND),YES)
    SLEPC_OPT := $(shell sed -n "s/$(SLEPC_INC_VAR) *= *//p" $(SLEPC_VARS))
    # Some additional external libraries might be defined in this file
    -include ${SLEPC_DIR}/${PETSC_ARCH}/lib/slepc/conf/slepcvariables
-   SLEPC_LIB := $(shell sed -n "s/$(SLEPC_LIB_VAR) *= *//p" $(SLEPC_VARS))
-   SLEPC_LIB := -Wl,-rpath,$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib\
-      -L$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib -lslepc $(SLEPC_LIB)
+   SLEPC_DEP := $(shell sed -n "s/$(SLEPC_LIB_VAR) *= *//p" $(SLEPC_VARS))
+   SLEPC_LIB = $(XLINKER)-rpath,$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib\
+      -L$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib -lslepc\
+      $(subst $(CXX_XLINKER),$(XLINKER),$(SLEPC_DEP))
 endif
 
 # MPFR library configuration
@@ -339,7 +343,7 @@ MPFR_LIB = -lmpfr
 CONDUIT_DIR = @MFEM_DIR@/../conduit
 CONDUIT_OPT = -I$(CONDUIT_DIR)/include/conduit
 CONDUIT_LIB = \
-   -Wl,-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
+   $(XLINKER)-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
    -lconduit -lconduit_relay -lconduit_blueprint  -ldl
 
 # Check if Conduit was built with hdf5 support, by looking
@@ -347,7 +351,7 @@ CONDUIT_LIB = \
 CONDUIT_HDF5_HEADER=$(CONDUIT_DIR)/include/conduit/conduit_relay_hdf5.hpp
 ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
    CONDUIT_OPT += -I$(HDF5_DIR)/include
-   CONDUIT_LIB += -Wl,-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
+   CONDUIT_LIB += $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
                   -lhdf5 $(ZLIB_LIB)
 endif
 
@@ -357,9 +361,9 @@ SIDRE_DIR = @MFEM_DIR@/../axom
 SIDRE_OPT = -I$(SIDRE_DIR)/include -I$(CONDUIT_DIR)/include/conduit\
  -I$(HDF5_DIR)/include
 SIDRE_LIB = \
-   -Wl,-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
-   -Wl,-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
-   -Wl,-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
+   $(XLINKER)-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
+   $(XLINKER)-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
+   $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
    -laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
 
 # PUMI
@@ -415,9 +419,9 @@ MKL_CPARDISO_DIR ?=
 MKL_MPI_WRAPPER ?= mkl_blacs_mpich_lp64
 MKL_LIBRARY_SUBDIR ?= lib
 MKL_CPARDISO_OPT = -I$(MKL_CPARDISO_DIR)/include
-MKL_CPARDISO_LIB = -Wl,-rpath,$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR)\
-                   -L$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR) -l$(MKL_MPI_WRAPPER)\
-                   -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
+MKL_CPARDISO_LIB = $(XLINKER)-rpath,$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR)\
+   -L$(MKL_CPARDISO_DIR)/$(MKL_LIBRARY_SUBDIR) -l$(MKL_MPI_WRAPPER)\
+   -lmkl_intel_lp64 -lmkl_sequential -lmkl_core
 
 # If YES, enable some informational messages
 VERBOSE = NO

examples/ex10p.cpp

@@ -178,6 +178,7 @@ int main(int argc, char *argv[])
    double visc = 1e-2;
    double mu = 0.25;
    double K = 5.0;
+   bool adaptive_lin_rtol = true;
    bool visualization = true;
    int vis_steps = 1;
 
@@ -206,6 +207,9 @@ int main(int argc, char *argv[])
                   "Shear modulus in the Neo-Hookean hyperelastic model.");
    args.AddOption(&K, "-K", "--bulk-modulus",
                   "Bulk modulus in the Neo-Hookean hyperelastic model.");
+   args.AddOption(&adaptive_lin_rtol, "-alrtol", "--adaptive-lin-rtol",
+                  "-no-alrtol", "--no-adaptive-lin-rtol",
+                  "Enable or disable adaptive linear solver rtol.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -574,6 +578,7 @@ HyperelasticOperator::HyperelasticOperator(ParFiniteElementSpace &f,
    newton_solver.SetPrintLevel(1); // print Newton iterations
    newton_solver.SetRelTol(rel_tol);
    newton_solver.SetAbsTol(0.0);
+   newton_solver.SetAdaptiveLinRtol(2, 0.5, 0.9);
    newton_solver.SetMaxIter(10);
 }
 

examples/ex16p.cpp

@@ -65,7 +65,7 @@ protected:
    HypreSmoother M_prec; // Preconditioner for the mass matrix M
 
    CGSolver T_solver;    // Implicit solver for T = M + dt K
-   HypreSmoother T_prec; // Preconditioner for the implicit solver
+   HypreBoomerAMG T_prec; // Preconditioner for the implicit solver
 
    double alpha, kappa;
 
@@ -353,7 +353,6 @@ int main(int argc, char *argv[])
          }
 #endif
       }
-      oper.SetParameters(u);
    }
 
 #ifdef MFEM_USE_ADIOS2
@@ -414,6 +413,7 @@ ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
    T_solver.SetPrintLevel(0);
    T_solver.SetPreconditioner(T_prec);
 
+   T_prec.SetPrintLevel(0);
    SetParameters(u);
 }
 
@@ -430,19 +430,53 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 void ConductionOperator::ImplicitSolve(const double dt,
                                        const Vector &u, Vector &du_dt)
 {
-   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt
-   if (!T)
-   {
-      T = Add(1.0, Mmat, dt, Kmat);
-      current_dt = dt;
+   // Here we use Picard iterations to solve a nonlinear equation
+   // for the Runge-Kutta stage vector k,
+   //
+   //    M*k = N(u+dt*k)                     (1)
+   //
+   // for nonlinear operator N. We assume N can be written as
+   //
+   //    N(u+dt*k) := L[u+dt*k](u+dt*k) + f(t)
+   //
+   // where L is a matrix-valued operator evaluated at u+dt*k and f(t)
+   // a (potentially zero) time-dependent forcing vector. (1) can be
+   // rewritten as a fixed-point equation
+   //
+   //    x = (M - dt*L[x])^{-1} (Mu + f)     (2)
+   //
+   // where x := u + dt*k, which can be solved using a Picard iteration,
+   // where a function G(x) = x is solved via iteraitons x_{k+1} = G(x_k).
+
+   double tol = 1e-6;
+   int maxiter = 100;
+
+   // Right-hand side for nonlinear iteration
+   Mmat.Mult(u, z);  // Add forcing function if one exists
+   du_dt = u;        // Set u as initial guess for x (2)
+   Vector temp(u);   // Vector to measure error
+   temp = u;
+   double error = 1;
+   int iter = 0;
+   while (error > tol) {
+      iter ++;
+      this->SetParameters(du_dt);   // Update nonlinear operator L[x]
+      T = Add(1.0, Mmat, dt, Kmat); // Form matrix (M - dt*L[x])
       T_solver.SetOperator(*T);
+      T_solver.Mult(z, du_dt);      // Apply (M - dt*L[x])^{-1}
+      temp -= du_dt;                // Measure error
+      error = std::sqrt(InnerProduct(MPI_COMM_WORLD, temp, temp));
+      temp = du_dt;
+      if (iter >= maxiter) {
+         mfem_warning("Nonlinear iteration did not converge!");
+         break;
+      }
    }
-   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Above we solved for x = u + dt*k, where k is the desired update
+   // Map du_dt -> k.
+   du_dt -= u;
+   du_dt /= dt;
 }
 
 void ConductionOperator::SetParameters(const Vector &u)
@@ -483,4 +517,4 @@ double InitialTemperature(const Vector &x)
    {
       return 1.0;
    }
-}
+}

examples/ex1p.cpp

@@ -231,10 +231,17 @@ int main(int argc, char *argv[])
    else
    {
       prec = new HypreBoomerAMG;
+      prec->SetOperator(*A);
    }
-   CGSolver cg(MPI_COMM_WORLD);
+   // CGSolver cg(MPI_COMM_WORLD);
+   AndersonAcceleration cg(MPI_COMM_WORLD);
+   cg.SetKDim(10);
+   cg.SetRestart(true);      // WORKS
+   cg.SetAAStart(0);          // WORKS
+   cg.SetWeight(1.0);         // Not robust but seems to work
+
    cg.SetRelTol(1e-12);
-   cg.SetMaxIter(2000);
+   cg.SetMaxIter(50);
    cg.SetPrintLevel(1);
    if (prec) { cg.SetPreconditioner(*prec); }
    cg.SetOperator(*A);

fem/bilininteg.cpp

@@ -3515,6 +3515,47 @@ VectorScalarProductInterpolator::AssembleElementMatrix2(
 }
 
 
+void
+ScalarCrossProductInterpolator::AssembleElementMatrix2(
+   const FiniteElement &dom_fe,
+   const FiniteElement &ran_fe,
+   ElementTransformation &Trans,
+   DenseMatrix &elmat)
+{
+   // Vector coefficient product with vector shape functions
+   struct VCrossVShapeCoefficient : public VectorCoefficient
+   {
+      VectorCoefficient &VQ;
+      const FiniteElement &fe;
+      DenseMatrix vshape;
+      Vector vc;
+
+      VCrossVShapeCoefficient(VectorCoefficient &vq, const FiniteElement &fe_)
+         : VectorCoefficient(fe_.GetDof()), VQ(vq), fe(fe_),
+           vshape(vdim, vq.GetVDim()), vc(vq.GetVDim()) { }
+
+      virtual void Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+      {
+         V.SetSize(vdim);
+         VQ.Eval(vc, T, ip);
+         fe.CalcPhysVShape(T, vshape);
+         for (int k = 0; k < vdim; k++)
+         {
+            V(k) = vc(0) * vshape(k,1) - vc(1) * vshape(k,0);
+         }
+      }
+   };
+
+   VCrossVShapeCoefficient dom_shape_coeff(*VQ, dom_fe);
+
+   elmat.SetSize(ran_fe.GetDof(),dom_fe.GetDof());
+
+   Vector elmat_as_vec(elmat.Data(), elmat.Height()*elmat.Width());
+
+   ran_fe.Project(dom_shape_coeff, Trans, elmat_as_vec);
+}
+
 void
 VectorCrossProductInterpolator::AssembleElementMatrix2(
    const FiniteElement &dom_fe,

fem/bilininteg.hpp

@@ -2671,11 +2671,28 @@ public:
 /** Integrator for the DG form:
     alpha < rho_u (u.n) {v},[w] > + beta < rho_u |u.n| [v],[w] >,
     where v and w are the trial and test variables, respectively, and rho/u are
-    given scalar/vector coefficients. The vector coefficient, u, is assumed to
-    be continuous across the faces and when given the scalar coefficient, rho,
-    is assumed to be discontinuous. The integrator uses the upwind value of rho,
-    rho_u, which is value from the side into which the vector coefficient, u,
-    points. */
+    given scalar/vector coefficients. {v} represents the average value of v on
+    the face and [v] is the jump such that {v}=(v1+v2)/2 and [v]=(v1-v2) for the
+    face between elements 1 and 2. For boundary elements, v2=0. The vector
+    coefficient, u, is assumed to be continuous across the faces and when given
+    the scalar coefficient, rho, is assumed to be discontinuous. The integrator
+    uses the upwind value of rho, rho_u, which is value from the side into which
+    the vector coefficient, u, points.
+
+    One use case for this integrator is to discretize the operator -u.grad(v)
+    with a DG formulation. The resulting formulation uses the
+    ConvectionIntegrator (with coefficient u, and parameter alpha = -1) and the
+    transpose of the DGTraceIntegrator (with coefficient u, and parameters
+    alpha = 1, beta = -1/2 to use the upwind face flux). This discretization and
+    the handling of the inflow and outflow boundaries is illustrated in Example
+    9/9p.
+
+    Another use case for this integrator is to discretize the operator -div(u v)
+    with a DG formulation. The resulting formulation is conservative and
+    consists of the transpose of the ConvectionIntegrator (with coefficient u,
+    and parameter alpha = 1) plus the DGTraceIntegrator (with coefficient u, and
+    parameters alpha = -1, beta = -1/2 to use the upwind face flux).
+    */
 class DGTraceIntegrator : public BilinearFormIntegrator
 {
 protected:
@@ -3090,6 +3107,22 @@ protected:
    VectorCoefficient *VQ;
 };
 
+/** Interpolator of the 2D cross product between a vector coefficient and an
+    H(curl)-conforming field onto an L2-conforming field. */
+class ScalarCrossProductInterpolator : public DiscreteInterpolator
+{
+public:
+   ScalarCrossProductInterpolator(VectorCoefficient & vc)
+      : VQ(&vc) { }
+
+   virtual void AssembleElementMatrix2(const FiniteElement &nd_fe,
+                                       const FiniteElement &l2_fe,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+protected:
+   VectorCoefficient *VQ;
+};
+
 /** Interpolator of the cross product between a vector coefficient and an
     H(curl)-conforming field onto an H(div)-conforming field. The range space
     can also be vector L2. */

fem/fe.cpp

@@ -596,17 +596,22 @@ void NodalFiniteElement::ProjectCurl_2D(
    const FiniteElement &fe, ElementTransformation &Trans,
    DenseMatrix &curl) const
 {
-   MFEM_ASSERT(GetMapType() == FiniteElement::INTEGRAL, "");
-
    DenseMatrix curl_shape(fe.GetDof(), 1);
 
    curl.SetSize(dof, fe.GetDof());
    for (int i = 0; i < dof; i++)
    {
       fe.CalcCurlShape(Nodes.IntPoint(i), curl_shape);
+
+      double w = 1.0;
+      if (GetMapType() == FiniteElement::VALUE)
+      {
+         Trans.SetIntPoint(&Nodes.IntPoint(i));
+         w /= Trans.Weight();
+      }
       for (int j = 0; j < fe.GetDof(); j++)
       {
-         curl(i,j) = curl_shape(j,0);
+         curl(i,j) = w * curl_shape(j,0);
       }
    }
 }
@@ -722,17 +727,34 @@ void NodalFiniteElement::Project(
 {
    if (fe.GetRangeType() == SCALAR)
    {
-      MFEM_ASSERT(map_type == fe.GetMapType(), "");
-
       Vector shape(fe.GetDof());
 
       I.SetSize(dof, fe.GetDof());
-      for (int k = 0; k < dof; k++)
+      if (map_type == fe.GetMapType())
       {
-         fe.CalcShape(Nodes.IntPoint(k), shape);
-         for (int j = 0; j < shape.Size(); j++)
+         for (int k = 0; k < dof; k++)
          {
-            I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
+            fe.CalcShape(Nodes.IntPoint(k), shape);
+            for (int j = 0; j < shape.Size(); j++)
+            {
+               I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
+            }
+         }
+      }
+      else
+      {
+         for (int k = 0; k < dof; k++)
+         {
+            Trans.SetIntPoint(&Nodes.IntPoint(k));
+            fe.CalcPhysShape(Trans, shape);
+            if (map_type == INTEGRAL)
+            {
+               shape *= Trans.Weight();
+            }
+            for (int j = 0; j < shape.Size(); j++)
+            {
+               I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
+            }
          }
       }
    }
@@ -1031,6 +1053,8 @@ void VectorFiniteElement::Project_RT(
 
          fe.CalcShape(ip, shape);
          Trans.SetIntPoint(&ip);
+         // Transform RT face normals from reference to physical space
+         // vk = adj(J)^T nk
          Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, vk);
          if (fe.GetMapType() == INTEGRAL)
          {
@@ -1048,6 +1072,8 @@ void VectorFiniteElement::Project_RT(
             {
                s = 0.0;
             }
+            // Project scalar basis function multiplied by each coordinate
+            // direction onto the transformed face normals
             for (int d = 0; d < sdim; d++)
             {
                I(k,j+d*shape.Size()) = s*vk[d];
@@ -1057,7 +1083,31 @@ void VectorFiniteElement::Project_RT(
    }
    else
    {
-      mfem_error("VectorFiniteElement::Project_RT (fe version)");
+      int sdim = Trans.GetSpaceDim();
+      double vk[Geometry::MaxDim];
+      DenseMatrix vshape(fe.GetDof(), sdim);
+      Vector vshapenk(fe.GetDof());
+      const bool square_J = (dim == sdim);
+
+      I.SetSize(dof, fe.GetDof());
+      for (int k = 0; k < dof; k++)
+      {
+         const IntegrationPoint &ip = Nodes.IntPoint(k);
+
+         Trans.SetIntPoint(&ip);
+         // Transform RT face normals from reference to physical space
+         // vk = adj(J)^T nk
+         Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*dim, vk);
+         // Compute fe basis functions in physical space
+         fe.CalcVShape(Trans, vshape);
+         // Project fe basis functions onto transformed face normals
+         vshape.Mult(vk, vshapenk);
+         if (!square_J) { vshapenk /= Trans.Weight(); }
+         for (int j=0; j<vshapenk.Size(); j++)
+         {
+            I(k,j) = vshapenk(j);
+         }
+      }
    }
 }
 
@@ -1218,6 +1268,8 @@ void VectorFiniteElement::Project_ND(
 
          fe.CalcShape(ip, shape);
          Trans.SetIntPoint(&ip);
+         // Transform ND edge tengents from reference to physical space
+         // vk = J tk
          Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
          if (fe.GetMapType() == INTEGRAL)
          {
@@ -1235,6 +1287,8 @@ void VectorFiniteElement::Project_ND(
             {
                s = 0.0;
             }
+            // Project scalar basis function multiplied by each coordinate
+            // direction onto the transformed edge tangents
             for (int d = 0; d < sdim; d++)
             {
                I(k, j + d*shape.Size()) = s*vk[d];
@@ -1244,7 +1298,29 @@ void VectorFiniteElement::Project_ND(
    }
    else
    {
-      mfem_error("VectorFiniteElement::Project_ND (fe version)");
+      int sdim = Trans.GetSpaceDim();
+      double vk[Geometry::MaxDim];
+      DenseMatrix vshape(fe.GetDof(), sdim);
+      Vector vshapetk(fe.GetDof());
+
+      I.SetSize(dof, fe.GetDof());
+      for (int k = 0; k < dof; k++)
+      {
+         const IntegrationPoint &ip = Nodes.IntPoint(k);
+
+         Trans.SetIntPoint(&ip);
+         // Transform ND edge tangents from reference to physical space
+         // vk = J tk
+         Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
+         // Compute fe basis functions in physical space
+         fe.CalcVShape(Trans, vshape);
+         // Project fe basis functions onto transformed edge tangents
+         vshape.Mult(vk, vshapetk);
+         for (int j=0; j<vshapetk.Size(); j++)
+         {
+            I(k, j) = vshapetk(j);
+         }
+      }
    }
 }
 

fem/fe.hpp

@@ -819,11 +819,25 @@ protected:
    void CalcVShape_ND(ElementTransformation &Trans,
                       DenseMatrix &shape) const;
 
+   /** @brief Project a vector coefficient onto the RT basis functions
+       @param nk    Face normal vectors for this element type
+       @param d2n   Offset into nk for each degree of freedom
+       @param vc    Vector coefficient to be projected
+       @param Trans Transformation from reference to physical coordinates
+       @param dofs  Expansion coefficients for the approximation of vc
+   */
    void Project_RT(const double *nk, const Array<int> &d2n,
                    VectorCoefficient &vc, ElementTransformation &Trans,
                    Vector &dofs) const;
 
    /// Projects the vector of values given at FE nodes to RT space
+   /** Project vector values onto the RT basis functions
+       @param nk    Face normal vectors for this element type
+       @param d2n   Offset into nk for each degree of freedom
+       @param vc    Vector values at each interpolation point
+       @param Trans Transformation from reference to physical coordinates
+       @param dofs  Expansion coefficients for the approximation of vc
+   */
    void Project_RT(const double *nk, const Array<int> &d2n,
                    Vector &vc, ElementTransformation &Trans,
                    Vector &dofs) const;
@@ -833,6 +847,19 @@ protected:
       const double *nk, const Array<int> &d2n,
       MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
 
+   /** @brief Project vector-valued basis functions onto the RT basis functions
+       @param nk    Face normal vectors for this element type
+       @param d2n   Offset into nk for each degree of freedom
+       @param fe    Vector-valued finite element basis
+       @param Trans Transformation from reference to physical coordinates
+       @param I     Expansion coefficients for the approximation of each basis
+                    function
+
+       Note: If the FiniteElement, fe, is scalar-valued the projection will
+             assume that a FiniteElementSpace is being used to define a vector
+             field using the scalar basis functions for each component of the
+             vector field.
+   */
    void Project_RT(const double *nk, const Array<int> &d2n,
                    const FiniteElement &fe, ElementTransformation &Trans,
                    DenseMatrix &I) const;
@@ -852,11 +879,25 @@ protected:
                        const FiniteElement &fe, ElementTransformation &Trans,
                        DenseMatrix &curl) const;
 
+   /** @brief Project a vector coefficient onto the ND basis functions
+       @param tk    Edge tangent vectors for this element type
+       @param d2t   Offset into tk for each degree of freedom
+       @param vc    Vector coefficient to be projected
+       @param Trans Transformation from reference to physical coordinates
+       @param dofs  Expansion coefficients for the approximation of vc
+   */
    void Project_ND(const double *tk, const Array<int> &d2t,
                    VectorCoefficient &vc, ElementTransformation &Trans,
                    Vector &dofs) const;
 
    /// Projects the vector of values given at FE nodes to ND space
+   /** Project vector values onto the ND basis functions
+       @param tk    Edge tangent vectors for this element type
+       @param d2t   Offset into tk for each degree of freedom
+       @param vc    Vector values at each interpolation point
+       @param Trans Transformation from reference to physical coordinates
+       @param dofs  Expansion coefficients for the approximation of vc
+   */
    void Project_ND(const double *tk, const Array<int> &d2t,
                    Vector &vc, ElementTransformation &Trans,
                    Vector &dofs) const;
@@ -866,6 +907,19 @@ protected:
       const double *tk, const Array<int> &d2t,
       MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const;
 
+   /** @brief Project vector-valued basis functions onto the ND basis functions
+       @param tk    Edge tangent vectors for this element type
+       @param d2t   Offset into tk for each degree of freedom
+       @param fe    Vector-valued finite element basis
+       @param Trans Transformation from reference to physical coordinates
+       @param I     Expansion coefficients for the approximation of each basis
+                    function
+
+       Note: If the FiniteElement, fe, is scalar-valued the projection will
+             assume that a FiniteElementSpace is being used to define a vector
+             field using the scalar basis functions for each component of the
+             vector field.
+   */
    void Project_ND(const double *tk, const Array<int> &d2t,
                    const FiniteElement &fe, ElementTransformation &Trans,
                    DenseMatrix &I) const;

fem/gridfunc.cpp

@@ -467,15 +467,26 @@ const
    fes->GetElementDofs(i, dofs);
    fes->DofsToVDofs(vdim-1, dofs);
    const FiniteElement *FElem = fes->GetFE(i);
-   MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
-               "invalid FE map type");
    int dof = FElem->GetDof();
    Vector DofVal(dof), loc_data(dof);
    GetSubVector(dofs, loc_data);
-   for (int k = 0; k < n; k++)
+   if (FElem->GetMapType() == FiniteElement::VALUE)
    {
-      FElem->CalcShape(ir.IntPoint(k), DofVal);
-      vals(k) = DofVal * loc_data;
+      for (int k = 0; k < n; k++)
+      {
+         FElem->CalcShape(ir.IntPoint(k), DofVal);
+         vals(k) = DofVal * loc_data;
+      }
+   }
+   else
+   {
+      ElementTransformation *Tr = fes->GetElementTransformation(i);
+      for (int k = 0; k < n; k++)
+      {
+         Tr->SetIntPoint(&ir.IntPoint(k));
+         FElem->CalcPhysShape(*Tr, DofVal);
+         vals(k) = DofVal * loc_data;
+      }
    }
 }
 
@@ -989,15 +1000,14 @@ void GridFunction::GetVectorValues(ElementTransformation &T,
 
    if (FElem->GetRangeType() == FiniteElement::SCALAR)
    {
-      MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
-                  "invalid FE map type");
       Vector shape(dof);
       int vdim = fes->GetVDim();
       vals.SetSize(vdim, nip);
       for (int j = 0; j < nip; j++)
       {
          const IntegrationPoint &ip = ir.IntPoint(j);
-         FElem->CalcShape(ip, shape);
+         T.SetIntPoint(&ip);
+         FElem->CalcPhysShape(T, shape);
 
          for (int k = 0; k < vdim; k++)
          {

linalg/hypre.cpp

@@ -2164,6 +2164,7 @@ HypreSmoother::HypreSmoother() : Solver()
    B = X = V = Z = NULL;
    X0 = X1 = NULL;
    fir_coeffs = NULL;
+   A_is_symmetric = false;
 }
 
 HypreSmoother::HypreSmoother(HypreParMatrix &_A, int _type,
@@ -2183,6 +2184,7 @@ HypreSmoother::HypreSmoother(HypreParMatrix &_A, int _type,
    B = X = V = Z = NULL;
    X0 = X1 = NULL;
    fir_coeffs = NULL;
+   A_is_symmetric = false;
 
    SetOperator(_A);
 }
@@ -2470,6 +2472,16 @@ void HypreSmoother::Mult(const Vector &b, Vector &x) const
    Mult(*B, *X);
 }
 
+void HypreSmoother::MultTranspose(const Vector &b, Vector &x) const
+{
+   if (A_is_symmetric || type == 0 || type == 1 || type == 5)
+   {
+      Mult(b, x);
+      return;
+   }
+   mfem_error("HypreSmoother::MultTranspose (...) : undefined!\n");
+}
+
 HypreSmoother::~HypreSmoother()
 {
    if (B) { delete B; }

linalg/hypre.hpp

@@ -638,6 +638,9 @@ protected:
    /// Combined coefficients for windowing and Chebyshev polynomials.
    double* fir_coeffs;
 
+   /// A flag that indicates whether the linear system matrix A is symmetric
+   bool A_is_symmetric;
+
 public:
    /** Hypre smoother types:
        0    = Jacobi
@@ -684,6 +687,12 @@ public:
        entries in the associated matrix. */
    void SetPositiveDiagonal(bool pos = true) { pos_l1_norms = pos; }
 
+   /** Explicitly indicate whether the linear system matrix A is symmetric. If A
+       is symmetric, the smoother will also be symmetric. In this case, calling
+       MultTranspose will be redirected to Mult. (This is also done if the
+       smoother is diagonal.) By default, A is assumed to be nonsymmetric. */
+   void SetOperatorSymmetry(bool is_sym) { A_is_symmetric = is_sym; }
+
    /** Set/update the associated operator. Must be called after setting the
        HypreSmoother type and options. */
    virtual void SetOperator(const Operator &op);
@@ -692,6 +701,9 @@ public:
    virtual void Mult(const HypreParVector &b, HypreParVector &x) const;
    virtual void Mult(const Vector &b, Vector &x) const;
 
+   /// Apply transpose of the smoother to relax the linear system Ax=b
+   virtual void MultTranspose(const Vector &b, Vector &x) const;
+
    virtual ~HypreSmoother();
 };
 

linalg/operator.hpp

@@ -651,7 +651,6 @@ public:
    virtual void SetOperator(const Operator &op) = 0;
 };
 
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {

linalg/solvers.cpp

@@ -1556,6 +1556,7 @@ void NewtonSolver::SetOperator(const Operator &op)
    width = op.Width();
    MFEM_ASSERT(height == width, "square Operator is required.");
 
+   xcur.SetSize(width);
    r.SetSize(width);
    c.SetSize(width);
 }
@@ -1615,9 +1616,20 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
          break;
       }
 
-      prec->SetOperator(oper->GetGradient(x));
+      grad = &oper->GetGradient(x);
+      prec->SetOperator(*grad);
 
-      prec->Mult(r, c);  // c = [DF(x_i)]^{-1} [F(x_i)-b]
+      if (lin_rtol_type)
+      {
+         AdaptiveLinRtolPreSolve(x, it, norm);
+      }
+
+      prec->Mult(r, c); // c = [DF(x_i)]^{-1} [F(x_i)-b]
+
+      if (lin_rtol_type)
+      {
+         AdaptiveLinRtolPostSolve(c, r, it, norm);
+      }
 
       const double c_scale = ComputeScalingFactor(x, b);
       if (c_scale == 0.0)
@@ -1641,6 +1653,403 @@ void NewtonSolver::Mult(const Vector &b, Vector &x) const
    final_norm = norm;
 }
 
+void NewtonSolver::SetAdaptiveLinRtol(const int type,
+                                      const double rtol0,
+                                      const double rtol_max,
+                                      const double alpha,
+                                      const double gamma)
+{
+   lin_rtol_type = type;
+   lin_rtol0 = rtol0;
+   lin_rtol_max = rtol_max;
+   this->alpha = alpha;
+   this->gamma = gamma;
+}
+
+void NewtonSolver::AdaptiveLinRtolPreSolve(const Vector &x,
+                                           const int it,
+                                           const double fnorm) const
+{
+   // Assume that when adaptive linear solver relative tolerance is activated,
+   // we are working with an iterative solver.
+   auto iterative_solver = static_cast<IterativeSolver *>(prec);
+   // Adaptive linear solver relative tolerance
+   double eta;
+   // Safeguard threshold
+   double sg_threshold = 0.1;
+
+   if (it == 0)
+   {
+      eta = lin_rtol0;
+   }
+   else
+   {
+      if (lin_rtol_type == 1)
+      {
+         // eta = gamma * abs(||F(x1)|| - ||F(x0) + DF(x0) s0||) / ||F(x0)||
+         eta = gamma * abs(fnorm - lnorm_last) / fnorm_last;
+      }
+      else if (lin_rtol_type == 2)
+      {
+         // eta = gamma * (||F(x1)|| / ||F(x0)||)^alpha
+         eta = gamma * pow(fnorm / fnorm_last, alpha);
+      }
+      else
+      {
+         MFEM_ABORT("Unknown adaptive linear solver rtol version");
+      }
+
+      // Safeguard rtol from "oversolving" ?!
+      const double sg_eta = gamma * pow(eta_last, alpha);
+      if (sg_eta > sg_threshold) { eta = std::max(eta, sg_eta); }
+   }
+
+   eta = std::min(eta, lin_rtol_max);
+   iterative_solver->SetRelTol(eta);
+   eta_last = eta;
+   if (print_level >= 0)
+   {
+      mfem::out << "Eisenstat-Walker rtol = " << eta << "\n";
+   }
+}
+
+void NewtonSolver::AdaptiveLinRtolPostSolve(const Vector &x,
+                                            const Vector &b,
+                                            const int it,
+                                            const double fnorm) const
+{
+   fnorm_last = fnorm;
+
+   // If version 1 is chosen, the true linear residual norm has to be computed
+   // and in most cases we can only retrieve the preconditioned linear residual
+   // norm.
+   if (lin_rtol_type == 1)
+   {
+      // lnorm_last = ||F(x0) + DF(x0) s0||
+      Vector linres(x.Size());
+      grad->Mult(x, linres);
+      linres -= b;
+      lnorm_last = Norm(linres);
+   }
+}
+
+void AndersonAcceleration::QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const
+{
+   Vector temp(Q[0]->Size());
+   for (int i=0; i<(maxVecs-1); i++) {
+      double d = sqrt( R(i,i+1)*R(i,i+1) + R(i+1,i+1)*R(i+1,i+1) );
+      double c = R(i,i+1) / d;
+      double s = R(i+1,i+1) / d;
+      R(i,i+1) = d;
+      R(i+1,i+1) = 0;
+
+      if (i < (maxVecs-2)) {
+         for (int j=(i+2); j<maxVecs; j++) {
+            d = c*R(i,j) + s*R(i+1,j);
+            R(i+1,j) = -s*R(i,j) + c*R(i+1,j);
+            R(i,j) = d;
+         }
+      }
+      // temp = c*Q[i] + s*Q[i+1];
+      add(c, *(Q[i]), s, *(Q[i+1]), temp);
+      // Q[i+1] = -s*Q[i] + c*Q[i+1];
+      *(Q[i+1]) *= c;
+      Q[i+1] -> Add(-s, *(Q[i]));
+      *(Q[i]) = temp;
+   }
+   
+   // Shift Q <- Q[:,0:(m-2)], i.e., delete last column of Q
+   delete Q.back();
+   Q.pop_back();
+
+   // Shift columns of R to the left by one, R = R[0:(m−2), 1:(m-1)]
+   for (int j=1; j<maxVecs; j++) {
+      for (int i=0; i<maxVecs; i++) {
+         R(i,j-1) = R(i,j);
+      }
+   }
+}
+
+void AndersonAcceleration::SetOperator(const Operator &op)
+{
+   oper = &op;
+   height = op.Height();
+   width = op.Width();
+   MFEM_ASSERT(height == width, "square Operator is required.");
+}
+
+void AndersonAcceleration::FixedPointMult(const Vector &b,
+   const Vector &x, Vector &y, Vector &r) const
+{
+   // Assume Operator represents a fixed-point operator *with
+   // right-hand side*, so we iterate x_{k+1} = M^{-1}G(x_k)
+   if (isFixedPointOp) {
+      oper->Mult(x, y);
+      if (prec) {
+         prec->Mult(y,r);
+         y = r;
+         r -= x;
+      }
+      else {
+         r = y;
+         r -= x;
+      }
+   }
+   // Otherwise add x to y and not to residual r = y - x
+   else {
+      oper->Mult(x, r);
+      r *= -1;
+      r += b;
+      if (prec) {
+         prec->Mult(r,y);
+         r = y;
+         y += x;
+      }
+      else {
+         y = r;
+         y += x;         
+      }
+   }
+}
+
+void AndersonAcceleration::Mult(const Vector &b, Vector &x) const
+{
+   MFEM_ASSERT(oper != NULL, "the Operator is not set (use SetOperator).");
+  
+   int n = width;
+   int numVecs = 0;
+   double resid, norm_df;
+   double min_diag = 1e-13;
+   final_norm = -1;
+
+   // Check vector is initialized, set to zero for
+   // iterative_mode = false
+   if (x.Size() != n) {
+      x.SetSize(n);
+      x = 0.0;
+   }
+   else if (!iterative_mode) {
+      x = 0.0;
+   }
+
+   // Storage containers for acceleration
+   std::deque<Vector *> G;
+   std::deque<Vector *> Q;
+   DenseMatrix R(maxVecs);
+   R = 0.0;
+   Vector g_old(n);
+   Vector g_current(n);
+   Vector f_old(n);
+   Vector f_current(n);
+   Vector gamma(maxVecs);
+   Vector rhs(maxVecs);
+   Vector correction;
+   if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+      correction.SetSize(n);
+   }
+   Vector *dg;
+   Vector *df;
+
+   // Loop over AA iterations
+   int k;
+   for (k=0; k<max_iter; k++) {
+
+      // Compute g_current = G(x), f_current = G(x) - x
+      this->FixedPointMult(b, x, g_current, f_current); 
+
+      // Check norm of current approximation to fixed point G(u) = u
+      resid = Norm(f_current);
+      MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+      if (print_level == 1)
+      {
+         mfem::out << "  Iteration : " << setw(3) << k
+                 << "  ||G(u) - u|| = " << resid << endl;
+      }
+      
+      // Set stopping tolerance on first iteration.
+      if (final_norm < 0) {
+         final_norm = std::max(rel_tol*resid, abs_tol);
+      }
+
+      // Check for convergence
+      if (resid <= final_norm)
+      {
+         final_norm = resid;
+         final_iter = k;
+         converged = 1;
+         goto finish;
+      }
+
+      // Start Anderson Acceleration after AAstart FP iterations
+      if (k > AAstart) {
+         // df = f_current - f_old; 
+         df = new Vector(n);
+         add(1.0, f_current, -1.0, f_old, *df);
+         // dg = g_current - g_old;
+         dg = new Vector(n);
+         add(1.0, g_current, -1.0, g_old, *dg);
+      
+         if (numVecs < maxVecs) {
+            G.push_back(dg);
+         }
+         else {
+            delete G[0];
+            G.pop_front();
+            G.push_back(dg);
+         }
+         numVecs++;
+         dg = NULL;
+      }
+      
+      f_old = f_current;
+      g_old = g_current;
+      
+      // First iteration or initial fixed-point iterations
+      if (numVecs == 0) {
+         x = g_current;
+         continue;
+      }
+
+      // All later iterations: orthogonalize and find best approximation
+      if (numVecs == 1) {
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(0,0) = norm_df;
+         df = NULL;
+      }
+      else {
+         // Remove first column in basis F and R, reorthogonalize
+         if (numVecs > maxVecs) {
+            this->QRdelete(Q, R);
+            numVecs--;
+         }
+         // Compute last column of R
+         for (int i=0; i<(numVecs-1); i++) {
+            R(i,numVecs-1) = Dot(*(Q[i]), *df);
+            // df -= R(i,numVecs-1) * Q[i]
+            df -> Add(-R(i,numVecs-1), *(Q[i]));
+         }
+         norm_df = Norm(*df);
+         MFEM_ASSERT(IsFinite(norm_df), "norm_df = " << norm_df);
+         (*df) /= norm_df;
+         Q.push_back(df);
+         R(numVecs-1, numVecs-1) = norm_df;
+         df = NULL;
+      }
+
+      // Back solve for new weights, R\gamma = Q^T * f_current
+      rhs = 0.0;
+      gamma = 0.0;
+      for (int i=0; i<numVecs; i++) {
+         rhs(i) = Dot(*(Q[i]), f_current);   // Form right hand side
+      }
+      for (int i=(numVecs-1); i>=0; i--) {
+         double temp = rhs(i);
+         for (int j=(i+1); j<numVecs; j++) {
+            temp -= R(i,j)*gamma(j);
+         }
+         if (std::abs(R(i,i)) < min_diag) {
+            gamma(i) = 0.0;
+            std::cout << "Diagonal of R -- " << R(i,i) << " ~ 0.\n";
+         }
+         else {
+            gamma(i) = temp / R(i,i);            
+         }
+      }
+
+      /// DEBUG --> test backsolve
+      Vector test(numVecs);
+      for (int i=0; i<numVecs; i++) {
+         test(i) = 0;
+         for (int j=i; j<numVecs; j++) {
+            test(i) += R(i,j) * gamma(j);
+         }
+         if (std::abs(rhs(i) - test(i)) > 1e-10) {
+            std::cout << "Bad solve! Err = " << rhs(i) - test(i) << "\n";
+         }
+      }
+
+      // Compute updated solution x = g_current − G*\gamma 
+      x = g_current;
+      for (int i=0; i<numVecs; i++) {
+         // x -= gamma(i)*G[i]
+         x.Add(-gamma(i), *(G[i]));
+      }
+
+      // Apply damped iteration for \omega \in (0,1),
+      //   x -= (1−omega) * (f_current − Q*R*gamma);
+      if (omega > 0 && std::abs(omega -  1) > 1e-14) {
+         // Redefine rhs = R*gamma
+         for(int i=0; i<numVecs; i++) {
+            rhs(i) = 0;
+            for (int j=i; j<numVecs; j++) {
+               rhs(i) += R(i,j)*gamma(j);
+            }
+         }
+         correction = f_current;
+         for (int i=0; i<numVecs; i++) {
+            // correction -= rhs(i)*Q[i]
+            correction.Add(-rhs(i), *(Q[i]));
+         }
+         // x -= (1 - omega) * correction;
+         x.Add( -(1 - omega), correction);
+      }
+
+      // Restart AA minimization by eliminating all vectors but the most recent
+      if (restart && (numVecs == maxVecs)) {
+         for (int i=0; i<(maxVecs-1); i++) {
+            delete G[0];
+            G.pop_front();          
+            delete Q[0];
+            Q.pop_front();
+         }
+         R = 0.0;
+         R(0,0) = norm_df;
+         numVecs = 1;
+         if (print_level == 1)
+         {
+            mfem::out << "Restarting..." << '\n';
+         }
+      }
+   }
+
+   // Compute final residual, save counts for solve
+   this->FixedPointMult(b, x, g_current, f_current); 
+   resid = Norm(f_current);
+   MFEM_ASSERT(IsFinite(resid), "||G(u) - u|| = " << resid);
+   final_norm = resid;
+   final_iter = max_iter;
+   if (resid <= final_norm) converged = 1;
+   else converged = 0;
+
+finish:
+   if (print_level == 3)
+   {
+      mfem::out << "  Iteration : " << setw(3) << k
+              << "  ||G(u) - u|| = " << resid << endl;
+   }
+   else if (print_level == 2)
+   {
+      mfem::out << "Anderson Acceleration: Number of iterations: " << final_iter << '\n';
+   }
+   if (print_level >= 0 && !converged)
+   {
+      mfem::out << "Anderson Acceleration: No convergence!\n";
+   }
+
+   // Cleanup pointers
+   for (int i=0; i<numVecs; i++) {
+      delete G[0];
+      G.pop_front();          
+      delete Q[0];
+      Q.pop_front();
+   }
+   delete dg;
+   delete df;
+}
+
 void LBFGSSolver::Mult(const Vector &b, Vector &x) const
 {
    MFEM_VERIFY(oper != NULL, "the Operator is not set (use SetOperator).");

linalg/solvers.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "densemat.hpp"
+#include <deque>
 
 #ifdef MFEM_USE_MPI
 #include <mpi.h>
@@ -406,7 +407,40 @@ void MINRES(const Operator &A, Solver &B, const Vector &b, Vector &x,
 class NewtonSolver : public IterativeSolver
 {
 protected:
-   mutable Vector r, c;
+   mutable Vector xcur, r, c;
+   mutable Operator *grad;
+
+   // Adaptive linear solver rtol variables
+
+   // Method to determine rtol, 0 means the adaptive algorithm is deactivated.
+   int lin_rtol_type = 0;
+   // rtol to use in first iteration
+   double lin_rtol0;
+   // Maximum rtol
+   double lin_rtol_max;
+   // Function norm ||F(x)|| of the previous iterate
+   mutable double fnorm_last = 0.0;
+   // Linear residual norm of the previous iterate
+   mutable double lnorm_last = 0.0;
+   // Forcing term (linear residual rtol) from the previous iterate
+   mutable double eta_last = 0.0;
+   // Eisenstat-Walker factor gamma
+   double gamma;
+   // Eisenstat-Walker factor alpha
+   double alpha;
+
+   /** @brief Method for the adaptive linear solver rtol invoked before the
+       linear solve. */
+   void AdaptiveLinRtolPreSolve(const Vector &x,
+                                const int it,
+                                const double fnorm) const;
+
+   /** @brief Method for the adaptive linear solver rtol invoked after the
+       linear solve. */
+   void AdaptiveLinRtolPostSolve(const Vector &x,
+                                 const Vector &b,
+                                 const int it,
+                                 const double fnorm) const;
 
 public:
    NewtonSolver() { }
@@ -434,8 +468,86 @@ public:
    /** @brief This method can be overloaded in derived classes to perform
        computations that need knowledge of the newest Newton state. */
    virtual void ProcessNewState(const Vector &x) const { }
+
+   const Vector &GetCurrentResidual() const { return r; }
+   const Vector &GetCurrentIterate() const { return xcur; }
+
+   /// Enable adaptive linear solver relative tolerance algorithm.
+   /** Compute a relative tolerance for the Krylov method after each nonlinear
+    iteration, based on the algorithm presented in [1].
+
+    The maximum linear solver relative tolerance @a rtol_max should be < 1. For
+    @a type 1 the parameters @a alpha and @a gamma are ignored. For @a type 2
+    @a alpha has to be between 0 and 1 and @a gamma between 1 and 2.
+
+    [1] Eisenstat, Stanley C., and Homer F. Walker. "Choosing the forcing terms
+    in an inexact Newton method."
+    */
+   void SetAdaptiveLinRtol(const int type = 2,
+                           const double rtol0 = 0.5,
+                           const double rtol_max = 0.9,
+                           const double alpha = 0.5 * (1.0 + sqrt(5.0)),
+                           const double gamma = 1.0);
 };
 
+
+/// Nonlinear Anderson Acceleration
+class AndersonAcceleration : public IterativeSolver
+{
+protected:
+   int maxVecs; // see SetKDim()
+   int AAstart;
+   bool restart;
+   bool isFixedPointOp;
+   double omega;
+
+   /// Apply fixed-point Mult and compute residual.
+   //  For !isFixedPointOp:
+   //    y <-- G(x) = x + A(x) - b and r <-- G(x) - x = A(x) - b
+   //  For isFixedPointOp:
+   //    y <-- G(x) = A(x) - b and r <-- G(x) - x = A(x) - b - x
+   void FixedPointMult(const Vector &b, const Vector &x, Vector &y, Vector &r) const;
+
+   // Helper function for Anderson Acceleration
+   void QRdelete(std::deque<Vector *> &Q, DenseMatrix &R) const;
+
+public:
+   AndersonAcceleration() : maxVecs(25), AAstart(0), omega(1),
+      restart(false), isFixedPointOp(false) { }
+
+#ifdef MFEM_USE_MPI
+   AndersonAcceleration(MPI_Comm _comm) : IterativeSolver(_comm),
+      maxVecs(25), AAstart(0), omega(1), restart(false),
+      isFixedPointOp(false) { }
+#endif
+
+   /// Boolean describing whether the action of the operator is such
+   //  that we want to solve G(x) = x (true) or G(x) = 0 (false) for
+   //  zero right-hand side vector b passed into Mult(). Default is
+   //  false in construction of class.
+   void IsFixedPointOperator(bool isFixedPointOp_)
+   { isFixedPointOp = isFixedPointOp_; }
+
+   /// Maximum number of vectors to store in Krylov-like space
+   void SetKDim(int dim) { maxVecs = dim; }
+
+   /// Number of fixed-point iterations to do before starting AA
+   void SetAAStart(int start_) { AAstart = start_; }
+
+   /// Boolean to restart, that is, erase entire space after maxVecs
+   //  are stored (AAstart=true) or use a sliding space (AAstart=false)
+   //  where one vector is deleted to make room for a new one.
+   void SetRestart(bool restart_) { restart = restart_; }
+
+   /// Set relaxation weight
+   void SetWeight(double omega_) { omega = omega_; }
+
+   virtual void SetOperator(const Operator &op);
+
+   virtual void Mult(const Vector &b, Vector &x) const;
+};
+
+
 /** L-BFGS method for solving F(x)=b for a given operator F, by minimizing
     the norm of F(x) - b. Requires only the action of the operator F. */
 class LBFGSSolver : public NewtonSolver

makefile

@@ -536,7 +536,7 @@ clean: $(addsuffix /clean,$(EM_DIRS) $(TEST_DIRS))
 distclean: clean config/clean doc/clean
 	rm -rf mfem/
 
-INSTALL_SHARED_LIB = $(MFEM_CXX) $(MFEM_BUILD_FLAGS) $(INSTALL_SOFLAGS)\
+INSTALL_SHARED_LIB = $(MFEM_CXX) $(MFEM_LINK_FLAGS) $(INSTALL_SOFLAGS)\
    $(OBJECT_FILES) $(EXT_LIBS) -o $(PREFIX_LIB)/libmfem.$(SO_VER) && \
    cd $(PREFIX_LIB) && ln -sf libmfem.$(SO_VER) libmfem.$(SO_EXT)
 

miniapps/meshing/mesh-optimizer.cpp

@@ -61,15 +61,15 @@
 //   * mesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
 //   Blade shape:
-//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8
+//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
 //   Blade shape with FD-based solver:
-//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8 -fd
+//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
 //   Blade limited shape:
-//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8 -lc 5000
+//     mesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
 //   ICF shape and equal size:
-//     mesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -qo 5
+//     mesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -art 2 -qo 5
 //   ICF shape and initial size:
-//     mesh-optimizer -o 3 -mid 9 -tid 3 -ni 100 -bnd -qt 1 -qo 8
+//     mesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
 //   ICF shape:
 //     mesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
 //   ICF limited shape:
@@ -114,6 +114,7 @@ int main(int argc, char *argv[])
    int solver_type       = 0;
    int solver_iter       = 20;
    double solver_rtol    = 1e-10;
+   int solver_art_type   = 0;
    int lin_solver        = 2;
    int max_lin_iter      = 100;
    bool move_bnd         = true;
@@ -192,12 +193,17 @@ int main(int argc, char *argv[])
                   "Maximum number of Newton iterations.");
    args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
                   "Relative tolerance for the Newton solver.");
+   args.AddOption(&solver_art_type, "-art", "--adaptive-rel-tol",
+                  "Type of adaptive relative linear solver tolerance:\n\t"
+                  "0: None (default)\n\t"
+                  "1: Eisenstat-Walker type 1\n\t"
+                  "2: Eisenstat-Walker type 2");
    args.AddOption(&lin_solver, "-ls", "--lin-solver",
                   "Linear solver:\n\t"
                   "0: l1-Jacobi\n\t"
                   "1: CG\n\t"
                   "2: MINRES\n\t"
-                  "3: MINRES + Jacobi preconditioner"
+                  "3: MINRES + Jacobi preconditioner\n\t"
                   "4: MINRES + l1-Jacobi preconditioner");
    args.AddOption(&max_lin_iter, "-li", "--lin-iter",
                   "Maximum number of iterations in the linear solve.");
@@ -205,7 +211,7 @@ int main(int argc, char *argv[])
                   "--fix-boundary",
                   "Enable motion along horizontal and vertical boundaries.");
    args.AddOption(&combomet, "-cmb", "--combo-type",
-                  "Combination of metrics options:"
+                  "Combination of metrics options:\n\t"
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
@@ -845,6 +851,10 @@ int main(int argc, char *argv[])
    solver.SetMaxIter(solver_iter);
    solver.SetRelTol(solver_rtol);
    solver.SetAbsTol(0.0);
+   if (solver_art_type > 0)
+   {
+      solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
+   }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
    solver.SetOperator(a);
    solver.Mult(b, x.GetTrueVector());

miniapps/meshing/pmesh-optimizer.cpp

@@ -61,15 +61,15 @@
 //   * mpirun -np 4 pmesh-optimizer -m stretched2D.mesh -o 2 -mid 2 -tid 1 -ni 50 -qo 5 -nor -vl 1 -alc 0.5 -fd -ae 1
 //
 //   Blade shape:
-//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8
+//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 3 -art 1 -bnd -qt 1 -qo 8
 //   Blade shape with FD-based solver:
-//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8 -fd
+//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 30 -ls 4 -bnd -qt 1 -qo 8 -fd
 //   Blade limited shape:
-//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -ni 200 -bnd -qt 1 -qo 8 -lc 5000
+//     mpirun -np 4 pmesh-optimizer -m blade.mesh -o 4 -mid 2 -tid 1 -bnd -qt 1 -qo 8 -lc 5000
 //   ICF shape and equal size:
-//     mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -qo 5
+//     mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 2 -ni 25 -ls 3 -art 2 -qo 5
 //   ICF shape and initial size:
-//     mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 3 -ni 100 -bnd -qt 1 -qo 8
+//     mpirun -np 4 pmesh-optimizer -o 3 -mid 9 -tid 3 -ni 30 -ls 3 -bnd -qt 1 -qo 8
 //   ICF shape:
 //     mpirun -np 4 pmesh-optimizer -o 3 -mid 1 -tid 1 -ni 100 -bnd -qt 1 -qo 8
 //   ICF limited shape:
@@ -121,6 +121,7 @@ int main (int argc, char *argv[])
    int solver_type       = 0;
    int solver_iter       = 20;
    double solver_rtol    = 1e-10;
+   int solver_art_type   = 0;
    int lin_solver        = 2;
    int max_lin_iter      = 100;
    bool move_bnd         = true;
@@ -201,12 +202,17 @@ int main (int argc, char *argv[])
                   "Maximum number of Newton iterations.");
    args.AddOption(&solver_rtol, "-rtol", "--newton-rel-tolerance",
                   "Relative tolerance for the Newton solver.");
+   args.AddOption(&solver_art_type, "-art", "--adaptive-rel-tol",
+                  "Type of adaptive relative linear solver tolerance:\n\t"
+                  "0: None (default)\n\t"
+                  "1: Eisenstat-Walker type 1\n\t"
+                  "2: Eisenstat-Walker type 2");
    args.AddOption(&lin_solver, "-ls", "--lin-solver",
                   "Linear solver:\n\t"
                   "0: l1-Jacobi\n\t"
                   "1: CG\n\t"
                   "2: MINRES\n\t"
-                  "3: MINRES + Jacobi preconditioner"
+                  "3: MINRES + Jacobi preconditioner\n\t"
                   "4: MINRES + l1-Jacobi preconditioner");
    args.AddOption(&max_lin_iter, "-li", "--lin-iter",
                   "Maximum number of iterations in the linear solve.");
@@ -214,7 +220,7 @@ int main (int argc, char *argv[])
                   "--fix-boundary",
                   "Enable motion along horizontal and vertical boundaries.");
    args.AddOption(&combomet, "-cmb", "--combo-type",
-                  "Combination of metrics options:"
+                  "Combination of metrics options:\n\t"
                   "0: Use single metric\n\t"
                   "1: Shape + space-dependent size given analytically\n\t"
                   "2: Shape + adapted size given discretely; shared target");
@@ -895,6 +901,10 @@ int main (int argc, char *argv[])
    solver.SetMaxIter(solver_iter);
    solver.SetRelTol(solver_rtol);
    solver.SetAbsTol(0.0);
+   if (solver_art_type > 0)
+   {
+      solver.SetAdaptiveLinRtol(solver_art_type, 0.5, 0.9);
+   }
    solver.SetPrintLevel(verbosity_level >= 1 ? 1 : -1);
    solver.SetOperator(a);
    solver.Mult(b, x.GetTrueVector());