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 _equation_system->setReleaseNodalForces(_r_neq, process_id);
121
123
126}
127
129 std::vector<GlobalVector*>& x,
130 std::vector<GlobalVector*> const& x_prev,
131 std::function<void(int, bool, std::vector<GlobalVector*> const&)> const&
132 postIterationCallback,
133 int const process_id)
134{
135 namespace LinAlg = MathLib::LinAlg;
136 auto& sys = *_equation_system;
137
140
141 std::vector<GlobalVector*> x_new{x};
142 x_new[process_id] =
144 LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
145
146 bool error_norms_met = false;
147
148 _convergence_criterion->preFirstIteration();
149
150 int iteration = 1;
151 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
152 {
153 BaseLib::RunTime timer_dirichlet;
154 double time_dirichlet = 0.0;
155
156 BaseLib::RunTime time_iteration;
157 time_iteration.start();
158
159 INFO("Iteration #{:d} started.", iteration);
160 timer_dirichlet.start();
161 auto& x_new_process = *x_new[process_id];
163 sys.computeKnownSolutions(x_new_process, process_id);
164 sys.applyKnownSolutions(x_new_process);
165 time_dirichlet += timer_dirichlet.elapsed();
166
167 sys.preIteration(iteration, x_new_process);
168
169 BaseLib::RunTime time_assembly;
170 time_assembly.start();
171 sys.assemble(x_new, x_prev, process_id);
172 sys.getA(A);
173 sys.getRhs(*x_prev[process_id], rhs);
174
175 // Normalize the linear equation system, if required
176 if (sys.requiresNormalization() &&
177 !_linear_solver.canSolveRectangular())
178 {
179 sys.getAandRhsNormalized(A, rhs);
180 WARN(
181 "The equation system is rectangular, but the current linear "
182 "solver only supports square systems. "
183 "The system will be normalized, which lead to a squared "
184 "condition number and potential numerical issues. "
185 "It is recommended to use a solver that supports rectangular "
186 "equation systems for better numerical stability.");
187 }
188
189 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
190
191 // Subtract non-equilibrium initial residuum if set
192 if (_r_neq != nullptr)
193 {
194 LinAlg::axpy(rhs, -1, *_r_neq);
195 }
196
197 timer_dirichlet.start();
198 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
199 time_dirichlet += timer_dirichlet.elapsed();
200 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
201
202 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
203 {
204 GlobalVector res;
205 LinAlg::matMult(A, x_new_process, res); // res = A * x_new
206 LinAlg::axpy(res, -1.0, rhs); // res -= rhs
207 _convergence_criterion->checkResidual(res);
208 }
209
210 bool iteration_succeeded =
211 detail::solvePicard(_linear_solver, A, rhs, x_new_process,
212 sys.linearSolverNeedsToCompute());
213
214 if (iteration_succeeded)
215 {
216 if (postIterationCallback)
217 {
218 postIterationCallback(iteration, error_norms_met, x_new);
219 }
220
221 switch (sys.postIteration(x_new_process))
222 {
224 // Don't copy here. The old x might still be used further
225 // below. Although currently it is not.
226 break;
228 ERR("Picard: The postIteration() hook reported a "
229 "non-recoverable error.");
230 iteration_succeeded = false;
231 // Copy new solution to x.
232 // Thereby the failed solution can be used by the caller for
233 // debugging purposes.
234 LinAlg::copy(x_new_process, *x[process_id]);
235 break;
237 INFO(
238 "Picard: The postIteration() hook decided that this "
239 "iteration has to be repeated.");
241 *x[process_id],
242 x_new_process); // throw the iteration result away
243 continue;
244 }
245 }
246
247 if (!iteration_succeeded)
248 {
249 // Don't compute error norms, break here.
250 error_norms_met = false;
251 break;
252 }
253
254 if (sys.isLinear())
255 {
256 error_norms_met = true;
257 }
258 else
259 {
260 if (_convergence_criterion->hasDeltaXCheck())
261 {
262 GlobalVector minus_delta_x(*x[process_id]);
263 LinAlg::axpy(minus_delta_x, -1.0,
264 x_new_process); // minus_delta_x = x - x_new
265 _convergence_criterion->checkDeltaX(minus_delta_x,
266 x_new_process);
267 }
268
269 error_norms_met = _convergence_criterion->isSatisfied();
270 }
271
272 // Update x s.t. in the next iteration we will compute the right delta x
273 LinAlg::copy(x_new_process, *x[process_id]);
274
275 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
276 time_iteration.elapsed());
277
278 if (error_norms_met)
279 {
280 break;
281 }
282
283 // Avoid increment of the 'iteration' if the error norms are not met,
284 // but maximum number of iterations is reached.
285 if (iteration >= _maxiter)
286 {
287 break;
288 }
289 }
290
291 if (iteration > _maxiter)
292 {
293 ERR("Picard: Could not solve the given nonlinear system within {:d} "
294 "iterations",
295 _maxiter);
296 }
297
301
302 return {error_norms_met, iteration};
303}
304
307 std::vector<GlobalVector*> const& x,
308 std::vector<GlobalVector*> const& x_prev, int const process_id)
309{
310 if (!_compensate_non_equilibrium_initial_residuum)
311 {
312 return;
313 }
314
315 INFO("Calculate non-equilibrium initial residuum.");
316
317 _equation_system->assemble(x, x_prev, process_id);
319 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
320
321 // Set the values of the selected entries of _r_neq, which are associated
322 // with the equations that do not need initial residual compensation, to
323 // zero.
324 const auto selected_global_indices =
325 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
326 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
327
328 _r_neq->set(selected_global_indices, zero_entries);
329 _equation_system->setReleaseNodalForces(_r_neq, process_id);
330
332}
333
335 std::vector<GlobalVector*>& x,
336 std::vector<GlobalVector*> const& x_prev,
337 std::function<void(int, bool, std::vector<GlobalVector*> const&)> const&
338 postIterationCallback,
339 int const process_id)
340{
341 namespace LinAlg = MathLib::LinAlg;
342 auto& sys = *_equation_system;
343
345 auto& minus_delta_x =
348
349 bool error_norms_met = false;
350
351 // TODO be more efficient
352 // init minus_delta_x to the right size
353 LinAlg::copy(*x[process_id], minus_delta_x);
354
355 _convergence_criterion->preFirstIteration();
356
357 int iteration = 1;
358#if !defined(USE_PETSC) && !defined(USE_LIS)
359 int next_iteration_inv_jacobian_recompute = 1;
360#endif
361 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
362 {
363 BaseLib::RunTime timer_dirichlet;
364 double time_dirichlet = 0.0;
365
366 BaseLib::RunTime time_iteration;
367 INFO("Iteration #{:d} started.", iteration);
368 time_iteration.start();
369
370 timer_dirichlet.start();
371 sys.computeKnownSolutions(*x[process_id], process_id);
372 time_dirichlet += timer_dirichlet.elapsed();
373
374 sys.preIteration(iteration, *x[process_id]);
375
376 BaseLib::RunTime time_assembly;
377 time_assembly.start();
378 bool mpi_rank_assembly_ok = true;
379 try
380 {
381 sys.assemble(x, x_prev, process_id);
382 }
383 catch (AssemblyException const& e)
384 {
385 ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
386 e.what());
387 error_norms_met = false;
388 iteration = _maxiter;
389 mpi_rank_assembly_ok = false;
390 }
391 if (BaseLib::MPI::anyOf(!mpi_rank_assembly_ok))
392 {
393 break;
394 }
395 sys.getResidual(*x[process_id], *x_prev[process_id], res);
396 sys.getJacobian(J);
397 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
398
399 // Subtract non-equilibrium initial residuum if set
400 if (_r_neq != nullptr)
401 LinAlg::axpy(res, -1, *_r_neq);
402
403 minus_delta_x.setZero();
404
405 timer_dirichlet.start();
406 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
407 time_dirichlet += timer_dirichlet.elapsed();
408 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
409
410 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
411 {
412 _convergence_criterion->checkResidual(res);
413 }
414
415 BaseLib::RunTime time_linear_solver;
416 time_linear_solver.start();
417#if !defined(USE_PETSC) && !defined(USE_LIS)
418 auto linear_solver_behaviour = MathLib::LinearSolverBehaviour::REUSE;
419 if (iteration == next_iteration_inv_jacobian_recompute)
420 {
421 linear_solver_behaviour =
423 next_iteration_inv_jacobian_recompute =
424 next_iteration_inv_jacobian_recompute + _recompute_jacobian;
425 }
426
427 bool iteration_succeeded = false;
428 if (!_linear_solver.compute(J, linear_solver_behaviour))
429 {
430 ERR("Newton: The linear solver failed in the compute() step.");
431 }
432 else
433 {
434 iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
435 }
436#else
437 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
438#endif
439 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
440
441 if (!iteration_succeeded)
442 {
443 ERR("Newton: The linear solver failed.");
444 }
445 else
446 {
447 // TODO could be solved in a better way
448 // cf.
449 // https://petsc.org/release/manualpages/Vec/VecWAXPY
450
451 // Copy pointers, replace the one for the given process id.
452 std::vector<GlobalVector*> x_new{x};
453 x_new[process_id] =
455 *x[process_id], _x_new_id);
456 double damping = _damping;
457 if (_convergence_criterion->hasNonNegativeDamping())
458 {
459 damping = _convergence_criterion->getDampingFactor(
460 minus_delta_x, *x[process_id], _damping);
461 }
462 if (_damping_reduction)
463 {
464 damping =
465 damping +
466 (1 - damping) *
467 std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
468 }
469 LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
470
471 if (postIterationCallback)
472 {
473 postIterationCallback(iteration, error_norms_met, x_new);
474 }
475
476 switch (sys.postIteration(*x_new[process_id]))
477 {
479 break;
481 ERR("Newton: The postIteration() hook reported a "
482 "non-recoverable error.");
483 iteration_succeeded = false;
484 break;
486 INFO(
487 "Newton: The postIteration() hook decided that this "
488 "iteration has to be repeated.");
489 // TODO introduce some onDestroy hook.
491 *x_new[process_id]);
492 continue; // That throws the iteration result away.
493 }
494
495 LinAlg::copy(*x_new[process_id],
496 *x[process_id]); // copy new solution to x
498 *x_new[process_id]);
499 }
500
501 if (!iteration_succeeded)
502 {
503 // Don't compute further error norms, but break here.
504 error_norms_met = false;
505 break;
506 }
507
508 if (sys.isLinear())
509 {
510 error_norms_met = true;
511 }
512 else
513 {
514 if (_convergence_criterion->hasDeltaXCheck())
515 {
516 // Note: x contains the new solution!
517 _convergence_criterion->checkDeltaX(minus_delta_x,
518 *x[process_id]);
519 }
520
521 error_norms_met = _convergence_criterion->isSatisfied();
522 }
523
524 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
525 time_iteration.elapsed());
526
527 if (error_norms_met)
528 {
529 break;
530 }
531
532 // Avoid increment of the 'iteration' if the error norms are not met,
533 // but maximum number of iterations is reached.
534 if (iteration >= _maxiter)
535 {
536 break;
537 }
538 }
539
540 if (iteration > _maxiter)
541 {
542 ERR("Newton: Could not solve the given nonlinear system within {:d} "
543 "iterations",
544 _maxiter);
545 }
546
550
551 return {error_norms_met, iteration};
552}
553
554std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
556 BaseLib::ConfigTree const& config)
557{
559 auto const type = config.getConfigParameter<std::string>("type");
561 auto const max_iter = config.getConfigParameter<int>("max_iter");
562
563 if (type == "Picard")
564 {
565 auto const tag = NonlinearSolverTag::Picard;
566 using ConcreteNLS = NonlinearSolver<tag>;
567 return std::make_pair(
568 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
569 }
570 if (type == "Newton")
571 {
572 auto const recompute_jacobian =
574 config.getConfigParameter<int>("recompute_jacobian", 1);
575
577 auto const damping = config.getConfigParameter<double>("damping", 1.0);
578 if (damping <= 0)
579 {
580 OGS_FATAL(
581 "The damping factor for the Newon method must be positive, got "
582 "{:g}.",
583 damping);
584 }
585 auto const damping_reduction =
587 config.getConfigParameterOptional<double>("damping_reduction");
588 auto const tag = NonlinearSolverTag::Newton;
589 using ConcreteNLS = NonlinearSolver<tag>;
590 return std::make_pair(std::make_unique<ConcreteNLS>(
591 linear_solver, max_iter, recompute_jacobian,
592 damping, damping_reduction),
593 tag);
594 }
595#ifdef USE_PETSC
596 if (boost::iequals(type, "PETScSNES"))
597 {
598 auto prefix =
600 config.getConfigParameter<std::string>("prefix", "");
601 auto const tag = NonlinearSolverTag::Newton;
602 using ConcreteNLS = PETScNonlinearSolver;
603 return std::make_pair(std::make_unique<ConcreteNLS>(
604 linear_solver, max_iter, std::move(prefix)),
605 tag);
606 }
607
608#endif
609 OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
610}
611
613{
614 if (_r_neq != nullptr)
615 {
617 }
618}
619
627
628} // namespace NumLib
#define OGS_FATAL(...)
Definition Error.h:26
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:36
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:48
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:42
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:26
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{OGS_COMM_WORLD})
Definition MPI.h:191
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.