version 0.4.4
newtonraphson.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
2// SPDX-License-Identifier: LGPL-3.0-or-later
3
9#pragma once
10
20
21namespace Ikarus {
22
23template <typename F, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault,
24 typename IDBCF = utils::IDBCForceDefault>
25class NewtonRaphson;
26
28{
29 double tol{1e-8};
30 int maxIter{20};
31 int minIter{0};
32};
33
38template <typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault,
39 typename IDBCF = utils::IDBCForceDefault>
41{
42 using LinearSolver = LS;
43 using UpdateFunction = UF;
44 using IDBCForceFunction = IDBCF;
49
50 template <typename UF2>
53 .linearSolver = linearSolver,
54 .updateFunction = std::forward<UF2>(updateFunction),
55 .idbcForceFunction = idbcForceFunction};
56 return settings;
57 }
58
59 template <typename IDBCF2>
62 .linearSolver = linearSolver,
63 .updateFunction = updateFunction,
64 .idbcForceFunction = std::forward<IDBCF2>(idbcForceFunction)};
65 return settings;
66 }
67
68 template <typename F>
70};
71
72// THE CTAD is broken for designated initializers in clang 16, when we drop support this can be simplified
73#ifndef DOXYGEN
77
78template <typename LS>
80
81template <typename UF>
84
85template <typename IDBCF>
88
89template <typename LS, typename UF>
91
92template <typename LS, typename UF, typename IDBCF>
94#endif
95
104template <typename F, typename NRConfig>
105requires traits::isSpecialization<NewtonRaphsonConfig, std::remove_cvref_t<NRConfig>>::value
106auto createNonlinearSolver(NRConfig&& config, F&& f) {
108 using UF = std::remove_cvref_t<NRConfig>::UpdateFunction;
109 using IDBCF = std::remove_cvref_t<NRConfig>::IDBCForceFunction;
110 auto solverFactory = []<class F2, class LS2, class UF2, class IDBCF2>(F2&& f2, LS2&& ls, UF2&& uf, IDBCF2&& if_) {
111 return std::make_shared<NewtonRaphson<std::remove_cvref_t<F2>, std::remove_cvref_t<LS2>, std::remove_cvref_t<UF2>,
112 std::remove_cvref_t<IDBCF2>>>(
113 f2, std::forward<LS2>(ls), std::forward<UF2>(uf), std::forward<IDBCF2>(if_));
114 };
115
116 if constexpr (std::remove_cvref_t<F>::nDerivatives == 2) {
117 auto solver =
118 solverFactory(derivative(f), std::forward<NRConfig>(config).linearSolver,
119 std::forward<NRConfig>(config).updateFunction, std::forward<NRConfig>(config).idbcForceFunction);
120 solver->setup(config.parameters);
121 return solver;
122 } else {
123 static_assert(std::remove_cvref_t<F>::nDerivatives >= 1,
124 "The number of derivatives in the differentiable function have to be more than 0");
125 auto solver =
126 solverFactory(f, std::forward<NRConfig>(config).linearSolver, std::forward<NRConfig>(config).updateFunction,
127 std::forward<NRConfig>(config).idbcForceFunction);
128 solver->setup(std::forward<NRConfig>(config).parameters);
129 return solver;
130 }
131}
132
143template <typename F, typename LS, typename UF, typename IDBCF>
145{
146public:
148 using SignatureTraits = typename F::Traits;
149 using CorrectionType = typename SignatureTraits::template Range<0>;
150 using JacobianType = typename SignatureTraits::template Range<1>;
153
154 using Domain = typename SignatureTraits::Domain;
155
156 using UpdateFunction = UF;
158 using IDBCForceFunction = IDBCF;
159
167 template <typename LS2 = LS, typename UF2 = UF, typename IDBCF2 = IDBCF>
168 explicit NewtonRaphson(const DifferentiableFunction& residual, LS2&& linearSolver = {}, UF2&& updateFunction = {},
169 IDBCF2&& idbcForceFunction = {})
170 : residualFunction_{residual},
171 jacobianFunction_{derivative(residualFunction_)},
172 linearSolver_{std::forward<LS2>(linearSolver)},
173 updateFunction_{std::forward<UF2>(updateFunction)},
174 idbcForceFunction_{std::forward<IDBCF2>(idbcForceFunction)} {}
175
180 void setup(const Settings& settings) {
181 if (settings.minIter > settings.maxIter)
182 DUNE_THROW(Dune::InvalidStateException,
183 "Minimum number of iterations cannot be greater than maximum number of iterations");
184 settings_ = settings;
185 }
186
193 [[nodiscard(
194 "The solve method returns information of the solution process. You should store this information and check if "
195 "it was successful")]] NonLinearSolverInformation
196 solve(Domain& x, double stepSize = 0.0) {
197 using enum NonLinearSolverMessages;
198
199 NonLinearSolverInformation solverInformation{};
200 auto state = typename NewtonRaphson::State(x, correction_, solverInformation);
201 this->notify(INIT, state);
202 solverInformation.success = true;
203
204 decltype(auto) rx = residualFunction_(x);
205 decltype(auto) Ax = jacobianFunction_(x);
206 auto rNorm = floatingPointNorm(rx);
207 solverInformation.residualNorm = static_cast<double>(rNorm);
208
209 decltype(rNorm) dNorm;
210 int iter{0};
211 if constexpr (isLinearSolver)
212 linearSolver_.analyzePattern(Ax);
213
214 if constexpr (not std::same_as<IDBCForceFunction, utils::IDBCForceDefault>) {
215 rx += idbcForceFunction_(x) * stepSize;
216 rNorm = floatingPointNorm(rx);
217 }
218
219 while ((rNorm > settings_.tol && iter < settings_.maxIter) or iter < settings_.minIter) {
220 this->notify(ITERATION_STARTED, state);
221 if constexpr (isLinearSolver) {
222 linearSolver_.factorize(Ax);
223 linearSolver_.solve(correction_, -rx);
224 } else {
225 correction_ = -linearSolver_(rx, Ax);
226 }
227 dNorm = floatingPointNorm(correction_);
228 solverInformation.correctionNorm = static_cast<double>(dNorm);
229
230 this->notify(CORRECTION_UPDATED, state);
231
232 if constexpr (requires { x.parameter(); })
233 updateFunction_(x.globalSolution(), correction_);
234 else
235 updateFunction_(x, correction_);
236
237 if constexpr (not std::same_as<IDBCForceFunction, utils::IDBCForceDefault>)
238 x.syncParameterAndGlobalSolution(updateFunction_);
239
240 this->notify(CORRECTIONNORM_UPDATED, state);
241 this->notify(SOLUTION_CHANGED, state);
242 rx = residualFunction_(x);
243 Ax = jacobianFunction_(x);
244 rNorm = floatingPointNorm(rx);
245 solverInformation.residualNorm = static_cast<double>(rNorm);
246 this->notify(RESIDUALNORM_UPDATED, state);
247 ++iter;
248 solverInformation.iterations = iter;
249 this->notify(ITERATION_ENDED, state);
250 }
251 if (iter == settings_.maxIter)
252 solverInformation.success = false;
253 solverInformation.iterations = iter;
254 if (solverInformation.success)
255 this->notify(FINISHED_SUCESSFULLY, state);
256 return solverInformation;
257 }
258
263 auto& residual() { return residualFunction_; }
264
269 const auto& residual() const { return residualFunction_; }
270
275 const UpdateFunction& updateFunction() const { return updateFunction_; }
276
281 const IDBCForceFunction& idbcForceFunction() const { return idbcForceFunction_; }
282
283private:
284 DifferentiableFunction residualFunction_;
285 typename DifferentiableFunction::Derivative jacobianFunction_;
286 std::remove_cvref_t<CorrectionType> correction_;
287 LS linearSolver_;
288 UpdateFunction updateFunction_;
289 IDBCForceFunction idbcForceFunction_;
290 Settings settings_;
291};
292
303template <typename F, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault,
304 typename IDBCF = utils::IDBCForceDefault>
305auto makeNewtonRaphson(const F& f, LS&& linearSolver = {}, UF&& updateFunction = {}, IDBCF&& idbcForceFunction = {}) {
306 return std::make_shared<NewtonRaphson<F, LS, UF, IDBCF>>(f, std::forward<LS>(linearSolver), std::move(updateFunction),
307 std::move(idbcForceFunction));
308}
309
310template <typename F, typename LS = utils::SolverDefault, typename UF = utils::UpdateDefault,
311 typename IDBCF = utils::IDBCForceDefault>
312NewtonRaphson(const F& f, LS&& linearSolver = {}, UF&& updateFunction = {}, IDBCF&& idbcForceFunction = {})
313 -> NewtonRaphson<F, std::remove_cvref_t<LS>, std::remove_cvref_t<UF>, std::remove_cvref_t<IDBCF>>;
314
315} // namespace Ikarus
Helper for the autodiff library.
Collection of fallback default functions.
Implementation of the observer design pattern with broadcasters.
Enums for observer messages.
Type-erased linear solver with templated scalar type.
Implementation of the solver information returned by the nonlinear solvers.
Base for all nonlinear solvers.
State for all nonlinear solvers.
auto floatingPointNorm(const Eigen::MatrixBase< Derived > &v)
Adding free floatingPointNorm function to Eigen types this is an indirection since otherwise norm fai...
Definition: linearalgebrahelper.hh:272
Definition: assemblermanipulatorbuildingblocks.hh:22
LinearSolverTemplate< double > LinearSolver
Definition: linearsolver.hh:235
auto makeNewtonRaphson(const F &f, LS &&linearSolver={}, UF &&updateFunction={}, IDBCF &&idbcForceFunction={})
Function to create a NewtonRaphson solver instance.
Definition: newtonraphson.hh:305
::value auto createNonlinearSolver(NRConfig &&config, F &&f)
Function to create a NewtonRaphson solver instance.
Definition: newtonraphson.hh:106
NonLinearSolverMessages
Enum class defining non-linear solver-related messages.
Definition: broadcastermessages.hh:22
NewtonRaphson(const F &f, LS &&linearSolver={}, UF &&updateFunction={}, IDBCF &&idbcForceFunction={}) -> NewtonRaphson< F, std::remove_cvref_t< LS >, std::remove_cvref_t< UF >, std::remove_cvref_t< IDBCF > >
Implementation of the Newton-Raphson method for solving nonlinear equations.
Definition: newtonraphson.hh:145
F DifferentiableFunction
Type of the non-linear operator.
Definition: newtonraphson.hh:157
IDBCF IDBCForceFunction
Type representing the force function to handle inhomogeneous Dirichlet BCs.
Definition: newtonraphson.hh:158
void setup(const Settings &settings)
Set up the solver with the given settings.
Definition: newtonraphson.hh:180
typename F::Traits SignatureTraits
Definition: newtonraphson.hh:148
typename SignatureTraits::Domain Domain
Type representing the parameter vector of the function.
Definition: newtonraphson.hh:154
UF UpdateFunction
Type representing the update function.
Definition: newtonraphson.hh:156
const UpdateFunction & updateFunction() const
Access the update function.
Definition: newtonraphson.hh:275
typename SignatureTraits::template Range< 1 > JacobianType
Compile-time boolean indicating if the linear solver satisfies the non-linear solver concept.
Definition: newtonraphson.hh:151
const IDBCForceFunction & idbcForceFunction() const
Access the force function calculating internal forces due to inhomogeneous Dirichlet BCs.
Definition: newtonraphson.hh:281
NonLinearSolverInformation solve(Domain &x, double stepSize=0.0)
Solve the nonlinear system.
Definition: newtonraphson.hh:196
NRSettings Settings
Definition: newtonraphson.hh:147
const auto & residual() const
Access the function.
Definition: newtonraphson.hh:269
static constexpr bool isLinearSolver
Definition: newtonraphson.hh:152
auto & residual()
Access the function.
Definition: newtonraphson.hh:263
NewtonRaphson(const DifferentiableFunction &residual, LS2 &&linearSolver={}, UF2 &&updateFunction={}, IDBCF2 &&idbcForceFunction={})
Constructor for NewtonRaphson.
Definition: newtonraphson.hh:168
typename SignatureTraits::template Range< 0 > CorrectionType
Type of the correction of x += deltaX.
Definition: newtonraphson.hh:149
Definition: newtonraphson.hh:28
int maxIter
Definition: newtonraphson.hh:30
int minIter
Definition: newtonraphson.hh:31
double tol
Definition: newtonraphson.hh:29
Config for the Newton-Raphson solver.
Definition: newtonraphson.hh:41
UF updateFunction
Definition: newtonraphson.hh:47
auto rebindUpdateFunction(UF2 &&updateFunction) const
Definition: newtonraphson.hh:51
IDBCF IDBCForceFunction
Definition: newtonraphson.hh:44
LS LinearSolver
Definition: newtonraphson.hh:42
IDBCF idbcForceFunction
Definition: newtonraphson.hh:48
NRSettings parameters
Definition: newtonraphson.hh:45
LS linearSolver
Definition: newtonraphson.hh:46
auto rebindIDBCForceFunction(IDBCF2 &&idbcForceFunction) const
Definition: newtonraphson.hh:60
UF UpdateFunction
Definition: newtonraphson.hh:43
Base for all nonlinear solvers. Defines the message interface that can be broadcasted to listeners.
Definition: nonlinearsolverbase.hh:24
NonlinearSolverStateType< F > State
Definition: nonlinearsolverbase.hh:25
Information about the result of a non-linear solver.
Definition: solverinfos.hh:21
void notify(NonLinearSolverMessages message, const State &data)
This calls all the registered functions.
Definition: broadcaster.hh:61
Default functor for solving operations.
Definition: defaultfunctions.hh:28
Default functor for updating operations.
Definition: defaultfunctions.hh:65
Concept to check if a linear solver implements all the needed functions for given vector and matrix t...
Definition: utils/concepts.hh:226
Several concepts.