OGS 6.2.1-720-gf9530c21a.dirty.20191210150148
NumLib::NonlinearSolver< NonlinearSolverTag::Newton > Class Template Referencefinal

Detailed Description

template<>
class NumLib::NonlinearSolver< NonlinearSolverTag::Newton >

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

Definition at line 86 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 assemble (std::vector< GlobalVector *> const &x, int const process_id) const override
 
void calculateNonEquilibriumInitialResiduum (std::vector< GlobalVector *> const &x, int const process_id) override
 
NonlinearSolverStatus solve (std::vector< GlobalVector *> &x, 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...
 
bool _compensate_non_equilibrium_initial_residuum = false
 

Member Typedef Documentation

◆ System

Type of the nonlinear equation system to be solved.

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

104  : _linear_solver(linear_solver), _maxiter(maxiter), _damping(damping)
105  {
106  }
int const _maxiter
maximum number of iterations

Member Function Documentation

◆ assemble()

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

Only assemble the equation system.

Note
This method is needed to preload CrankNicolson time discretization scheme. It is not used for the general solver steps; in those only the solve() method is needed.
Parameters
xthe state at which the equation system will be assembled.
process_idthe process' id which will be assembled.

Implements NumLib::NonlinearSolverBase.

Definition at line 214 of file NonlinearSolver.cpp.

References ProcessLib::calculateNonEquilibriumInitialResiduum().

216 {
217  _equation_system->assemble(x, process_id);
218  // TODO if the equation system would be reset to nullptr after each
219  // assemble() or solve() call, the user would be forced to set the
220  // equation every time and could not forget it.
221 }
virtual void assemble(std::vector< GlobalVector *> const &x, int const process_id)=0

◆ calculateNonEquilibriumInitialResiduum()

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

Implements NumLib::NonlinearSolverBase.

Definition at line 224 of file NonlinearSolver.cpp.

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

226 {
228  {
229  return;
230  }
231 
232  _equation_system->assemble(x, process_id);
234  _equation_system->getResidual(*x[process_id], *_r_neq);
235 }
GlobalVector * _r_neq
non-equilibrium initial residuum.
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
static NUMLIB_EXPORT VectorProvider & provider
virtual void assemble(std::vector< GlobalVector *> const &x, int const process_id)=0
virtual void getResidual(GlobalVector const &x, GlobalVector &res) const =0

◆ compensateNonEquilibriumInitialResiduum()

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

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

References ProcessLib::calculateNonEquilibriumInitialResiduum().

◆ solve()

NonlinearSolverStatus NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::solve ( std::vector< GlobalVector *> &  x,
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.
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 237 of file NonlinearSolver.cpp.

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

242 {
243  namespace LinAlg = MathLib::LinAlg;
244  auto& sys = *_equation_system;
245 
247  _res_id);
248  auto& minus_delta_x =
251  auto& J =
253 
254  bool error_norms_met = false;
255 
256  // TODO be more efficient
257  // init minus_delta_x to the right size
258  LinAlg::copy(*x[process_id], minus_delta_x);
259 
261 
262  int iteration = 1;
263  for (; iteration <= _maxiter;
264  ++iteration, _convergence_criterion->reset())
265  {
266  BaseLib::RunTime timer_dirichlet;
267  double time_dirichlet = 0.0;
268 
269  BaseLib::RunTime time_iteration;
270  time_iteration.start();
271 
272  timer_dirichlet.start();
273  sys.computeKnownSolutions(*x[process_id], process_id);
274  sys.applyKnownSolutions(*x[process_id]);
275  time_dirichlet += timer_dirichlet.elapsed();
276 
277  sys.preIteration(iteration, *x[process_id]);
278 
279  BaseLib::RunTime time_assembly;
280  time_assembly.start();
281  try
282  {
283  sys.assemble(x, process_id);
284  }
285  catch (AssemblyException const& e)
286  {
287  ERR("Abort nonlinear iteration. Repeating timestep. Reason: %s",
288  e.what());
289  error_norms_met = false;
290  iteration = _maxiter;
291  break;
292  }
293  sys.getResidual(*x[process_id], res);
294  sys.getJacobian(J);
295  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
296 
297  // Subract non-equilibrium initial residuum if set
298  if (_r_neq != nullptr)
299  LinAlg::axpy(res, -1, *_r_neq);
300 
301  minus_delta_x.setZero();
302 
303  timer_dirichlet.start();
304  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
305  time_dirichlet += timer_dirichlet.elapsed();
306  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
307 
308  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
309  {
311  }
312 
313  BaseLib::RunTime time_linear_solver;
314  time_linear_solver.start();
315  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
316  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
317 
318  if (!iteration_succeeded)
319  {
320  ERR("Newton: The linear solver failed.");
321  }
322  else
323  {
324  // TODO could be solved in a better way
325  // cf.
326  // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
327 
328  // Copy pointers, replace the one for the given process id.
329  std::vector<GlobalVector*> x_new{x};
330  x_new[process_id] =
332  *x[process_id], _x_new_id);
333  LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
334 
335  if (postIterationCallback)
336  {
337  postIterationCallback(iteration, x_new);
338  }
339 
340  switch (sys.postIteration(*x_new[process_id]))
341  {
343  break;
345  ERR("Newton: The postIteration() hook reported a "
346  "non-recoverable error.");
347  iteration_succeeded = false;
348  break;
350  INFO(
351  "Newton: The postIteration() hook decided that this "
352  "iteration"
353  " has to be repeated.");
354  // TODO introduce some onDestroy hook.
356  *x_new[process_id]);
357  continue; // That throws the iteration result away.
358  }
359 
360  LinAlg::copy(*x_new[process_id],
361  *x[process_id]); // copy new solution to x
363  *x_new[process_id]);
364  }
365 
366  if (!iteration_succeeded)
367  {
368  // Don't compute further error norms, but break here.
369  error_norms_met = false;
370  break;
371  }
372 
373  if (sys.isLinear()) {
374  error_norms_met = true;
375  } else {
377  // Note: x contains the new solution!
378  _convergence_criterion->checkDeltaX(minus_delta_x,
379  *x[process_id]);
380  }
381 
382  error_norms_met = _convergence_criterion->isSatisfied();
383  }
384 
385  INFO("[time] Iteration #%u took %g s.", iteration,
386  time_iteration.elapsed());
387 
388  if (error_norms_met)
389  {
390  break;
391  }
392 
393  // Avoid increment of the 'iteration' if the error norms are not met,
394  // but maximum number of iterations is reached.
395  if (iteration >= _maxiter)
396  {
397  break;
398  }
399  }
400 
401  if (iteration > _maxiter)
402  {
403  ERR("Newton: Could not solve the given nonlinear system within %u "
404  "iterations",
405  _maxiter);
406  }
407 
411  minus_delta_x);
412 
413  return {error_norms_met, iteration};
414 }
GlobalVector * _r_neq
non-equilibrium initial residuum.
std::size_t _J_id
ID of the Jacobian matrix.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
static NUMLIB_EXPORT MatrixProvider & provider
virtual bool hasResidualCheck() const =0
static NUMLIB_EXPORT VectorProvider & provider
std::size_t _x_new_id
ID of the vector storing .
int const _maxiter
maximum number of iterations
virtual void checkDeltaX(GlobalVector const &minus_delta_x, GlobalVector const &x)=0
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
Count the running time.
Definition: RunTime.h:31
void start()
Start the timer.
Definition: RunTime.h:35
std::size_t _res_id
ID of the residual vector.
virtual bool hasDeltaXCheck() const =0
virtual void releaseVector(GlobalVector const &x)=0
virtual bool isSatisfied() const
Tell if the convergence criterion is satisfied.
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
virtual void checkResidual(GlobalVector const &residual)=0
Check if the residual satisfies the convergence criterion.
double elapsed() const
Get the elapsed time after started.
Definition: RunTime.h:51

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

◆ _convergence_criterion

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

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

◆ _equation_system

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

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

◆ _linear_solver

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

Definition at line 134 of file NonlinearSolver.h.

◆ _maxiter

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

maximum number of iterations

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

◆ _r_neq

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

non-equilibrium initial residuum.

Definition at line 147 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 148 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 151 of file NonlinearSolver.h.


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