Loading [MathJax]/jax/output/HTML-CSS/config.js
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
425 bool iteration_succeeded = false;
426 if (!_linear_solver.compute(J, linear_solver_behaviour))
427 {
428 ERR("Newton: The linear solver failed in the compute() step.");
429 }
430 else
431 {
432 iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
433 }
434#else
435 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
436#endif
437 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
438
439 if (!iteration_succeeded)
440 {
441 ERR("Newton: The linear solver failed.");
442 }
443 else
444 {
445 // TODO could be solved in a better way
446 // cf.
447 // https://petsc.org/release/manualpages/Vec/VecWAXPY
448
449 // Copy pointers, replace the one for the given process id.
450 std::vector<GlobalVector*> x_new{x};
451 x_new[process_id] =
453 *x[process_id], _x_new_id);
454 double damping = _damping;
455 if (_convergence_criterion->hasNonNegativeDamping())
456 {
457 damping = _convergence_criterion->getDampingFactor(
458 minus_delta_x, *x[process_id], _damping);
459 }
460 if (_damping_reduction)
461 {
462 damping =
463 damping +
464 (1 - damping) *
465 std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
466 }
467 LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
468
469 if (postIterationCallback)
470 {
471 postIterationCallback(iteration, error_norms_met, x_new);
472 }
473
474 switch (sys.postIteration(*x_new[process_id]))
475 {
477 break;
479 ERR("Newton: The postIteration() hook reported a "
480 "non-recoverable error.");
481 iteration_succeeded = false;
482 break;
484 INFO(
485 "Newton: The postIteration() hook decided that this "
486 "iteration has to be repeated.");
487 // TODO introduce some onDestroy hook.
489 *x_new[process_id]);
490 continue; // That throws the iteration result away.
491 }
492
493 LinAlg::copy(*x_new[process_id],
494 *x[process_id]); // copy new solution to x
496 *x_new[process_id]);
497 }
498
499 if (!iteration_succeeded)
500 {
501 // Don't compute further error norms, but break here.
502 error_norms_met = false;
503 break;
504 }
505
506 if (sys.isLinear())
507 {
508 error_norms_met = true;
509 }
510 else
511 {
512 if (_convergence_criterion->hasDeltaXCheck())
513 {
514 // Note: x contains the new solution!
515 _convergence_criterion->checkDeltaX(minus_delta_x,
516 *x[process_id]);
517 }
518
519 error_norms_met = _convergence_criterion->isSatisfied();
520 }
521
522 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
523 time_iteration.elapsed());
524
525 if (error_norms_met)
526 {
527 break;
528 }
529
530 // Avoid increment of the 'iteration' if the error norms are not met,
531 // but maximum number of iterations is reached.
532 if (iteration >= _maxiter)
533 {
534 break;
535 }
536 }
537
538 if (iteration > _maxiter)
539 {
540 ERR("Newton: Could not solve the given nonlinear system within {:d} "
541 "iterations",
542 _maxiter);
543 }
544
548
549 return {error_norms_met, iteration};
550}
551
552std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
554 BaseLib::ConfigTree const& config)
555{
557 auto const type = config.getConfigParameter<std::string>("type");
559 auto const max_iter = config.getConfigParameter<int>("max_iter");
560
561 if (type == "Picard")
562 {
563 auto const tag = NonlinearSolverTag::Picard;
564 using ConcreteNLS = NonlinearSolver<tag>;
565 return std::make_pair(
566 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
567 }
568 if (type == "Newton")
569 {
570 auto const recompute_jacobian =
572 config.getConfigParameter<int>("recompute_jacobian", 1);
573
575 auto const damping = config.getConfigParameter<double>("damping", 1.0);
576 if (damping <= 0)
577 {
578 OGS_FATAL(
579 "The damping factor for the Newon method must be positive, got "
580 "{:g}.",
581 damping);
582 }
583 auto const damping_reduction =
585 config.getConfigParameterOptional<double>("damping_reduction");
586 auto const tag = NonlinearSolverTag::Newton;
587 using ConcreteNLS = NonlinearSolver<tag>;
588 return std::make_pair(std::make_unique<ConcreteNLS>(
589 linear_solver, max_iter, recompute_jacobian,
590 damping, damping_reduction),
591 tag);
592 }
593#ifdef USE_PETSC
594 if (boost::iequals(type, "PETScSNES"))
595 {
596 auto prefix =
598 config.getConfigParameter<std::string>("prefix", "");
599 auto const tag = NonlinearSolverTag::Newton;
600 using ConcreteNLS = PETScNonlinearSolver;
601 return std::make_pair(std::make_unique<ConcreteNLS>(
602 linear_solver, max_iter, std::move(prefix)),
603 tag);
604 }
605
606#endif
607 OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
608}
609
611{
612 if (_r_neq != nullptr)
613 {
615 }
616}
617
625
626} // 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: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.