From f012c46d7ab9f360b6cc6eb212914b806e575fde Mon Sep 17 00:00:00 2001
From: Peter Munch <peterrmuench@gmail.com>
Date: Thu, 29 Dec 2022 20:48:45 +0100
Subject: [PATCH] Apply suggestions

---
 cmake/configure/configure_20_trilinos.cmake |   2 +-
 include/deal.II/base/config.h.in            |   1 +
 include/deal.II/trilinos/nox.h              | 110 ++++++++++++-----
 include/deal.II/trilinos/nox.templates.h    | 127 +++++++++++++-------
 source/trilinos/CMakeLists.txt              |  12 +-
 source/trilinos/nox.cc                      |   2 +-
 6 files changed, 167 insertions(+), 87 deletions(-)

diff --git a/cmake/configure/configure_20_trilinos.cmake b/cmake/configure/configure_20_trilinos.cmake
index ea317beb7a..86d50873ed 100644
--- a/cmake/configure/configure_20_trilinos.cmake
+++ b/cmake/configure/configure_20_trilinos.cmake
@@ -150,7 +150,7 @@ macro(feature_trilinos_find_external var)
       #
       # Check for modules.
       #
-      foreach(_optional_module Belos EpetraExt Kokkos MueLu ROL Sacado SEACAS Tpetra Zoltan)
+      foreach(_optional_module Belos EpetraExt Kokkos MueLu NOX ROL Sacado SEACAS Tpetra Zoltan)
         item_matches(_module_found ${_optional_module} ${Trilinos_PACKAGE_LIST})
         if(_module_found)
           message(STATUS "  Found ${_optional_module}")
diff --git a/include/deal.II/base/config.h.in b/include/deal.II/base/config.h.in
index 47e0644ff1..933e46b873 100644
--- a/include/deal.II/base/config.h.in
+++ b/include/deal.II/base/config.h.in
@@ -242,6 +242,7 @@
 #cmakedefine DEAL_II_TRILINOS_WITH_BELOS
 #cmakedefine DEAL_II_TRILINOS_WITH_EPETRAEXT
 #cmakedefine DEAL_II_TRILINOS_WITH_MUELU
+#cmakedefine DEAL_II_TRILINOS_WITH_NOX
 #cmakedefine DEAL_II_TRILINOS_WITH_ROL
 #cmakedefine DEAL_II_TRILINOS_WITH_SACADO
 #cmakedefine DEAL_II_TRILINOS_WITH_SEACAS
diff --git a/include/deal.II/trilinos/nox.h b/include/deal.II/trilinos/nox.h
index 16b323de49..ff1ac236cb 100644
--- a/include/deal.II/trilinos/nox.h
+++ b/include/deal.II/trilinos/nox.h
@@ -18,7 +18,7 @@
 
 #include <deal.II/base/config.h>
 
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
 
 #  include <deal.II/base/exceptions.h>
 
@@ -26,6 +26,8 @@
 
 #  include <Teuchos_ParameterList.hpp>
 
+#  include <functional>
+
 DEAL_II_NAMESPACE_OPEN
 
 namespace TrilinosWrappers
@@ -35,9 +37,31 @@ namespace TrilinosWrappers
 
 
   /**
-   * Wrapper around the non-linear solver from the NOX
+   * Wrapper around the nonlinear solver from the NOX
    * packge (https://docs.trilinos.org/dev/packages/nox/doc/html/index.html),
    * targeting deal.II data structures.
+   *
+   * The following code shows the steps how to use this class:
+   * @code
+   * // set configuration
+   * TrilinosWrappers::NOXSolver<VectorType>::AdditionalData additional_data;
+   *
+   * // create nonlinear solver
+   * TrilinosWrappers::NOXSolver<VectorType> solver(additional_data);
+   *
+   * // set user functions to compute residual,
+   * // set up Jacobian, and to apply the inverse of
+   * // Jacobian; note that there are more functions that
+   * // can be set
+   * solver.residual = [](const auto &src, auto &dst) {...};
+   * solver.setup_jacobian = [](const auto &src) {...};
+   * solver.solve_with_jacobian =
+   *   [](const auto &src, auto &dst, const auto) {...};
+   *
+   * // solver nonlinear system with solution containing the intial guess and
+   * // the final solution
+   * solver.solve(solution);
+   * @endcode
    */
   template <typename VectorType>
   class NOXSolver
