--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+# include <deal.II/base/exceptions.h>
+
+# include <deal.II/lac/solver_control.h>
+
+# include <Teuchos_ParameterList.hpp>
+
+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 <typename VectorType>
+ 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<Teuchos::ParameterList> ¶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<int(const VectorType &x, VectorType &f)> residual;
+
+ /**
+ * User function that sets up the Jacobian.
+ *
+ * @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.
+ *
+ * @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<int(const VectorType &x)> 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<int(const VectorType &x, VectorType &v)> 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<SolverControl::State(const unsigned int i,
+ const double norm_f,
+ const VectorType & x,
+ const VectorType & f)>
+ 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<bool()> 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<Teuchos::ParameterList> 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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/config.h>
+
+#include <deal.II/trilinos/nox.h>
+
+#ifdef DEAL_II_WITH_TRILINOS
+
+# include <NOX_Abstract_Group.H>
+# include <NOX_Abstract_Vector.H>
+# include <NOX_Solver_Factory.H>
+# include <NOX_Solver_Generic.H>
+# include <NOX_StatusTest_Combo.H>
+# include <NOX_StatusTest_MaxIters.H>
+# include <NOX_StatusTest_NormF.H>
+# include <NOX_StatusTest_RelativeNormF.H>
+
+DEAL_II_NAMESPACE_OPEN
+
+# ifndef DOXYGEN
+
+namespace TrilinosWrappers
+{
+ namespace internal
+ {
+ namespace NOXWrappers
+ {
+ template <typename VectorType>
+ 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 <typename VectorType>
+ 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<VectorType>();
+
+ const auto y_ = dynamic_cast<const Vector<VectorType> *>(&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<const Vector<VectorType> *>(&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<const Vector<VectorType> *>(&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<const Vector<VectorType> *>(&a);
+ const auto b_ = dynamic_cast<const Vector<VectorType> *>(&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<NOX::Abstract::Vector>
+ clone(NOX::CopyType copy_type) const override
+ {
+ auto new_vector = Teuchos::rcp(new Vector<VectorType>());
+ new_vector->vector = std::make_shared<VectorType>();
+ 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<const Vector<VectorType> *>(&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<VectorType> vector;
+
+ friend Group<VectorType>;
+ };
+
+ /**
+ * 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 <typename VectorType>
+ 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<int(const VectorType &, VectorType &)> &residual,
+ const std::function<int(const VectorType &)> &setup_jacobian,
+ const std::function<int(const VectorType &, VectorType &)>
+ & apply_jacobian,
+ const std::function<int(const VectorType &,
+ VectorType &,
+ const double)> &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<const Group<VectorType> *>(&source);
+
+ Assert(other, ExcInternalError());
+
+ if (other->x.vector)
+ {
+ if (this->x.vector == nullptr)
+ this->x.vector = std::make_shared<VectorType>();
+
+ *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<VectorType>();
+
+ *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<VectorType>();
+
+ *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<VectorType>();
+
+ *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<const Group *>(&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<VectorType>();
+ 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<const NOX::Abstract::Vector>
+ getXPtr() const override
+ {
+ AssertThrow(false, ExcNotImplemented());
+ return {};
+ }
+
+ /**
+ * Return RCP to F(x).
+ */
+ Teuchos::RCP<const NOX::Abstract::Vector>
+ getFPtr() const override
+ {
+ AssertThrow(false, ExcNotImplemented());
+ return {};
+ }
+
+ /**
+ * Return RCP to gradient.
+ */
+ Teuchos::RCP<const NOX::Abstract::Vector>
+ getGradientPtr() const override
+ {
+ AssertThrow(false, ExcNotImplemented());
+ return {};
+ }
+
+ /**
+ * Return RCP to Newton direction.
+ */
+ Teuchos::RCP<const NOX::Abstract::Vector>
+ 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<NOX::Abstract::Group>
+ clone(NOX::CopyType copy_type) const override
+ {
+ auto new_group =
+ Teuchos::rcp(new Group<VectorType>(*x.vector,
+ residual,
+ setup_jacobian,
+ apply_jacobian,
+ solve_with_jacobian));
+
+ if (x.vector)
+ {
+ new_group->x.vector = std::make_shared<VectorType>();
+ new_group->x.vector->reinit(*x.vector);
+ }
+
+ if (f.vector)
+ {
+ new_group->f.vector = std::make_shared<VectorType>();
+ new_group->f.vector->reinit(*f.vector);
+ }
+
+ if (gradient.vector)
+ {
+ new_group->gradient.vector = std::make_shared<VectorType>();
+ new_group->gradient.vector->reinit(*gradient.vector);
+ }
+
+ if (newton.vector)
+ {
+ new_group->newton.vector = std::make_shared<VectorType>();
+ 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<VectorType>();
+
+ newton.vector->reinit(*f.vector, false);
+
+ const double tolerance = p.get<double>("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<const Vector<VectorType> *>(&input);
+ const auto *result_ =
+ dynamic_cast<const Vector<VectorType> *>(&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<VectorType> x, f, gradient, newton;
+
+ // 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)>
+ 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 <typename VectorType>
+ class NOXCheck : public NOX::StatusTest::Generic
+ {
+ public:
+ /**
+ * Constructor.
+ */
+ NOXCheck(const std::function<SolverControl::State(const unsigned int,
+ const double,
+ const VectorType &,
+ const VectorType &)>
+ 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<VectorType> *>(
+ &problem.getSolutionGroup().getX());
+ const VectorType &f__ = *dynamic_cast<
+ const internal::NOXWrappers::Vector<VectorType> *>(
+ &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<SolverControl::State(const unsigned int i,
+ const double f_norm,
+ const VectorType & x,
+ const VectorType & f)>
+ check_iteration_status;
+
+ /**
+ * Last retured value of checkStatus(), which is used for
+ * getStatus() and print().
+ */
+ NOX::StatusTest::StatusType status;
+ };
+ } // namespace NOXWrappers
+ } // namespace internal
+
+
+
+ template <typename VectorType>
+ NOXSolver<VectorType>::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 <typename VectorType>
+ NOXSolver<VectorType>::NOXSolver(
+ AdditionalData & additional_data,
+ const Teuchos::RCP<Teuchos::ParameterList> ¶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 <typename VectorType>
+ void
+ NOXSolver<VectorType>::clear()
+ {
+ // clear interal counters
+ n_residual_evaluations = 0;
+ n_jacobian_applications = 0;
+ n_nonlinear_iterations = 0;
+ n_last_linear_iterations = 0;
+ }
+
+
+
+ template <typename VectorType>
+ unsigned int
+ NOXSolver<VectorType>::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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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 <deal.II/base/mpi.h>
+
+#include <deal.II/lac/la_parallel_vector.h>
+
+#include <deal.II/trilinos/nox.h>
+
+// as reference solution
+#include <deal.II/lac/trilinos_precondition.h>
+#include <deal.II/lac/trilinos_solver.h>
+#include <deal.II/lac/trilinos_sparse_matrix.h>
+#include <deal.II/lac/trilinos_sparsity_pattern.h>
+
+#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<Epetra_CrsMatrix *>(&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<Number>;
+
+ // 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<VectorType>::AdditionalData additional_data(
+ n_max_iterations, abs_tolerance, rel_tolerance);
+
+ // set up parameters
+ Teuchos::RCP<Teuchos::ParameterList> 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<VectorType> 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<NoxInterface> interface = Teuchos::rcp(new NoxInterface());
+
+ Teuchos::RCP<NOX::Epetra::Interface::Required> iReq = interface;
+ Teuchos::RCP<NOX::Epetra::Interface::Jacobian> iJac = interface;
+
+ Teuchos::RCP<NOX::Epetra::LinearSystemAztecOO> linSys =
+ Teuchos::rcp(new NOX::Epetra::LinearSystemAztecOO(
+ printParams, lsParams, iReq, iJac, A, InitialGuess));
+
+ NOX::Epetra::Vector noxInitGuess(InitialGuess, NOX::DeepCopy);
+ Teuchos::RCP<NOX::Epetra::Group> 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<NOX::Solver::Generic> 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;
+ }
+}
projects / dealii.git / commitdiff