OGS 6.3.0-179-g962fdcd4e.dirty.20200403132553
NonlinearSolver.cpp
Go to the documentation of this file.
1
11 #include "NonlinearSolver.h"
12
13 #include <logog/include/logog.hpp>
14
15 #include "BaseLib/ConfigTree.h"
16 #include "BaseLib/Error.h"
17 #include "BaseLib/RunTime.h"
18 #include "ConvergenceCriterion.h"
19 #include "MathLib/LinAlg/LinAlg.h"
21 #include "NumLib/Exceptions.h"
22
23 namespace NumLib
24 {
26  std::vector<GlobalVector*> const& x, int const process_id) const
27 {
28  _equation_system->assemble(x, process_id);
29 }
30
32  calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
33  int const process_id)
34 {
35  if (!_compensate_non_equilibrium_initial_residuum)
36  {
37  return;
38  }
39
42  _equation_system->assemble(x, process_id);
43  _equation_system->getA(A);
44  _equation_system->getRhs(rhs);
45
46  // r_neq = A * x - rhs
48  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
49  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
50 }
51
53  std::vector<GlobalVector*>& x,
54  std::function<void(int, std::vector<GlobalVector*> const&)> const&
55  postIterationCallback,
56  int const process_id)
57 {
58  namespace LinAlg = MathLib::LinAlg;
59  auto& sys = *_equation_system;
60
61  auto& A =
64  _rhs_id);
65
66  std::vector<GlobalVector*> x_new{x};
67  x_new[process_id] =
69  LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
70
71  bool error_norms_met = false;
72
73  _convergence_criterion->preFirstIteration();
74
75  int iteration = 1;
76  for (; iteration <= _maxiter;
77  ++iteration, _convergence_criterion->reset())
78  {
79  BaseLib::RunTime timer_dirichlet;
80  double time_dirichlet = 0.0;
81
82  BaseLib::RunTime time_iteration;
83  time_iteration.start();
84
85  timer_dirichlet.start();
86  sys.computeKnownSolutions(*x_new[process_id], process_id);
87  sys.applyKnownSolutions(*x_new[process_id]);
88  time_dirichlet += timer_dirichlet.elapsed();
89
90  sys.preIteration(iteration, *x_new[process_id]);
91
92  BaseLib::RunTime time_assembly;
93  time_assembly.start();
94  sys.assemble(x_new, process_id);
95  sys.getA(A);
96  sys.getRhs(rhs);
97  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
98
99  // Subract non-equilibrium initial residuum if set
100  if (_r_neq != nullptr)
101  {
102  LinAlg::axpy(rhs, -1, *_r_neq);
103  }
104
105  timer_dirichlet.start();
106  sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
107  time_dirichlet += timer_dirichlet.elapsed();
108  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
109
110  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
111  GlobalVector res;
112  LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
113  LinAlg::axpy(res, -1.0, rhs); // res -= rhs
114  _convergence_criterion->checkResidual(res);
115  }
116
117  BaseLib::RunTime time_linear_solver;
118  time_linear_solver.start();
119  bool iteration_succeeded =
120  _linear_solver.solve(A, rhs, *x_new[process_id]);
121  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
122
123  if (!iteration_succeeded)
124  {
125  ERR("Picard: The linear solver failed.");
126  }
127  else
128  {
129  if (postIterationCallback)
130  {
131  postIterationCallback(iteration, x_new);
132  }
133
134  switch (sys.postIteration(*x_new[process_id]))
135  {
137  // Don't copy here. The old x might still be used further
138  // below. Although currently it is not.
139  break;
141  ERR("Picard: The postIteration() hook reported a "
142  "non-recoverable error.");
143  iteration_succeeded = false;
144  // Copy new solution to x.
145  // Thereby the failed solution can be used by the caller for
146  // debugging purposes.
147  LinAlg::copy(*x_new[process_id], *x[process_id]);
148  break;
150  INFO(
151  "Picard: The postIteration() hook decided that this "
152  "iteration has to be repeated.");
153  LinAlg::copy(
154  *x[process_id],
155  *x_new[process_id]); // throw the iteration result away
156  continue;
157  }
158  }
159
160  if (!iteration_succeeded)
161  {
162  // Don't compute error norms, break here.
163  error_norms_met = false;
164  break;
165  }
166
167  if (sys.isLinear()) {
168  error_norms_met = true;
169  } else {
170  if (_convergence_criterion->hasDeltaXCheck()) {
171  GlobalVector minus_delta_x(*x[process_id]);
172  LinAlg::axpy(minus_delta_x, -1.0,
173  *x_new[process_id]); // minus_delta_x = x - x_new
174  _convergence_criterion->checkDeltaX(minus_delta_x,
175  *x_new[process_id]);
176  }
177
178  error_norms_met = _convergence_criterion->isSatisfied();
179  }
180
181  // Update x s.t. in the next iteration we will compute the right delta x
182  LinAlg::copy(*x_new[process_id], *x[process_id]);
183
184  INFO("[time] Iteration #%u took %g s.", iteration,
185  time_iteration.elapsed());
186
187  if (error_norms_met)
188  {
189  break;
190  }
191
192  // Avoid increment of the 'iteration' if the error norms are not met,
193  // but maximum number of iterations is reached.
194  if (iteration >= _maxiter)
195  {
196  break;
197  }
198  }
199
200  if (iteration > _maxiter)
201  {
202  ERR("Picard: Could not solve the given nonlinear system within %u "
203  "iterations",
204  _maxiter);
205  }
206
210
211  return {error_norms_met, iteration};
212 }
213
215  std::vector<GlobalVector*> const& x, int const process_id) const
216 {
217  _equation_system->assemble(x, process_id);
218  // TODO if the equation system would be reset to nullptr after each
219  // assemble() or solve() call, the user would be forced to set the
220  // equation every time and could not forget it.
221 }
222
224  calculateNonEquilibriumInitialResiduum(std::vector<GlobalVector*> const& x,
225  int const process_id)
226 {
227  if (!_compensate_non_equilibrium_initial_residuum)
228  {
229  return;
230  }
231
232  _equation_system->assemble(x, process_id);
234  _equation_system->getResidual(*x[process_id], *_r_neq);
235 }
236
238  std::vector<GlobalVector*>& x,
239  std::function<void(int, std::vector<GlobalVector*> const&)> const&
240  postIterationCallback,
241  int const process_id)
242 {
243  namespace LinAlg = MathLib::LinAlg;
244  auto& sys = *_equation_system;
245
247  _res_id);
248  auto& minus_delta_x =
250  _minus_delta_x_id);
251  auto& J =
253
254  bool error_norms_met = false;
255
256  // TODO be more efficient
257  // init minus_delta_x to the right size
258  LinAlg::copy(*x[process_id], minus_delta_x);
259
260  _convergence_criterion->preFirstIteration();
261
262  int iteration = 1;
263  for (; iteration <= _maxiter;
264  ++iteration, _convergence_criterion->reset())
265  {
266  BaseLib::RunTime timer_dirichlet;
267  double time_dirichlet = 0.0;
268
269  BaseLib::RunTime time_iteration;
270  time_iteration.start();
271
272  timer_dirichlet.start();
273  sys.computeKnownSolutions(*x[process_id], process_id);
274  sys.applyKnownSolutions(*x[process_id]);
275  time_dirichlet += timer_dirichlet.elapsed();
276
277  sys.preIteration(iteration, *x[process_id]);
278
279  BaseLib::RunTime time_assembly;
280  time_assembly.start();
281  try
282  {
283  sys.assemble(x, process_id);
284  }
285  catch (AssemblyException const& e)
286  {
287  ERR("Abort nonlinear iteration. Repeating timestep. Reason: %s",
288  e.what());
289  error_norms_met = false;
290  iteration = _maxiter;
291  break;
292  }
293  sys.getResidual(*x[process_id], res);
294  sys.getJacobian(J);
295  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
296
297  // Subract non-equilibrium initial residuum if set
298  if (_r_neq != nullptr)
299  LinAlg::axpy(res, -1, *_r_neq);
300
301  minus_delta_x.setZero();
302
303  timer_dirichlet.start();
304  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
305  time_dirichlet += timer_dirichlet.elapsed();
306  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
307
308  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
309  {
310  _convergence_criterion->checkResidual(res);
311  }
312
313  BaseLib::RunTime time_linear_solver;
314  time_linear_solver.start();
315  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
316  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
317
318  if (!iteration_succeeded)
319  {
320  ERR("Newton: The linear solver failed.");
321  }
322  else
323  {
324  // TODO could be solved in a better way
325  // cf.
326  // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
327
328  // Copy pointers, replace the one for the given process id.
329  std::vector<GlobalVector*> x_new{x};
330  x_new[process_id] =
332  *x[process_id], _x_new_id);
333  LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
334
335  if (postIterationCallback)
336  {
337  postIterationCallback(iteration, x_new);
338  }
339
340  switch (sys.postIteration(*x_new[process_id]))
341  {
343  break;
345  ERR("Newton: The postIteration() hook reported a "
346  "non-recoverable error.");
347  iteration_succeeded = false;
348  break;
350  INFO(
351  "Newton: The postIteration() hook decided that this "
352  "iteration"
353  " has to be repeated.");
354  // TODO introduce some onDestroy hook.
356  *x_new[process_id]);
357  continue; // That throws the iteration result away.
358  }
359
360  LinAlg::copy(*x_new[process_id],
361  *x[process_id]); // copy new solution to x
363  *x_new[process_id]);
364  }
365
366  if (!iteration_succeeded)
367  {
368  // Don't compute further error norms, but break here.
369  error_norms_met = false;
370  break;
371  }
372
373  if (sys.isLinear()) {
374  error_norms_met = true;
375  } else {
376  if (_convergence_criterion->hasDeltaXCheck()) {
377  // Note: x contains the new solution!
378  _convergence_criterion->checkDeltaX(minus_delta_x,
379  *x[process_id]);
380  }
381
382  error_norms_met = _convergence_criterion->isSatisfied();
383  }
384
385  INFO("[time] Iteration #%u took %g s.", iteration,
386  time_iteration.elapsed());
387
388  if (error_norms_met)
389  {
390  break;
391  }
392
393  // Avoid increment of the 'iteration' if the error norms are not met,
394  // but maximum number of iterations is reached.
395  if (iteration >= _maxiter)
396  {
397  break;
398  }
399  }
400
401  if (iteration > _maxiter)
402  {
403  ERR("Newton: Could not solve the given nonlinear system within %u "
404  "iterations",
405  _maxiter);
406  }
407
411  minus_delta_x);
412
413  return {error_norms_met, iteration};
414 }
415
416 std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
417 createNonlinearSolver(GlobalLinearSolver& linear_solver,
418  BaseLib::ConfigTree const& config)
419 {
421  auto const type = config.getConfigParameter<std::string>("type");
423  auto const max_iter = config.getConfigParameter<int>("max_iter");
424
425  if (type == "Picard") {
426  auto const tag = NonlinearSolverTag::Picard;
427  using ConcreteNLS = NonlinearSolver<tag>;
428  return std::make_pair(
429  std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
430  }
431  if (type == "Newton")
432  {
434  auto const damping = config.getConfigParameter<double>("damping", 1.0);
435  if (damping <= 0)
436  {
437  OGS_FATAL(
438  "The damping factor for the Newon method must be positive, got "
439  "%g.",
440  damping);
441  }
442  auto const tag = NonlinearSolverTag::Newton;
443  using ConcreteNLS = NonlinearSolver<tag>;
444  return std::make_pair(
445  std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
446  tag);
447  }
448  OGS_FATAL("Unsupported nonlinear solver type");
449 }
450 } // namespace NumLib
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
virtual void releaseMatrix(GlobalMatrix const &A)=0
std::pair< std::unique_ptr< NonlinearSolverBase >, NonlinearSolverTag > createNonlinearSolver(GlobalLinearSolver &linear_solver, BaseLib::ConfigTree const &config)
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
T getConfigParameter(std::string const &param) const
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Definition of the RunTime class.
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
Status of the non-linear solver.
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
Count the running time.
Definition: RunTime.h:27
void start()
Start the timer.
Definition: RunTime.h:31
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector *> process_solutions)
Definition: TimeLoop.cpp:170
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
virtual void releaseVector(GlobalVector const &x)=0
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:41