13#include <boost/algorithm/string.hpp>
29#if !defined(USE_PETSC) && !defined(USE_LIS)
32 bool const compute_necessary)
35 time_linear_solver.
start();
37 if (compute_necessary)
39 if (!linear_solver.compute(A))
41 ERR(
"Picard: The linear solver failed in the compute() step.");
46 bool const iteration_succeeded = linear_solver.
solve(rhs, x);
48 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
50 if (iteration_succeeded)
55 ERR(
"Picard: The linear solver failed in the solve() step.");
61 bool const compute_necessary)
63 if (!compute_necessary)
66 "The performance optimization to skip the linear solver compute() "
67 "step is not implemented for PETSc or LIS linear solvers.");
71 time_linear_solver.
start();
73 bool const iteration_succeeded = linear_solver.
solve(A, rhs, x);
75 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
77 if (iteration_succeeded)
82 ERR(
"Picard: The linear solver failed in the solve() step.");
90 std::vector<GlobalVector*>
const& x,
91 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
93 if (!_compensate_non_equilibrium_initial_residuum)
98 INFO(
"Calculate non-equilibrium initial residuum.");
102 _equation_system->assemble(x, x_prev, process_id);
103 _equation_system->getA(A);
104 _equation_system->getRhs(*x_prev[process_id], rhs);
114 const auto selected_global_indices =
115 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
117 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
118 _r_neq->set(selected_global_indices, zero_entries);
127 std::vector<GlobalVector*>& x,
128 std::vector<GlobalVector*>
const& x_prev,
129 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
130 postIterationCallback,
131 int const process_id)
134 auto& sys = *_equation_system;
139 std::vector<GlobalVector*> x_new{x};
144 bool error_norms_met =
false;
146 _convergence_criterion->preFirstIteration();
149 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
152 double time_dirichlet = 0.0;
155 time_iteration.
start();
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();
164 sys.preIteration(iteration, x_new_process);
167 time_assembly.
start();
168 sys.assemble(x_new, x_prev, process_id);
170 sys.getRhs(*x_prev[process_id], rhs);
171 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
174 if (_r_neq !=
nullptr)
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);
184 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
189 _convergence_criterion->checkResidual(res);
192 bool iteration_succeeded =
194 sys.linearSolverNeedsToCompute());
196 if (iteration_succeeded)
198 if (postIterationCallback)
200 postIterationCallback(iteration, x_new);
203 switch (sys.postIteration(x_new_process))
210 ERR(
"Picard: The postIteration() hook reported a "
211 "non-recoverable error.");
212 iteration_succeeded =
false;
220 "Picard: The postIteration() hook decided that this "
221 "iteration has to be repeated.");
229 if (!iteration_succeeded)
232 error_norms_met =
false;
238 error_norms_met =
true;
242 if (_convergence_criterion->hasDeltaXCheck())
247 _convergence_criterion->checkDeltaX(minus_delta_x,
251 error_norms_met = _convergence_criterion->isSatisfied();
257 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
267 if (iteration >= _maxiter)
273 if (iteration > _maxiter)
275 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
284 return {error_norms_met, iteration};
289 std::vector<GlobalVector*>
const& x,
290 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
292 if (!_compensate_non_equilibrium_initial_residuum)
297 INFO(
"Calculate non-equilibrium initial residuum.");
299 _equation_system->assemble(x, x_prev, process_id);
301 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
306 const auto selected_global_indices =
307 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
308 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
310 _r_neq->set(selected_global_indices, zero_entries);
316 std::vector<GlobalVector*>& x,
317 std::vector<GlobalVector*>
const& x_prev,
318 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
319 postIterationCallback,
320 int const process_id)
323 auto& sys = *_equation_system;
326 auto& minus_delta_x =
330 bool error_norms_met =
false;
336 _convergence_criterion->preFirstIteration();
339 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
342 double time_dirichlet = 0.0;
345 time_iteration.
start();
347 timer_dirichlet.
start();
348 sys.computeKnownSolutions(*x[process_id], process_id);
349 sys.applyKnownSolutions(*x[process_id]);
350 time_dirichlet += timer_dirichlet.
elapsed();
352 sys.preIteration(iteration, *x[process_id]);
355 time_assembly.
start();
358 sys.assemble(x, x_prev, process_id);
362 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
364 error_norms_met =
false;
365 iteration = _maxiter;
368 sys.getResidual(*x[process_id], *x_prev[process_id], res);
370 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
373 if (_r_neq !=
nullptr)
376 minus_delta_x.setZero();
378 timer_dirichlet.
start();
379 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
380 time_dirichlet += timer_dirichlet.
elapsed();
381 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
383 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
385 _convergence_criterion->checkResidual(res);
389 time_linear_solver.
start();
390 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
391 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
393 if (!iteration_succeeded)
395 ERR(
"Newton: The linear solver failed.");
404 std::vector<GlobalVector*> x_new{x};
407 *x[process_id], _x_new_id);
408 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
410 if (postIterationCallback)
412 postIterationCallback(iteration, x_new);
415 switch (sys.postIteration(*x_new[process_id]))
420 ERR(
"Newton: The postIteration() hook reported a "
421 "non-recoverable error.");
422 iteration_succeeded =
false;
426 "Newton: The postIteration() hook decided that this "
427 "iteration has to be repeated.");
440 if (!iteration_succeeded)
443 error_norms_met =
false;
449 error_norms_met =
true;
453 if (_convergence_criterion->hasDeltaXCheck())
456 _convergence_criterion->checkDeltaX(minus_delta_x,
460 error_norms_met = _convergence_criterion->isSatisfied();
463 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
473 if (iteration >= _maxiter)
479 if (iteration > _maxiter)
481 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
490 return {error_norms_met, iteration};
502 if (type ==
"Picard")
506 return std::make_pair(
507 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
509 if (type ==
"Newton")
516 "The damping factor for the Newon method must be positive, got "
522 return std::make_pair(
523 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
527 if (boost::iequals(type,
"PETScSNES"))
534 return std::make_pair(std::make_unique<ConcreteNLS>(
535 linear_solver, max_iter, std::move(prefix)),
540 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());
545 if (_r_neq !=
nullptr)
553 if (_r_neq !=
nullptr)
void INFO(fmt::format_string< Args... > fmt, Args &&... args)
void ERR(fmt::format_string< Args... > fmt, Args &&... args)
void WARN(fmt::format_string< Args... > fmt, Args &&... args)
Definition of the RunTime class.
T getConfigParameter(std::string const ¶m) const
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
bool solve(EigenMatrix &A, EigenVector &b, EigenVector &x)
Global vector based on Eigen vector.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
std::pair< std::unique_ptr< NonlinearSolverBase >, NonlinearSolverTag > createNonlinearSolver(GlobalLinearSolver &linear_solver, BaseLib::ConfigTree const &config)
void finalizeAssembly(PETScMatrix &A)
void copy(PETScVector const &x, PETScVector &y)
void setLocalAccessibleVector(PETScVector const &x)
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
void axpy(PETScVector &y, PetscScalar const a, PETScVector const &x)
bool solvePicard(GlobalLinearSolver &linear_solver, GlobalMatrix &A, GlobalVector &rhs, GlobalVector &x, bool const compute_necessary)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.