@@ -62,22 +86,28 @@ namespace TrilinosWrappers
                      const bool         reuse_solver                   = false);
 
       /**
-       * Max number of non-linear iterations.
+       * Max number of nonlinear iterations.
        */
       unsigned int max_iter;
 
       /**
-       * Absolute l2 tolerance to be reached.
+       * Absolute l2 tolerance of the residual to be reached.
+       *
+       * @note Solver terminates successfully if either the absolut or
+       * the relative tolerance has been reached.
        */
       double abs_tol;
 
       /**
-       * Relative l2 tolerance to be reached.
+       * Relative l2 tolerance of the residual to be reached.
+       *
+       * @note Solver terminates successfully if either the absolut or
+       * the relative tolerance has been reached.
        */
       double rel_tol;
 
       /**
-       * Number of non-linear iterations after which the preconditioner
+       * Number of nonlinear iterations after which the preconditioner
        * should be updated.
        */
       unsigned int threshold_nonlinear_iterations;
@@ -90,7 +120,10 @@ namespace TrilinosWrappers
       unsigned int threshold_n_linear_iterations;
 
       /**
-       * Reuse NOX solver instance in the next non-linear solution.
+       * Reuse nonlinear solver the next time solve() is called. In particular,
+       * this parameter allows to reuse the preconditioner from the last
+       * solution step, enabling preconditioner lagging over multiple nonlinear
+       * solutions.
        */
       bool reuse_solver;
     };
@@ -98,7 +131,7 @@ namespace TrilinosWrappers
     /**
      * Constructor.
      *
-     * If @p parameters is not filled, a Newton solver with full step is used.
+     * If @p parameters is not filled, a Newton solver with a full step is used.
      * An overview of possible parameters is given at
      * https://docs.trilinos.org/dev/packages/nox/doc/html/parameters.html.
      */
