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

Detailed Description

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

Definition at line 75 of file NonlinearSolver.h.

#include <NonlinearSolver.h>

Inheritance diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:
Collaboration diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:

Public Types

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

Public Member Functions

 NonlinearSolver (GlobalLinearSolver &linear_solver, int const maxiter, double const damping=1.0)
 
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
 
const int _maxiter
 maximum number of iterations More...
 
const double _damping
 
GlobalVector * _r_neq = nullptr
 non-equilibrium initial residuum. More...
 
std::size_t _res_id = 0u
 ID of the residual vector. More...
 
std::size_t _J_id = 0u
 ID of the Jacobian matrix. More...
 
std::size_t _minus_delta_x_id = 0u
 ID of the \( -\Delta x\) vector. More...
 
std::size_t _x_new_id
 ID of the vector storing \( x - (-\Delta x) \). More...
 
bool _compensate_non_equilibrium_initial_residuum = false
 

Member Typedef Documentation

◆ System

Type of the nonlinear equation system to be solved.

Definition at line 80 of file NonlinearSolver.h.

Constructor & Destructor Documentation

◆ NonlinearSolver()

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

Constructs a new instance.

Parameters
linear_solverthe linear solver used by this nonlinear solver.
maxiterthe maximum number of iterations used to solve the equation.
dampingA positive damping factor.
See also
_damping

Definition at line 90 of file NonlinearSolver.h.

93  : _linear_solver(linear_solver), _maxiter(maxiter), _damping(damping)
94  {
95  }

Member Function Documentation

◆ 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 213 of file NonlinearSolver.cpp.

216 {
218  {
219  return;
220  }
221 
222  _equation_system->assemble(x, x_prev, process_id);
224  _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
225 }

References NumLib::VectorProvider::getVector(), and NumLib::GlobalVectorProvider::provider.

◆ compensateNonEquilibriumInitialResiduum()

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

Definition at line 117 of file NonlinearSolver.h.

118  {
120  }

◆ 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 99 of file NonlinearSolver.h.

100  {
101  _equation_system = &eq;
102  _convergence_criterion = &conv_crit;
103  }

◆ 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
xin: the initial guess, out: the solution.
x_prevprevious time step solution.
postIterationCallbackcalled after each iteration if set.
process_idusually used in staggered schemes.
Return values
trueif the equation system could be solved
falseotherwise

Implements NumLib::NonlinearSolverBase.

Definition at line 227 of file NonlinearSolver.cpp.

