95 std::vector<GlobalVector*>
const& x,
96 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
98 if (!_compensate_non_equilibrium_initial_residuum)
103 INFO(
"Calculate non-equilibrium initial residuum.");
107 _equation_system->assemble(x, x_prev, process_id);
108 _equation_system->getA(A);
109 _equation_system->getRhs(*x_prev[process_id], rhs);
119 const auto selected_global_indices =
120 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
122 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
123 _r_neq->set(selected_global_indices, zero_entries);
132 std::vector<GlobalVector*>& x,
133 std::vector<GlobalVector*>
const& x_prev,
134 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
135 postIterationCallback,
136 int const process_id)
139 auto& sys = *_equation_system;
144 std::vector<GlobalVector*> x_new{x};
149 bool error_norms_met =
false;
151 _convergence_criterion->preFirstIteration();
154 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
157 double time_dirichlet = 0.0;
160 time_iteration.
start();
162 timer_dirichlet.
start();
163 auto& x_new_process = *x_new[process_id];
165 sys.computeKnownSolutions(x_new_process, process_id);
166 sys.applyKnownSolutions(x_new_process);
167 time_dirichlet += timer_dirichlet.
elapsed();
169 sys.preIteration(iteration, x_new_process);
172 time_assembly.
start();
173 sys.assemble(x_new, x_prev, process_id);
175 sys.getRhs(*x_prev[process_id], rhs);
176 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
179 if (_r_neq !=
nullptr)
184 timer_dirichlet.
start();
185 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
186 time_dirichlet += timer_dirichlet.
elapsed();
187 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
189 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
194 _convergence_criterion->checkResidual(res);
197 bool iteration_succeeded =
199 sys.linearSolverNeedsToCompute());
201 if (iteration_succeeded)
203 if (postIterationCallback)
205 postIterationCallback(iteration, x_new);
208 switch (sys.postIteration(x_new_process))
215 ERR(
"Picard: The postIteration() hook reported a "
216 "non-recoverable error.");
217 iteration_succeeded =
false;
225 "Picard: The postIteration() hook decided that this "
226 "iteration has to be repeated.");
234 if (!iteration_succeeded)
237 error_norms_met =
false;
243 error_norms_met =
true;
247 if (_convergence_criterion->hasDeltaXCheck())
252 _convergence_criterion->checkDeltaX(minus_delta_x,
256 error_norms_met = _convergence_criterion->isSatisfied();
262 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
272 if (iteration >= _maxiter)
278 if (iteration > _maxiter)
280 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
289 return {error_norms_met, iteration};
294 std::vector<GlobalVector*>
const& x,
295 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
297 if (!_compensate_non_equilibrium_initial_residuum)
302 INFO(
"Calculate non-equilibrium initial residuum.");
304 _equation_system->assemble(x, x_prev, process_id);
306 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
311 const auto selected_global_indices =
312 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
313 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
315 _r_neq->set(selected_global_indices, zero_entries);
321 std::vector<GlobalVector*>& x,
322 std::vector<GlobalVector*>
const& x_prev,
323 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
324 postIterationCallback,
325 int const process_id)
328 auto& sys = *_equation_system;
331 auto& minus_delta_x =
335 bool error_norms_met =
false;
341 _convergence_criterion->preFirstIteration();
344 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
347 double time_dirichlet = 0.0;
350 time_iteration.
start();
352 timer_dirichlet.
start();
353 sys.computeKnownSolutions(*x[process_id], process_id);
354 sys.applyKnownSolutions(*x[process_id]);
355 time_dirichlet += timer_dirichlet.
elapsed();
357 sys.preIteration(iteration, *x[process_id]);
360 time_assembly.
start();
363 sys.assemble(x, x_prev, process_id);
367 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
369 error_norms_met =
false;
370 iteration = _maxiter;
373 sys.getResidual(*x[process_id], *x_prev[process_id], res);
375 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
378 if (_r_neq !=
nullptr)
381 minus_delta_x.setZero();
383 timer_dirichlet.
start();
384 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
385 time_dirichlet += timer_dirichlet.
elapsed();
386 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
388 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
390 _convergence_criterion->checkResidual(res);
394 time_linear_solver.
start();
395 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
396 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
398 if (!iteration_succeeded)
400 ERR(
"Newton: The linear solver failed.");
409 std::vector<GlobalVector*> x_new{x};
412 *x[process_id], _x_new_id);
413 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
415 if (postIterationCallback)
417 postIterationCallback(iteration, x_new);
420 switch (sys.postIteration(*x_new[process_id]))
425 ERR(
"Newton: The postIteration() hook reported a "
426 "non-recoverable error.");
427 iteration_succeeded =
false;
431 "Newton: The postIteration() hook decided that this "
432 "iteration has to be repeated.");
445 if (!iteration_succeeded)
448 error_norms_met =
false;
454 error_norms_met =
true;
458 if (_convergence_criterion->hasDeltaXCheck())
461 _convergence_criterion->checkDeltaX(minus_delta_x,
465 error_norms_met = _convergence_criterion->isSatisfied();
468 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
478 if (iteration >= _maxiter)
484 if (iteration > _maxiter)
486 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
495 return {error_norms_met, iteration};
507 if (type ==
"Picard")
511 return std::make_pair(
512 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
514 if (type ==
"Newton")
521 "The damping factor for the Newon method must be positive, got "
527 return std::make_pair(
528 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
532 if (boost::iequals(type,
"PETScSNES"))
539 return std::make_pair(std::make_unique<ConcreteNLS>(
540 linear_solver, max_iter, std::move(prefix)),
545 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());