#
# 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}")
#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
#include <deal.II/base/config.h>
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
# include <deal.II/base/exceptions.h>
# include <Teuchos_ParameterList.hpp>
+# include <functional>
+
DEAL_II_NAMESPACE_OPEN
namespace TrilinosWrappers
/**
- * Wrapper around the non-linear solver from the NOX
+ * Wrapper around the nonlinear solver from the NOX
* packge (https://docs.trilinos.org/dev/packages/nox/doc/html/index.html),
* targeting deal.II data structures.
+ *
+ * The following code shows the steps how to use this class:
+ * @code
+ * // set configuration
+ * TrilinosWrappers::NOXSolver<VectorType>::AdditionalData additional_data;
+ *
+ * // create nonlinear solver
+ * TrilinosWrappers::NOXSolver<VectorType> solver(additional_data);
+ *
+ * // set user functions to compute residual,
+ * // set up Jacobian, and to apply the inverse of
+ * // Jacobian; note that there are more functions that
+ * // can be set
+ * solver.residual = [](const auto &src, auto &dst) {...};
+ * solver.setup_jacobian = [](const auto &src) {...};
+ * solver.solve_with_jacobian =
+ * [](const auto &src, auto &dst, const auto) {...};
+ *
+ * // solver nonlinear system with solution containing the intial guess and
+ * // the final solution
+ * solver.solve(solution);
+ * @endcode
*/
template <typename VectorType>
class NOXSolver
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;
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;
};
/**
* 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.
*/
Teuchos::rcp(new Teuchos::ParameterList));
/**
- * Clear internal state.
+ * Clear the internal state.
*/
void
clear();
/**
- * Solve non-linear problem and return number of iterations.
+ * Solve the nonlinear problem and return the number of nonlinear
+ * iterations executed.
*/
unsigned int
solve(VectorType &solution);
/**
- * User function that computes the residual.
+ * A user function that computes the residual @p f based on the
+ * current solution @p x.
*
* @note This function should return 0 in the case of success.
*/
std::function<int(const VectorType &x, VectorType &f)> residual;
/**
- * User function that sets up the Jacobian.
+ * A user function that sets up the Jacobian, based on the
+ * current solution @p x.
*
* @note This function should return 0 in the case of success.
*/
std::function<int(const VectorType &x)> setup_jacobian;
/**
- * User function that sets up the preconditioner for inverting
- * the Jacobian.
+ * A user function that sets up the preconditioner for inverting
+ * the Jacobian, based on the current solution @p x.
*
- * @note The function is optional and is used when setup_jacobian is
+ * @note This function is optional and is used when setup_jacobian is
* called and the preconditioner needs to be updated (see
* update_preconditioner_predicate and
* AdditionalData::threshold_nonlinear_iterations).
std::function<int(const VectorType &x)> setup_preconditioner;
/**
- * User function that applies the Jacobian.
+ * A user function that applies the Jacobian to @p x and writes
+ * the result in @p v.
*
- * @note The function is optional and is used in the case of certain
+ * @note This function is optional and is used in the case of certain
* configurations. For instance, this function is required if the
* polynomial line search (@p NOX::LineSearch::Polynomial) is
* chosen, whereas for the full step case (@p NOX::LineSearch::FullStep)
std::function<int(const VectorType &x, VectorType &v)> apply_jacobian;
/**
- * User function that applies the inverse of the Jacobian.
+ * A user function that applies the inverse of the Jacobian to
+ * @p x and writes the result in @p x. The parameter @p tolerance
+ * species the error reduction if a interative solver is used.
*
- * @note The function is optional and is used in the case of certain
+ * @note This function is optional and is used in the case of certain
* configurations.
*
* @note This function should return 0 in the case of success.
solve_with_jacobian;
/**
- * User function that applies the inverse of the Jacobian and
- * returns the numer of linear iterations the linear solver needed.
+ * A user function that applies the inverse of the Jacobian to
+ * @p x, writes the result in @p x and returns the numer of
+ * linear iterations the linear solver needed.
+ * The parameter @p tolerance species the error reduction if a
+ * interative solver is used.
*
- * @note This function should return -1 in the case of failure.
+ * @note This function is optional and is used in the case of certain
+ * configurations.
*/
std::function<
int(const VectorType &f, VectorType &x, const double tolerance)>
solve_with_jacobian_and_track_n_linear_iterations;
/**
- * User function that allows to check convergence in addition to
+ * A user function that allows to check convergence in addition to
* ones checking the l2-norm and the number of iterations (see
- * AdditionalData). It is run after each non-linear iteration.
+ * AdditionalData). It is run after each nonlinear iteration.
*
* The input are the current iteration number @p i, the l2-norm
* @p norm_f of the residual vector, the current solution @p x,
* and the current residual vector @p f.
*
- * @note The function is optional.
+ * @note This function is optional.
*/
std::function<SolverControl::State(const unsigned int i,
const double norm_f,
check_iteration_status;
/**
- * Function that allows to force to update the preconditioner in
- * addition to AdditionalData::threshold_nonlinear_iterations. A reason
+ * A user function that, in addition to
+ * AdditionalData::threshold_nonlinear_iterations,
+ * allows to force to update the preconditioner. A reason
* for wanting to update the preconditioner is when the expected number
- * of linear iterations exceeds.
+ * of linear iterations is exceeded.
*
- * @note The function is optional. If no function is attached, this
+ * @note This function is optional. If no function is attached, this
* means implicitly a return value of false.
*/
std::function<bool()> update_preconditioner_predicate;
const Teuchos::RCP<Teuchos::ParameterList> parameters;
/**
- * Counter for number of (accumulated) residual evaluations.
+ * A counter for the number of (accumulated) residual evaluations.
*/
unsigned int n_residual_evaluations;
/**
- * Counter for number of (accumulated) Jacobi applications.
+ * A counter for the number of (accumulated) Jacobi applications.
*/
unsigned int n_jacobian_applications;
/**
- * Counter for number of (accumulated) non-linear iterations.
+ * A counter for the number of (accumulated) nonlinear iterations.
*/
unsigned int n_nonlinear_iterations;
/**
- * Number of linear iterations of the last Jacobian solve.
+ * The number of linear iterations of the last Jacobian solve.
*/
unsigned int n_last_linear_iterations;
};
#include <deal.II/trilinos/nox.h>
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
# include <NOX_Abstract_Group.H>
# include <NOX_Abstract_Vector.H>
# include <NOX_StatusTest_NormF.H>
# include <NOX_StatusTest_RelativeNormF.H>
+# include <memory>
+
DEAL_II_NAMESPACE_OPEN
# ifndef DOXYGEN
NOX::Abstract::Vector &
init(double gamma) override
{
- *vector = gamma;
+ *vector = static_cast<typename VectorType::value_type>(gamma);
return *this;
}
NOX::Abstract::Vector &
scale(double gamma) override
{
- *vector *= gamma;
+ *vector *= static_cast<typename VectorType::value_type>(gamma);
return *this;
}
Assert(a_, ExcInternalError());
- vector->sadd(gamma, alpha, *a_->vector);
+ vector->sadd(static_cast<typename VectorType::value_type>(gamma),
+ static_cast<typename VectorType::value_type>(alpha),
+ *a_->vector);
return *this;
}
Assert(a_, ExcInternalError());
Assert(b_, ExcInternalError());
- vector->operator*=(gamma);
- vector->add(alpha, *a_->vector, beta, *b_->vector);
+ vector->operator*=(
+ static_cast<typename VectorType::value_type>(gamma));
+ vector->add(static_cast<typename VectorType::value_type>(alpha),
+ *a_->vector,
+ static_cast<typename VectorType::value_type>(beta),
+ *b_->vector);
return *this;
}
}
/**
- * Norm.
+ * Return the vector norm.
*/
double
norm(NOX::Abstract::Vector::NormType type =
}
/**
- * Weighted 2-Norm.
+ * Return the vector's weighted 2-Norm.
*/
double
norm(const NOX::Abstract::Vector &weights) const override
}
/**
- * Inner product with y.
+ * Return the inner product of the vector with @p y.
*/
double
innerProduct(const NOX::Abstract::Vector &y) const override
private:
/**
- * Underlying deal.II vector.
+ * The underlying deal.II vector.
*/
std::shared_ptr<VectorType> vector;
}
/**
- * Set the solution vector x to y.
+ * Set the solution vector `x` to @y.
*/
void
setX(const NOX::Abstract::Vector &y) override
}
/**
- * Compute x = grp.x + step * d.
+ * Compute the solution update `x = grp.x + step * d`.
*/
void
computeX(const NOX::Abstract::Group & grp,
}
/**
- * Compute and store F(x).
+ * Compute and store the residual update `F(x)`.
*/
NOX::Abstract::Group::ReturnType
computeF() override
}
/**
- * Return true if F is valid.
+ * Return true if the residual vector `F` is valid.
*/
bool
isF() const override
}
/**
- * Compute and store Jacobian.
+ * Compute and store the Jacobian.
*/
NOX::Abstract::Group::ReturnType
computeJacobian() override
}
/**
- * Return solution vector.
+ * Return the (total) solution vector.
*/
const NOX::Abstract::Vector &
getX() const override
}
/**
- * Return F(x).
+ * Return the residual `F(x)`.
*/
const NOX::Abstract::Vector &
getF() const override
}
/**
- * Return 2-norm of F(x)
+ * Return the 2-norm of `F(x)`.
*/
double
getNormF() const override
}
/**
- * Return gradient.
+ * Return the gradient.
*/
const NOX::Abstract::Vector &
getGradient() const override
}
/**
- * Return Newton direction.
+ * Return the Newton descent direction.
*/
const NOX::Abstract::Vector &
getNewton() const override
}
/**
- * Return RCP to solution vector.
+ * Return an `RCP` to solution vector.
*/
Teuchos::RCP<const NOX::Abstract::Vector>
getXPtr() const override
}
/**
- * Return RCP to F(x).
+ * Return an RCP to the residual `F(x)`.
*/
Teuchos::RCP<const NOX::Abstract::Vector>
getFPtr() const override
}
/**
- * Return RCP to gradient.
+ * Return an `RCP` to gradient.
*/
Teuchos::RCP<const NOX::Abstract::Vector>
getGradientPtr() const override
}
/**
- * Return RCP to Newton direction.
+ * Return an `RCP` to the Newton descent direction.
*/
Teuchos::RCP<const NOX::Abstract::Vector>
getNewtonPtr() const override
/**
* Create a new Group of the same derived type as this one by
- * cloning this one, and return a ref count pointer to the new group.
+ * cloning this one, and return a reference counting
+ * pointer to the new group.
*/
Teuchos::RCP<NOX::Abstract::Group>
clone(NOX::CopyType copy_type) const override
}
/**
- * 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
}
/**
- * 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,
is_valid_j = false;
}
- // internal vectors
- Vector<VectorType> x, f, gradient, newton;
+ // internal vector for the current solution
+ Vector<VectorType> x;
+
+ // internal vector for the residual
+ Vector<VectorType> f;
+
+ // internal vector for the solution gradient
+ Vector<VectorType> gradient;
+
+ // internal vector for the newton step
+ Vector<VectorType> newton;
+
+ // helper function to compute residual
+ std::function<int(const VectorType &x, VectorType &f)> residual;
+
+ // helper function to setup Jacobian
+ std::function<int(const VectorType &x)> setup_jacobian;
- // helper functions to compute residual, to setup jacobian, and
- // solve jacobian
- std::function<int(const VectorType &, VectorType &)> residual;
- std::function<int(const VectorType &)> setup_jacobian;
- std::function<int(const VectorType &, VectorType &)> apply_jacobian;
- std::function<int(const VectorType &, VectorType &, const double)>
+ // helper function to apply Jacobian
+ std::function<int(const VectorType &x, VectorType &v)> apply_jacobian;
+
+ // helper function to solve jacobian
+ std::function<
+ int(const VectorType &f, VectorType &x, const double tolerance)>
solve_with_jacobian;
- // internal state (are residuum and jacobian computed?)
- bool is_valid_f, is_valid_j;
+ // internal state (is residual computed?)
+ bool is_valid_f;
+
+ // internal state (is Jacobian computed?)
+ bool is_valid_j;
};
{}
/**
- * Check status.
+ * Check the status of the nonlinear solver.
*/
NOX::StatusTest::StatusType
checkStatus(const NOX::Solver::Generic &problem,
}
/**
- * Return last return value of checkStatus().
+ * Return the last value that was given by `checkStatus()`.
*/
NOX::StatusTest::StatusType
getStatus() const override
}
/**
- * 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
private:
/**
- * User function that allows to check convergence.
+ * The user function that allows the solver to check for
+ * convergence.
*/
const std::function<SolverControl::State(const unsigned int i,
const double f_norm,
check_iteration_status;
/**
- * Last retured value of checkStatus(), which is used for
- * getStatus() and print().
+ * The last returned value of `checkStatus()`.
*/
NOX::StatusTest::StatusType status;
};
// invert Jacobian
if (solve_with_jacobian)
{
+ Assert(
+ !solve_with_jacobian_and_track_n_linear_iterations,
+ ExcMessage(
+ "It does not make sense to provide both solve_with_jacobian and "
+ "solve_with_jacobian_and_track_n_linear_iterations!"));
+
// without tracking of linear iterations
return solve_with_jacobian(f, x, tolerance);
}
if (this->check_iteration_status)
{
- const auto info = Teuchos::rcp(
- new internal::NOXWrappers::NOXCheck(this->check_iteration_status));
+ const auto info =
+ Teuchos::rcp(new internal::NOXWrappers::NOXCheck<VectorType>(
+ this->check_iteration_status));
check->addStatusTest(info);
}
##
## ---------------------------------------------------------------------
-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}")
#include <deal.II/trilinos/nox.templates.h>
-#ifdef DEAL_II_WITH_TRILINOS
+#ifdef DEAL_II_TRILINOS_WITH_NOX
DEAL_II_NAMESPACE_OPEN
projects / dealii.git / commitdiff