black_box_newton_solver.cc
Go to the documentation of this file.
1// LIC// ====================================================================
2// LIC// This file forms part of oomph-lib, the object-oriented,
3// LIC// multi-physics finite-element library, available
4// LIC// at http://www.oomph-lib.org.
5// LIC//
6// LIC// Copyright (C) 2006-2022 Matthias Heil and Andrew Hazel
7// LIC//
8// LIC// This library is free software; you can redistribute it and/or
9// LIC// modify it under the terms of the GNU Lesser General Public
10// LIC// License as published by the Free Software Foundation; either
11// LIC// version 2.1 of the License, or (at your option) any later version.
12// LIC//
13// LIC// This library is distributed in the hope that it will be useful,
14// LIC// but WITHOUT ANY WARRANTY; without even the implied warranty of
15// LIC// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
16// LIC// Lesser General Public License for more details.
17// LIC//
18// LIC// You should have received a copy of the GNU Lesser General Public
19// LIC// License along with this library; if not, write to the Free Software
20// LIC// Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21// LIC// 02110-1301 USA.
22// LIC//
23// LIC// The authors may be contacted at oomph-lib@maths.man.ac.uk.
24// LIC//
25// LIC//====================================================================
26
27
28// Config header generated by autoconfig
29#ifdef HAVE_CONFIG_H
30#include <oomph-lib-config.h>
31#endif
32
34
35
36namespace oomph
37{
38 //======================================================================
39 /// Namespace for black-box FD Newton solver.
40 //======================================================================
41 namespace BlackBoxFDNewtonSolver
42 {
43 /// Max. # of Newton iterations
44 unsigned Max_iter = 20;
45
46 /// Number of Newton iterations taken in most recent invocation
47 unsigned N_iter_taken = 0;
48
49 /// Flag to indicate if progress of Newton iteration is to be
50 /// documented (defaults to false)
51 bool Doc_Progress = false;
52
53 /// FD step
54 double FD_step = 1.0e-8;
55
56 /// Tolerance
57 double Tol = 1.0e-8;
58
59 /// Use steplength control do make globally convergent (default false)
61
62 /// Black-box FD Newton solver:
63 /// Calling sequence for residual function is
64 /// \code residual_fct(parameters,unknowns,residuals) \endcode
65 /// where all arguments are double Vectors.
66 /// unknowns.size() = residuals.size()
68 const Vector<double>& params,
69 Vector<double>& unknowns,
70 JacobianFctPt jacobian_fct)
71 {
72 // Jacobian, current and advanced residual Vectors
73 unsigned ndof = unknowns.size();
74 DenseDoubleMatrix jacobian(ndof);
75 Vector<double> residuals(ndof);
76 Vector<double> residuals_pls(ndof);
77 Vector<double> dx(ndof);
78 Vector<double> gradient(ndof);
79 Vector<double> newton_direction(ndof);
80
81 double half_residual_squared = 0.0;
82 double max_step = 0.0;
83
84 /// Reset number of Newton iterations taken in most recent invocation
85 N_iter_taken = 0;
86
87 // Newton iterations
88 for (unsigned iloop = 0; iloop < Max_iter; iloop++)
89 {
90 // Evaluate current residuals
91 residual_fct(params, unknowns, residuals);
92
93 // Get half of squared residual and find maximum step length
94 // for step length control
96 {
97 half_residual_squared = 0.0;
98 double sum = 0.0;
99 for (unsigned i = 0; i < ndof; i++)
100 {
101 sum += unknowns[i] * unknowns[i];
102 half_residual_squared += residuals[i] * residuals[i];
103 }
104 half_residual_squared *= 0.5;
105 max_step = 100.0 * std::max(sqrt(sum), double(ndof));
106 }
107
108
109 // Check max. residuals
110 double max_res = std::fabs(*std::max_element(
111 residuals.begin(), residuals.end(), AbsCmp<double>()));
112
113
114 // Doc progress?
115 if (Doc_Progress)
116 {
117 oomph_info << "\nNewton iteration iter=" << iloop
118 << "\ni residual[i] unknown[i] " << std::endl;
119 for (unsigned i = 0; i < ndof; i++)
120 {
121 oomph_info << i << " " << residuals[i] << " " << unknowns[i]
122 << std::endl;
123 }
124 }
125
126
127 // Converged?
128 if (max_res < Tol)
129 {
130 return;
131 }
132
133 // Next iteration...
134 N_iter_taken++;
135
136 // ...and how would Sir like his Jacobian?
137 if (jacobian_fct == 0)
138 {
139 // FD loop for Jacobian
140 for (unsigned i = 0; i < ndof; i++)
141 {
142 double backup = unknowns[i];
143 unknowns[i] += FD_step;
144
145 // Evaluate advanced residuals
146 residual_fct(params, unknowns, residuals_pls);
147
148 // Do FD
149 for (unsigned j = 0; j < ndof; j++)
150 {
151 jacobian(j, i) = (residuals_pls[j] - residuals[j]) / FD_step;
152 }
153
154 // Reset fd step
155 unknowns[i] = backup;
156 }
157 }
158 // Analytical Jacobian
159 else
160 {
161 jacobian_fct(params, unknowns, jacobian);
162 }
163
164
165 if (Doc_Progress)
166 {
167 oomph_info << "\n\nJacobian: " << std::endl;
169 oomph_info << std::endl;
170 }
171
172 // Get gradient
174 {
175 for (unsigned i = 0; i < ndof; i++)
176 {
177 double sum = 0.0;
178 for (unsigned j = 0; j < ndof; j++)
179 {
180 sum += jacobian(j, i) * residuals[j];
181 }
182 gradient[i] = sum;
183 }
184 }
185
186 // Solve
187 jacobian.solve(residuals, newton_direction);
188
189 // Update
191 {
192 for (unsigned i = 0; i < ndof; i++)
193 {
194 newton_direction[i] *= -1.0;
195 }
196 // Update with steplength control
197 Vector<double> unknowns_old(unknowns);
198 double half_residual_squared_old = half_residual_squared;
199 line_search(unknowns_old,
200 half_residual_squared_old,
201 gradient,
202 residual_fct,
203 params,
204 newton_direction,
205 unknowns,
206 half_residual_squared,
207 max_step);
208 }
209 else
210 {
211 // Direct Newton update:
212 for (unsigned i = 0; i < ndof; i++)
213 {
214 unknowns[i] -= newton_direction[i];
215 }
216 }
217 }
218
219
220 // Failed to converge
221 std::ostringstream error_stream;
222 error_stream << "Newton solver did not converge in " << Max_iter
223 << " steps " << std::endl;
224
225 throw OomphLibError(
226 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
227 }
228
229
230 //=======================================================================
231 /// Line search helper for globally convergent Newton method
232 //=======================================================================
233 void line_search(const Vector<double>& x_old,
234 const double half_residual_squared_old,
235 const Vector<double>& gradient,
236 ResidualFctPt residual_fct,
237 const Vector<double>& params,
238 Vector<double>& newton_dir,
240 double& half_residual_squared,
241 const double& stpmax)
242 {
243 const double min_fct_decrease = 1.0e-4;
244 double convergence_tol_on_x = 1.0e-16;
245 double f_aux = 0.0;
246 double lambda_aux = 0.0;
247 double proposed_lambda = 0.0;
248 unsigned n = x_old.size();
249 double sum = 0.0;
250 for (unsigned i = 0; i < n; i++)
251 {
252 sum += newton_dir[i] * newton_dir[i];
253 }
254 sum = sqrt(sum);
255 if (sum > stpmax)
256 {
257 for (unsigned i = 0; i < n; i++)
258 {
259 newton_dir[i] *= stpmax / sum;
260 }
261 }
262 double slope = 0.0;
263 for (unsigned i = 0; i < n; i++)
264 {
265 slope += gradient[i] * newton_dir[i];
266 }
267 if (slope >= 0.0)
268 {
269 std::ostringstream error_stream;
270 error_stream << "Roundoff problem in lnsrch: slope=" << slope << "\n";
271 throw OomphLibError(
272 error_stream.str(), OOMPH_CURRENT_FUNCTION, OOMPH_EXCEPTION_LOCATION);
273 }
274 double test = 0.0;
275 for (unsigned i = 0; i < n; i++)
276 {
277 double temp =
278 std::fabs(newton_dir[i]) / std::max(std::fabs(x_old[i]), 1.0);
279 if (temp > test) test = temp;
280 }
281 double lambda_min = convergence_tol_on_x / test;
282 double lambda = 1.0;
283 while (true)
284 {
285 for (unsigned i = 0; i < n; i++)
286 {
287 x[i] = x_old[i] + lambda * newton_dir[i];
288 }
289
290 // Evaluate current residuals
291 Vector<double> residuals(n);
292 residual_fct(params, x, residuals);
293 half_residual_squared = 0.0;
294 for (unsigned i = 0; i < n; i++)
295 {
296 half_residual_squared += residuals[i] * residuals[i];
297 }
298 half_residual_squared *= 0.5;
299
300 if (lambda < lambda_min)
301 {
302 for (unsigned i = 0; i < n; i++) x[i] = x_old[i];
303
304 // Create an Oomph Lib warning
305 OomphLibWarning("Warning: Line search converged on x only!",
306 "BlackBoxFDNewtonSolver::line_search()",
307 OOMPH_EXCEPTION_LOCATION);
308 return;
309 }
310 else if (half_residual_squared <=
311 half_residual_squared_old + min_fct_decrease * lambda * slope)
312 {
313 return;
314 }
315 else
316 {
317 if (lambda == 1.0)
318 {
319 proposed_lambda =
320 -slope / (2.0 * (half_residual_squared -
321 half_residual_squared_old - slope));
322 }
323 else
324 {
325 double r1 = half_residual_squared - half_residual_squared_old -
326 lambda * slope;
327 double r2 = f_aux - half_residual_squared_old - lambda_aux * slope;
328 double a_poly =
329 (r1 / (lambda * lambda) - r2 / (lambda_aux * lambda_aux)) /
330 (lambda - lambda_aux);
331 double b_poly = (-lambda_aux * r1 / (lambda * lambda) +
332 lambda * r2 / (lambda_aux * lambda_aux)) /
333 (lambda - lambda_aux);
334 if (a_poly == 0.0)
335 {
336 proposed_lambda = -slope / (2.0 * b_poly);
337 }
338 else
339 {
340 double discriminant = b_poly * b_poly - 3.0 * a_poly * slope;
341 if (discriminant < 0.0)
342 {
343 proposed_lambda = 0.5 * lambda;
344 }
345 else if (b_poly <= 0.0)
346 {
347 proposed_lambda =
348 (-b_poly + sqrt(discriminant)) / (3.0 * a_poly);
349 }
350 else
351 {
352 proposed_lambda = -slope / (b_poly + sqrt(discriminant));
353 }
354 }
355 if (proposed_lambda > 0.5 * lambda)
356 {
357 proposed_lambda = 0.5 * lambda;
358 }
359 }
360 }
361 lambda_aux = lambda;
362 f_aux = half_residual_squared;
363 lambda = std::max(proposed_lambda, 0.1 * lambda);
364 }
365 }
366 } // namespace BlackBoxFDNewtonSolver
367
368
369} // namespace oomph
cstr elem_len * i
Definition: cfortran.h:603
Function-type-object to perform absolute comparison of objects. Apparently this inlines better.
Class of matrices containing doubles, and stored as a DenseMatrix<double>, but with solving functiona...
Definition: matrices.h:1271
void solve(DoubleVector &rhs)
Complete LU solve (replaces matrix by its LU decomposition and overwrites RHS with solution)....
Definition: matrices.cc:50
void sparse_indexed_output(std::ostream &outfile, const unsigned &precision=0, const bool &output_bottom_right_zero=false) const
Indexed output function to print a matrix to the stream outfile as i,j,a(i,j) for a(i,...
Definition: matrices.h:182
std::ostream *& stream_pt()
Access function for the stream pointer.
An OomphLibError object which should be thrown when an run-time error is encountered....
An OomphLibWarning object which should be created as a temporary object to issue a warning....
unsigned Max_iter
Max. # of Newton iterations.
unsigned N_iter_taken
Number of Newton iterations taken in most recent invocation.
void line_search(const Vector< double > &x_old, const double half_residual_squared_old, const Vector< double > &gradient, ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &newton_dir, Vector< double > &x, double &half_residual_squared, const double &stpmax)
Line search helper for globally convergent Newton method.
void(* JacobianFctPt)(const Vector< double > &, const Vector< double > &, DenseDoubleMatrix &jacobian)
bool Doc_Progress
Flag to indicate if progress of Newton iteration is to be documented (defaults to false)
void black_box_fd_newton_solve(ResidualFctPt residual_fct, const Vector< double > &params, Vector< double > &unknowns, JacobianFctPt jacobian_fct)
Black-box FD Newton solver: Calling sequence for residual function is.
void(* ResidualFctPt)(const Vector< double > &parameters, const Vector< double > &unknowns, Vector< double > &residuals)
bool Use_step_length_control
Use steplength control do make globally convergent (default false)
//////////////////////////////////////////////////////////////////// ////////////////////////////////...
OomphInfo oomph_info
Single (global) instantiation of the OomphInfo object – this is used throughout the library as a "rep...