OGS  v6.4.0
NumLib::NonlinearSolver< NonlinearSolverTag::Picard > Class Referencefinal

## Detailed Description

Find a solution to a nonlinear equation using the Picard fixpoint iteration method.

Definition at line 160 of file NonlinearSolver.h.

#include <NonlinearSolver.h>

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

## Public Types

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

## Public Member Functions

NonlinearSolver (GlobalLinearSolver &linear_solver, const int maxiter)

~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

const int _maxiter
maximum number of iterations More...

GlobalVector_r_neq = nullptr
non-equilibrium initial residuum. More...

std::size_t _A_id = 0u
ID of the $$A$$ matrix. More...

std::size_t _rhs_id = 0u
ID of the right-hand side vector. More...

std::size_t _x_new_id = 0u

std::size_t _r_neq_id = 0u

bool _compensate_non_equilibrium_initial_residuum = false

## ◆ System

Type of the nonlinear equation system to be solved.

Definition at line 165 of file NonlinearSolver.h.

## ◆ NonlinearSolver()

 NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::NonlinearSolver ( GlobalLinearSolver & linear_solver, const int maxiter )
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.

Definition at line 173 of file NonlinearSolver.h.

175  : _linear_solver(linear_solver), _maxiter(maxiter)
176  {
177  }
const int _maxiter
maximum number of iterations

## ◆ ~NonlinearSolver()

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

Definition at line 456 of file NonlinearSolver.cpp.

457 {
458  if (_r_neq != nullptr)
459  {
461  }
462 }
GlobalVector * _r_neq
non-equilibrium initial residuum.
virtual void releaseVector(GlobalVector const &x)=0
static NUMLIB_EXPORT VectorProvider & provider

## ◆ calculateNonEquilibriumInitialResiduum()

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

Implements NumLib::NonlinearSolverBase.

Definition at line 27 of file NonlinearSolver.cpp.

31 {
33  {
34  return;
35  }
36
39  _equation_system->assemble(x, x_prev, process_id);
41  _equation_system->getRhs(*x_prev[process_id], rhs);
42
43  // r_neq = A * x - rhs
45  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
46  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
47 }
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
std::size_t _rhs_id
ID of the right-hand side vector.
virtual void getA(GlobalMatrix &A) const =0
virtual void assemble(std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id)=0
virtual void getRhs(GlobalVector const &x_prev, GlobalVector &rhs) const =0
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
void axpy(PETScVector &y, double const a, PETScVector const &x)
Definition: LinAlg.cpp:55
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:143
static NUMLIB_EXPORT MatrixProvider & provider

## ◆ compensateNonEquilibriumInitialResiduum()

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

Definition at line 201 of file NonlinearSolver.h.

202  {
204  }

## ◆ setEquationSystem()

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

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

Definition at line 183 of file NonlinearSolver.h.

184  {
185  _equation_system = &eq;
186  _convergence_criterion = &conv_crit;
187  }

## ◆ solve()

 NonlinearSolverStatus NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::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 49 of file NonlinearSolver.cpp.

