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"
22#include "NumLib/Exceptions.h"
24
25namespace NumLib
26{
27namespace detail
28{
29#if !defined(USE_PETSC) && !defined(USE_LIS)
30bool solvePicard(GlobalLinearSolver& linear_solver, GlobalMatrix& A,
32 bool const compute_necessary)
33{
34 BaseLib::RunTime time_linear_solver;
35 time_linear_solver.start();
36
37 if (compute_necessary)
38 {
39 if (!linear_solver.compute(A))
40 {
41 ERR("Picard: The linear solver failed in the compute() step.");
42 return false;
43 }
44 }
45
46 bool const iteration_succeeded = linear_solver.solve(rhs, x);
47
48 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
49
50 if (iteration_succeeded)
51 {
52 return true;
53 }
54
55 ERR("Picard: The linear solver failed in the solve() step.");
56 return false;
57}
58#else
61 bool const compute_necessary)
62{
63 if (!compute_necessary)
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 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
172
173 // Subtract non-equilibrium initial residuum if set
174 if (_r_neq != nullptr)
175 {
176 LinAlg::axpy(rhs, -1, *_r_neq);
177 }
178
179 timer_dirichlet.start();
180 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
181 time_dirichlet += timer_dirichlet.elapsed();
182 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
183
184 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
185 {
186 GlobalVector res;
187 LinAlg::matMult(A, x_new_process, res); // res = A * x_new
188 LinAlg::axpy(res, -1.0, rhs); // res -= rhs
189 _convergence_criterion->checkResidual(res);
190 }
191
192 bool iteration_succeeded =
193 detail::solvePicard(_linear_solver, A, rhs, x_new_process,
194 sys.linearSolverNeedsToCompute());
195
196 if (iteration_succeeded)
197 {
198 if (postIterationCallback)
199 {
200 postIterationCallback(iteration, x_new);
201 }
202
203 switch (sys.postIteration(x_new_process))
204 {
206 // Don't copy here. The old x might still be used further
207 // below. Although currently it is not.
208 break;
210 ERR("Picard: The postIteration() hook reported a "
211 "non-recoverable error.");
212 iteration_succeeded = false;
213 // Copy new solution to x.
214 // Thereby the failed solution can be used by the caller for
215 // debugging purposes.
216 LinAlg::copy(x_new_process, *x[process_id]);
217 break;
219 INFO(
220 "Picard: The postIteration() hook decided that this "
221 "iteration has to be repeated.");
223 *x[process_id],
224 x_new_process); // throw the iteration result away
225 continue;
226 }
227 }
228
229 if (!iteration_succeeded)
230 {
231 // Don't compute error norms, break here.
232 error_norms_met = false;
233 break;
234 }
235
236 if (sys.isLinear())
237 {
238 error_norms_met = true;
239 }
240 else
241 {
242 if (_convergence_criterion->hasDeltaXCheck())
243 {
244 GlobalVector minus_delta_x(*x[process_id]);
245 LinAlg::axpy(minus_delta_x, -1.0,
246 x_new_process); // minus_delta_x = x - x_new
247 _convergence_criterion->checkDeltaX(minus_delta_x,
248 x_new_process);
249 }
250
251 error_norms_met = _convergence_criterion->isSatisfied();
252 }
253
254 // Update x s.t. in the next iteration we will compute the right delta x
255 LinAlg::copy(x_new_process, *x[process_id]);
256
257 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
258 time_iteration.elapsed());
259
260 if (error_norms_met)
261 {
262 break;
263 }
264
265 // Avoid increment of the 'iteration' if the error norms are not met,
266 // but maximum number of iterations is reached.
267 if (iteration >= _maxiter)
268 {
269 break;
270 }
271 }
272
273 if (iteration > _maxiter)
274 {
275 ERR("Picard: Could not solve the given nonlinear system within {:d} "
276 "iterations",
277 _maxiter);
278 }
279
283
284 return {error_norms_met, iteration};
285}
286
289 std::vector<GlobalVector*> const& x,
290 std::vector<GlobalVector*> const& x_prev, int const process_id)
291{
292 if (!_compensate_non_equilibrium_initial_residuum)
293 {
294 return;
295 }
296
297 INFO("Calculate non-equilibrium initial residuum.");
298
299 _equation_system->assemble(x, x_prev, process_id);
301 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
302
303 // Set the values of the selected entries of _r_neq, which are associated
304 // with the equations that do not need initial residual compensation, to
305 // zero.
306 const auto selected_global_indices =
307 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
308 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
309
310 _r_neq->set(selected_global_indices, zero_entries);
311
313}
314
316 std::vector<GlobalVector*>& x,
317 std::vector<GlobalVector*> const& x_prev,
318 std::function<void(int, std::vector<GlobalVector*> const&)> const&
319 postIterationCallback,
320 int const process_id)
321{
322 namespace LinAlg = MathLib::LinAlg;
323 auto& sys = *_equation_system;
324
326 auto& minus_delta_x =
329
330 bool error_norms_met = false;
331
332 // TODO be more efficient
333 // init minus_delta_x to the right size
334 LinAlg::copy(*x[process_id], minus_delta_x);
335
336 _convergence_criterion->preFirstIteration();
337
338 int iteration = 1;
339 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
340 {
341 BaseLib::RunTime timer_dirichlet;
342 double time_dirichlet = 0.0;
343
344 BaseLib::RunTime time_iteration;
345 time_iteration.start();
346
347 timer_dirichlet.start();
348 sys.computeKnownSolutions(*x[process_id], process_id);
349 sys.applyKnownSolutions(*x[process_id]);
350 time_dirichlet += timer_dirichlet.elapsed();
351
352 sys.preIteration(iteration, *x[process_id]);
353
354 BaseLib::RunTime time_assembly;
355 time_assembly.start();
356 try
357 {
358 sys.assemble(x, x_prev, process_id);
359 }
360 catch (AssemblyException const& e)
361 {
362 ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
363 e.what());
364 error_norms_met = false;
365 iteration = _maxiter;
366 break;
367 }
368 sys.getResidual(*x[process_id], *x_prev[process_id], res);
369 sys.getJacobian(J);
370 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
371
372 // Subtract non-equilibrium initial residuum if set
373 if (_r_neq != nullptr)
374 LinAlg::axpy(res, -1, *_r_neq);
375
376 minus_delta_x.setZero();
377
378 timer_dirichlet.start();
379 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
380 time_dirichlet += timer_dirichlet.elapsed();
381 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
382
383 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
384 {
385 _convergence_criterion->checkResidual(res);
386 }
387
388 BaseLib::RunTime time_linear_solver;
389 time_linear_solver.start();
390 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
391 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
392
393 if (!iteration_succeeded)
394 {
395 ERR("Newton: The linear solver failed.");
396 }
397 else
398 {
399 // TODO could be solved in a better way
400 // cf.
401 // https://www.petsc.org/release/docs/manualpages/Vec/VecWAXPY.html
402
403 // Copy pointers, replace the one for the given process id.
404 std::vector<GlobalVector*> x_new{x};
405 x_new[process_id] =
407 *x[process_id], _x_new_id);
408 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
409
410 if (postIterationCallback)
411 {
412 postIterationCallback(iteration, x_new);
413 }
414
415 switch (sys.postIteration(*x_new[process_id]))
416 {
418 break;
420 ERR("Newton: The postIteration() hook reported a "
421 "non-recoverable error.");
422 iteration_succeeded = false;
423 break;
425 INFO(
426 "Newton: The postIteration() hook decided that this "
427 "iteration has to be repeated.");
428 // TODO introduce some onDestroy hook.
430 *x_new[process_id]);
431 continue; // That throws the iteration result away.
432 }
433
434 LinAlg::copy(*x_new[process_id],
435 *x[process_id]); // copy new solution to x
437 *x_new[process_id]);
438 }
439
440 if (!iteration_succeeded)
441 {
442 // Don't compute further error norms, but break here.
443 error_norms_met = false;
444 break;
445 }
446
447 if (sys.isLinear())
448 {
449 error_norms_met = true;
450 }
451 else
452 {
453 if (_convergence_criterion->hasDeltaXCheck())
454 {
455 // Note: x contains the new solution!
456 _convergence_criterion->checkDeltaX(minus_delta_x,
457 *x[process_id]);
458 }
459
460 error_norms_met = _convergence_criterion->isSatisfied();
461 }
462
463 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
464 time_iteration.elapsed());
465
466 if (error_norms_met)
467 {
468 break;
469 }
470
471 // Avoid increment of the 'iteration' if the error norms are not met,
472 // but maximum number of iterations is reached.
473 if (iteration >= _maxiter)
474 {
475 break;
476 }
477 }
478
479 if (iteration > _maxiter)
480 {
481 ERR("Newton: Could not solve the given nonlinear system within {:d} "
482 "iterations",
483 _maxiter);
484 }
485
489
490 return {error_norms_met, iteration};
491}
492
493std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
495 BaseLib::ConfigTree const& config)
496{
498 auto const type = config.getConfigParameter<std::string>("type");
500 auto const max_iter = config.getConfigParameter<int>("max_iter");
501
502 if (type == "Picard")
503 {
504 auto const tag = NonlinearSolverTag::Picard;
505 using ConcreteNLS = NonlinearSolver<tag>;
506 return std::make_pair(
507 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
508 }
509 if (type == "Newton")
510 {
512 auto const damping = config.getConfigParameter<double>("damping", 1.0);
513 if (damping <= 0)
514 {
515 OGS_FATAL(
516 "The damping factor for the Newon method must be positive, got "
517 "{:g}.",
518 damping);
519 }
520 auto const tag = NonlinearSolverTag::Newton;
521 using ConcreteNLS = NonlinearSolver<tag>;
522 return std::make_pair(
523 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
524 tag);
525 }
526#ifdef USE_PETSC
527 if (boost::iequals(type, "PETScSNES"))
528 {
529 auto prefix =
531 config.getConfigParameter<std::string>("prefix", "");
532 auto const tag = NonlinearSolverTag::Newton;
533 using ConcreteNLS = PETScNonlinearSolver;
534 return std::make_pair(std::make_unique<ConcreteNLS>(
535 linear_solver, max_iter, std::move(prefix)),
536 tag);
537 }
538
539#endif
540 OGS_FATAL("Unsupported nonlinear solver type '{:s}'.", type.c_str());
541}
542
544{
545 if (_r_neq != nullptr)
546 {
548 }
549}
550
552{
553 if (_r_neq != nullptr)
554 {
556 }
557}
558
559} // 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)
void finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:163
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:142
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, bool const compute_necessary)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.