OGS
NonlinearSolver.cpp
Go to the documentation of this file.
1
11#include "NonlinearSolver.h"
12
13#include <algorithm>
14#include <boost/algorithm/string.hpp>
15
16#include "BaseLib/ConfigTree.h"
17#include "BaseLib/Error.h"
18#include "BaseLib/Logging.h"
19#include "BaseLib/MPI.h"
20#include "BaseLib/RunTime.h"
24#include "NumLib/Exceptions.h"
26
27namespace NumLib
28{
29namespace detail
30{
31#if !defined(USE_PETSC) && !defined(USE_LIS)
32bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
34 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
35{
36 BaseLib::RunTime time_linear_solver;
37 time_linear_solver.start();
38
39 if (!linear_solver.compute(A, linear_solver_behaviour))
40 {
41 ERR("Picard: The linear solver failed in the compute() step.");
42 return false;
43 }
44
45 bool const iteration_succeeded = linear_solver.solve(rhs, x);
46
47 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
48
49 if (iteration_succeeded)
50 {
51 return true;
52 }
53
54 ERR("Picard: The linear solver failed in the solve() step.");
55 return false;
56}
57#else
60 MathLib::LinearSolverBehaviour const linear_solver_behaviour)
61{
62 if (linear_solver_behaviour ==
64 linear_solver_behaviour == MathLib::LinearSolverBehaviour::REUSE)
65 {
66 WARN(
67 "The performance optimization to skip the linear solver compute() "
68 "step is not implemented for PETSc or LIS linear solvers.");
69 }
70
71 BaseLib::RunTime time_linear_solver;
72 time_linear_solver.start();
73
74 bool const iteration_succeeded = linear_solver.solve(A, rhs, x);
75
76 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
77
78 if (iteration_succeeded)
79 {
80 return true;
81 }
82
83 ERR("Picard: The linear solver failed in the solve() step.");
84 return false;
85}
86#endif
87} // namespace detail
88
91 std::vector<GlobalVector*> const& x,
92 std::vector<GlobalVector*> const& x_prev, int const process_id)
93{
94 if (!_compensate_non_equilibrium_initial_residuum)
95 {
96 return;
97 }
98
99 INFO("Calculate non-equilibrium initial residuum.");
100
103 _equation_system->assemble(x, x_prev, process_id);
104 _equation_system->getA(A);
105 _equation_system->getRhs(*x_prev[process_id], rhs);
106
107 // r_neq = A * x - rhs
109 MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
110 MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
111
112 // Set the values of the selected entries of _r_neq, which are associated
113 // with the equations that do not need initial residual compensation, to
114 // zero.
115 const auto selected_global_indices =
116 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
117
118 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
119 _r_neq->set(selected_global_indices, zero_entries);
120
122
125}
126
128 std::vector<GlobalVector*>& x,
129 std::vector<GlobalVector*> const& x_prev,
130 std::function<void(int, bool, std::vector<GlobalVector*> const&)> const&
131 postIterationCallback,
132 int const process_id)
133{
134 namespace LinAlg = MathLib::LinAlg;
135 auto& sys = *_equation_system;
136
139
140 std::vector<GlobalVector*> x_new{x};
141 x_new[process_id] =
143 LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
144
145 bool error_norms_met = false;
146
147 _convergence_criterion->preFirstIteration();
148
149 int iteration = 1;
150 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
151 {
152 BaseLib::RunTime timer_dirichlet;
153 double time_dirichlet = 0.0;
154
155 BaseLib::RunTime time_iteration;
156 time_iteration.start();
157
158 INFO("Iteration #{:d} started.", iteration);
159 timer_dirichlet.start();
160 auto& x_new_process = *x_new[process_id];
162 sys.computeKnownSolutions(x_new_process, process_id);
163 sys.applyKnownSolutions(x_new_process);
164 time_dirichlet += timer_dirichlet.elapsed();
165
166 sys.preIteration(iteration, x_new_process);
167
168 BaseLib::RunTime time_assembly;
169 time_assembly.start();
170 sys.assemble(x_new, x_prev, process_id);
171 sys.getA(A);
172 sys.getRhs(*x_prev[process_id], rhs);
173
174 // Normalize the linear equation system, if required
175 if (sys.requiresNormalization() &&
176 !_linear_solver.canSolveRectangular())
177 {
178 sys.getAandRhsNormalized(A, rhs);
179 WARN(
180 "The equation system is rectangular, but the current linear "
181 "solver only supports square systems. "
182 "The system will be normalized, which lead to a squared "
183 "condition number and potential numerical issues. "
184 "It is recommended to use a solver that supports rectangular "
185 "equation systems for better numerical stability.");
186 }
187
188 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
189
190 // Subtract non-equilibrium initial residuum if set
191 if (_r_neq != nullptr)
192 {
193 LinAlg::axpy(rhs, -1, *_r_neq);
194 }
195
196 timer_dirichlet.start();
197 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
198 time_dirichlet += timer_dirichlet.elapsed();
199 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
200
201 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
202 {
203 GlobalVector res;
204 LinAlg::matMult(A, x_new_process, res); // res = A * x_new
205 LinAlg::axpy(res, -1.0, rhs); // res -= rhs
206 _convergence_criterion->checkResidual(res);
207 }
208
209 bool iteration_succeeded =
210 detail::solvePicard(_linear_solver, A, rhs, x_new_process,
211 sys.linearSolverNeedsToCompute());
212
213 if (iteration_succeeded)
214 {
215 if (postIterationCallback)
216 {
217 postIterationCallback(iteration, error_norms_met, x_new);
218 }
219
220 switch (sys.postIteration(x_new_process))
221 {
223 // Don't copy here. The old x might still be used further
224 // below. Although currently it is not.
225 break;
227 ERR("Picard: The postIteration() hook reported a "
228 "non-recoverable error.");
229 iteration_succeeded = false;
230 // Copy new solution to x.
231 // Thereby the failed solution can be used by the caller for
232 // debugging purposes.
233 LinAlg::copy(x_new_process, *x[process_id]);
234 break;
236 INFO(
237 "Picard: The postIteration() hook decided that this "
238 "iteration has to be repeated.");
240 *x[process_id],
241 x_new_process); // throw the iteration result away
242 continue;
243 }
244 }
245
246 if (!iteration_succeeded)
247 {
248 // Don't compute error norms, break here.
249 error_norms_met = false;
250 break;
251 }
252
253 if (sys.isLinear())
254 {
255 error_norms_met = true;
256 }
257 else
258 {
259 if (_convergence_criterion->hasDeltaXCheck())
260 {
261 GlobalVector minus_delta_x(*x[process_id]);
262 LinAlg::axpy(minus_delta_x, -1.0,
263 x_new_process); // minus_delta_x = x - x_new
264 _convergence_criterion->checkDeltaX(minus_delta_x,
265 x_new_process);
266 }
267
268 error_norms_met = _convergence_criterion->isSatisfied();
269 }
270
271 // Update x s.t. in the next iteration we will compute the right delta x
272 LinAlg::copy(x_new_process, *x[process_id]);
273
274 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
275 time_iteration.elapsed());
276
277 if (error_norms_met)
278 {
279 break;
280 }
281
282 // Avoid increment of the 'iteration' if the error norms are not met,
283 // but maximum number of iterations is reached.
284 if (iteration >= _maxiter)
285 {
286 break;
287 }
288 }
289
290 if (iteration > _maxiter)
291 {
292 ERR("Picard: Could not solve the given nonlinear system within {:d} "
293 "iterations",
294 _maxiter);
295 }
296
300
301 return {error_norms_met, iteration};
302}
303
306 std::vector<GlobalVector*> const& x,
307 std::vector<GlobalVector*> const& x_prev, int const process_id)
308{
309 if (!_compensate_non_equilibrium_initial_residuum)
310 {
311 return;
312 }
313
314 INFO("Calculate non-equilibrium initial residuum.");
315
316 _equation_system->assemble(x, x_prev, process_id);
318 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
319
320 // Set the values of the selected entries of _r_neq, which are associated
321 // with the equations that do not need initial residual compensation, to
322 // zero.
323 const auto selected_global_indices =
324 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
325 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
326
327 _r_neq->set(selected_global_indices, zero_entries);
328
330}
331
333 std::vector<GlobalVector*>& x,
334 std::vector<GlobalVector*> const& x_prev,
335 std::function<void(int, bool, std::vector<GlobalVector*> const&)> const&
336 postIterationCallback,
337 int const process_id)
338{
339 namespace LinAlg = MathLib::LinAlg;
340 auto& sys = *_equation_system;
341
343 auto& minus_delta_x =
346
347 bool error_norms_met = false;
348
349 // TODO be more efficient
350 // init minus_delta_x to the right size
351 LinAlg::copy(*x[process_id], minus_delta_x);
352
353 _convergence_criterion->preFirstIteration();
354
355 int iteration = 1;
356#if !defined(USE_PETSC) && !defined(USE_LIS)
357 int next_iteration_inv_jacobian_recompute = 1;
358#endif
359 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
360 {
361 BaseLib::RunTime timer_dirichlet;
362 double time_dirichlet = 0.0;
363
364 BaseLib::RunTime time_iteration;
365 INFO("Iteration #{:d} started.", iteration);
366 time_iteration.start();
367
368 timer_dirichlet.start();
369 sys.computeKnownSolutions(*x[process_id], process_id);
370 time_dirichlet += timer_dirichlet.elapsed();
371
372 sys.preIteration(iteration, *x[process_id]);
373
374 BaseLib::RunTime time_assembly;
375 time_assembly.start();
376 bool mpi_rank_assembly_ok = true;
377 try
378 {
379 sys.assemble(x, x_prev, process_id);
380 }
381 catch (AssemblyException const& e)
382 {
383 ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
384 e.what());
385 error_norms_met = false;
386 iteration = _maxiter;
387 mpi_rank_assembly_ok = false;
388 }
389 if (BaseLib::MPI::anyOf(!mpi_rank_assembly_ok))
390 {
391 break;
392 }
393 sys.getResidual(*x[process_id], *x_prev[process_id], res);
394 sys.getJacobian(J);
395 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
396
397 // Subtract non-equilibrium initial residuum if set
398 if (_r_neq != nullptr)
399 LinAlg::axpy(res, -1, *_r_neq);
400
401 minus_delta_x.setZero();
402
403 timer_dirichlet.start();
404 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
405 time_dirichlet += timer_dirichlet.elapsed();
406 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
407
408 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
409 {
410 _convergence_criterion->checkResidual(res);
411 }
412
413 BaseLib::RunTime time_linear_solver;
414 time_linear_solver.start();
415#if !defined(USE_PETSC) && !defined(USE_LIS)
416 auto linear_solver_behaviour = MathLib::LinearSolverBehaviour::REUSE;
417 if (iteration == next_iteration_inv_jacobian_recompute)
418 {
419 linear_solver_behaviour =
421 next_iteration_inv_jacobian_recompute =
422 next_iteration_inv_jacobian_recompute + _recompute_jacobian;
423 }
424 if (!_linear_solver.compute(J, linear_solver_behaviour))
425 {
426 ERR("Newton: The linear solver failed in the compute() step.");
427 }
428 bool iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
429#else
430 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
431#endif
432 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
433
434 if (!iteration_succeeded)
435 {
436 ERR("Newton: The linear solver failed.");
437 }
438 else
439 {
440 // TODO could be solved in a better way
441 // cf.
442 // https://petsc.org/release/manualpages/Vec/VecWAXPY
443
444 // Copy pointers, replace the one for the given process id.
445 std::vector<GlobalVector*> x_new{x};
446 x_new[process_id] =
448 *x[process_id], _x_new_id);
449 double damping = _damping;
450 if (_convergence_criterion->hasNonNegativeDamping())
451 {
452 damping = _convergence_criterion->getDampingFactor(
453 minus_delta_x, *x[process_id], _damping);
454 }
455 if (_damping_reduction)
456 {
457 damping =
458 damping +
459 (1 - damping) *
460 std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
461 }
462 LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
463
464 if (postIterationCallback)
465 {
466 postIterationCallback(iteration, error_norms_met, x_new);
467 }
468
469 switch (sys.postIteration(*x_new[process_id]))
470 {
472 break;
474 ERR("Newton: The postIteration() hook reported a "
475 "non-recoverable error.");
476 iteration_succeeded = false;
477 break;
479 INFO(
480 "Newton: The postIteration() hook decided that this "
481 "iteration has to be repeated.");
482 // TODO introduce some onDestroy hook.
484 *x_new[process_id]);
485 continue; // That throws the iteration result away.
486 }
487
488 LinAlg::copy(*x_new[process_id],
489 *x[process_id]); // copy new solution to x
491 *x_new[process_id]);
492 }
493
494 if (!iteration_succeeded)
495 {
496 // Don't compute further error norms, but break here.
497 error_norms_met = false;
498 break;
499 }
500
501 if (sys.isLinear())
502 {
503 error_norms_met = true;
504 }
505 else
506 {
507 if (_convergence_criterion->hasDeltaXCheck())
508 {
509 // Note: x contains the new solution!
510 _convergence_criterion->checkDeltaX(minus_delta_x,
511 *x[process_id]);
512 }
513
514 error_norms_met = _convergence_criterion->isSatisfied();
515 }
516
517 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
518 time_iteration.elapsed());
519
520 if (error_norms_met)
521 {
522 break;
523 }
524
525 // Avoid increment of the 'iteration' if the error norms are not met,
526 // but maximum number of iterations is reached.
527 if (iteration >= _maxiter)
528 {
529 break;
530 }
531 }
532
533 if (iteration > _maxiter)
534 {
535 ERR("Newton: Could not solve the given nonlinear system within {:d} "
536 "iterations",
537 _maxiter);
538 }
539
543
544 return {error_norms_met, iteration};
545}
546
547std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
549 BaseLib::ConfigTree const& config)
550{
552 auto const type = config.getConfigParameter<std::string>("type");
554 auto const max_iter = config.getConfigParameter<int>("max_iter");
555
556 if (type == "Picard")
557 {
558 auto const tag = NonlinearSolverTag::Picard;
559 using ConcreteNLS = NonlinearSolver<tag>;
560 return std::make_pair(
561 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
562 }
563 if (type == "Newton")
564 {
565 auto const recompute_jacobian =
567 config.getConfigParameter<int>("recompute_jacobian", 1);
568
570 auto const damping = config.getConfigParameter<double>("damping", 1.0);
571 if (damping <= 0)
572 {
573 OGS_FATAL(
574 "The damping factor for the Newon method must be positive, got "
575 "{:g}.",
576 damping);
577 }
578 auto const damping_reduction =
580 config.getConfigParameterOptional<double>("damping_reduction");
581 auto const tag = NonlinearSolverTag::Newton;
582 using ConcreteNLS = NonlinearSolver<tag>;
583 return std::make_pair(std::make_unique<ConcreteNLS>(
584 linear_solver, max_iter, recompute_jacobian,
585 damping, damping_reduction),
586 tag);
587 }
588#ifdef USE_PETSC
589 if (boost::iequals(type, "PETScSNES"))
590 {
591 auto prefix =
593 config.getConfigParameter<std::string>("prefix", "");
594 auto const tag = NonlinearSolverTag::Newton;
595 using ConcreteNLS = PETScNonlinearSolver;
596 return std::make_pair(std::make_unique<ConcreteNLS>(
597 linear_solver, max_iter, std::move(prefix)),
598 tag);
599 }
600
601#endif
602 OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
603}
604
606{
607 if (_r_neq != nullptr)
608 {
610 }
611}
612
620
621} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:40
Definition of the RunTime class.
std::optional< T > getConfigParameterOptional(std::string const &param) const
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
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
Global vector based on Eigen vector.
Definition EigenVector.h:25
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
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)
static bool anyOf(bool const val, Mpi const &mpi=Mpi{MPI_COMM_WORLD})
Definition MPI.h:187
void finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:198
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void setLocalAccessibleVector(PETScVector const &x)
Definition LinAlg.cpp:27
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:149
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition LinAlg.cpp:57
bool solvePicard(GlobalLinearSolver &linear_solver, GlobalMatrix &A, GlobalVector &rhs, GlobalVector &x, MathLib::LinearSolverBehaviour const linear_solver_behaviour)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.