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

Detailed Description

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

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

Definition at line 167 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)
 
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
 
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
 
bool _compensate_non_equilibrium_initial_residuum = false
 

Member Typedef Documentation

◆ System

Type of the nonlinear equation system to be solved.

Definition at line 172 of file NonlinearSolver.h.

Constructor & Destructor Documentation

◆ NonlinearSolver()

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

Definition at line 180 of file NonlinearSolver.h.

182  : _linear_solver(linear_solver), _maxiter(maxiter)
183  {
184  }
const int _maxiter
maximum number of iterations

Member Function Documentation

◆ assemble()

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

References ProcessLib::calculateNonEquilibriumInitialResiduum().

27 {
28  _equation_system->assemble(x, process_id);
29 }
virtual void assemble(std::vector< GlobalVector *> const &x, int const process_id)=0

◆ calculateNonEquilibriumInitialResiduum()

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

Implements NumLib::NonlinearSolverBase.

Definition at line 32 of file NonlinearSolver.cpp.

References MathLib::LinAlg::axpy(), NumLib::MatrixProvider::getMatrix(), NumLib::VectorProvider::getVector(), MathLib::LinAlg::matMult(), NumLib::GlobalVectorProvider::provider, and NumLib::GlobalMatrixProvider::provider.

34 {
36  {
37  return;
38  }
39 
42  _equation_system->assemble(x, process_id);
45 
46  // r_neq = A * x - rhs
48  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
49  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
50 }
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
std::size_t _rhs_id
ID of the right-hand side vector.
GlobalVector * _r_neq
non-equilibrium initial residuum.
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
virtual void assemble(std::vector< GlobalVector *> const &x, int const process_id)=0
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
virtual void getA(GlobalMatrix &A) const =0
virtual void getRhs(GlobalVector &rhs) const =0

◆ compensateNonEquilibriumInitialResiduum()

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

Definition at line 206 of file NonlinearSolver.h.

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

References ProcessLib::calculateNonEquilibriumInitialResiduum().

◆ solve()

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

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

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

◆ _A_id

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

ID of the $ A $ matrix.

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

◆ _convergence_criterion

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

Definition at line 216 of file NonlinearSolver.h.

◆ _equation_system

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

Definition at line 213 of file NonlinearSolver.h.

◆ _linear_solver

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

Definition at line 212 of file NonlinearSolver.h.

◆ _maxiter

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

maximum number of iterations

Definition at line 217 of file NonlinearSolver.h.

◆ _r_neq

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

non-equilibrium initial residuum.

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


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