@@ -107,36 +140,39 @@ namespace TrilinosWrappers
                 Teuchos::rcp(new Teuchos::ParameterList));
 
     /**
-     * Clear internal state.
+     * Clear the internal state.
      */
     void
     clear();
 
     /**
-     * Solve non-linear problem and return number of iterations.
+     * Solve the nonlinear problem and return the number of nonlinear
+     * iterations executed.
      */
     unsigned int
     solve(VectorType &solution);
 
     /**
-     * User function that computes the residual.
+     * A user function that computes the residual @p f based on the
+     * current solution @p x.
      *
      * @note This function should return 0 in the case of success.
      */
     std::function<int(const VectorType &x, VectorType &f)> residual;
 
     /**
-     * User function that sets up the Jacobian.
+     * A user function that sets up the Jacobian, based on the
+     * current solution @p x.
      *
      * @note This function should return 0 in the case of success.
      */
     std::function<int(const VectorType &x)> setup_jacobian;
 
     /**
-     * User function that sets up the preconditioner for inverting
-     * the Jacobian.
+     * A user function that sets up the preconditioner for inverting
+     * the Jacobian, based on the current solution @p x.
      *
-     * @note The function is optional and is used when setup_jacobian is
+     * @note This function is optional and is used when setup_jacobian is
      * called and the preconditioner needs to be updated (see
      * update_preconditioner_predicate and
      * AdditionalData::threshold_nonlinear_iterations).
@@ -146,9 +182,10 @@ namespace TrilinosWrappers
     std::function<int(const VectorType &x)> setup_preconditioner;
 
     /**
-     * User function that applies the Jacobian.
+     * A user function that applies the Jacobian to @p x and writes
+     * the result in @p v.
      *
-     * @note The function is optional and is used in the case of certain
+     * @note This function is optional and is used in the case of certain
      * configurations. For instance, this function is required if the
      * polynomial line search (@p NOX::LineSearch::Polynomial) is
      * chosen, whereas for the full step case (@p NOX::LineSearch::FullStep)
@@ -159,9 +196,11 @@ namespace TrilinosWrappers
     std::function<int(const VectorType &x, VectorType &v)> apply_jacobian;
 
     /**
-     * User function that applies the inverse of the Jacobian.
+     * A user function that applies the inverse of the Jacobian to
+     * @p x and writes the result in @p x. The parameter @p tolerance
+     * species the error reduction if a interative solver is used.
      *
-     * @note The function is optional and is used in the case of certain
+     * @note This function is optional and is used in the case of certain
      * configurations.
      *
      * @note This function should return 0 in the case of success.
@@ -171,25 +210,29 @@ namespace TrilinosWrappers
       solve_with_jacobian;
 
     /**
-     * User function that applies the inverse of the Jacobian and
-     * returns the numer of linear iterations the linear solver needed.
+     * A user function that applies the inverse of the Jacobian to
+     * @p x, writes the result in @p x and returns the numer of
+     * linear iterations the linear solver needed.
+     * The parameter @p tolerance species the error reduction if a
+     * interative solver is used.
      *
-     * @note This function should return -1 in the case of failure.
+     * @note This function is optional and is used in the case of certain
+     * configurations.
      */
     std::function<
       int(const VectorType &f, VectorType &x, const double tolerance)>
       solve_with_jacobian_and_track_n_linear_iterations;
 
     /**
-     * User function that allows to check convergence in addition to
+     * A user function that allows to check convergence in addition to
      * ones checking the l2-norm and the number of iterations (see
-     * AdditionalData). It is run after each non-linear iteration.
+     * AdditionalData). It is run after each nonlinear iteration.
      *
      * The input are the current iteration number @p i, the l2-norm
      * @p norm_f of the residual vector, the current solution @p x,
      * and the current residual vector @p f.
      *
-     * @note The function is optional.
+     * @note This function is optional.
      */
     std::function<SolverControl::State(const unsigned int i,
                                        const double       norm_f,
@@ -198,12 +241,13 @@ namespace TrilinosWrappers
       check_iteration_status;
 
     /**
-     * Function that allows to force to update the preconditioner in
-     * addition to AdditionalData::threshold_nonlinear_iterations. A reason
+     * A user function that, in addition to
+     * AdditionalData::threshold_nonlinear_iterations,
+     * allows to force to update the preconditioner. A reason
      * for wanting to update the preconditioner is when the expected number
-     * of linear iterations exceeds.
+     * of linear iterations is exceeded.
      *
-     * @note The function is optional. If no function is attached, this
+     * @note This function is optional. If no function is attached, this
      * means implicitly a return value of false.
      */
     std::function<bool()> update_preconditioner_predicate;
@@ -222,22 +266,22 @@ namespace TrilinosWrappers
     const Teuchos::RCP<Teuchos::ParameterList> parameters;
 
     /**
-     * Counter for number of (accumulated) residual evaluations.
+     * A counter for the number of (accumulated) residual evaluations.
      */
     unsigned int n_residual_evaluations;
 
     /**
-     * Counter for number of (accumulated) Jacobi applications.
+     * A counter for the number of (accumulated) Jacobi applications.
      */
     unsigned int n_jacobian_applications;
 
     /**
-     * Counter for number of (accumulated) non-linear iterations.
+     * A counter for the number of (accumulated) nonlinear iterations.
      */
     unsigned int n_nonlinear_iterations;
 
     /**
-     * Number of linear iterations of the last Jacobian solve.
+     * The number of linear iterations of the last Jacobian solve.
      */
     unsigned int n_last_linear_iterations;
   };
