OGS
NewtonRaphson.h
Go to the documentation of this file.
1
11#pragma once
12
13#include <Eigen/Dense>
14#include <optional>
15
16#include "BaseLib/Logging.h"
17
18namespace NumLib
19{
21{
25};
26
31template <typename LinearSolver, typename JacobianMatrix,
32 typename JacobianMatrixUpdate, typename ResidualVector,
33 typename ResidualUpdate, typename SolutionUpdate>
34class NewtonRaphson final
35{
36public:
37 NewtonRaphson(LinearSolver& linear_solver,
38 JacobianMatrixUpdate jacobian_update,
39 ResidualUpdate residual_update,
40 SolutionUpdate solution_update,
41 NewtonRaphsonSolverParameters const& solver_parameters)
42 : _linear_solver(linear_solver),
43 _jacobian_update(jacobian_update),
44 _residual_update(residual_update),
45 _solution_update(solution_update),
46 _maximum_iterations(solver_parameters.maximum_iterations),
47 _residuum_tolerance_squared(solver_parameters.residuum_tolerance *
48 solver_parameters.residuum_tolerance),
49 _increment_tolerance_squared(solver_parameters.increment_tolerance *
50 solver_parameters.increment_tolerance)
51 {
52 }
53
56 std::optional<int> solve(JacobianMatrix& jacobian) const
57 {
58 int iteration = 0;
59 ResidualVector increment;
60 ResidualVector residual;
61 do
62 {
63 // The jacobian and the residual are updated simultaneously to keep
64 // consistency. The jacobian is used after the non-linear solver
65 // onward.
66 _jacobian_update(jacobian);
67 _residual_update(residual);
68
69 if (residual.squaredNorm() < _residuum_tolerance_squared)
70 {
71 break; // convergence criteria fulfilled.
72 }
73
74 increment.noalias() =
75 _linear_solver.compute(jacobian).solve(-residual);
76 // DBUG("Local linear solver accuracy |J dx - r| = {:g}",
77 // (jacobian * increment + residual).norm());
78
79 _solution_update(increment);
80
81 if (increment.squaredNorm() < _increment_tolerance_squared)
82 {
83 break; // increment to small.
84 }
85
86 // DBUG("Local Newton: Iteration #{:d} |dx| = {:g}, |r| = {:g}",
87 // iteration, increment.norm(), residual.norm());
88 // fmt::print("Local Newton: Increment {}\n",
89 // fmt::join(increment.data(),
90 // increment.data() + increment.size(), ", "));
91 // fmt::print("Local Newton: Residuum {}\n",
92 // fmt::join(residual.data(),
93 // residual.data() + residual.size(), ", "));
94 } while (iteration++ < _maximum_iterations);
95
96 if (iteration > _maximum_iterations)
97 {
98 ERR("The local Newton method did not converge within the given "
99 "number of iterations. Iteration: {:d}, increment {:g}, "
100 "residual: "
101 "{:g}",
102 iteration - 1, increment.norm(), residual.norm());
103 return std::nullopt;
104 }
105
106 return iteration;
107 };
108
109private:
110 LinearSolver& _linear_solver;
111 JacobianMatrixUpdate _jacobian_update;
112 ResidualUpdate _residual_update;
113 SolutionUpdate _solution_update;
115 const double
117 const double
119};
120} // namespace NumLib
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:44
SolutionUpdate _solution_update
NewtonRaphson(LinearSolver &linear_solver, JacobianMatrixUpdate jacobian_update, ResidualUpdate residual_update, SolutionUpdate solution_update, NewtonRaphsonSolverParameters const &solver_parameters)
Definition: NewtonRaphson.h:37
const double _increment_tolerance_squared
Error tolerance for the increment.
std::optional< int > solve(JacobianMatrix &jacobian) const
Definition: NewtonRaphson.h:56
const int _maximum_iterations
Maximum number of iterations.
const double _residuum_tolerance_squared
Error tolerance for the residuum.
LinearSolver & _linear_solver
JacobianMatrixUpdate _jacobian_update
ResidualUpdate _residual_update