From f012c46d7ab9f360b6cc6eb212914b806e575fde Mon Sep 17 00:00:00 2001 From: Peter Munch 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 -#ifdef DEAL_II_WITH_TRILINOS +#ifdef DEAL_II_TRILINOS_WITH_NOX # include @@ -26,6 +26,8 @@ # include +# include + 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::AdditionalData additional_data; + * + * // create nonlinear solver + * TrilinosWrappers::NOXSolver 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 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 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 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 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 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 update_preconditioner_predicate; @@ -222,22 +266,22 @@ namespace TrilinosWrappers const Teuchos::RCP 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 -#ifdef DEAL_II_WITH_TRILINOS +#ifdef DEAL_II_TRILINOS_WITH_NOX # include # include @@ -31,6 +31,8 @@ # include # include +# include + DEAL_II_NAMESPACE_OPEN # ifndef DOXYGEN @@ -73,7 +75,7 @@ namespace TrilinosWrappers NOX::Abstract::Vector & init(double gamma) override { - *vector = gamma; + *vector = static_cast(gamma); return *this; } @@ -143,7 +145,7 @@ namespace TrilinosWrappers NOX::Abstract::Vector & scale(double gamma) override { - *vector *= gamma; + *vector *= static_cast(gamma); return *this; } @@ -175,7 +177,9 @@ namespace TrilinosWrappers Assert(a_, ExcInternalError()); - vector->sadd(gamma, alpha, *a_->vector); + vector->sadd(static_cast(gamma), + static_cast(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(gamma)); + vector->add(static_cast(alpha), + *a_->vector, + static_cast(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 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 getXPtr() const override @@ -540,7 +548,7 @@ namespace TrilinosWrappers } /** - * Return RCP to F(x). + * Return an RCP to the residual `F(x)`. */ Teuchos::RCP getFPtr() const override @@ -550,7 +558,7 @@ namespace TrilinosWrappers } /** - * Return RCP to gradient. + * Return an `RCP` to gradient. */ Teuchos::RCP getGradientPtr() const override @@ -560,7 +568,7 @@ namespace TrilinosWrappers } /** - * Return RCP to Newton direction. + * Return an `RCP` to the Newton descent direction. */ Teuchos::RCP 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 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 x, f, gradient, newton; + // internal vector for the current solution + Vector x; + + // internal vector for the residual + Vector f; + + // internal vector for the solution gradient + Vector gradient; + + // internal vector for the newton step + Vector newton; + + // helper function to compute residual + std::function residual; + + // helper function to setup Jacobian + std::function setup_jacobian; - // helper functions to compute residual, to setup jacobian, and - // solve jacobian - std::function residual; - std::function setup_jacobian; - std::function apply_jacobian; - std::function + // helper function to apply Jacobian + std::function 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::functioncheck_iteration_status) { - const auto info = Teuchos::rcp( - new internal::NOXWrappers::NOXCheck(this->check_iteration_status)); + const auto info = + Teuchos::rcp(new internal::NOXWrappers::NOXCheck( + 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 -#ifdef DEAL_II_WITH_TRILINOS +#ifdef DEAL_II_TRILINOS_WITH_NOX DEAL_II_NAMESPACE_OPEN -- 2.39.5