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 161 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.
 

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 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
 
const int _maxiter
 maximum number of iterations
 
GlobalVector_r_neq = nullptr
 non-equilibrium initial residuum.
 
std::size_t _A_id = 0u
 ID of the \( A \) matrix.
 
std::size_t _rhs_id = 0u
 ID of the right-hand side vector.
 
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 166 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 174 of file NonlinearSolver.h.

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

◆ ~NonlinearSolver()

Definition at line 543 of file NonlinearSolver.cpp.

544{
545 if (_r_neq != nullptr)
546 {
548 }
549}
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 88 of file NonlinearSolver.cpp.

92{
94 {
95 return;
96 }
97
98 INFO("Calculate non-equilibrium initial residuum.");
99
102 _equation_system->assemble(x, x_prev, process_id);
104 _equation_system->getRhs(*x_prev[process_id], rhs);
105
106 // r_neq = A * x - rhs
108 MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
109 MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs); // res -= rhs
110
111 // Set the values of the selected entries of _r_neq, which are associated
112 // with the equations that do not need initial residual compensation, to
113 // zero.
114 const auto selected_global_indices =
116
117 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
118 _r_neq->set(selected_global_indices, zero_entries);
119
121
124}
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 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 _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 202 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 184 of file NonlinearSolver.h.

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

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

132{
133 namespace LinAlg = MathLib::LinAlg;
134 auto& sys = *_equation_system;
135
138
139 std::vector<GlobalVector*> x_new{x};
140 x_new[process_id] =
142 LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
143
144 bool error_norms_met = false;
145
147
148 int iteration = 1;
149 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
150 {
151 BaseLib::RunTime timer_dirichlet;
152 double time_dirichlet = 0.0;
153
154 BaseLib::RunTime time_iteration;
155 time_iteration.start();
156
157 timer_dirichlet.start();
158 auto& x_new_process = *x_new[process_id];
160 sys.computeKnownSolutions(x_new_process, process_id);
161 sys.applyKnownSolutions(x_new_process);
162 time_dirichlet += timer_dirichlet.elapsed();
163
164 sys.preIteration(iteration, x_new_process);
165
166 BaseLib::RunTime time_assembly;
167 time_assembly.start();
168 sys.assemble(x_new, x_prev, process_id);
169 sys.getA(A);
170 sys.getRhs(*x_prev[process_id], rhs);
171 INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
172
173 // Subtract non-equilibrium initial residuum if set
174 if (_r_neq != nullptr)
175 {
176 LinAlg::axpy(rhs, -1, *_r_neq);
177 }
178
179 timer_dirichlet.start();
180 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
181 time_dirichlet += timer_dirichlet.elapsed();
182 INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
183
184 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
185 {
186 GlobalVector res;
187 LinAlg::matMult(A, x_new_process, res); // res = A * x_new
188 LinAlg::axpy(res, -1.0, rhs); // res -= rhs
190 }
191
192 bool iteration_succeeded =
193 detail::solvePicard(_linear_solver, A, rhs, x_new_process,
194 sys.linearSolverNeedsToCompute());
195
196 if (iteration_succeeded)
197 {
198 if (postIterationCallback)
199 {
200 postIterationCallback(iteration, x_new);
201 }
202
203 switch (sys.postIteration(x_new_process))
204 {
206 // Don't copy here. The old x might still be used further
207 // below. Although currently it is not.
208 break;
210 ERR("Picard: The postIteration() hook reported a "
211 "non-recoverable error.");
212 iteration_succeeded = false;
213 // Copy new solution to x.
214 // Thereby the failed solution can be used by the caller for
215 // debugging purposes.
216 LinAlg::copy(x_new_process, *x[process_id]);
217 break;
219 INFO(
220 "Picard: The postIteration() hook decided that this "
221 "iteration has to be repeated.");
223 *x[process_id],
224 x_new_process); // throw the iteration result away
225 continue;
226 }
227 }
228
229 if (!iteration_succeeded)
230 {
231 // Don't compute error norms, break here.
232 error_norms_met = false;
233 break;
234 }
235
236 if (sys.isLinear())
237 {
238 error_norms_met = true;
239 }
240 else
241 {
243 {
244 GlobalVector minus_delta_x(*x[process_id]);
245 LinAlg::axpy(minus_delta_x, -1.0,
246 x_new_process); // minus_delta_x = x - x_new
247 _convergence_criterion->checkDeltaX(minus_delta_x,
248 x_new_process);
249 }
250
251 error_norms_met = _convergence_criterion->isSatisfied();
252 }
253
254 // Update x s.t. in the next iteration we will compute the right delta x
255 LinAlg::copy(x_new_process, *x[process_id]);
256
257 INFO("[time] Iteration #{:d} took {:g} s.", iteration,
258 time_iteration.elapsed());
259
260 if (error_norms_met)
261 {
262 break;
263 }
264
265 // Avoid increment of the 'iteration' if the error norms are not met,
266 // but maximum number of iterations is reached.
267 if (iteration >= _maxiter)
268 {
269 break;
270 }
271 }
272
273 if (iteration > _maxiter)
274 {
275 ERR("Picard: Could not solve the given nonlinear system within {:d} "
276 "iterations",
277 _maxiter);
278 }
279
283
284 return {error_norms_met, iteration};
285}
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
Global vector based on Eigen vector.
Definition: EigenVector.h:25
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
void setLocalAccessibleVector(PETScVector const &x)
Definition: LinAlg.cpp:27
bool solvePicard(GlobalLinearSolver &linear_solver, GlobalMatrix &A, GlobalVector &rhs, GlobalVector &x, bool const compute_necessary)

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, MathLib::LinAlg::setLocalAccessibleVector(), NumLib::detail::solvePicard(), 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 216 of file NonlinearSolver.h.

◆ _compensate_non_equilibrium_initial_residuum

bool NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_compensate_non_equilibrium_initial_residuum = false
private

Definition at line 225 of file NonlinearSolver.h.

◆ _convergence_criterion

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

Definition at line 212 of file NonlinearSolver.h.

◆ _equation_system

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

Definition at line 209 of file NonlinearSolver.h.

◆ _linear_solver

Definition at line 208 of file NonlinearSolver.h.

◆ _maxiter

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

maximum number of iterations

Definition at line 213 of file NonlinearSolver.h.

◆ _r_neq

non-equilibrium initial residuum.

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


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