OGS
NewtonRaphson.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <spdlog/fmt/bundled/ranges.h>
14
15#include <Eigen/Core>
16#include <optional>
17
18#include "BaseLib/Logging.h"
19
20namespace NumLib
21{
28
33template <typename LinearSolver, typename JacobianMatrixUpdate,
34 typename ResidualUpdate, typename SolutionUpdate>
35class NewtonRaphson final
36{
37private:
38 template <typename T>
40 {
41 static_assert(
42 !std::is_same_v<T, T>,
43 "For the type T in the expression !std::is_same_v<T, T>: Could not "
44 "find overload for FirstArgument<T>. Expected T to be a callable "
45 "with one argument. Could be a missing overload of FirstArgument "
46 "for the type T.");
47 };
48
49 // Specialization for lambdas and functors with const operator()
50 template <typename C, typename R, typename Arg, typename... Args>
51 struct FirstArgument<R (C::*)(Arg, Args...) const>
52 {
53 using type = std::remove_reference_t<Arg>;
54 };
55
56 template <typename F>
58 decltype(&std::remove_reference_t<F>::operator())>::type;
59
62
63public:
64 NewtonRaphson(LinearSolver& linear_solver,
65 JacobianMatrixUpdate jacobian_update,
66 ResidualUpdate residual_update,
67 SolutionUpdate solution_update,
68 NewtonRaphsonSolverParameters const& solver_parameters)
69 : _linear_solver(linear_solver),
70 _jacobian_update(jacobian_update),
71 _residual_update(residual_update),
72 _solution_update(solution_update),
73 _maximum_iterations(solver_parameters.maximum_iterations),
74 _residuum_tolerance_squared(solver_parameters.residuum_tolerance *
75 solver_parameters.residuum_tolerance),
76 _increment_tolerance_squared(solver_parameters.increment_tolerance *
77 solver_parameters.increment_tolerance)
78 {
79 }
80
83 std::optional<int> solve(JacobianMatrix& jacobian) const
84 {
85 int iteration = 0;
86 ResidualVector increment;
87 ResidualVector residual;
88 do
89 {
90 // The jacobian and the residual are updated simultaneously to keep
91 // consistency. The jacobian is used after the non-linear solver
92 // onward.
93 _jacobian_update(jacobian);
94 _residual_update(residual);
95
96 if (residual.squaredNorm() < _residuum_tolerance_squared)
97 {
98 break; // convergence criteria fulfilled.
99 }
100
101 increment.noalias() =
102 _linear_solver.compute(jacobian).solve(-residual);
103 // DBUG("Local linear solver accuracy |J dx - r| = {:g}",
104 // (jacobian * increment + residual).norm());
105
106 _solution_update(increment);
107
108 if (increment.squaredNorm() < _increment_tolerance_squared)
109 {
110 break; // increment to small.
111 }
112
113 // DBUG("Local Newton: Iteration #{:d} |dx| = {:g}, |r| = {:g}",
114 // iteration, increment.norm(), residual.norm());
115 // fmt::print("Local Newton: Increment {}\n",
116 // fmt::join(increment.data(),
117 // increment.data() + increment.size(), ", "));
118 // fmt::print("Local Newton: Residuum {}\n",
119 // fmt::join(residual.data(),
120 // residual.data() + residual.size(), ", "));
121 } while (iteration++ < _maximum_iterations);
122
123 if (iteration > _maximum_iterations)
124 {
125 ERR("The local Newton method did not converge within the given "
126 "number of iterations. Iteration: {:d}, increment {:g}, "
127 "residual: "
128 "{:g}",
129 iteration - 1, increment.norm(), residual.norm());
130 return std::nullopt;
131 }
132
133 return iteration;
134 };
135
136private:
137 LinearSolver& _linear_solver;
138 JacobianMatrixUpdate _jacobian_update;
139 ResidualUpdate _residual_update;
140 SolutionUpdate _solution_update;
142 const double
144 const double
146};
147} // namespace NumLib
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
ResidualUpdate _residual_update
JacobianMatrixUpdate _jacobian_update
const int _maximum_iterations
Maximum number of iterations.
FirstArgumentType< JacobianMatrixUpdate > JacobianMatrix
NewtonRaphson(LinearSolver &linear_solver, JacobianMatrixUpdate jacobian_update, ResidualUpdate residual_update, SolutionUpdate solution_update, NewtonRaphsonSolverParameters const &solver_parameters)
LinearSolver & _linear_solver
const double _residuum_tolerance_squared
Error tolerance for the residuum.
SolutionUpdate _solution_update
FirstArgumentType< ResidualUpdate > ResidualVector
typename FirstArgument< decltype(&std::remove_reference_t< F >::operator())>::type FirstArgumentType
std::optional< int > solve(JacobianMatrix &jacobian) const
const double _increment_tolerance_squared
Error tolerance for the increment.