From ecd063bff53fa0d4fb30836b13adc1779d4c430c Mon Sep 17 00:00:00 2001 From: Peter Munch Date: Sun, 2 Oct 2022 20:30:32 +0200 Subject: [PATCH] Add NOX wrapper --- doc/news/changes/major/20221002Munch | 4 + include/deal.II/trilinos/nox.h | 250 ++++ include/deal.II/trilinos/nox.templates.h | 1027 +++++++++++++++++ source/CMakeLists.txt | 1 + source/trilinos/CMakeLists.txt | 31 + source/trilinos/nox.cc | 43 + source/trilinos/nox.inst.in | 19 + tests/trilinos/nox_solver_01.cc | 222 ++++ .../nox_solver_01.with_trilinos=true.output | 3 + 9 files changed, 1600 insertions(+) create mode 100644 doc/news/changes/major/20221002Munch create mode 100644 include/deal.II/trilinos/nox.h create mode 100644 include/deal.II/trilinos/nox.templates.h create mode 100644 source/trilinos/CMakeLists.txt create mode 100644 source/trilinos/nox.cc create mode 100644 source/trilinos/nox.inst.in create mode 100644 tests/trilinos/nox_solver_01.cc create mode 100644 tests/trilinos/nox_solver_01.with_trilinos=true.output diff --git a/doc/news/changes/major/20221002Munch b/doc/news/changes/major/20221002Munch new file mode 100644 index 0000000000..194171c10c --- /dev/null +++ b/doc/news/changes/major/20221002Munch @@ -0,0 +1,4 @@ +New: Add a wrapper to the non-linear solver NOX from the +Trilinos library. +
+(Peter Munch, Vladimir Ivannikov, 2022/10/02) diff --git a/include/deal.II/trilinos/nox.h b/include/deal.II/trilinos/nox.h new file mode 100644 index 0000000000..16b323de49 --- /dev/null +++ b/include/deal.II/trilinos/nox.h @@ -0,0 +1,250 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_nox +#define dealii_trilinos_nox + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +# include + +# include + +# include + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ + // Indicate that NOXSolver has not converged. + DeclException0(ExcNOXNoConvergence); + + + /** + * Wrapper around the non-linear solver from the NOX + * packge (https://docs.trilinos.org/dev/packages/nox/doc/html/index.html), + * targeting deal.II data structures. + */ + template + class NOXSolver + { + public: + /** + * Struct that helps to configure NOXSolver. More advanced + * parameters are passed to the constructor NOXSolver + * directly via a Teuchos::ParameterList. + */ + struct AdditionalData + { + public: + /** + * Constructor. + */ + AdditionalData(const unsigned int max_iter = 10, + const double abs_tol = 1.e-20, + const double rel_tol = 1.e-5, + const unsigned int threshold_nonlinear_iterations = 1, + const unsigned int threshold_n_linear_iterations = 0, + const bool reuse_solver = false); + + /** + * Max number of non-linear iterations. + */ + unsigned int max_iter; + + /** + * Absolute l2 tolerance to be reached. + */ + double abs_tol; + + /** + * Relative l2 tolerance to be reached. + */ + double rel_tol; + + /** + * Number of non-linear iterations after which the preconditioner + * should be updated. + */ + unsigned int threshold_nonlinear_iterations; + + /** + * Max number of linear iterations after which the preconditioner + * should be updated. This is only used if a lambda is attached to + * solve_with_jacobian_and_track_n_linear_iterations. + */ + unsigned int threshold_n_linear_iterations; + + /** + * Reuse NOX solver instance in the next non-linear solution. + */ + bool reuse_solver; + }; + + /** + * Constructor. + * + * If @p parameters is not filled, a Newton solver with full step is used. + * An overview of possible parameters is given at + * https://docs.trilinos.org/dev/packages/nox/doc/html/parameters.html. + */ + NOXSolver(AdditionalData & additional_data, + const Teuchos::RCP ¶meters = + Teuchos::rcp(new Teuchos::ParameterList)); + + /** + * Clear internal state. + */ + void + clear(); + + /** + * Solve non-linear problem and return number of iterations. + */ + unsigned int + solve(VectorType &solution); + + /** + * User function that computes the residual. + * + * @note This function should return 0 in the case of success. + */ + std::function residual; + + /** + * User function that sets up the Jacobian. + * + * @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. + * + * @note The 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). + * + * @note This function should return 0 in the case of success. + */ + std::function setup_preconditioner; + + /** + * User function that applies the Jacobian. + * + * @note The 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) + * it won't be called. + * + * @note This function should return 0 in the case of success. + */ + std::function apply_jacobian; + + /** + * User function that applies the inverse of the Jacobian. + * + * @note The function is optional and is used in the case of certain + * configurations. + * + * @note This function should return 0 in the case of success. + */ + std::function< + int(const VectorType &f, VectorType &x, const double tolerance)> + solve_with_jacobian; + + /** + * User function that applies the inverse of the Jacobian and + * returns the numer of linear iterations the linear solver needed. + * + * @note This function should return -1 in the case of failure. + */ + 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 + * ones checking the l2-norm and the number of iterations (see + * AdditionalData). It is run after each non-linear 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. + */ + std::function + check_iteration_status; + + /** + * Function that allows to force to update the preconditioner in + * addition to AdditionalData::threshold_nonlinear_iterations. A reason + * for wanting to update the preconditioner is when the expected number + * of linear iterations exceeds. + * + * @note The function is optional. If no function is attached, this + * means implicitly a return value of false. + */ + std::function update_preconditioner_predicate; + + private: + /** + * Additional data with basic settings. + */ + AdditionalData additional_data; + + /** + * Additional data with advanced settings. An overview of + * possible parameters is given at + * https://docs.trilinos.org/dev/packages/nox/doc/html/parameters.html. + */ + const Teuchos::RCP parameters; + + /** + * Counter for number of (accumulated) residual evaluations. + */ + unsigned int n_residual_evaluations; + + /** + * Counter for number of (accumulated) Jacobi applications. + */ + unsigned int n_jacobian_applications; + + /** + * Counter for number of (accumulated) non-linear iterations. + */ + unsigned int n_nonlinear_iterations; + + /** + * Number of linear iterations of the last Jacobian solve. + */ + unsigned int n_last_linear_iterations; + }; +} // namespace TrilinosWrappers + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/include/deal.II/trilinos/nox.templates.h b/include/deal.II/trilinos/nox.templates.h new file mode 100644 index 0000000000..20c23eadfb --- /dev/null +++ b/include/deal.II/trilinos/nox.templates.h @@ -0,0 +1,1027 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +#ifndef dealii_trilinos_nox_templates +#define dealii_trilinos_nox_templates + +#include + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +# include +# include +# include +# include +# include +# include +# include +# include + +DEAL_II_NAMESPACE_OPEN + +# ifndef DOXYGEN + +namespace TrilinosWrappers +{ + namespace internal + { + namespace NOXWrappers + { + template + class Group; + + /** + * Implementation of the abstract interface + * NOX::Abstract::Vector for deal.II vectors. For details, + * see + * https://docs.trilinos.org/dev/packages/nox/doc/html/classNOX_1_1Abstract_1_1Vector.html. + */ + template + class Vector : public NOX::Abstract::Vector + { + public: + /** + * Create empty vector. + */ + Vector() = default; + + /** + * Wrap an existing vector. The ownership is not transferred. + */ + Vector(VectorType &vector) + { + this->vector.reset(&vector, [](auto *) { /*nothing to do*/ }); + } + + /** + * Initialize every element of this vector with gamma. + */ + NOX::Abstract::Vector & + init(double gamma) override + { + *vector = gamma; + return *this; + } + + /** + * Initialize each element of this vector with a random value. + */ + NOX::Abstract::Vector & + random(bool useSeed = false, int seed = 1) override + { + AssertThrow(false, ExcNotImplemented()); + + (void)useSeed; + (void)seed; + + return *this; + } + + /** + * Put element-wise absolute values of source vector y into this vector. + */ + NOX::Abstract::Vector & + abs(const NOX::Abstract::Vector &y) override + { + AssertThrow(false, ExcNotImplemented()); + + (void)y; + + return *this; + } + + /** + * Copy source vector y into this vector. + */ + NOX::Abstract::Vector & + operator=(const NOX::Abstract::Vector &y) override + { + if (vector == nullptr) + vector = std::shared_ptr(); + + const auto y_ = dynamic_cast *>(&y); + + Assert(y_, ExcInternalError()); + + vector->reinit(*y_->vector); + + *vector = *y_->vector; + + return *this; + } + + /** + * Put element-wise reciprocal of source vector y into this vector. + */ + NOX::Abstract::Vector & + reciprocal(const NOX::Abstract::Vector &y) override + { + AssertThrow(false, ExcNotImplemented()); + + (void)y; + + return *this; + } + + /** + * Scale each element of this vector by gamma. + */ + NOX::Abstract::Vector & + scale(double gamma) override + { + *vector *= gamma; + + return *this; + } + + /** + * Scale this vector element-by-element by the vector a. + */ + NOX::Abstract::Vector & + scale(const NOX::Abstract::Vector &a) override + { + const auto a_ = dynamic_cast *>(&a); + + Assert(a_, ExcInternalError()); + + vector->scale(*a_->vector); + + return *this; + } + + /** + * Compute x = (alpha * a) + (gamma * x) where x is this vector. + */ + NOX::Abstract::Vector & + update(double alpha, + const NOX::Abstract::Vector &a, + double gamma = 0.0) override + { + const auto a_ = dynamic_cast *>(&a); + + Assert(a_, ExcInternalError()); + + vector->sadd(gamma, alpha, *a_->vector); + + return *this; + } + + /** + * Compute x = (alpha * a) + (beta * b) + (gamma * x) where x is this + * vector. + */ + NOX::Abstract::Vector & + update(double alpha, + const NOX::Abstract::Vector &a, + double beta, + const NOX::Abstract::Vector &b, + double gamma = 0.0) override + { + const auto a_ = dynamic_cast *>(&a); + const auto b_ = dynamic_cast *>(&b); + + Assert(a_, ExcInternalError()); + Assert(b_, ExcInternalError()); + + vector->operator*=(gamma); + vector->add(alpha, *a_->vector, beta, *b_->vector); + + return *this; + } + + /** + * Create a new Vector of the same underlying type by cloning "this", + * and return a pointer to the new vector. + */ + Teuchos::RCP + clone(NOX::CopyType copy_type) const override + { + auto new_vector = Teuchos::rcp(new Vector()); + new_vector->vector = std::make_shared(); + new_vector->vector->reinit(*this->vector); + + if (copy_type == NOX::CopyType::DeepCopy) + *new_vector->vector = *this->vector; + else + Assert(copy_type == NOX::CopyType::ShapeCopy, ExcInternalError()); + + return new_vector; + } + + /** + * Norm. + */ + double + norm(NOX::Abstract::Vector::NormType type = + NOX::Abstract::Vector::TwoNorm) const override + { + if (type == NOX::Abstract::Vector::NormType::TwoNorm) + return vector->l2_norm(); + if (type == NOX::Abstract::Vector::NormType::OneNorm) + return vector->l1_norm(); + if (type == NOX::Abstract::Vector::NormType::MaxNorm) + return vector->linfty_norm(); + + Assert(false, ExcInternalError()); + + return 0.0; + } + + /** + * Weighted 2-Norm. + */ + double + norm(const NOX::Abstract::Vector &weights) const override + { + AssertThrow(false, ExcNotImplemented()); + + (void)weights; + + return 0.0; + } + + /** + * Inner product with y. + */ + double + innerProduct(const NOX::Abstract::Vector &y) const override + { + const auto y_ = dynamic_cast *>(&y); + + Assert(y_, ExcInternalError()); + + return (*vector) * (*y_->vector); + } + + /** + * Return the length of vector. + */ + NOX::size_type + length() const override + { + return vector->size(); + } + + /** + * Return underlying vector. + */ + operator VectorType &() const + { + AssertThrow(vector, ExcInternalError()); + + return *vector; + } + + private: + /** + * Underlying deal.II vector. + */ + std::shared_ptr vector; + + friend Group; + }; + + /** + * Implementation of the abstract interface + * NOX::Abstract::Group for deal.II vectors and deal.II solvers. For + * details, see + * https://docs.trilinos.org/dev/packages/nox/doc/html/classNOX_1_1Abstract_1_1Group.html. + */ + template + class Group : public NOX::Abstract::Group + { + public: + /** + * Constructor. The class is intialized by the solution vector and + * functions to compute the residual, to setup the jacobian, and + * to solve the Jacobian. + */ + Group( + VectorType & solution, + const std::function &residual, + const std::function &setup_jacobian, + const std::function + & apply_jacobian, + const std::function &solve_with_jacobian) + : x(solution) + , residual(residual) + , setup_jacobian(setup_jacobian) + , apply_jacobian(apply_jacobian) + , solve_with_jacobian(solve_with_jacobian) + , is_valid_f(false) + , is_valid_j(false) + {} + + /** + * Copies the source group into this group. + */ + NOX::Abstract::Group & + operator=(const NOX::Abstract::Group &source) override + { + if (this != &source) + { + const auto other = + dynamic_cast *>(&source); + + Assert(other, ExcInternalError()); + + if (other->x.vector) + { + if (this->x.vector == nullptr) + this->x.vector = std::make_shared(); + + *this->x.vector = *other->x.vector; + } + else + { + this->x.vector = {}; + } + + if (other->f.vector) + { + if (this->f.vector == nullptr) + this->f.vector = std::make_shared(); + + *this->f.vector = *other->f.vector; + } + else + { + this->f.vector = {}; + } + + if (other->gradient.vector) + { + if (this->gradient.vector == nullptr) + this->gradient.vector = std::make_shared(); + + *this->gradient.vector = *other->gradient.vector; + } + else + { + this->gradient.vector = {}; + } + + if (other->newton.vector) + { + if (this->newton.vector == nullptr) + this->newton.vector = std::make_shared(); + + *this->newton.vector = *other->newton.vector; + } + else + { + this->newton.vector = {}; + } + + this->residual = other->residual; + this->setup_jacobian = other->setup_jacobian; + this->apply_jacobian = other->apply_jacobian; + this->solve_with_jacobian = other->solve_with_jacobian; + + this->is_valid_f = other->is_valid_f; + this->is_valid_j = other->is_valid_j; + } + + return *this; + } + + /** + * Set the solution vector x to y. + */ + void + setX(const NOX::Abstract::Vector &y) override + { + reset(); + + x = y; + } + + /** + * Compute x = grp.x + step * d. + */ + void + computeX(const NOX::Abstract::Group & grp, + const NOX::Abstract::Vector &d, + double step) override + { + reset(); + + const auto grp_ = dynamic_cast(&grp); + + Assert(grp_, ExcInternalError()); + + x.update(1.0, grp_->x, step, d); + } + + /** + * Compute and store F(x). + */ + NOX::Abstract::Group::ReturnType + computeF() override + { + if (isF() == false) + { + f.vector = std::make_shared(); + f.vector->reinit(*x.vector); + + if (residual(*x.vector, *f.vector) != 0) + return NOX::Abstract::Group::Failed; + + is_valid_f = true; + } + + return NOX::Abstract::Group::Ok; + } + + /** + * Return true if F is valid. + */ + bool + isF() const override + { + return is_valid_f; + } + + /** + * Compute and store Jacobian. + */ + NOX::Abstract::Group::ReturnType + computeJacobian() override + { + if (isJacobian() == false) + { + if (setup_jacobian(*x.vector) != 0) + return NOX::Abstract::Group::Failed; + + is_valid_j = true; + } + + return NOX::Abstract::Group::Ok; + } + + /** + * Return true if the Jacobian is valid. + */ + bool + isJacobian() const override + { + return is_valid_j; + } + + /** + * Return solution vector. + */ + const NOX::Abstract::Vector & + getX() const override + { + return x; + } + + /** + * Return F(x). + */ + const NOX::Abstract::Vector & + getF() const override + { + return f; + } + + /** + * Return 2-norm of F(x) + */ + double + getNormF() const override + { + return f.norm(); + } + + /** + * Return gradient. + */ + const NOX::Abstract::Vector & + getGradient() const override + { + return gradient; + } + + /** + * Return Newton direction. + */ + const NOX::Abstract::Vector & + getNewton() const override + { + return newton; + } + + /** + * Return RCP to solution vector. + */ + Teuchos::RCP + getXPtr() const override + { + AssertThrow(false, ExcNotImplemented()); + return {}; + } + + /** + * Return RCP to F(x). + */ + Teuchos::RCP + getFPtr() const override + { + AssertThrow(false, ExcNotImplemented()); + return {}; + } + + /** + * Return RCP to gradient. + */ + Teuchos::RCP + getGradientPtr() const override + { + AssertThrow(false, ExcNotImplemented()); + return {}; + } + + /** + * Return RCP to Newton direction. + */ + Teuchos::RCP + getNewtonPtr() const override + { + AssertThrow(false, ExcNotImplemented()); + return {}; + } + + /** + * 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. + */ + Teuchos::RCP + clone(NOX::CopyType copy_type) const override + { + auto new_group = + Teuchos::rcp(new Group(*x.vector, + residual, + setup_jacobian, + apply_jacobian, + solve_with_jacobian)); + + if (x.vector) + { + new_group->x.vector = std::make_shared(); + new_group->x.vector->reinit(*x.vector); + } + + if (f.vector) + { + new_group->f.vector = std::make_shared(); + new_group->f.vector->reinit(*f.vector); + } + + if (gradient.vector) + { + new_group->gradient.vector = std::make_shared(); + new_group->gradient.vector->reinit(*gradient.vector); + } + + if (newton.vector) + { + new_group->newton.vector = std::make_shared(); + new_group->newton.vector->reinit(*newton.vector); + } + + if (copy_type == NOX::CopyType::DeepCopy) + { + if (x.vector) + *new_group->x.vector = *x.vector; + + if (f.vector) + *new_group->f.vector = *f.vector; + + if (gradient.vector) + *new_group->gradient.vector = *gradient.vector; + + if (newton.vector) + *new_group->newton.vector = *newton.vector; + + new_group->is_valid_f = is_valid_f; + new_group->is_valid_j = is_valid_j; + } + else + Assert(copy_type == NOX::CopyType::ShapeCopy, ExcInternalError()); + + return new_group; + } + + /** + * Compute the Newton direction, using parameters for the linear solve. + */ + NOX::Abstract::Group::ReturnType + computeNewton(Teuchos::ParameterList &p) override + { + if (isNewton()) + return NOX::Abstract::Group::Ok; + + if (isF() == false || isJacobian() == false) + return NOX::Abstract::Group::BadDependency; + + if (newton.vector == nullptr) + newton.vector = std::make_shared(); + + newton.vector->reinit(*f.vector, false); + + const double tolerance = p.get("Tolerance"); + + if (solve_with_jacobian(*f.vector, *newton.vector, tolerance) != 0) + return NOX::Abstract::Group::NotConverged; + + newton.scale(-1.0); + + return NOX::Abstract::Group::Ok; + } + + /** + * Applies Jacobian to the given input vector and puts + * the answer in the result. + */ + NOX::Abstract::Group::ReturnType + applyJacobian(const NOX::Abstract::Vector &input, + NOX::Abstract::Vector & result) const override + { + if (apply_jacobian == nullptr) + return NOX::Abstract::Group::NotDefined; + + if (!isJacobian()) + return NOX::Abstract::Group::BadDependency; + + const auto *input_ = dynamic_cast *>(&input); + const auto *result_ = + dynamic_cast *>(&result); + + if (apply_jacobian(*input_->vector, *result_->vector) != 0) + return NOX::Abstract::Group::Failed; + + return NOX::Abstract::Group::Ok; + } + + private: + /** + * Reset state. + */ + void + reset() + { + is_valid_f = false; + is_valid_j = false; + } + + // internal vectors + Vector x, f, gradient, newton; + + // helper functions to compute residual, to setup jacobian, and + // solve jacobian + std::function residual; + std::function setup_jacobian; + std::function apply_jacobian; + std::function + solve_with_jacobian; + + // internal state (are residuum and jacobian computed?) + bool is_valid_f, is_valid_j; + }; + + + /** + * Wrapper class around the user function that allows to check + * convergence. + */ + template + class NOXCheck : public NOX::StatusTest::Generic + { + public: + /** + * Constructor. + */ + NOXCheck(const std::function + check_iteration_status) + : check_iteration_status(check_iteration_status) + , status(NOX::StatusTest::Unevaluated) + {} + + /** + * Check status. + */ + NOX::StatusTest::StatusType + checkStatus(const NOX::Solver::Generic &problem, + NOX::StatusTest::CheckType checkType) override + { + if (checkType == NOX::StatusTest::None) + { + status = NOX::StatusTest::Unevaluated; + } + else + { + if (check_iteration_status == nullptr) + { + status = NOX::StatusTest::Converged; + } + else + { + // unwrap the various vectors + const VectorType &x__ = *dynamic_cast< + const internal::NOXWrappers::Vector *>( + &problem.getSolutionGroup().getX()); + const VectorType &f__ = *dynamic_cast< + const internal::NOXWrappers::Vector *>( + &problem.getSolutionGroup().getF()); + + // forward to the user-provided function and checks + // convergence + const auto state = this->check_iteration_status( + problem.getNumIterations(), f__.l2_norm(), x__, f__); + + // translate the returned value back to Trilinos data + // structure + switch (state) + { + case SolverControl::iterate: + status = NOX::StatusTest::Unconverged; + break; + case SolverControl::failure: + status = NOX::StatusTest::Failed; + break; + case SolverControl::success: + status = NOX::StatusTest::Converged; + break; + default: + AssertThrow(false, ExcNotImplemented()); + } + } + } + + return status; + } + + /** + * Return last return value of checkStatus(). + */ + NOX::StatusTest::StatusType + getStatus() const override + { + return status; + } + + /** + * Print last return value of print(). + */ + virtual std::ostream & + print(std::ostream &stream, int indent = 0) const override + { + for (int j = 0; j < indent; ++j) + stream << ' '; + stream << status << std::endl; + return stream; + } + + private: + /** + * User function that allows to check convergence. + */ + const std::function + check_iteration_status; + + /** + * Last retured value of checkStatus(), which is used for + * getStatus() and print(). + */ + NOX::StatusTest::StatusType status; + }; + } // namespace NOXWrappers + } // namespace internal + + + + template + NOXSolver::AdditionalData::AdditionalData( + const unsigned int max_iter, + const double abs_tol, + const double rel_tol, + const unsigned int threshold_nonlinear_iterations, + const unsigned int threshold_n_linear_iterations, + const bool reuse_solver) + : max_iter(max_iter) + , abs_tol(abs_tol) + , rel_tol(rel_tol) + , threshold_nonlinear_iterations(threshold_nonlinear_iterations) + , threshold_n_linear_iterations(threshold_n_linear_iterations) + , reuse_solver(reuse_solver) + {} + + + + template + NOXSolver::NOXSolver( + AdditionalData & additional_data, + const Teuchos::RCP ¶meters) + : additional_data(additional_data) + , parameters(parameters) + , n_residual_evaluations(0) + , n_jacobian_applications(0) + , n_nonlinear_iterations(0) + , n_last_linear_iterations(0) + {} + + + + template + void + NOXSolver::clear() + { + // clear interal counters + n_residual_evaluations = 0; + n_jacobian_applications = 0; + n_nonlinear_iterations = 0; + n_last_linear_iterations = 0; + } + + + + template + unsigned int + NOXSolver::solve(VectorType &solution) + { + if (additional_data.reuse_solver == false) + clear(); // clear state + + // create group + const auto group = Teuchos::rcp(new internal::NOXWrappers::Group< + VectorType>( + solution, + [&](const VectorType &x, VectorType &f) -> int { + Assert( + residual, + ExcMessage( + "No residual function has been attached to the NOXSolver object.")); + + n_residual_evaluations++; + + // evalute residual + return residual(x, f); + }, + [&](const VectorType &x) -> int { + Assert( + setup_jacobian, + ExcMessage( + "No setup_jacobian function has been attached to the NOXSolver object.")); + + // setup Jacobian + int flag = setup_jacobian(x); + + if (flag != 0) + return flag; + + if (setup_preconditioner) + { + // check if preconditioner needs to be updated + bool update_preconditioner = + ((additional_data.threshold_nonlinear_iterations > 0) && + ((n_nonlinear_iterations % + additional_data.threshold_nonlinear_iterations) == 0)) || + (solve_with_jacobian_and_track_n_linear_iterations && + (n_last_linear_iterations > + additional_data.threshold_n_linear_iterations)); + + if ((update_preconditioner == false) && + (update_preconditioner_predicate != nullptr)) + update_preconditioner = update_preconditioner_predicate(); + + if (update_preconditioner) // update preconditioner + flag = setup_preconditioner(x); + } + + return flag; + }, + [&](const VectorType &x, VectorType &v) -> int { + Assert( + apply_jacobian, + ExcMessage( + "No apply_jacobian function has been attached to the NOXSolver object.")); + + n_jacobian_applications++; + + // apply Jacobian + return apply_jacobian(x, v); + }, + [&](const VectorType &f, VectorType &x, const double tolerance) -> int { + n_nonlinear_iterations++; + + // invert Jacobian + if (solve_with_jacobian) + { + // without tracking of linear iterations + return solve_with_jacobian(f, x, tolerance); + } + else if (solve_with_jacobian_and_track_n_linear_iterations) + { + // with tracking of linear iterations + const int n_linear_iterations = + solve_with_jacobian_and_track_n_linear_iterations(f, + x, + tolerance); + + if (n_linear_iterations == -1) + return 1; + + this->n_last_linear_iterations = n_linear_iterations; + + return 0; + } + else + { + Assert( + false, + ExcMessage( + "Neither a solve_with_jacobian or a " + "solve_with_jacobian_and_track_n_linear_iterations function " + "has been attached to the NOXSolver object.")); + + Assert(false, ExcNotImplemented()); + return 1; + } + })); + + // setup solver control + auto check = + Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR)); + + if (this->check_iteration_status) + { + const auto info = Teuchos::rcp( + new internal::NOXWrappers::NOXCheck(this->check_iteration_status)); + check->addStatusTest(info); + } + + if (additional_data.abs_tol > 0.0) + { + const auto additional_data_norm_f_abs = + Teuchos::rcp(new NOX::StatusTest::NormF(additional_data.abs_tol)); + check->addStatusTest(additional_data_norm_f_abs); + } + + if (additional_data.rel_tol > 0.0) + { + const auto additional_data_norm_f_rel = Teuchos::rcp( + new NOX::StatusTest::RelativeNormF(additional_data.rel_tol)); + check->addStatusTest(additional_data_norm_f_rel); + } + + if (additional_data.max_iter > 0) + { + const auto additional_data_max_iterations = + Teuchos::rcp(new NOX::StatusTest::MaxIters(additional_data.max_iter)); + check->addStatusTest(additional_data_max_iterations); + } + + // create non-linear solver + const auto solver = NOX::Solver::buildSolver(group, check, parameters); + + // solve + const auto status = solver->solve(); + + AssertThrow(status == NOX::StatusTest::Converged, ExcNOXNoConvergence()); + + return solver->getNumIterations(); + } + +} // namespace TrilinosWrappers + +# endif + +DEAL_II_NAMESPACE_CLOSE + +#endif + +#endif diff --git a/source/CMakeLists.txt b/source/CMakeLists.txt index 83a53336b6..ce9b1ff232 100644 --- a/source/CMakeLists.txt +++ b/source/CMakeLists.txt @@ -73,6 +73,7 @@ add_subdirectory(physics) add_subdirectory(optimization/rol) add_subdirectory(non_matching) add_subdirectory(sundials) +add_subdirectory(trilinos) add_subdirectory(arborx) foreach(build ${DEAL_II_BUILD_TYPES}) diff --git a/source/trilinos/CMakeLists.txt b/source/trilinos/CMakeLists.txt new file mode 100644 index 0000000000..428dfa69f0 --- /dev/null +++ b/source/trilinos/CMakeLists.txt @@ -0,0 +1,31 @@ +## --------------------------------------------------------------------- +## +## Copyright (C) 2012 - 2021 by the deal.II authors +## +## This file is part of the deal.II library. +## +## The deal.II library is free software; you can use it, redistribute +## it, and/or modify it under the terms of the GNU Lesser General +## Public License as published by the Free Software Foundation; either +## version 2.1 of the License, or (at your option) any later version. +## The full text of the license can be found in the file LICENSE.md at +## the top level directory of deal.II. +## +## --------------------------------------------------------------------- + +INCLUDE_DIRECTORIES(BEFORE ${CMAKE_CURRENT_BINARY_DIR}) + +SET(_src + nox.cc + ) + +SET(_inst + nox.inst.in + ) + +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}") diff --git a/source/trilinos/nox.cc b/source/trilinos/nox.cc new file mode 100644 index 0000000000..52ac18e771 --- /dev/null +++ b/source/trilinos/nox.cc @@ -0,0 +1,43 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#ifdef DEAL_II_WITH_TRILINOS + +DEAL_II_NAMESPACE_OPEN + +namespace TrilinosWrappers +{ +# include "nox.inst" +} + +DEAL_II_NAMESPACE_CLOSE + +#endif diff --git a/source/trilinos/nox.inst.in b/source/trilinos/nox.inst.in new file mode 100644 index 0000000000..3509448dcf --- /dev/null +++ b/source/trilinos/nox.inst.in @@ -0,0 +1,19 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +for (VEC : REAL_VECTOR_TYPES) + { + template class NOXSolver; + } diff --git a/tests/trilinos/nox_solver_01.cc b/tests/trilinos/nox_solver_01.cc new file mode 100644 index 0000000000..2c5e598198 --- /dev/null +++ b/tests/trilinos/nox_solver_01.cc @@ -0,0 +1,222 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2022 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + + +// Check TrilinosWrappers::NOXSolver by solving f(x) = x^2 with initial +// condition x=2. + +#include + +#include + +#include + +// as reference solution +#include +#include +#include +#include + +#include "../tests.h" + +#include "NOX_Epetra_Group.H" +#include "NOX_Epetra_Interface_Jacobian.H" +#include "NOX_Epetra_Interface_Required.H" +#include "NOX_Epetra_LinearSystem_AztecOO.H" + + +class NoxInterface : public NOX::Epetra::Interface::Required, + public NOX::Epetra::Interface::Jacobian +{ +public: + bool + computeF(const Epetra_Vector & x, + Epetra_Vector & f, + NOX::Epetra::Interface::Required::FillType F) override + { + (void)F; + + f[0] = x[0] * x[0]; + + return true; + } + + bool + computeJacobian(const Epetra_Vector &x, Epetra_Operator &Jac) override + { + auto jac = dynamic_cast(&Jac); + + AssertThrow(jac, ExcNotImplemented()); + + jac->PutScalar(2.0 * x[0]); + + return true; + } +}; + + + +int +main(int argc, char **argv) +{ + initlog(); + + Utilities::MPI::MPI_InitFinalize mpi_init(argc, argv, 1); + + using Number = double; + using VectorType = LinearAlgebra::distributed::Vector; + + // set up solver control + const unsigned int n_max_iterations = 100; + const double abs_tolerance = 1e-9; + const double rel_tolerance = 1e-5; + const double lin_rel_tolerance = 1e-3; + + TrilinosWrappers::NOXSolver::AdditionalData additional_data( + n_max_iterations, abs_tolerance, rel_tolerance); + + // set up parameters + Teuchos::RCP non_linear_parameters = + Teuchos::rcp(new Teuchos::ParameterList); + + non_linear_parameters->set("Nonlinear Solver", "Line Search Based"); + non_linear_parameters->sublist("Printing").set("Output Information", 15); + non_linear_parameters->sublist("Direction").set("Method", "Newton"); + non_linear_parameters->sublist("Direction") + .sublist("Newton") + .sublist("Linear Solver") + .set("Tolerance", lin_rel_tolerance); + non_linear_parameters->sublist("Line Search").set("Method", "Polynomial"); + + if (true) + { + // set up solver + TrilinosWrappers::NOXSolver solver(additional_data, + non_linear_parameters); + + // ... helper functions + double J = 0.0; + + solver.residual = [](const auto &src, auto &dst) { + // compute residual + dst[0] = src[0] * src[0]; + return 0; + }; + + solver.setup_jacobian = [&](const auto &src) { + // compute Jacobian + J = 2.0 * src[0]; + return 0; + }; + + solver.apply_jacobian = [&](const auto &src, auto &dst) { + // solve with Jacobian + dst[0] = src[0] * J; + return 0; + }; + + solver.solve_with_jacobian = [&](const auto &src, auto &dst, const auto) { + // solve with Jacobian + dst[0] = src[0] / J; + return 0; + }; + + // initial guess + VectorType solution(1); + solution[0] = 2.0; + + // solve with the given initial guess + solver.solve(solution); + + deallog << "The solution is: " << solution[0] << std::endl; + } + + if (true) + { + // convert data structures to Epetra structures + IndexSet is(1); + is.add_index(0); + + TrilinosWrappers::MPI::Vector solution(is, MPI_COMM_WORLD); + solution[0] = 2.0; + + TrilinosWrappers::SparsityPattern dsp(is, MPI_COMM_WORLD); + dsp.add(0, 0); + dsp.compress(); + + TrilinosWrappers::SparseMatrix system_matrix; + system_matrix.reinit(dsp); + + + // setup linear system and group + + Epetra_Vector InitialGuess(Copy, solution.trilinos_vector(), 0); + + auto A = + Teuchos::rcp(new Epetra_CrsMatrix(system_matrix.trilinos_matrix())); + + Teuchos::ParameterList &printParams = + non_linear_parameters->sublist("Printing"); + Teuchos::ParameterList &lsParams = + non_linear_parameters->sublist("Newton").sublist("Linear Solver"); + + Teuchos::RCP interface = Teuchos::rcp(new NoxInterface()); + + Teuchos::RCP iReq = interface; + Teuchos::RCP iJac = interface; + + Teuchos::RCP linSys = + Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO( + printParams, lsParams, iReq, iJac, A, InitialGuess)); + + NOX::Epetra::Vector noxInitGuess(InitialGuess, NOX::DeepCopy); + Teuchos::RCP group = Teuchos::rcp( + new NOX::Epetra::Group(printParams, iReq, noxInitGuess, linSys)); + + // setup solver control + auto check = + Teuchos::rcp(new NOX::StatusTest::Combo(NOX::StatusTest::Combo::OR)); + + if (additional_data.abs_tol > 0.0) + { + const auto additional_data_norm_f_abs = + Teuchos::rcp(new NOX::StatusTest::NormF(additional_data.abs_tol)); + check->addStatusTest(additional_data_norm_f_abs); + } + + if (additional_data.rel_tol > 0.0) + { + const auto additional_data_norm_f_rel = Teuchos::rcp( + new NOX::StatusTest::RelativeNormF(additional_data.rel_tol)); + check->addStatusTest(additional_data_norm_f_rel); + } + + if (additional_data.max_iter > 0) + { + const auto additional_data_max_iterations = Teuchos::rcp( + new NOX::StatusTest::MaxIters(additional_data.max_iter)); + check->addStatusTest(additional_data_max_iterations); + } + + Teuchos::RCP solver = + NOX::Solver::buildSolver(group, check, non_linear_parameters); + auto status = solver->solve(); + + AssertThrow(status == NOX::StatusTest::Converged, ExcInternalError()); + + deallog << "The solution is: " << group->getX().norm() << std::endl; + } +} diff --git a/tests/trilinos/nox_solver_01.with_trilinos=true.output b/tests/trilinos/nox_solver_01.with_trilinos=true.output new file mode 100644 index 0000000000..b068dbf325 --- /dev/null +++ b/tests/trilinos/nox_solver_01.with_trilinos=true.output @@ -0,0 +1,3 @@ + +DEAL::The solution is: 0.00390625 +DEAL::The solution is: 0.00390625 -- 2.39.5