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.
 

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 void calculateNonEquilibriumInitialResiduum (std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)=0
 
virtual 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)=0
 
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
 

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

497{
498 if (_r_neq != nullptr)
499 {
501 }
502}
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 232 of file NonlinearSolver.cpp.

236{
238 {
239 return;
240 }
241
242 INFO("Calculate non-equilibrium initial residuum.");
243
244 _equation_system->assemble(x, x_prev, process_id);
246 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
247
248 // Set the values of the selected entries of _r_neq, which are associated
249 // with the equations that do not need initial residual compensation, to
250 // zero.
251 const auto selected_global_indices =
253 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
254
255 _r_neq->set(selected_global_indices, zero_entries);
256
258}
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:163

References MathLib::LinAlg::finalizeAssembly(), NumLib::VectorProvider::getVector(), INFO(), and NumLib::GlobalVectorProvider::provider.

◆ compensateNonEquilibriumInitialResiduum()

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

Definition at line 119 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 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 260 of file NonlinearSolver.cpp.

266{
267 namespace LinAlg = MathLib::LinAlg;
268 auto& sys = *_equation_system;
269
271 auto& minus_delta_x =
274
275 bool error_norms_met = false;
276
277 // TODO be more efficient
278 // init minus_delta_x to the right size
279 LinAlg::copy(*x[process_id], minus_delta_x);
280
282
283 int iteration = 1;
284 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
285 {
286 BaseLib::RunTime timer_dirichlet;
287 double time_dirichlet = 0.0;
288
289 BaseLib::RunTime time_iteration;
290 time_iteration.start();
291
292 timer_dirichlet.start();
293 sys.computeKnownSolutions(*x[process_id], process_id);
294 sys.applyKnownSolutions(*x[process_id]);
295 time_dirichlet += timer_dirichlet.elapsed();
296
297 sys.preIteration(iteration, *x[process_id]);
298
299 BaseLib::RunTime time_assembly;
300 time_assembly.start();
301 try
302 {
303 sys.assemble(x, x_prev, process_id);
304 }
305 catch (AssemblyException const& e)
306 {
307 ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
308 e.what());
309 error_norms_met = false;
310 iteration = _maxiter;
311 break;
312 }
313 sys.getResidual(*x[process_id], *x_prev[process_id], res);
314 sys.getJacobian(J);
315 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
316
317 // Subtract non-equilibrium initial residuum if set
318 if (_r_neq != nullptr)
319 LinAlg::axpy(res, -1, *_r_neq);
320
321 minus_delta_x.setZero();
322
323 timer_dirichlet.start();
324 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
325 time_dirichlet += timer_dirichlet.elapsed();
326 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
327
328 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
329 {
331 }
332
333 BaseLib::RunTime time_linear_solver;
334 time_linear_solver.start();
335 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
336 INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
337
338 if (!iteration_succeeded)
339 {
340 ERR("Newton: The linear solver failed.");
341 }
342 else
343 {
344 // TODO could be solved in a better way
345 // cf.
346 // https://www.petsc.org/release/docs/manualpages/Vec/VecWAXPY.html
347
348 // Copy pointers, replace the one for the given process id.
349 std::vector<GlobalVector*> x_new{x};
350 x_new[process_id] =
352 *x[process_id], _x_new_id);
353 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
354
355 if (postIterationCallback)
356 {
357 postIterationCallback(iteration, x_new);
358 }
359
360 switch (sys.postIteration(*x_new[process_id]))
361 {
363 break;
365 ERR("Newton: The postIteration() hook reported a "
366 "non-recoverable error.");
367 iteration_succeeded = false;
368 break;
370 INFO(
371 "Newton: The postIteration() hook decided that this "
372 "iteration has to be repeated.");
373 // TODO introduce some onDestroy hook.
375 *x_new[process_id]);
376 continue; // That throws the iteration result away.
377 }
378
379 LinAlg::copy(*x_new[process_id],
380 *x[process_id]); // copy new solution to x
382 *x_new[process_id]);
383 }
384
385 if (!iteration_succeeded)
386 {
387 // Don't compute further error norms, but break here.
388 error_norms_met = false;
389 break;
390 }
391
392 if (sys.isLinear())
393 {
394 error_norms_met = true;
395 }
396 else
397 {
399 {
400 // Note: x contains the new solution!
401 _convergence_criterion->checkDeltaX(minus_delta_x,
402 *x[process_id]);
403 }
404
405 error_norms_met = _convergence_criterion->isSatisfied();
406 }
407
408 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
409 time_iteration.elapsed());
410
411 if (error_norms_met)
412 {
413 break;
414 }
415
416 // Avoid increment of the 'iteration' if the error norms are not met,
417 // but maximum number of iterations is reached.
418 if (iteration >= _maxiter)
419 {
420 break;
421 }
422 }
423
424 if (iteration > _maxiter)
425 {
426 ERR("Newton: Could not solve the given nonlinear system within {:d} "
427 "iterations",
428 _maxiter);
429 }
430
434
435 return {error_norms_met, iteration};
436}
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

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: