OGS
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 
37  INFO("Calculate non-equilibrium initial residuum.");
38 
41  _equation_system->assemble(x, x_prev, process_id);
42  _equation_system->getA(A);
43  _equation_system->getRhs(*x_prev[process_id], rhs);
44 
45  // r_neq = A * x - rhs
47  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
48  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
49 
50  // Set the values of the selected entries of _r_neq, which are associated
51  // with the equations that do not need initial residual compensation, to
52  // zero.
53  const auto selected_global_indices =
54  _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
55 
56  std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
57  _r_neq->set(selected_global_indices, zero_entries);
58 
60 
63 }
64 
66  std::vector<GlobalVector*>& x,
67  std::vector<GlobalVector*> const& x_prev,
68  std::function<void(int, std::vector<GlobalVector*> const&)> const&
69  postIterationCallback,
70  int const process_id)
71 {
72  namespace LinAlg = MathLib::LinAlg;
73  auto& sys = *_equation_system;
74 
77 
78  std::vector<GlobalVector*> x_new{x};
79  x_new[process_id] =
81  LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
82 
83  bool error_norms_met = false;
84 
85  _convergence_criterion->preFirstIteration();
86 
87  int iteration = 1;
88  for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
89  {
90  BaseLib::RunTime timer_dirichlet;
91  double time_dirichlet = 0.0;
92 
93  BaseLib::RunTime time_iteration;
94  time_iteration.start();
95 
96  timer_dirichlet.start();
97  sys.computeKnownSolutions(*x_new[process_id], process_id);
98  sys.applyKnownSolutions(*x_new[process_id]);
99  time_dirichlet += timer_dirichlet.elapsed();
100 
101  sys.preIteration(iteration, *x_new[process_id]);
102 
103  BaseLib::RunTime time_assembly;
104  time_assembly.start();
105  sys.assemble(x_new, x_prev, process_id);
106  sys.getA(A);
107  sys.getRhs(*x_prev[process_id], rhs);
108  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
109 
110  // Subtract non-equilibrium initial residuum if set
111  if (_r_neq != nullptr)
112  {
113  LinAlg::axpy(rhs, -1, *_r_neq);
114  }
115 
116  timer_dirichlet.start();
117  sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
118  time_dirichlet += timer_dirichlet.elapsed();
119  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
120 
121  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
122  {
123  GlobalVector res;
124  LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
125  LinAlg::axpy(res, -1.0, rhs); // res -= rhs
126  _convergence_criterion->checkResidual(res);
127  }
128 
129  BaseLib::RunTime time_linear_solver;
130  time_linear_solver.start();
131  bool iteration_succeeded =
132  _linear_solver.solve(A, rhs, *x_new[process_id]);
133  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
134 
135  if (!iteration_succeeded)
136  {
137  ERR("Picard: The linear solver failed.");
138  }
139  else
140  {
141  if (postIterationCallback)
142  {
143  postIterationCallback(iteration, x_new);
144  }
145 
146  switch (sys.postIteration(*x_new[process_id]))
147  {
149  // Don't copy here. The old x might still be used further
150  // below. Although currently it is not.
151  break;
153  ERR("Picard: The postIteration() hook reported a "
154  "non-recoverable error.");
155  iteration_succeeded = false;
156  // Copy new solution to x.
157  // Thereby the failed solution can be used by the caller for
158  // debugging purposes.
159  LinAlg::copy(*x_new[process_id], *x[process_id]);
160  break;
162  INFO(
163  "Picard: The postIteration() hook decided that this "
164  "iteration has to be repeated.");
165  LinAlg::copy(
166  *x[process_id],
167  *x_new[process_id]); // throw the iteration result away
168  continue;
169  }
170  }
171 
172  if (!iteration_succeeded)
173  {
174  // Don't compute error norms, break here.
175  error_norms_met = false;
176  break;
177  }
178 
179  if (sys.isLinear())
180  {
181  error_norms_met = true;
182  }
183  else
184  {
185  if (_convergence_criterion->hasDeltaXCheck())
186  {
187  GlobalVector minus_delta_x(*x[process_id]);
188  LinAlg::axpy(minus_delta_x, -1.0,
189  *x_new[process_id]); // minus_delta_x = x - x_new
190  _convergence_criterion->checkDeltaX(minus_delta_x,
191  *x_new[process_id]);
192  }
193 
194  error_norms_met = _convergence_criterion->isSatisfied();
195  }
196 
197  // Update x s.t. in the next iteration we will compute the right delta x
198  LinAlg::copy(*x_new[process_id], *x[process_id]);
199 
200  INFO("[time] Iteration #{:d} took {:g} s.", iteration,
201  time_iteration.elapsed());
202 
203  if (error_norms_met)
204  {
205  break;
206  }
207 
208  // Avoid increment of the 'iteration' if the error norms are not met,
209  // but maximum number of iterations is reached.
210  if (iteration >= _maxiter)
211  {
212  break;
213  }
214  }
215 
216  if (iteration > _maxiter)
217  {
218  ERR("Picard: Could not solve the given nonlinear system within {:d} "
219  "iterations",
220  _maxiter);
221  }
222 
226 
227  return {error_norms_met, iteration};
228 }
229 
232  std::vector<GlobalVector*> const& x,
233  std::vector<GlobalVector*> const& x_prev, int const process_id)
234 {
235  if (!_compensate_non_equilibrium_initial_residuum)
236  {
237  return;
238  }
239 
240  INFO("Calculate non-equilibrium initial residuum.");
241 
242  _equation_system->assemble(x, x_prev, process_id);
244  _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
245 
246  // Set the values of the selected entries of _r_neq, which are associated
247  // with the equations that do not need initial residual compensation, to
248  // zero.
249  const auto selected_global_indices =
250  _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
251  std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
252 
253  _r_neq->set(selected_global_indices, zero_entries);
254 
256 }
257 
259  std::vector<GlobalVector*>& x,
260  std::vector<GlobalVector*> const& x_prev,
261  std::function<void(int, std::vector<GlobalVector*> const&)> const&
262  postIterationCallback,
263  int const process_id)
264 {
265  namespace LinAlg = MathLib::LinAlg;
266  auto& sys = *_equation_system;
267 
269  auto& minus_delta_x =
272 
273  bool error_norms_met = false;
274 
275  // TODO be more efficient
276  // init minus_delta_x to the right size
277  LinAlg::copy(*x[process_id], minus_delta_x);
278 
279  _convergence_criterion->preFirstIteration();
280 
281  int iteration = 1;
282  for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
283  {
284  BaseLib::RunTime timer_dirichlet;
285  double time_dirichlet = 0.0;
286 
287  BaseLib::RunTime time_iteration;
288  time_iteration.start();
289 
290  timer_dirichlet.start();
291  sys.computeKnownSolutions(*x[process_id], process_id);
292  sys.applyKnownSolutions(*x[process_id]);
293  time_dirichlet += timer_dirichlet.elapsed();
294 
295  sys.preIteration(iteration, *x[process_id]);
296 
297  BaseLib::RunTime time_assembly;
298  time_assembly.start();
299  try
300  {
301  sys.assemble(x, x_prev, process_id);
302  }
303  catch (AssemblyException const& e)
304  {
305  ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
306  e.what());
307  error_norms_met = false;
308  iteration = _maxiter;
309  break;
310  }
311  sys.getResidual(*x[process_id], *x_prev[process_id], res);
312  sys.getJacobian(J);
313  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
314 
315  // Subtract non-equilibrium initial residuum if set
316  if (_r_neq != nullptr)
317  LinAlg::axpy(res, -1, *_r_neq);
318 
319  minus_delta_x.setZero();
320 
321  timer_dirichlet.start();
322  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
323  time_dirichlet += timer_dirichlet.elapsed();
324  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
325 
326  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
327  {
328  _convergence_criterion->checkResidual(res);
329  }
330 
331  BaseLib::RunTime time_linear_solver;
332  time_linear_solver.start();
333  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
334  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
335 
336  if (!iteration_succeeded)
337  {
338  ERR("Newton: The linear solver failed.");
339  }
340  else
341  {
342  // TODO could be solved in a better way
343  // cf.
344  // https://www.petsc.org/release/docs/manualpages/Vec/VecWAXPY.html
345 
346  // Copy pointers, replace the one for the given process id.
347  std::vector<GlobalVector*> x_new{x};
348  x_new[process_id] =
350  *x[process_id], _x_new_id);
351  LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
352 
353  if (postIterationCallback)
354  {
355  postIterationCallback(iteration, x_new);
356  }
357 
358  switch (sys.postIteration(*x_new[process_id]))
359  {
361  break;
363  ERR("Newton: The postIteration() hook reported a "
364  "non-recoverable error.");
365  iteration_succeeded = false;
366  break;
368  INFO(
369  "Newton: The postIteration() hook decided that this "
370  "iteration has to be repeated.");
371  // TODO introduce some onDestroy hook.
373  *x_new[process_id]);
374  continue; // That throws the iteration result away.
375  }
376 
377  LinAlg::copy(*x_new[process_id],
378  *x[process_id]); // copy new solution to x
380  *x_new[process_id]);
381  }
382 
383  if (!iteration_succeeded)
384  {
385  // Don't compute further error norms, but break here.
386  error_norms_met = false;
387  break;
388  }
389 
390  if (sys.isLinear())
391  {
392  error_norms_met = true;
393  }
394  else
395  {
396  if (_convergence_criterion->hasDeltaXCheck())
397  {
398  // Note: x contains the new solution!
399  _convergence_criterion->checkDeltaX(minus_delta_x,
400  *x[process_id]);
401  }
402 
403  error_norms_met = _convergence_criterion->isSatisfied();
404  }
405 
406  INFO("[time] Iteration #{:d} took {:g} s.", iteration,
407  time_iteration.elapsed());
408 
409  if (error_norms_met)
410  {
411  break;
412  }
413 
414  // Avoid increment of the 'iteration' if the error norms are not met,
415  // but maximum number of iterations is reached.
416  if (iteration >= _maxiter)
417  {
418  break;
419  }
420  }
421 
422  if (iteration > _maxiter)
423  {
424  ERR("Newton: Could not solve the given nonlinear system within {:d} "
425  "iterations",
426  _maxiter);
427  }
428 
432 
433  return {error_norms_met, iteration};
434 }
435 
436 std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
438  BaseLib::ConfigTree const& config)
439 {
441  auto const type = config.getConfigParameter<std::string>("type");
443  auto const max_iter = config.getConfigParameter<int>("max_iter");
444 
445  if (type == "Picard")
446  {
447  auto const tag = NonlinearSolverTag::Picard;
448  using ConcreteNLS = NonlinearSolver<tag>;
449  return std::make_pair(
450  std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
451  }
452  if (type == "Newton")
453  {
455  auto const damping = config.getConfigParameter<double>("damping", 1.0);
456  if (damping <= 0)
457  {
458  OGS_FATAL(
459  "The damping factor for the Newon method must be positive, got "
460  "{:g}.",
461  damping);
462  }
463  auto const tag = NonlinearSolverTag::Newton;
464  using ConcreteNLS = NonlinearSolver<tag>;
465  return std::make_pair(
466  std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
467  tag);
468  }
469 #ifdef USE_PETSC
470  if (boost::iequals(type, "PETScSNES"))
471  {
472  auto prefix =
474  config.getConfigParameter<std::string>("prefix", "");
475  auto const tag = NonlinearSolverTag::Newton;
476  using ConcreteNLS = PETScNonlinearSolver;
477  return std::make_pair(std::make_unique<ConcreteNLS>(
478  linear_solver, max_iter, std::move(prefix)),
479  tag);
480  }
481 
482 #endif
483  OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
484 }
485 
487 {
488  if (_r_neq != nullptr)
489  {
491  }
492 }
493 
495 {
496  if (_r_neq != nullptr)
497  {
499  }
500 }
501 
502 } // namespace NumLib
#define OGS_FATAL(...)
Definition: Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:34
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:44
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:28
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 finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:163
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:37
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:142
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition: LinAlg.cpp:57
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:206
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.