OGS
NumLib::NonlinearSolver< NonlinearSolverTag::Newton > Class Referencefinal

## Detailed Description

Find a solution to a nonlinear equation using the Newton-Raphson method.

Definition at line 76 of file NonlinearSolver.h.

#include <NonlinearSolver.h>

Inheritance diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:
[legend]
Collaboration diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:
[legend]

## Public Types

using System = NonlinearSystem<NonlinearSolverTag::Newton>
Type of the nonlinear equation system to be solved.

## Public Member Functions

NonlinearSolver (GlobalLinearSolver &linear_solver, int const maxiter, double const damping=1.0)

~NonlinearSolver ()

void setEquationSystem (System &eq, ConvergenceCriterion &conv_crit)

void calculateNonEquilibriumInitialResiduum (std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id) override

NonlinearSolverStatus solve (std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::function< void(int, std::vector< GlobalVector * > const &)> const &postIterationCallback, int const process_id) override

void compensateNonEquilibriumInitialResiduum (bool const value)

Public Member Functions inherited from NumLib::NonlinearSolverBase
virtual ~NonlinearSolverBase ()=default

## Private Attributes

GlobalLinearSolver_linear_solver

System_equation_system = nullptr

ConvergenceCriterion_convergence_criterion = nullptr

int const _maxiter
maximum number of iterations

double const _damping

GlobalVector_r_neq = nullptr
non-equilibrium initial residuum.

std::size_t _res_id = 0u
ID of the residual vector.

std::size_t _J_id = 0u
ID of the Jacobian matrix.

std::size_t _minus_delta_x_id = 0u
ID of the $$-\Delta x$$ vector.

std::size_t _x_new_id
ID of the vector storing $$x - (-\Delta x)$$.

std::size_t _r_neq_id = 0u

bool _compensate_non_equilibrium_initial_residuum = false

## ◆ System

 using NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::System = NonlinearSystem

Type of the nonlinear equation system to be solved.

Definition at line 81 of file NonlinearSolver.h.

## ◆ NonlinearSolver()

 NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::NonlinearSolver ( GlobalLinearSolver & linear_solver, int const maxiter, double const damping = 1.0 )
inlineexplicit

Constructs a new instance.

Parameters
 linear_solver the linear solver used by this nonlinear solver. maxiter the maximum number of iterations used to solve the equation. damping A positive damping factor.
_damping

Definition at line 91 of file NonlinearSolver.h.

94 : _linear_solver(linear_solver), _maxiter(maxiter), _damping(damping)
95 {
96 }
int const _maxiter
maximum number of iterations

## ◆ ~NonlinearSolver()

 NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::~NonlinearSolver ( )

Definition at line 554 of file NonlinearSolver.cpp.

555{
556 if (_r_neq != nullptr)
557 {
559 }
560}
GlobalVector * _r_neq
non-equilibrium initial residuum.
virtual void releaseVector(GlobalVector const &x)=0
static NUMLIB_EXPORT VectorProvider & provider

## ◆ calculateNonEquilibriumInitialResiduum()

 void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::calculateNonEquilibriumInitialResiduum ( std::vector< GlobalVector * > const & x, std::vector< GlobalVector * > const & x_prev, int const process_id )
overridevirtual

Implements NumLib::NonlinearSolverBase.

Definition at line 291 of file NonlinearSolver.cpp.

295{
297 {
298 return;
299 }
300
301 INFO("Calculate non-equilibrium initial residuum.");
302
303 _equation_system->assemble(x, x_prev, process_id);
305 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
306
307 // Set the values of the selected entries of _r_neq, which are associated
308 // with the equations that do not need initial residual compensation, to
309 // zero.
310 const auto selected_global_indices =
312 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
313
314 _r_neq->set(selected_global_indices, zero_entries);
315
317}
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:35
void set(IndexType rowId, double v)
set entry
Definition EigenVector.h:73
virtual std::vector< GlobalIndexType > getIndicesOfResiduumWithoutInitialCompensation() const =0
virtual void assemble(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)=0
virtual void getResidual(GlobalVector const &x, GlobalVector const &x_prev, GlobalVector &res) const =0
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
void finalizeAssembly(PETScMatrix &A)
Definition LinAlg.cpp:198