233 {
234  namespace LinAlg = MathLib::LinAlg;
235  auto& sys = *_equation_system;
236 
238  _res_id);
239  auto& minus_delta_x =
242  auto& J =
244 
245  bool error_norms_met = false;
246 
247  // TODO be more efficient
248  // init minus_delta_x to the right size
249  LinAlg::copy(*x[process_id], minus_delta_x);
250 
252 
253  int iteration = 1;
254  for (; iteration <= _maxiter;
255  ++iteration, _convergence_criterion->reset())
256  {
257  BaseLib::RunTime timer_dirichlet;
258  double time_dirichlet = 0.0;
259 
260  BaseLib::RunTime time_iteration;
261  time_iteration.start();
262 
263  timer_dirichlet.start();
264  sys.computeKnownSolutions(*x[process_id], process_id);
265  sys.applyKnownSolutions(*x[process_id]);
266  time_dirichlet += timer_dirichlet.elapsed();
267 
268  sys.preIteration(iteration, *x[process_id]);
269 
270  BaseLib::RunTime time_assembly;
271  time_assembly.start();
272  try
273  {
274  sys.assemble(x, x_prev, process_id);
275  }
276  catch (AssemblyException const& e)
277  {
278  ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
279  e.what());
280  error_norms_met = false;
281  iteration = _maxiter;
282  break;
283  }
284  sys.getResidual(*x[process_id], *x_prev[process_id], res);
285  sys.getJacobian(J);
286  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
287 
288  // Subract non-equilibrium initial residuum if set
289  if (_r_neq != nullptr)
290  LinAlg::axpy(res, -1, *_r_neq);
291 
292  minus_delta_x.setZero();
293 
294  timer_dirichlet.start();
295  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
296  time_dirichlet += timer_dirichlet.elapsed();
297  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
298 
299  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
300  {
302  }
303 
304  BaseLib::RunTime time_linear_solver;
305  time_linear_solver.start();
306  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
307  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
308 
309  if (!iteration_succeeded)
310  {
311  ERR("Newton: The linear solver failed.");
312  }
313  else
314  {
315  // TODO could be solved in a better way
316  // cf.
317  // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
318 
319  // Copy pointers, replace the one for the given process id.
320  std::vector<GlobalVector*> x_new{x};
321  x_new[process_id] =
323  *x[process_id], _x_new_id);
324  LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
325 
326  if (postIterationCallback)
327  {
328  postIterationCallback(iteration, x_new);
329  }
330 
331  switch (sys.postIteration(*x_new[process_id]))
332  {
334  break;
336  ERR("Newton: The postIteration() hook reported a "
337  "non-recoverable error.");
338  iteration_succeeded = false;
339  break;
341  INFO(
342  "Newton: The postIteration() hook decided that this "
343  "iteration"
344  " has to be repeated.");
345  // TODO introduce some onDestroy hook.
347  *x_new[process_id]);
348  continue; // That throws the iteration result away.
349  }
350 
351  LinAlg::copy(*x_new[process_id],
352  *x[process_id]); // copy new solution to x
354  *x_new[process_id]);
355  }
356 
357  if (!iteration_succeeded)
358  {
359  // Don't compute further error norms, but break here.
360  error_norms_met = false;
361  break;
362  }
363 
364  if (sys.isLinear()) {
365  error_norms_met = true;
366  } else {
368  // Note: x contains the new solution!
369  _convergence_criterion->checkDeltaX(minus_delta_x,
370  *x[process_id]);
371  }
372 
373  error_norms_met = _convergence_criterion->isSatisfied();
374  }
375 
376  INFO("[time] Iteration #{:d} took {:g} s.", iteration,
377  time_iteration.elapsed());
378 
379  if (error_norms_met)
380  {
381  break;
382  }
383 
384  // Avoid increment of the 'iteration' if the error norms are not met,
385  // but maximum number of iterations is reached.
386  if (iteration >= _maxiter)
387  {
388  break;
389  }
390  }
391 
392  if (iteration > _maxiter)
393  {
394  ERR("Newton: Could not solve the given nonlinear system within {:d} "
395  "iterations",
396  _maxiter);
397  }
398 
402  minus_delta_x);
403 
404  return {error_norms_met, iteration};
405 }

References MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), BaseLib::RunTime::elapsed(), ERR(), NumLib::FAILURE, NumLib::MatrixProvider::getMatrix(), NumLib::VectorProvider::getVector(), INFO(), NumLib::GlobalVectorProvider::provider, NumLib::GlobalMatrixProvider::provider, NumLib::MatrixProvider::releaseMatrix(), NumLib::VectorProvider::releaseVector(), NumLib::REPEAT_ITERATION, BaseLib::RunTime::start(), and NumLib::SUCCESS.

Member Data Documentation

◆ _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 148 of file NonlinearSolver.h.

◆ _convergence_criterion

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

Definition at line 127 of file NonlinearSolver.h.

◆ _damping

const double 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 134 of file NonlinearSolver.h.

◆ _equation_system

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

Definition at line 124 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 138 of file NonlinearSolver.h.

◆ _linear_solver

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

Definition at line 123 of file NonlinearSolver.h.

◆ _maxiter

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

maximum number of iterations

Definition at line 128 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 139 of file NonlinearSolver.h.

◆ _r_neq

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

non-equilibrium initial residuum.