diff --git a/include/deal.II/trilinos/nox.templates.h b/include/deal.II/trilinos/nox.templates.h
index 20c23eadfb..42c9346298 100644
--- a/include/deal.II/trilinos/nox.templates.h
+++ b/include/deal.II/trilinos/nox.templates.h
@@ -20,7 +20,7 @@
 
 #include <deal.II/trilinos/nox.h>
 
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
 
 #  include <NOX_Abstract_Group.H>
 #  include <NOX_Abstract_Vector.H>
@@ -31,6 +31,8 @@
 #  include <NOX_StatusTest_NormF.H>
 #  include <NOX_StatusTest_RelativeNormF.H>
 
+#  include <memory>
+
 DEAL_II_NAMESPACE_OPEN
 
 #  ifndef DOXYGEN
@@ -73,7 +75,7 @@ namespace TrilinosWrappers
         NOX::Abstract::Vector &
         init(double gamma) override
         {
-          *vector = gamma;
+          *vector = static_cast<typename VectorType::value_type>(gamma);
           return *this;
         }
 
@@ -143,7 +145,7 @@ namespace TrilinosWrappers
         NOX::Abstract::Vector &
         scale(double gamma) override
         {
-          *vector *= gamma;
+          *vector *= static_cast<typename VectorType::value_type>(gamma);
 
           return *this;
         }
@@ -175,7 +177,9 @@ namespace TrilinosWrappers
 
           Assert(a_, ExcInternalError());
 
-          vector->sadd(gamma, alpha, *a_->vector);
+          vector->sadd(static_cast<typename VectorType::value_type>(gamma),
+                       static_cast<typename VectorType::value_type>(alpha),
+                       *a_->vector);
 
           return *this;
         }
@@ -197,8 +201,12 @@ namespace TrilinosWrappers
           Assert(a_, ExcInternalError());
           Assert(b_, ExcInternalError());
 
-          vector->operator*=(gamma);
-          vector->add(alpha, *a_->vector, beta, *b_->vector);
+          vector->operator*=(
+            static_cast<typename VectorType::value_type>(gamma));
+          vector->add(static_cast<typename VectorType::value_type>(alpha),
+                      *a_->vector,
+                      static_cast<typename VectorType::value_type>(beta),
+                      *b_->vector);
 
           return *this;
         }
