13 #include <boost/algorithm/string.hpp>
29 std::vector<GlobalVector*>
const& x,
30 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
32 if (!_compensate_non_equilibrium_initial_residuum)
39 _equation_system->assemble(x, x_prev, process_id);
40 _equation_system->getA(A);
41 _equation_system->getRhs(*x_prev[process_id], rhs);
50 std::vector<GlobalVector*>& x,
51 std::vector<GlobalVector*>
const& x_prev,
52 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
53 postIterationCallback,
57 auto& sys = *_equation_system;
64 std::vector<GlobalVector*> x_new{x};
69 bool error_norms_met =
false;
71 _convergence_criterion->preFirstIteration();
74 for (; iteration <= _maxiter;
75 ++iteration, _convergence_criterion->reset())
78 double time_dirichlet = 0.0;
81 time_iteration.
start();
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();
88 sys.preIteration(iteration, *x_new[process_id]);
91 time_assembly.
start();
92 sys.assemble(x_new, x_prev, process_id);
94 sys.getRhs(*x_prev[process_id], rhs);
95 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
98 if (_r_neq !=
nullptr)
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);
108 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
112 _convergence_criterion->checkResidual(res);
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());
121 if (!iteration_succeeded)
123 ERR(
"Picard: The linear solver failed.");
127 if (postIterationCallback)
129 postIterationCallback(iteration, x_new);
132 switch (sys.postIteration(*x_new[process_id]))
139 ERR(
"Picard: The postIteration() hook reported a "
140 "non-recoverable error.");
141 iteration_succeeded =
false;
149 "Picard: The postIteration() hook decided that this "
150 "iteration has to be repeated.");
158 if (!iteration_succeeded)
161 error_norms_met =
false;
165 if (sys.isLinear()) {
166 error_norms_met =
true;
168 if (_convergence_criterion->hasDeltaXCheck()) {
172 _convergence_criterion->checkDeltaX(minus_delta_x,
176 error_norms_met = _convergence_criterion->isSatisfied();
182 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
192 if (iteration >= _maxiter)
198 if (iteration > _maxiter)
200 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
209 return {error_norms_met, iteration};
214 std::vector<GlobalVector*>
const& x,
215 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
217 if (!_compensate_non_equilibrium_initial_residuum)
222 _equation_system->assemble(x, x_prev, process_id);
224 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
228 std::vector<GlobalVector*>& x,
229 std::vector<GlobalVector*>
const& x_prev,
230 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
231 postIterationCallback,
232 int const process_id)
235 auto& sys = *_equation_system;
239 auto& minus_delta_x =
245 bool error_norms_met =
false;
251 _convergence_criterion->preFirstIteration();
254 for (; iteration <= _maxiter;
255 ++iteration, _convergence_criterion->reset())
258 double time_dirichlet = 0.0;
261 time_iteration.
start();
263 timer_dirichlet.
start();
264 sys.computeKnownSolutions(*x[process_id], process_id);
265 sys.applyKnownSolutions(*x[process_id]);
266 time_dirichlet += timer_dirichlet.
elapsed();
268 sys.preIteration(iteration, *x[process_id]);
271 time_assembly.
start();
274 sys.assemble(x, x_prev, process_id);
278 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
280 error_norms_met =
false;
281 iteration = _maxiter;
284 sys.getResidual(*x[process_id], *x_prev[process_id], res);
286 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
289 if (_r_neq !=
nullptr)
292 minus_delta_x.setZero();
294 timer_dirichlet.
start();
295 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
296 time_dirichlet += timer_dirichlet.
elapsed();
297 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
299 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
301 _convergence_criterion->checkResidual(res);
305 time_linear_solver.
start();
306 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
307 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
309 if (!iteration_succeeded)
311 ERR(
"Newton: The linear solver failed.");
320 std::vector<GlobalVector*> x_new{x};
323 *x[process_id], _x_new_id);
324 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
326 if (postIterationCallback)
328 postIterationCallback(iteration, x_new);
331 switch (sys.postIteration(*x_new[process_id]))
336 ERR(
"Newton: The postIteration() hook reported a "
337 "non-recoverable error.");
338 iteration_succeeded =
false;
342 "Newton: The postIteration() hook decided that this "
344 " has to be repeated.");
357 if (!iteration_succeeded)
360 error_norms_met =
false;
364 if (sys.isLinear()) {
365 error_norms_met =
true;
367 if (_convergence_criterion->hasDeltaXCheck()) {
369 _convergence_criterion->checkDeltaX(minus_delta_x,
373 error_norms_met = _convergence_criterion->isSatisfied();
376 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
386 if (iteration >= _maxiter)
392 if (iteration > _maxiter)
394 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
404 return {error_norms_met, iteration};
416 if (type ==
"Picard") {
419 return std::make_pair(
420 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
422 if (type ==
"Newton")
429 "The damping factor for the Newon method must be positive, got "
435 return std::make_pair(
436 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
440 if (boost::iequals(type,
"PETScSNES"))
447 return std::make_pair(std::make_unique<ConcreteNLS>(
448 linear_solver, max_iter, std::move(prefix)),
453 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());
458 if (_r_neq !=
nullptr)
466 if (_r_neq !=
nullptr)
void INFO(char const *fmt, Args const &... args)
void ERR(char const *fmt, Args const &... 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.
Global vector based on Eigen vector.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
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 axpy(PETScVector &y, double const a, PETScVector const &x)
void copy(PETScVector const &x, PETScVector &y)
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.