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