@@ -223,7 +231,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Norm.
+         * Return the vector norm.
          */
         double
         norm(NOX::Abstract::Vector::NormType type =
@@ -242,7 +250,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Weighted 2-Norm.
+         * Return the vector's weighted 2-Norm.
          */
         double
         norm(const NOX::Abstract::Vector &weights) const override
@@ -255,7 +263,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Inner product with y.
+         * Return the inner product of the vector with @p y.
          */
         double
         innerProduct(const NOX::Abstract::Vector &y) const override
@@ -288,7 +296,7 @@ namespace TrilinosWrappers
 
       private:
         /**
-         * Underlying deal.II vector.
+         * The underlying deal.II vector.
          */
         std::shared_ptr<VectorType> vector;
 
@@ -402,7 +410,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Set the solution vector x to y.
+         * Set the solution vector `x` to @y.
          */
         void
         setX(const NOX::Abstract::Vector &y) override
@@ -413,7 +421,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Compute x = grp.x + step * d.
+         * Compute the solution update `x = grp.x + step * d`.
          */
         void
         computeX(const NOX::Abstract::Group & grp,
@@ -430,7 +438,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Compute and store F(x).
+         * Compute and store the residual update `F(x)`.
          */
         NOX::Abstract::Group::ReturnType
         computeF() override
@@ -450,7 +458,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return true if F is valid.
+         * Return true if the residual vector `F` is valid.
          */
         bool
         isF() const override
@@ -459,7 +467,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Compute and store Jacobian.
+         * Compute and store the Jacobian.
          */
         NOX::Abstract::Group::ReturnType
         computeJacobian() override
@@ -485,7 +493,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return solution vector.
+         * Return the (total) solution vector.
          */
         const NOX::Abstract::Vector &
         getX() const override
@@ -494,7 +502,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return F(x).
+         * Return the residual `F(x)`.
          */
         const NOX::Abstract::Vector &
         getF() const override
@@ -503,7 +511,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return 2-norm of F(x)
+         * Return the 2-norm of `F(x)`.
          */
         double
         getNormF() const override
@@ -512,7 +520,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return gradient.
+         * Return the gradient.
          */
         const NOX::Abstract::Vector &
         getGradient() const override
@@ -521,7 +529,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return Newton direction.
+         * Return the Newton descent direction.
          */
         const NOX::Abstract::Vector &
         getNewton() const override
@@ -530,7 +538,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return RCP to solution vector.
+         * Return an `RCP` to solution vector.
          */
         Teuchos::RCP<const NOX::Abstract::Vector>
         getXPtr() const override
@@ -540,7 +548,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return RCP to F(x).
+         * Return an RCP to the residual `F(x)`.
          */
         Teuchos::RCP<const NOX::Abstract::Vector>
         getFPtr() const override
@@ -550,7 +558,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return RCP to gradient.
+         * Return an `RCP` to gradient.
          */
         Teuchos::RCP<const NOX::Abstract::Vector>
         getGradientPtr() const override
@@ -560,7 +568,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return RCP to Newton direction.
+         * Return an `RCP` to the Newton descent direction.
          */
         Teuchos::RCP<const NOX::Abstract::Vector>
         getNewtonPtr() const override
@@ -571,7 +579,8 @@ namespace TrilinosWrappers
 
         /**
          * Create a new Group of the same derived type as this one by
-         * cloning this one, and return a ref count pointer to the new group.
+         * cloning this one, and return a reference counting
+         * pointer to the new group.
          */
         Teuchos::RCP<NOX::Abstract::Group>
         clone(NOX::CopyType copy_type) const override
@@ -631,7 +640,8 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Compute the Newton direction, using parameters for the linear solve.
+         * Compute the Newton direction, using the chosen
+         * parameters for the linear solve.
          */
         NOX::Abstract::Group::ReturnType
         computeNewton(Teuchos::ParameterList &p) override
@@ -658,8 +668,8 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Applies Jacobian to the given input vector and puts
-         * the answer in the result.
+         * Applies the Jacobian to the given input vector and assigns
+         * the output to the @p result.
          */
         NOX::Abstract::Group::ReturnType
         applyJacobian(const NOX::Abstract::Vector &input,
@@ -692,19 +702,37 @@ namespace TrilinosWrappers
           is_valid_j = false;
         }
 
-        // internal vectors
-        Vector<VectorType> x, f, gradient, newton;
+        // internal vector for the current solution
+        Vector<VectorType> x;
+
+        // internal vector for the residual
+        Vector<VectorType> f;
+
+        // internal vector for the solution gradient
+        Vector<VectorType> gradient;
+
+        // internal vector for the newton step
+        Vector<VectorType> newton;
+
+        // helper function to compute residual
+        std::function<int(const VectorType &x, VectorType &f)> residual;
+
+        // helper function to setup Jacobian
+        std::function<int(const VectorType &x)> setup_jacobian;
 
-        // helper functions to compute residual, to setup jacobian, and
-        // solve jacobian
-        std::function<int(const VectorType &, VectorType &)> residual;
-        std::function<int(const VectorType &)>               setup_jacobian;
-        std::function<int(const VectorType &, VectorType &)> apply_jacobian;
-        std::function<int(const VectorType &, VectorType &, const double)>
+        // helper function to apply Jacobian
+        std::function<int(const VectorType &x, VectorType &v)> apply_jacobian;
+
+        // helper function to solve jacobian
+        std::function<
+          int(const VectorType &f, VectorType &x, const double tolerance)>
           solve_with_jacobian;
 
-        // internal state (are residuum and jacobian computed?)
-        bool is_valid_f, is_valid_j;
+        // internal state (is residual computed?)
+        bool is_valid_f;
+
+        // internal state (is Jacobian computed?)
+        bool is_valid_j;
       };
 
 
@@ -729,7 +757,7 @@ namespace TrilinosWrappers
         {}
 
         /**
-         * Check status.
+         * Check the status of the nonlinear solver.
          */
         NOX::StatusTest::StatusType
         checkStatus(const NOX::Solver::Generic &problem,
@@ -783,7 +811,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Return last return value of checkStatus().
+         * Return the last value that was given by `checkStatus()`.
          */
         NOX::StatusTest::StatusType
         getStatus() const override
@@ -792,7 +820,7 @@ namespace TrilinosWrappers
         }
 
         /**
-         * Print last return value of print().
+         * Print the last value that was given by `checkStatus()`
          */
         virtual std::ostream &
         print(std::ostream &stream, int indent = 0) const override
@@ -805,7 +833,8 @@ namespace TrilinosWrappers
 
       private:
         /**
-         * User function that allows to check convergence.
+         * The user function that allows the solver to check for
+         * convergence.
          */
         const std::function<SolverControl::State(const unsigned int i,
                                                  const double       f_norm,
@@ -814,8 +843,7 @@ namespace TrilinosWrappers
           check_iteration_status;
 
         /**
-         *  Last retured value of checkStatus(), which is used for
-         * getStatus() and print().
+         * The last returned value of `checkStatus()`.
          */
         NOX::StatusTest::StatusType status;
       };
@@ -941,6 +969,12 @@ namespace TrilinosWrappers
         // invert Jacobian
         if (solve_with_jacobian)
           {
+            Assert(
+              !solve_with_jacobian_and_track_n_linear_iterations,
+              ExcMessage(
+                "It does not make sense to provide both solve_with_jacobian and "
+                "solve_with_jacobian_and_track_n_linear_iterations!"));
+
             // without tracking of linear iterations
             return solve_with_jacobian(f, x, tolerance);
           }
@@ -979,8 +1013,9 @@ namespace TrilinosWrappers
 
     if (this->check_iteration_status)
       {
-        const auto info = Teuchos::rcp(
-          new internal::NOXWrappers::NOXCheck(this->check_iteration_status));
+        const auto info =
+          Teuchos::rcp(new internal::NOXWrappers::NOXCheck<VectorType>(
+            this->check_iteration_status));
         check->addStatusTest(info);
       }
 
diff --git a/source/trilinos/CMakeLists.txt b/source/trilinos/CMakeLists.txt
index 428dfa69f0..021ecd072e 100644
--- a/source/trilinos/CMakeLists.txt
+++ b/source/trilinos/CMakeLists.txt
@@ -13,19 +13,19 @@
 ##
 ## ---------------------------------------------------------------------
 
-INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
+include_directories(BEFORE ${CMAKE_CURRENT_BINARY_DIR})
 
-SET(_src
+set(_src
   nox.cc
   )
 
-SET(_inst
+set(_inst
   nox.inst.in
   )
 
-FILE(GLOB _header
+file(GLOB _header
   ${CMAKE_SOURCE_DIR}/include/deal.II/trilinos/*.h
   )
 
-DEAL_II_ADD_LIBRARY(obj_trilinos OBJECT ${_src} ${_header} ${_inst})
-EXPAND_INSTANTIATIONS(obj_trilinos "${_inst}")
+deal_ii_add_library(obj_trilinos OBJECT ${_src} ${_header} ${_inst})
+expand_instantiations(obj_trilinos "${_inst}")
diff --git a/source/trilinos/nox.cc b/source/trilinos/nox.cc
index 52ac18e771..9a1db9e186 100644
--- a/source/trilinos/nox.cc
+++ b/source/trilinos/nox.cc
@@ -29,7 +29,7 @@
 
 #include <deal.II/trilinos/nox.templates.h>
 
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
 
 DEAL_II_NAMESPACE_OPEN
 
-- 
2.39.5