55 {
56  namespace LinAlg = MathLib::LinAlg;
57  auto& sys = *_equation_system;
58
59  auto& A =
62  _rhs_id);
63
64  std::vector<GlobalVector*> x_new{x};
65  x_new[process_id] =
67  LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
68
69  bool error_norms_met = false;
70
72
73  int iteration = 1;
74  for (; iteration <= _maxiter;
75  ++iteration, _convergence_criterion->reset())
76  {
77  BaseLib::RunTime timer_dirichlet;
78  double time_dirichlet = 0.0;
79
80  BaseLib::RunTime time_iteration;
81  time_iteration.start();
82
83  timer_dirichlet.start();
84  sys.computeKnownSolutions(*x_new[process_id], process_id);
85  sys.applyKnownSolutions(*x_new[process_id]);
86  time_dirichlet += timer_dirichlet.elapsed();
87
88  sys.preIteration(iteration, *x_new[process_id]);
89
90  BaseLib::RunTime time_assembly;
91  time_assembly.start();
92  sys.assemble(x_new, x_prev, process_id);
93  sys.getA(A);
94  sys.getRhs(*x_prev[process_id], rhs);
95  INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
96
97  // Subtract non-equilibrium initial residuum if set
98  if (_r_neq != nullptr)
99  {
100  LinAlg::axpy(rhs, -1, *_r_neq);
101  }
102
103  timer_dirichlet.start();
104  sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
105  time_dirichlet += timer_dirichlet.elapsed();
106  INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
107
108  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
109  GlobalVector res;
110  LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
111  LinAlg::axpy(res, -1.0, rhs); // res -= rhs
113  }
114
115  BaseLib::RunTime time_linear_solver;
116  time_linear_solver.start();
117  bool iteration_succeeded =
118  _linear_solver.solve(A, rhs, *x_new[process_id]);
119  INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
120
121  if (!iteration_succeeded)
122  {
123  ERR("Picard: The linear solver failed.");
124  }
125  else
126  {
127  if (postIterationCallback)
128  {
129  postIterationCallback(iteration, x_new);
130  }
131
132  switch (sys.postIteration(*x_new[process_id]))
133  {
135  // Don't copy here. The old x might still be used further
136  // below. Although currently it is not.
137  break;
139  ERR("Picard: The postIteration() hook reported a "
140  "non-recoverable error.");
141  iteration_succeeded = false;
142  // Copy new solution to x.
143  // Thereby the failed solution can be used by the caller for
144  // debugging purposes.
145  LinAlg::copy(*x_new[process_id], *x[process_id]);
146  break;
148  INFO(
149  "Picard: The postIteration() hook decided that this "
150  "iteration has to be repeated.");
151  LinAlg::copy(
152  *x[process_id],
153  *x_new[process_id]); // throw the iteration result away
154  continue;
155  }
156  }
157
158  if (!iteration_succeeded)
159  {
160  // Don't compute error norms, break here.
161  error_norms_met = false;
162  break;
163  }
164
165  if (sys.isLinear()) {
166  error_norms_met = true;
167  } else {
169  GlobalVector minus_delta_x(*x[process_id]);
170  LinAlg::axpy(minus_delta_x, -1.0,
171  *x_new[process_id]); // minus_delta_x = x - x_new
172  _convergence_criterion->checkDeltaX(minus_delta_x,
173  *x_new[process_id]);
174  }
175
176  error_norms_met = _convergence_criterion->isSatisfied();
177  }
178
179  // Update x s.t. in the next iteration we will compute the right delta x
180  LinAlg::copy(*x_new[process_id], *x[process_id]);
181
182  INFO("[time] Iteration #{:d} took {:g} s.", iteration,
183  time_iteration.elapsed());
184
185  if (error_norms_met)
186  {
187  break;
188  }
189
190  // Avoid increment of the 'iteration' if the error norms are not met,
191  // but maximum number of iterations is reached.
192  if (iteration >= _maxiter)
193  {
194  break;
195  }
196  }
197
198  if (iteration > _maxiter)
199  {
200  ERR("Picard: Could not solve the given nonlinear system within {:d} "
201  "iterations",
202  _maxiter);
203  }
204
208
209  return {error_norms_met, iteration};
210 }
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:27
double elapsed() const
Get the elapsed time in seconds.
Definition: RunTime.h:40
void start()
Start the timer.
Definition: RunTime.h:30
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
Global vector based on Eigen vector.
Definition: EigenVector.h:27
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
void copy(PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:36

## ◆ _A_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_A_id = 0u
private

ID of the $$A$$ matrix.

Definition at line 215 of file NonlinearSolver.h.

## ◆ _compensate_non_equilibrium_initial_residuum

 bool NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_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 224 of file NonlinearSolver.h.

## ◆ _convergence_criterion

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

Definition at line 211 of file NonlinearSolver.h.

## ◆ _equation_system

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

Definition at line 208 of file NonlinearSolver.h.

## ◆ _linear_solver

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

Definition at line 207 of file NonlinearSolver.h.

## ◆ _maxiter

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

maximum number of iterations

Definition at line 212 of file NonlinearSolver.h.

## ◆ _r_neq

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

non-equilibrium initial residuum.

Definition at line 214 of file NonlinearSolver.h.

## ◆ _r_neq_id

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

ID of the non-equilibrium initial residuum vector.

Definition at line 219 of file NonlinearSolver.h.

## ◆ _rhs_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_rhs_id = 0u
private

ID of the right-hand side vector.

Definition at line 216 of file NonlinearSolver.h.

## ◆ _x_new_id

 std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_x_new_id = 0u
private

ID of the vector storing the solution of the linearized equation.

Definition at line 217 of file NonlinearSolver.h.

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