OGS
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 >:
[legend]
Collaboration diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Picard >:
[legend]

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
 

Member Typedef Documentation

◆ System

Type of the nonlinear equation system to be solved.

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

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

◆ ~NonlinearSolver()

Definition at line 486 of file NonlinearSolver.cpp.

487 {
488  if (_r_neq != nullptr)
489  {
491  }
492 }
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::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 
37  INFO("Calculate non-equilibrium initial residuum.");
38 
41  _equation_system->assemble(x, x_prev, process_id);
43  _equation_system->getRhs(*x_prev[process_id], rhs);
44 
45  // r_neq = A * x - rhs
47  MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
48  MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
49 
50  // Set the values of the selected entries of _r_neq, which are associated
51  // with the equations that do not need initial residual compensation, to
52  // zero.
53  const auto selected_global_indices =
55 
56  std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
57  _r_neq->set(selected_global_indices, zero_entries);
58 
60 
63 }
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
Definition: Logging.h:34
void set(IndexType rowId, double v)
set entry
Definition: EigenVector.h:76
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 _rhs_id
ID of the right-hand side vector.
virtual void getA(GlobalMatrix &A) const =0
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 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 finalizeAssembly(PETScMatrix &A)
Definition: LinAlg.cpp:163
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
Definition: LinAlg.cpp:142
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::finalizeAssembly(), NumLib::MatrixProvider::getMatrix(), NumLib::VectorProvider::getVector(), INFO(), MathLib::LinAlg::matMult(), NumLib::GlobalVectorProvider::provider, NumLib::GlobalMatrixProvider::provider, NumLib::MatrixProvider::releaseMatrix(), and NumLib::VectorProvider::releaseVector().

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

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

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

Member Data Documentation

◆ _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

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

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: