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