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 >:
[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. More...
 

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 More...
 
double const _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...
 
std::size_t _r_neq_id = 0u
 
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  }
int const _maxiter
maximum number of iterations

◆ ~NonlinearSolver()

Definition at line 465 of file NonlinearSolver.cpp.

466 {
467  if (_r_neq != nullptr)
468  {
470  }
471 }
GlobalVector * _r_neq
non-equilibrium initial residuum.
virtual void releaseVector(GlobalVector const &x)=0
static NUMLIB_EXPORT VectorProvider & provider

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

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

218 {
220  {
221  return;
222  }
223 
224  _equation_system->assemble(x, x_prev, process_id);
226  _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
227 }
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.

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

◆ compensateNonEquilibriumInitialResiduum()

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

Definition at line 119 of file NonlinearSolver.h.

120  {
122  }

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

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

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

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

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

◆ _convergence_criterion

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

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

◆ _equation_system

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

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

◆ _linear_solver

Definition at line 125 of file NonlinearSolver.h.

◆ _maxiter

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

maximum number of iterations

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

◆ _r_neq

non-equilibrium initial residuum.

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


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