Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[WIP] Add functionality to have user-specific convergence criteria for nonlinear solvers #341

Draft
wants to merge 5 commits into
base: main
Choose a base branch
from
Draft
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Original file line number Diff line number Diff line change
Expand Up @@ -190,7 +190,7 @@ private:

int minIter = isAutoDiff ? 1 : 0;
// THE CTAD is broken for designated initializers in clang 16, when we drop support this can be simplified
NewtonRaphsonConfig<decltype(linearSolver), decltype(updateFunction)> nrs{
NewtonRaphsonConfig<ConvergenceCriterion::ResiduumNorm, decltype(linearSolver), decltype(updateFunction)> nrs{
.parameters = {.tol = tol_, .maxIter = 100, .minIter = minIter},
.linearSolver = linearSolver,
.updateFunction = updateFunction
Expand Down
2 changes: 1 addition & 1 deletion ikarus/solver/nonlinearsolver/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,6 @@

# install headers
install(FILES newtonraphson.hh trustregion.hh newtonraphsonwithscalarsubsidiaryfunction.hh
nonlinearsolverfactory.hh solverinfos.hh
nonlinearsolverfactory.hh solverinfos.hh convergencecriteria.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/ikarus/solver/nonlinearsolver
)
136 changes: 136 additions & 0 deletions ikarus/solver/nonlinearsolver/convergencecriteria.hh
Original file line number Diff line number Diff line change
@@ -0,0 +1,136 @@
// SPDX-FileCopyrightText: 2021-2024 The Ikarus Developers [email protected]
// SPDX-License-Identifier: LGPL-3.0-or-later

/**
* \file convergencecriteria.hh
* \brief A collection of convergence criteria used by the nonlinear solvers.
*/

#include <dune/common/float_cmp.hh>

#include <ikarus/controlroutines/pathfollowingfunctions.hh>
#include <ikarus/utils/defaultfunctions.hh>
#include <ikarus/utils/linearalgebrahelper.hh>
#include <ikarus/utils/makeenum.hh>

#pragma once

namespace Ikarus {

/**
* \enum PreConditioner
* \brief Enumeration of available preconditioners for the trust region solver.
*/
enum class PreConditioner
{
IncompleteCholesky,
IdentityPreconditioner,
DiagonalPreconditioner
};

MAKE_ENUM(ConvergenceCriterion, ResiduumNorm, CorrectionNorm, MaximumOfCorrectionVector);

template <typename NLO, ConvergenceCriterion CC = ConvergenceCriterion::ResiduumNorm,
typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
class NewtonRaphson;

template <typename NLO, ConvergenceCriterion CC = ConvergenceCriterion::ResiduumNorm,
typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
class NewtonRaphsonWithSubsidiaryFunction;

template <typename NLO, PreConditioner preConditioner = PreConditioner::IncompleteCholesky,
typename UF = utils::UpdateDefault>
class TrustRegion;

template <ConvergenceCriterion CC>
requires(CC != ConvergenceCriterion::BEGIN and CC != ConvergenceCriterion::END)
struct NRConvergenceCriteria
{
template <typename NLO, typename NSS, typename CV>
bool operator()(const NLO& nonLinearOperator, const NSS& solverSettings, const CV& correctionVector) {
if constexpr (CC == ConvergenceCriterion::ResiduumNorm) {
const auto& rx = nonLinearOperator.value();
auto rNorm = norm(rx);
return Dune::FloatCmp::le(static_cast<double>(rNorm), solverSettings.tol);
} else if constexpr (CC == ConvergenceCriterion::CorrectionNorm) {
auto dNorm = norm(correctionVector);
return Dune::FloatCmp::le(static_cast<double>(dNorm), solverSettings.tol);
} else if constexpr (CC == ConvergenceCriterion::MaximumOfCorrectionVector) {
auto maxCorr = max(correctionVector);
return Dune::FloatCmp::le(static_cast<double>(maxCorr), solverSettings.tol);
} else {
return false; // TODO Change to Dune::AlwaysFalse
}
}
};

template <ConvergenceCriterion CC>
requires(CC != ConvergenceCriterion::BEGIN and CC != ConvergenceCriterion::END)
struct NRWSFConvergenceCriteria
{
template <typename NLO, typename NSS, typename DD>
bool operator()(const NLO& nonLinearOperator, const NSS& solverSettings, const DD& deltaD, double deltaLambda,
const SubsidiaryArgs& subsidiaryArgs) {
if constexpr (CC == ConvergenceCriterion::ResiduumNorm) {
const auto& rx = nonLinearOperator.value();
Eigen::VectorXd totalResidual(rx.size() + 1);
totalResidual << rx, subsidiaryArgs.f;
auto rNorm = norm(totalResidual);
return Dune::FloatCmp::le(static_cast<double>(rNorm), solverSettings.tol);
} else if constexpr (CC == ConvergenceCriterion::CorrectionNorm) {
Eigen::VectorXd correctionVector(deltaD.size() + 1);
correctionVector << deltaD, deltaLambda;
auto dNorm = norm(correctionVector);
return Dune::FloatCmp::le(static_cast<double>(dNorm), solverSettings.tol);
} else if constexpr (CC == ConvergenceCriterion::MaximumOfCorrectionVector) {
Eigen::VectorXd correctionVector(deltaD.size() + 1);
correctionVector << deltaD, deltaLambda;
auto maxCorr = max(correctionVector);
return Dune::FloatCmp::le(static_cast<double>(maxCorr), solverSettings.tol);
} else {
return false; // TODO Change to Dune::AlwaysFalse
}
}
};

template <ConvergenceCriterion CC>
requires(CC != ConvergenceCriterion::BEGIN and CC != ConvergenceCriterion::END)
struct TRConvergenceCriteria
{
template <typename NLO, typename NSS, typename CV>
bool operator()(const NLO& nonLinearOperator, const NSS& solverSettings, const CV& correctionVector) {
if constexpr (CC == ConvergenceCriterion::ResiduumNorm) {
return false;
} else if constexpr (CC == ConvergenceCriterion::CorrectionNorm) {
return false;
} else if constexpr (CC == ConvergenceCriterion::MaximumOfCorrectionVector) {
return false;
} else {
return false; // TODO Change to Dune::AlwaysFalse
}
}
};

template <typename NLS, ConvergenceCriterion CC>
struct ConvergenceCriteria
{
using NR = NewtonRaphson<typename NLS::NonLinearOperator, NLS::criteraType, typename NLS::LinearSolverType,
typename NLS::UpdateFunctionType>;
using NRWSF = NewtonRaphsonWithSubsidiaryFunction<typename NLS::NonLinearOperator, NLS::criteraType,
typename NLS::LinearSolverType, typename NLS::UpdateFunctionType>;
using TR = TrustRegion<typename NLS::NonLinearOperator>;

using Criteria = std::conditional_t<
std::is_same_v<NLS, NR>, NRConvergenceCriteria<CC>,
std::conditional_t<std::is_same_v<NLS, NRWSF>, NRWSFConvergenceCriteria<CC>, TRConvergenceCriteria<CC>>>;

template <typename... Args>
bool operator()(Args&&... args) {
return criteria(std::forward<Args>(args)...);
}

private:
Criteria criteria{};
};

} // namespace Ikarus
55 changes: 34 additions & 21 deletions ikarus/solver/nonlinearsolver/newtonraphson.hh
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#pragma once

#include <ikarus/solver/linearsolver/linearsolver.hh>
#include <ikarus/solver/nonlinearsolver/convergencecriteria.hh>
#include <ikarus/solver/nonlinearsolver/solverinfos.hh>
#include <ikarus/utils/concepts.hh>
#include <ikarus/utils/defaultfunctions.hh>
Expand All @@ -18,9 +19,6 @@

namespace Ikarus {

template <typename NLO, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
class NewtonRaphson;

struct NRSettings
{
double tol{1e-8};
Expand All @@ -32,24 +30,28 @@ struct NRSettings
* \struct NewtonRaphsonConfig
* \brief Config for the Newton-Raphson solver.
*/
template <typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
template <ConvergenceCriterion CC = ConvergenceCriterion::ResiduumNorm, typename LS = utils::SolverDefault,
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

template<typename NLS>
struct ResidualNorm

template<>
struct ResidualNorm<NRwithsubsidary>

typename UF = utils::UpdateDefault>
struct NewtonRaphsonConfig
{
using LinearSolver = LS;
using UpdateFunction = UF;

static constexpr auto criteria = CC;

NRSettings parameters;
LS linearSolver;
UF updateFunction;

template <typename UF2>
auto rebindUpdateFunction(UF2&& updateFunction) const {
NewtonRaphsonConfig<LS, UF2> settings{
NewtonRaphsonConfig<CC, LS, UF2> settings{
.parameters = parameters, .linearSolver = linearSolver, .updateFunction = std::forward<UF2>(updateFunction)};
return settings;
}

template <typename NLO>
using Solver = NewtonRaphson<NLO, LS, UF>;
using Solver = NewtonRaphson<NLO, CC, LS, UF>;
};

/**
Expand All @@ -61,13 +63,14 @@ struct NewtonRaphsonConfig
* \return Shared pointer to the NewtonRaphson solver instance.
*/
template <typename NLO, typename NRConfig>
requires traits::isSpecialization<NewtonRaphsonConfig, std::remove_cvref_t<NRConfig>>::value
requires traits::isSpecializationNonTypeAndTypes<NewtonRaphsonConfig, std::remove_cvref_t<NRConfig>>::value
auto createNonlinearSolver(NRConfig&& config, NLO&& nonLinearOperator) {
using LS = std::remove_cvref_t<NRConfig>::LinearSolver;
using UF = std::remove_cvref_t<NRConfig>::UpdateFunction;
auto solverFactory = []<class NLO2, class LS2, class UF2>(NLO2&& nlo2, LS2&& ls, UF2&& uf) {
using LS = std::remove_cvref_t<NRConfig>::LinearSolver;
using UF = std::remove_cvref_t<NRConfig>::UpdateFunction;
static constexpr auto criteria = std::remove_cvref_t<NRConfig>::criteria;
auto solverFactory = []<class NLO2, class LS2, class UF2>(NLO2&& nlo2, LS2&& ls, UF2&& uf) {
return std::make_shared<
NewtonRaphson<std::remove_cvref_t<NLO2>, std::remove_cvref_t<LS2>, std::remove_cvref_t<UF2>>>(
NewtonRaphson<std::remove_cvref_t<NLO2>, criteria, std::remove_cvref_t<LS2>, std::remove_cvref_t<UF2>>>(
nlo2, std::forward<LS2>(ls), std::forward<UF2>(uf));
};

Expand Down Expand Up @@ -98,7 +101,7 @@ auto createNonlinearSolver(NRConfig&& config, NLO&& nonLinearOperator) {
* \relates makeNewtonRaphson
* \ingroup solvers
*/
template <typename NLO, typename LS, typename UF>
template <typename NLO, ConvergenceCriterion CC, typename LS, typename UF>
class NewtonRaphson : public IObservable<NonLinearSolverMessages>
{
public:
Expand All @@ -110,8 +113,10 @@ public:
///< Type representing the parameter vector of the nonlinear operator.
using ValueType = typename NLO::template ParameterValue<0>;

using UpdateFunction = UF; ///< Type representing the update function.
using NonLinearOperator = NLO; ///< Type of the non-linear operator
using LinearSolverType = LS;
using UpdateFunctionType = UF; ///< Type representing the update function.
using NonLinearOperator = NLO; ///< Type of the non-linear operator
static constexpr auto criteraType = CC;

/**
* \brief Constructor for NewtonRaphson.
Expand Down Expand Up @@ -170,7 +175,12 @@ public:
int iter{0};
if constexpr (isLinearSolver)
linearSolver_.analyzePattern(Ax);
while ((rNorm > settings_.tol && iter < settings_.maxIter) or iter < settings_.minIter) {

auto criteria = ConvergenceCriteria<NewtonRaphson<NLO, CC, LS, UF>, CC>{};

bool converged = criteria(nonLinearOperator(), settings_, correction_);

while ((not(converged) && iter < settings_.maxIter) or iter < settings_.minIter) {
this->notify(NonLinearSolverMessages::ITERATION_STARTED);
if constexpr (isLinearSolver) {
linearSolver_.factorize(Ax);
Expand All @@ -189,6 +199,7 @@ public:
this->notify(NonLinearSolverMessages::RESIDUALNORM_UPDATED, static_cast<double>(rNorm));
this->notify(NonLinearSolverMessages::ITERATION_ENDED);
++iter;
converged = criteria(nonLinearOperator(), settings_, correction_);
}
if (iter == settings_.maxIter)
solverInformation.success = false;
Expand All @@ -210,7 +221,7 @@ private:
NonLinearOperator nonLinearOperator_;
typename NonLinearOperator::ValueType correction_;
LS linearSolver_;
UpdateFunction updateFunction_;
UpdateFunctionType updateFunction_;
Settings settings_;
};

Expand All @@ -224,14 +235,16 @@ private:
* \param updateFunction Update function (default is UpdateDefault).
* \return Shared pointer to the NewtonRaphson solver instance.
*/
template <typename NLO, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
template <typename NLO, ConvergenceCriterion CC = ConvergenceCriterion::ResiduumNorm,
typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
auto makeNewtonRaphson(const NLO& nonLinearOperator, LS&& linearSolver = {}, UF&& updateFunction = {}) {
return std::make_shared<NewtonRaphson<NLO, LS, UF>>(nonLinearOperator, std::forward<LS>(linearSolver),
std::move(updateFunction));
return std::make_shared<NewtonRaphson<NLO, CC, LS, UF>>(nonLinearOperator, std::forward<LS>(linearSolver),
std::move(updateFunction));
}

template <typename NLO, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
template <typename NLO, ConvergenceCriterion CC = ConvergenceCriterion::ResiduumNorm,
typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault>
NewtonRaphson(const NLO& nonLinearOperator, LS&& linearSolver = {},
UF&& updateFunction = {}) -> NewtonRaphson<NLO, std::remove_cvref_t<LS>, std::remove_cvref_t<UF>>;
UF&& updateFunction = {}) -> NewtonRaphson<NLO, CC, std::remove_cvref_t<LS>, std::remove_cvref_t<UF>>;

} // namespace Ikarus
Loading
Loading