## ◆ compensateNonEquilibriumInitialResiduum()

 void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::compensateNonEquilibriumInitialResiduum ( bool const value )
inline

Definition at line 120 of file NonlinearSolver.h.

## ◆ setEquationSystem()

 void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::setEquationSystem ( System & eq, ConvergenceCriterion & conv_crit )
inline

Set the nonlinear equation system that will be solved. TODO doc

Definition at line 102 of file NonlinearSolver.h.

103 {
104 _equation_system = &eq;
105 _convergence_criterion = &conv_crit;
106 }

## ◆ solve()

 NonlinearSolverStatus NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::solve ( std::vector< GlobalVector * > & x, std::vector< GlobalVector * > const & x_prev, std::function< void(int, std::vector< GlobalVector * > const &)> const & postIterationCallback, int const process_id )
overridevirtual

Assemble and solve the equation system.

Parameters
 x in: the initial guess, out: the solution. x_prev previous time step solution. postIterationCallback called after each iteration if set. process_id usually used in staggered schemes.
Return values
 true if the equation system could be solved false otherwise

Implements NumLib::NonlinearSolverBase.

Definition at line 319 of file NonlinearSolver.cpp.

325{
326 namespace LinAlg = MathLib::LinAlg;
327 auto& sys = *_equation_system;
328
330 auto& minus_delta_x =
333
334 bool error_norms_met = false;
335
336 // TODO be more efficient
337 // init minus_delta_x to the right size
338 LinAlg::copy(*x[process_id], minus_delta_x);
339
341
342 int iteration = 1;
343 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
344 {
345 BaseLib::RunTime timer_dirichlet;
346 double time_dirichlet = 0.0;
347
348 BaseLib::RunTime time_iteration;
349 time_iteration.start();
350
351 timer_dirichlet.start();
352 sys.computeKnownSolutions(*x[process_id], process_id);
353 time_dirichlet += timer_dirichlet.elapsed();
354
355 sys.preIteration(iteration, *x[process_id]);
356
357 BaseLib::RunTime time_assembly;
358 time_assembly.start();
359 try
360 {
361 sys.assemble(x, x_prev, process_id);
362 }
363 catch (AssemblyException const& e)
364 {
365 ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
366 e.what());
367 error_norms_met = false;
368 iteration = _maxiter;
369 break;
370 }
371 sys.getResidual(*x[process_id], *x_prev[process_id], res);
372 sys.getJacobian(J);
373 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
374
375 // Subtract non-equilibrium initial residuum if set
376 if (_r_neq != nullptr)
377 LinAlg::axpy(res, -1, *_r_neq);
378
379 minus_delta_x.setZero();
380
381 timer_dirichlet.start();
382 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
383 time_dirichlet += timer_dirichlet.elapsed();
384 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
385
386 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
387 {
389 }
390
391 BaseLib::RunTime time_linear_solver;
392 time_linear_solver.start();
393 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
394 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
395
396 if (!iteration_succeeded)
397 {
398 ERR("Newton: The linear solver failed.");
399 }
400 else
401 {
402 // TODO could be solved in a better way
403 // cf.
404 // https://petsc.org/release/manualpages/Vec/VecWAXPY
405
406 // Copy pointers, replace the one for the given process id.
407 std::vector<GlobalVector*> x_new{x};
408 x_new[process_id] =
410 *x[process_id], _x_new_id);
411 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
412
413 if (postIterationCallback)
414 {
415 postIterationCallback(iteration, x_new);
416 }
417
418 switch (sys.postIteration(*x_new[process_id]))
419 {
421 break;
423 ERR("Newton: The postIteration() hook reported a "
424 "non-recoverable error.");
425 iteration_succeeded = false;
426 break;
428 INFO(
429 "Newton: The postIteration() hook decided that this "
430 "iteration has to be repeated.");
431 // TODO introduce some onDestroy hook.
433 *x_new[process_id]);
434 continue; // That throws the iteration result away.
435 }
436
437 LinAlg::copy(*x_new[process_id],
438 *x[process_id]); // copy new solution to x
440 *x_new[process_id]);
441 }
442
443 if (!iteration_succeeded)
444 {
445 // Don't compute further error norms, but break here.
446 error_norms_met = false;
447 break;
448 }
449
450 if (sys.isLinear())
451 {
452 error_norms_met = true;
453 }
454 else
455 {
457 {
458 // Note: x contains the new solution!
459 _convergence_criterion->checkDeltaX(minus_delta_x,
460 *x[process_id]);
461 }
462
463 error_norms_met = _convergence_criterion->isSatisfied();
464 }
465
466 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
467 time_iteration.elapsed());
468
469 if (error_norms_met)
470 {
471 break;
472 }
473
474 // Avoid increment of the 'iteration' if the error norms are not met,
475 // but maximum number of iterations is reached.
476 if (iteration >= _maxiter)
477 {
478 break;
479 }
480 }
481
482 if (iteration > _maxiter)
483 {
484 ERR("Newton: Could not solve the given nonlinear system within {:d} "
485 "iterations",
486 _maxiter);
487 }
488
492
493 return {error_norms_met, iteration};
494}
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
Definition Logging.h:45
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)
virtual void checkResidual(GlobalVector const &residual)=0
Check if the residual satisfies the convergence criterion.
virtual bool hasResidualCheck() const =0
virtual bool isSatisfied() const
Tell if the convergence criterion is satisfied.
virtual void checkDeltaX(GlobalVector const &minus_delta_x, GlobalVector const &x)=0
virtual bool hasDeltaXCheck() const =0
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
std::size_t _J_id
ID of the Jacobian matrix.
std::size_t _x_new_id
ID of the vector storing .
std::size_t _res_id
ID of the residual vector.
void copy(PETScVector const &x, PETScVector &y)
Definition LinAlg.cpp:37
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
Definition LinAlg.cpp:57
static NUMLIB_EXPORT MatrixProvider & provider

