OGS  master
NonlinearSolver.cpp
Go to the documentation of this file.
1 
11 #include "NonlinearSolver.h"
12 
13 #include <boost/algorithm/string.hpp>
14 
15 #include "BaseLib/ConfigTree.h"
16 #include "BaseLib/Error.h"
17 #include "BaseLib/Logging.h"
18 #include "BaseLib/RunTime.h"
19 #include "ConvergenceCriterion.h"
20 #include "MathLib/LinAlg/LinAlg.h"
22 #include "NumLib/Exceptions.h"
23 #include "PETScNonlinearSolver.h"
24 
25 namespace NumLib
26 {
29  std::vector<GlobalVector*> const& x,
30  std::vector<GlobalVector*> const& x_prev, int const process_id)
31 {
32  if (!_compensate_non_equilibrium_initial_residuum)
33  {
34  return;
35  }
36 
39  _equation_system->assemble(x, x_prev, process_id);
40  _equation_system->getA(A);
41  _equation_system->getRhs(*x_prev[process_id], rhs);
42 
43  // r_neq = A * x - rhs
45  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
46  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
47 }
48 
50  std::vector<GlobalVector*>& x,
51  std::vector<GlobalVector*> const& x_prev,
52  std::function<void(int, std::vector<GlobalVector*> const&)> const&
53  postIterationCallback,
54  int const process_id)
55 {
56  namespace LinAlg = MathLib::LinAlg;
57  auto& sys = *_equation_system;
58 
61 
62  std::vector<GlobalVector*> x_new{x};
63  x_new[process_id] =
65  LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
66 
67  bool error_norms_met = false;
68 
69  _convergence_criterion->preFirstIteration();
70 
71  int iteration = 1;
72  for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
73  {
74  BaseLib::RunTime timer_dirichlet;
75  double time_dirichlet = 0.0;
76 
77  BaseLib::RunTime time_iteration;
78  time_iteration.start();
79 
80  timer_dirichlet.start();
81  sys.computeKnownSolutions(*x_new[process_id], process_id);
82  sys.applyKnownSolutions(*x_new[process_id]);
83  time_dirichlet += timer_dirichlet.elapsed();
84 
85  sys.preIteration(iteration, *x_new[process_id]);
86 
87  BaseLib::RunTime time_assembly;
88  time_assembly.start();
89  sys.assemble(x_new, x_prev, process_id);
90  sys.getA(A);
91  sys.getRhs(*x_prev[process_id], rhs);
92  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
93 
94  // Subtract non-equilibrium initial residuum if set
95  if (_r_neq != nullptr)
96  {
97  LinAlg::axpy(rhs, -1, *_r_neq);
98  }
99 
100  timer_dirichlet.start();
101  sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
102  time_dirichlet += timer_dirichlet.elapsed();
103  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
104 
105  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
106  {
107  GlobalVector res;
108  LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
109  LinAlg::axpy(res, -1.0, rhs); // res -= rhs
110  _convergence_criterion->checkResidual(res);
111  }
112 
113  BaseLib::RunTime time_linear_solver;
114  time_linear_solver.start();
115  bool iteration_succeeded =
116  _linear_solver.solve(A, rhs, *x_new[process_id]);
117  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
118 
119  if (!iteration_succeeded)
120  {
121  ERR("Picard: The linear solver failed.");
122  }
123  else
124  {
125  if (postIterationCallback)
126  {
127  postIterationCallback(iteration, x_new);
128  }
129 
130  switch (sys.postIteration(*x_new[process_id]))
131  {
133  // Don't copy here. The old x might still be used further
134  // below. Although currently it is not.
135  break;
137  ERR("Picard: The postIteration() hook reported a "
138  "non-recoverable error.");
139  iteration_succeeded = false;
140  // Copy new solution to x.
141  // Thereby the failed solution can be used by the caller for
142  // debugging purposes.
143  LinAlg::copy(*x_new[process_id], *x[process_id]);
144  break;
146  INFO(
147  "Picard: The postIteration() hook decided that this "
148  "iteration has to be repeated.");
149  LinAlg::copy(
150  *x[process_id],
151  *x_new[process_id]); // throw the iteration result away
152  continue;
153  }
154  }
155 
156  if (!iteration_succeeded)
157  {
158  // Don't compute error norms, break here.
159  error_norms_met = false;
160  break;
161  }
162 
163  if (sys.isLinear())
164  {
165  error_norms_met = true;
166  }
167  else
168  {
169  if (_convergence_criterion->hasDeltaXCheck())
170  {
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 #{:d} 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 {:d} "
203  "iterations",
204  _maxiter);
205  }
206 
210 
211  return {error_norms_met, iteration};
212 }
213 
216  std::vector<GlobalVector*> const& x,
217  std::vector<GlobalVector*> const& x_prev, int const process_id)
218 {
219  if (!_compensate_non_equilibrium_initial_residuum)
220  {
221  return;
222  }
223 
224  _equation_system->assemble(x, x_prev, process_id);
226  _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
227 }
228 
230  std::vector<GlobalVector*>& x,
231  std::vector<GlobalVector*> const& x_prev,
232  std::function<void(int, std::vector<GlobalVector*> const&)> const&
233  postIterationCallback,
234  int const process_id)
235 {
236  namespace LinAlg = MathLib::LinAlg;
237  auto& sys = *_equation_system;
238 
240  auto& minus_delta_x =
243 
244  bool error_norms_met = false;
245 
246  // TODO be more efficient
247  // init minus_delta_x to the right size
248  LinAlg::copy(*x[process_id], minus_delta_x);
249 
250  _convergence_criterion->preFirstIteration();
251 
252  int iteration = 1;
253  for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
254  {
255  BaseLib::RunTime timer_dirichlet;
256  double time_dirichlet = 0.0;
257 
258  BaseLib::RunTime time_iteration;
259  time_iteration.start();
260 
261  timer_dirichlet.start();
262  sys.computeKnownSolutions(*x[process_id], process_id);
263  sys.applyKnownSolutions(*x[process_id]);
264  time_dirichlet += timer_dirichlet.elapsed();
265 
266  sys.preIteration(iteration, *x[process_id]);
267 
268  BaseLib::RunTime time_assembly;
269  time_assembly.start();
270  try
271  {
272  sys.assemble(x, x_prev, process_id);
273  }
274  catch (AssemblyException const& e)
275  {
276  ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
277  e.what());
278  error_norms_met = false;
279  iteration = _maxiter;
280  break;
281  }
282  sys.getResidual(*x[process_id], *x_prev[process_id], res);
283  sys.getJacobian(J);
284  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
285 
286  // Subtract non-equilibrium initial residuum if set
287  if (_r_neq != nullptr)
288  LinAlg::axpy(res, -1, *_r_neq);
289 
290  minus_delta_x.setZero();
291 
292  timer_dirichlet.start();
293  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
294  time_dirichlet += timer_dirichlet.elapsed();
295  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
296 
297  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
298  {
299  _convergence_criterion->checkResidual(res);
300  }
301 
302  BaseLib::RunTime time_linear_solver;
303  time_linear_solver.start();
304  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
305  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
306 
307  if (!iteration_succeeded)
308  {
309  ERR("Newton: The linear solver failed.");
310  }
311  else
312  {
313  // TODO could be solved in a better way
314  // cf.
315  // https://www.petsc.org/release/docs/manualpages/Vec/VecWAXPY.html
316 
317  // Copy pointers, replace the one for the given process id.
318  std::vector<GlobalVector*> x_new{x};
319  x_new[process_id] =
321  *x[process_id], _x_new_id);
322  LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
323 
324  if (postIterationCallback)
325  {
326  postIterationCallback(iteration, x_new);
327  }
328 
329  switch (sys.postIteration(*x_new[process_id]))
330  {
332  break;
334  ERR("Newton: The postIteration() hook reported a "
335  "non-recoverable error.");
336  iteration_succeeded = false;
337  break;
339  INFO(
340  "Newton: The postIteration() hook decided that this "
341  "iteration has to be repeated.");
342  // TODO introduce some onDestroy hook.
344  *x_new[process_id]);
345  continue; // That throws the iteration result away.
346  }
347 
348  LinAlg::copy(*x_new[process_id],
349  *x[process_id]); // copy new solution to x
351  *x_new[process_id]);
352  }
353 
354  if (!iteration_succeeded)
355  {
356  // Don't compute further error norms, but break here.
357  error_norms_met = false;
358  break;
359  }
360 
361  if (sys.isLinear())
362  {
363  error_norms_met = true;
364  }
365  else
366  {
367  if (_convergence_criterion->hasDeltaXCheck())
368  {
369  // Note: x contains the new solution!
370  _convergence_criterion->checkDeltaX(minus_delta_x,
371  *x[process_id]);
372  }
373 
374  error_norms_met = _convergence_criterion->isSatisfied();
375  }
376 
377  INFO("[time] Iteration #{:d} took {:g} s.", iteration,
378  time_iteration.elapsed());
379 
380  if (error_norms_met)
381  {
382  break;
383  }
384 
385  // Avoid increment of the 'iteration' if the error norms are not met,
386  // but maximum number of iterations is reached.
387  if (iteration >= _maxiter)
388  {
389  break;
390  }
391  }
392 
393  if (iteration > _maxiter)
394  {
395  ERR("Newton: Could not solve the given nonlinear system within {:d} "
396  "iterations",
397  _maxiter);
398  }
399 
403 
404  return {error_norms_met, iteration};
405 }
406 
407 std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
409  BaseLib::ConfigTree const& config)
410 {
412  auto const type = config.getConfigParameter<std::string>("type");
414  auto const max_iter = config.getConfigParameter<int>("max_iter");
415 
416  if (type == "Picard")
417  {
418  auto const tag = NonlinearSolverTag::Picard;
419  using ConcreteNLS = NonlinearSolver<tag>;
420  return std::make_pair(
421  std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
422  }
423  if (type == "Newton")
424  {
426  auto const damping = config.getConfigParameter<double>("damping", 1.0);
427  if (damping <= 0)
428  {
429  OGS_FATAL(
430  "The damping factor for the Newon method must be positive, got "
431  "{:g}.",
432  damping);
433  }
434  auto const tag = NonlinearSolverTag::Newton;
435  using ConcreteNLS = NonlinearSolver<tag>;
436  return std::make_pair(
437  std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
438  tag);
439  }
440 #ifdef USE_PETSC
441  if (boost::iequals(type, "PETScSNES"))
442  {
443  auto prefix =
445  config.getConfigParameter<std::string>("prefix", "");
446  auto const tag = NonlinearSolverTag::Newton;
447  using ConcreteNLS = PETScNonlinearSolver;
448  return std::make_pair(std::make_unique<ConcreteNLS>(
449  linear_solver, max_iter, std::move(prefix)),
450  tag);
451  }
452 
453 #endif
454  OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
455 }
456 
458 {
459  if (_r_neq != nullptr)
460  {
462  }
463 }
464 
466 {
467  if (_r_neq != nullptr)
468  {
470  }
471 }
472 
473 } // namespace NumLib
#define OGS_FATAL(...)
Definition: Error.h:25
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
Definition of the RunTime class.
T getConfigParameter(std::string const &param) const
Count the running time.
Definition: RunTime.h:29
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:42
void start()
Start the timer.
Definition: RunTime.h:32
Global vector based on Eigen vector.
Definition: EigenVector.h:27
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
std::pair< std::unique_ptr< NonlinearSolverBase >, NonlinearSolverTag > createNonlinearSolver(GlobalLinearSolver &linear_solver, BaseLib::ConfigTree const &config)
void axpy(PETScVector &y, double const a, PETScVector const &x)
Definition: LinAlg.cpp:57
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:141
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
Definition: TimeLoop.cpp:251
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.