Definition at line 136 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 137 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 140 of file NonlinearSolver.h.


The documentation for this class was generated from the following files:
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_r_neq
GlobalVector * _r_neq
non-equilibrium initial residuum.
Definition: NonlinearSolver.h:136
NumLib::MatrixProvider::releaseMatrix
virtual void releaseMatrix(GlobalMatrix const &A)=0
NumLib::VectorProvider::getVector
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_res_id
std::size_t _res_id
ID of the residual vector.
Definition: NonlinearSolver.h:137
NumLib::NonlinearSystem< NonlinearSolverTag::Newton >::getResidual
virtual void getResidual(GlobalVector const &x, GlobalVector const &x_prev, GlobalVector &res) const =0
NumLib::ConvergenceCriterion::reset
virtual void reset()
Definition: ConvergenceCriterion.h:81
INFO
void INFO(char const *fmt, Args const &... args)
Definition: Logging.h:32
NumLib::VectorProvider::releaseVector
virtual void releaseVector(GlobalVector const &x)=0
NumLib::ConvergenceCriterion::isSatisfied
virtual bool isSatisfied() const
Tell if the convergence criterion is satisfied.
Definition: ConvergenceCriterion.h:84
NumLib::ConvergenceCriterion::hasResidualCheck
virtual bool hasResidualCheck() const =0
NumLib::GlobalVectorProvider::provider
static NUMLIB_EXPORT VectorProvider & provider
Definition: GlobalMatrixProviders.h:21
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_equation_system
System * _equation_system
Definition: NonlinearSolver.h:124
MathLib::LinAlg::copy
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_convergence_criterion
ConvergenceCriterion * _convergence_criterion
Definition: NonlinearSolver.h:127
NumLib::IterationResult::FAILURE
@ FAILURE
MathLib::LinAlg::axpy
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_J_id
std::size_t _J_id
ID of the Jacobian matrix.
Definition: NonlinearSolver.h:138
NumLib::ConvergenceCriterion::checkDeltaX
virtual void checkDeltaX(GlobalVector const &minus_delta_x, GlobalVector const &x)=0
NumLib::IterationResult::REPEAT_ITERATION
@ REPEAT_ITERATION
NumLib::GlobalMatrixProvider::provider
static NUMLIB_EXPORT MatrixProvider & provider
Definition: GlobalMatrixProviders.h:26
ERR
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_x_new_id
std::size_t _x_new_id
ID of the vector storing .
Definition: NonlinearSolver.h:140
NumLib::IterationResult::SUCCESS
@ SUCCESS
BaseLib::RunTime
Count the running time.
Definition: RunTime.h:27
MathLib::LinAlg
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_minus_delta_x_id
std::size_t _minus_delta_x_id
ID of the vector.
Definition: NonlinearSolver.h:139
BaseLib::RunTime::elapsed
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:41
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_maxiter
const int _maxiter
maximum number of iterations
Definition: NonlinearSolver.h:128
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_damping
const double _damping
Definition: NonlinearSolver.h:134
NumLib::ConvergenceCriterion::checkResidual
virtual void checkResidual(GlobalVector const &residual)=0
Check if the residual satisfies the convergence criterion.
NumLib::ConvergenceCriterion::preFirstIteration
virtual void preFirstIteration()
Definition: ConvergenceCriterion.h:68
NumLib::NonlinearSystem< NonlinearSolverTag::Newton >::assemble
virtual void assemble(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)=0
NumLib::ConvergenceCriterion::hasDeltaXCheck
virtual bool hasDeltaXCheck() const =0
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_linear_solver
GlobalLinearSolver & _linear_solver
Definition: NonlinearSolver.h:123
BaseLib::RunTime::start
void start()
Start the timer.
Definition: RunTime.h:31
NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_compensate_non_equilibrium_initial_residuum
bool _compensate_non_equilibrium_initial_residuum
Definition: NonlinearSolver.h:148
NumLib::MatrixProvider::getMatrix
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.