1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 1998 - 2018 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15
16 #include <deal.II/base/logstream.h>
17 #include <deal.II/base/parameter_handler.h>
18 #include <deal.II/base/signaling_nan.h>
19
20 #include <deal.II/lac/solver_control.h>
21
22 #include <cmath>
23 #include <sstream>
24
25 DEAL_II_NAMESPACE_OPEN
26
27 /*----------------------- SolverControl ---------------------------------*/
28
29
SolverControl(const unsigned int maxiter,const double tolerance,const bool m_log_history,const bool m_log_result)30 SolverControl::SolverControl(const unsigned int maxiter,
31 const double tolerance,
32 const bool m_log_history,
33 const bool m_log_result)
34 : maxsteps(maxiter)
35 , tol(tolerance)
36 , lcheck(failure)
37 , initial_val(numbers::signaling_nan<double>())
38 , lvalue(numbers::signaling_nan<double>())
39 , lstep(numbers::invalid_unsigned_int)
40 , check_failure(false)
41 , relative_failure_residual(0)
42 , failure_residual(0)
43 , m_log_history(m_log_history)
44 , m_log_frequency(1)
45 , m_log_result(m_log_result)
46 , history_data_enabled(false)
47 {}
48
49
50
51 SolverControl::State
check(const unsigned int step,const double check_value)52 SolverControl::check(const unsigned int step, const double check_value)
53 {
54 // if this is the first time we
55 // come here, then store the
56 // residual for later comparisons
57 if (step == 0)
58 {
59 initial_val = check_value;
60 }
61
62 if (m_log_history && ((step % m_log_frequency) == 0))
63 deallog << "Check " << step << "\t" << check_value << std::endl;
64
65 lstep = step;
66 lvalue = check_value;
67
68 if (step == 0)
69 {
70 if (check_failure)
71 failure_residual = relative_failure_residual * check_value;
72
73 if (m_log_result)
74 deallog << "Starting value " << check_value << std::endl;
75 }
76
77 if (history_data_enabled)
78 history_data.push_back(check_value);
79
80 if (check_value <= tol)
81 {
82 if (m_log_result)
83 deallog << "Convergence step " << step << " value " << check_value
84 << std::endl;
85 lcheck = success;
86 return success;
87 }
88
89 if ((step >= maxsteps) || std::isnan(check_value) ||
90 (check_failure && (check_value > failure_residual)))
91 {
92 if (m_log_result)
93 deallog << "Failure step " << step << " value " << check_value
94 << std::endl;
95 lcheck = failure;
96 return failure;
97 }
98
99 lcheck = iterate;
100 return iterate;
101 }
102
103
104
105 SolverControl::State
last_check() const106 SolverControl::last_check() const
107 {
108 return lcheck;
109 }
110
111
112 double
initial_value() const113 SolverControl::initial_value() const
114 {
115 return initial_val;
116 }
117
118
119 double
last_value() const120 SolverControl::last_value() const
121 {
122 return lvalue;
123 }
124
125
126 unsigned int
last_step() const127 SolverControl::last_step() const
128 {
129 return lstep;
130 }
131
132
133 unsigned int
log_frequency(unsigned int f)134 SolverControl::log_frequency(unsigned int f)
135 {
136 if (f == 0)
137 f = 1;
138 unsigned int old = m_log_frequency;
139 m_log_frequency = f;
140 return old;
141 }
142
143
144 void
enable_history_data()145 SolverControl::enable_history_data()
146 {
147 history_data_enabled = true;
148 }
149
150
151
152 const std::vector<double> &
get_history_data() const153 SolverControl::get_history_data() const
154 {
155 Assert(history_data_enabled, ExcHistoryDataRequired());
156 Assert(
157 history_data.size() > 0,
158 ExcMessage(
159 "The SolverControl object was asked for the solver history "
160 "data, but there is no data. Possibly you requested the data before the "
161 "solver was run."));
162
163 return history_data;
164 }
165
166
167
168 double
average_reduction() const169 SolverControl::average_reduction() const
170 {
171 if (lstep == 0)
172 return 0.;
173
174 Assert(history_data_enabled, ExcHistoryDataRequired());
175 Assert(history_data.size() > lstep, ExcInternalError());
176 Assert(history_data[0] > 0., ExcInternalError());
177 Assert(history_data[lstep] > 0., ExcInternalError());
178
179 return std::pow(history_data[lstep] / history_data[0], 1. / lstep);
180 }
181
182
183
184 double
step_reduction(unsigned int step) const185 SolverControl::step_reduction(unsigned int step) const
186 {
187 Assert(history_data_enabled, ExcHistoryDataRequired());
188 Assert(history_data.size() > lstep, ExcInternalError());
189 Assert(step <= lstep, ExcIndexRange(step, 1, lstep + 1));
190 Assert(step > 0, ExcIndexRange(step, 1, lstep + 1));
191
192 return history_data[step] / history_data[step - 1];
193 }
194
195
196 double
final_reduction() const197 SolverControl::final_reduction() const
198 {
199 return step_reduction(lstep);
200 }
201
202
203 void
declare_parameters(ParameterHandler & param)204 SolverControl::declare_parameters(ParameterHandler ¶m)
205 {
206 param.declare_entry("Max steps", "100", Patterns::Integer());
207 param.declare_entry("Tolerance", "1.e-10", Patterns::Double());
208 param.declare_entry("Log history", "false", Patterns::Bool());
209 param.declare_entry("Log frequency", "1", Patterns::Integer());
210 param.declare_entry("Log result", "true", Patterns::Bool());
211 }
212
213
214 void
parse_parameters(ParameterHandler & param)215 SolverControl::parse_parameters(ParameterHandler ¶m)
216 {
217 set_max_steps(param.get_integer("Max steps"));
218 set_tolerance(param.get_double("Tolerance"));
219 log_history(param.get_bool("Log history"));
220 log_result(param.get_bool("Log result"));
221 log_frequency(param.get_integer("Log frequency"));
222 }
223
224 /*----------------------- ReductionControl ---------------------------------*/
225
226
ReductionControl(const unsigned int n,const double tol,const double red,const bool m_log_history,const bool m_log_result)227 ReductionControl::ReductionControl(const unsigned int n,
228 const double tol,
229 const double red,
230 const bool m_log_history,
231 const bool m_log_result)
232 : SolverControl(n, tol, m_log_history, m_log_result)
233 , reduce(red)
234 , reduced_tol(numbers::signaling_nan<double>())
235 {}
236
237
ReductionControl(const SolverControl & c)238 ReductionControl::ReductionControl(const SolverControl &c)
239 : SolverControl(c)
240 , reduce(numbers::signaling_nan<double>())
241 , reduced_tol(numbers::signaling_nan<double>())
242 {
243 set_reduction(0.);
244 }
245
246
247 ReductionControl &
operator =(const SolverControl & c)248 ReductionControl::operator=(const SolverControl &c)
249 {
250 SolverControl::operator=(c);
251 set_reduction(0.);
252 return *this;
253 }
254
255
256
257 SolverControl::State
check(const unsigned int step,const double check_value)258 ReductionControl::check(const unsigned int step, const double check_value)
259 {
260 // if this is the first time we
261 // come here, then store the
262 // residual for later comparisons
263 if (step == 0)
264 {
265 initial_val = check_value;
266 reduced_tol = check_value * reduce;
267 }
268
269 // check whether desired reduction
270 // has been achieved. also check
271 // for equality in case initial
272 // residual already was zero
273 if (check_value <= reduced_tol)
274 {
275 if (m_log_result)
276 deallog << "Convergence step " << step << " value " << check_value
277 << std::endl;
278 lstep = step;
279 lvalue = check_value;
280 lcheck = success;
281
282 if (history_data_enabled)
283 history_data.push_back(check_value);
284
285 return success;
286 }
287 else
288 return SolverControl::check(step, check_value);
289 }
290
291
292
293 void
declare_parameters(ParameterHandler & param)294 ReductionControl::declare_parameters(ParameterHandler ¶m)
295 {
296 SolverControl::declare_parameters(param);
297 param.declare_entry("Reduction", "1.e-2", Patterns::Double());
298 }
299
300
301 void
parse_parameters(ParameterHandler & param)302 ReductionControl::parse_parameters(ParameterHandler ¶m)
303 {
304 SolverControl::parse_parameters(param);
305 set_reduction(param.get_double("Reduction"));
306 }
307
308
309 /*---------------------- IterationNumberControl -----------------------------*/
310
311
IterationNumberControl(const unsigned int n,const double tolerance,const bool m_log_history,const bool m_log_result)312 IterationNumberControl::IterationNumberControl(const unsigned int n,
313 const double tolerance,
314 const bool m_log_history,
315 const bool m_log_result)
316 : SolverControl(n, tolerance, m_log_history, m_log_result)
317 {}
318
319
320
321 SolverControl::State
check(const unsigned int step,const double check_value)322 IterationNumberControl::check(const unsigned int step, const double check_value)
323 {
324 // check whether the given number of iterations was reached, and return
325 // success in that case. Otherwise, go on to the check of the base class.
326 if (step >= this->maxsteps)
327 {
328 if (m_log_result)
329 deallog << "Convergence step " << step << " value " << check_value
330 << std::endl;
331 lstep = step;
332 lvalue = check_value;
333
334 lcheck = success;
335 return success;
336 }
337 else
338 return SolverControl::check(step, check_value);
339 }
340
341 /*------------------------ ConsecutiveControl -------------------------------*/
342
343
ConsecutiveControl(const unsigned int n,const double tolerance,const unsigned int n_consecutive_iterations,const bool m_log_history,const bool m_log_result)344 ConsecutiveControl::ConsecutiveControl(
345 const unsigned int n,
346 const double tolerance,
347 const unsigned int n_consecutive_iterations,
348 const bool m_log_history,
349 const bool m_log_result)
350 : SolverControl(n, tolerance, m_log_history, m_log_result)
351 , n_consecutive_iterations(n_consecutive_iterations)
352 , n_converged_iterations(0)
353 {
354 AssertThrow(n_consecutive_iterations > 0,
355 ExcMessage("n_consecutive_iterations should be positive"));
356 }
357
358
359
ConsecutiveControl(const SolverControl & c)360 ConsecutiveControl::ConsecutiveControl(const SolverControl &c)
361 : SolverControl(c)
362 , n_consecutive_iterations(1)
363 , n_converged_iterations(0)
364 {}
365
366
367
368 ConsecutiveControl &
operator =(const SolverControl & c)369 ConsecutiveControl::operator=(const SolverControl &c)
370 {
371 SolverControl::operator =(c);
372 n_consecutive_iterations = 1;
373 n_converged_iterations = 0;
374 return *this;
375 }
376
377
378
379 SolverControl::State
check(const unsigned int step,const double check_value)380 ConsecutiveControl::check(const unsigned int step, const double check_value)
381 {
382 // reset the counter if ConsecutiveControl is being reused
383 if (step == 0)
384 n_converged_iterations = 0;
385 else
386 {
387 // check two things:
388 // (i) steps are ascending without repetitions
389 // (ii) user started from zero even when solver is being reused.
390 Assert(step - 1 == lstep,
391 ExcMessage("steps should be ascending integers."));
392 }
393
394 SolverControl::State state = SolverControl::check(step, check_value);
395 // check if we need to override the success:
396 if (state == success)
397 {
398 n_converged_iterations++;
399 if (n_converged_iterations == n_consecutive_iterations)
400 {
401 return success;
402 }
403 else
404 {
405 lcheck = iterate;
406 return iterate;
407 }
408 }
409 else
410 {
411 n_converged_iterations = 0;
412 return state;
413 }
414 }
415
416 DEAL_II_NAMESPACE_CLOSE
417