## ◆ _compensate_non_equilibrium_initial_residuum

 bool NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_compensate_non_equilibrium_initial_residuum = false
private

Enables computation of the non-equilibrium initial residuum $$r_{\rm neq}$$ before the first time step. The forces are zero if the external forces are in equilibrium with the initial state/initial conditions. During the simulation the new residuum reads $$\tilde r = r - r_{\rm neq}$$.

Definition at line 153 of file NonlinearSolver.h.

## ◆ _convergence_criterion

 ConvergenceCriterion* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_convergence_criterion = nullptr
private

Definition at line 130 of file NonlinearSolver.h.

## ◆ _damping

 double const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_damping
private

A positive damping factor. The default value 1.0 gives a non-damped Newton method. Common values are in the range 0.5 to 0.7 for somewhat conservative method and seldom become smaller than 0.2 for very conservative approach.

Definition at line 137 of file NonlinearSolver.h.

## ◆ _equation_system

 System* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_equation_system = nullptr
private

Definition at line 127 of file NonlinearSolver.h.

## ◆ _J_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_J_id = 0u
private

ID of the Jacobian matrix.

Definition at line 141 of file NonlinearSolver.h.

## ◆ _linear_solver

 GlobalLinearSolver& NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_linear_solver
private

Definition at line 126 of file NonlinearSolver.h.

## ◆ _maxiter

 int const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_maxiter
private

maximum number of iterations

Definition at line 131 of file NonlinearSolver.h.

## ◆ _minus_delta_x_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_minus_delta_x_id = 0u
private

ID of the $$-\Delta x$$ vector.

Definition at line 142 of file NonlinearSolver.h.

## ◆ _r_neq

 GlobalVector* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_r_neq = nullptr
private

non-equilibrium initial residuum.

Definition at line 139 of file NonlinearSolver.h.

## ◆ _r_neq_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_r_neq_id = 0u
private

ID of the non-equilibrium initial residuum vector.

Definition at line 145 of file NonlinearSolver.h.

## ◆ _res_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_res_id = 0u
private

ID of the residual vector.

Definition at line 140 of file NonlinearSolver.h.

## ◆ _x_new_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_x_new_id
private
Initial value:
=
0u

ID of the vector storing $$x - (-\Delta x)$$.

Definition at line 143 of file NonlinearSolver.h.

The documentation for this class was generated from the following files: