91 std::vector<GlobalVector*>
const& x,
92 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
94 if (!_compensate_non_equilibrium_initial_residuum)
99 INFO(
"Calculate non-equilibrium initial residuum.");
103 _equation_system->assemble(x, x_prev, process_id);
104 _equation_system->getA(A);
105 _equation_system->getRhs(*x_prev[process_id], rhs);
115 const auto selected_global_indices =
116 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
118 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
119 _r_neq->set(selected_global_indices, zero_entries);
120 _equation_system->setReleaseNodalForces(_r_neq, process_id);
129 std::vector<GlobalVector*>& x,
130 std::vector<GlobalVector*>
const& x_prev,
131 std::function<
void(
int,
bool, std::vector<GlobalVector*>
const&)>
const&
132 postIterationCallback,
133 int const process_id)
136 auto& sys = *_equation_system;
141 std::vector<GlobalVector*> x_new{x};
146 bool error_norms_met =
false;
148 _convergence_criterion->preFirstIteration();
151 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
154 double time_dirichlet = 0.0;
157 time_iteration.
start();
159 INFO(
"Iteration #{:d} started.", iteration);
160 timer_dirichlet.
start();
161 auto& x_new_process = *x_new[process_id];
163 sys.computeKnownSolutions(x_new_process, process_id);
164 sys.applyKnownSolutions(x_new_process);
165 time_dirichlet += timer_dirichlet.
elapsed();
167 sys.preIteration(iteration, x_new_process);
170 time_assembly.
start();
171 sys.assemble(x_new, x_prev, process_id);
173 sys.getRhs(*x_prev[process_id], rhs);
176 if (sys.requiresNormalization() &&
177 !_linear_solver.canSolveRectangular())
179 sys.getAandRhsNormalized(A, rhs);
181 "The equation system is rectangular, but the current linear "
182 "solver only supports square systems. "
183 "The system will be normalized, which lead to a squared "
184 "condition number and potential numerical issues. "
185 "It is recommended to use a solver that supports rectangular "
186 "equation systems for better numerical stability.");
189 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
192 if (_r_neq !=
nullptr)
197 timer_dirichlet.
start();
198 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
199 time_dirichlet += timer_dirichlet.
elapsed();
200 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
202 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
207 _convergence_criterion->checkResidual(res);
210 bool iteration_succeeded =
212 sys.linearSolverNeedsToCompute());
214 if (iteration_succeeded)
216 if (postIterationCallback)
218 postIterationCallback(iteration, error_norms_met, x_new);
221 switch (sys.postIteration(x_new_process))
228 ERR(
"Picard: The postIteration() hook reported a "
229 "non-recoverable error.");
230 iteration_succeeded =
false;
238 "Picard: The postIteration() hook decided that this "
239 "iteration has to be repeated.");
247 if (!iteration_succeeded)
250 error_norms_met =
false;
256 error_norms_met =
true;
260 if (_convergence_criterion->hasDeltaXCheck())
265 _convergence_criterion->checkDeltaX(minus_delta_x,
269 error_norms_met = _convergence_criterion->isSatisfied();
275 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
285 if (iteration >= _maxiter)
291 if (iteration > _maxiter)
293 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
302 return {error_norms_met, iteration};
307 std::vector<GlobalVector*>
const& x,
308 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
310 if (!_compensate_non_equilibrium_initial_residuum)
315 INFO(
"Calculate non-equilibrium initial residuum.");
317 _equation_system->assemble(x, x_prev, process_id);
319 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
324 const auto selected_global_indices =
325 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
326 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
328 _r_neq->set(selected_global_indices, zero_entries);
329 _equation_system->setReleaseNodalForces(_r_neq, process_id);
335 std::vector<GlobalVector*>& x,
336 std::vector<GlobalVector*>
const& x_prev,
337 std::function<
void(
int,
bool, std::vector<GlobalVector*>
const&)>
const&
338 postIterationCallback,
339 int const process_id)
342 auto& sys = *_equation_system;
345 auto& minus_delta_x =
349 bool error_norms_met =
false;
355 _convergence_criterion->preFirstIteration();
358#if !defined(USE_PETSC) && !defined(USE_LIS)
359 int next_iteration_inv_jacobian_recompute = 1;
361 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
364 double time_dirichlet = 0.0;
367 INFO(
"Iteration #{:d} started.", iteration);
368 time_iteration.
start();
370 timer_dirichlet.
start();
371 sys.computeKnownSolutions(*x[process_id], process_id);
372 time_dirichlet += timer_dirichlet.
elapsed();
374 sys.preIteration(iteration, *x[process_id]);
377 time_assembly.
start();
378 bool mpi_rank_assembly_ok =
true;
381 sys.assemble(x, x_prev, process_id);
385 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
387 error_norms_met =
false;
388 iteration = _maxiter;
389 mpi_rank_assembly_ok =
false;
395 sys.getResidual(*x[process_id], *x_prev[process_id], res);
397 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
400 if (_r_neq !=
nullptr)
403 minus_delta_x.setZero();
405 timer_dirichlet.
start();
406 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
407 time_dirichlet += timer_dirichlet.
elapsed();
408 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
410 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
412 _convergence_criterion->checkResidual(res);
416 time_linear_solver.
start();
417#if !defined(USE_PETSC) && !defined(USE_LIS)
419 if (iteration == next_iteration_inv_jacobian_recompute)
421 linear_solver_behaviour =
423 next_iteration_inv_jacobian_recompute =
424 next_iteration_inv_jacobian_recompute + _recompute_jacobian;
427 bool iteration_succeeded =
false;
428 if (!_linear_solver.compute(J, linear_solver_behaviour))
430 ERR(
"Newton: The linear solver failed in the compute() step.");
434 iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
437 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
439 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
441 if (!iteration_succeeded)
443 ERR(
"Newton: The linear solver failed.");
452 std::vector<GlobalVector*> x_new{x};
455 *x[process_id], _x_new_id);
456 double damping = _damping;
457 if (_convergence_criterion->hasNonNegativeDamping())
459 damping = _convergence_criterion->getDampingFactor(
460 minus_delta_x, *x[process_id], _damping);
462 if (_damping_reduction)
467 std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
469 LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
471 if (postIterationCallback)
473 postIterationCallback(iteration, error_norms_met, x_new);
476 switch (sys.postIteration(*x_new[process_id]))
481 ERR(
"Newton: The postIteration() hook reported a "
482 "non-recoverable error.");
483 iteration_succeeded =
false;
487 "Newton: The postIteration() hook decided that this "
488 "iteration has to be repeated.");
501 if (!iteration_succeeded)
504 error_norms_met =
false;
510 error_norms_met =
true;
514 if (_convergence_criterion->hasDeltaXCheck())
517 _convergence_criterion->checkDeltaX(minus_delta_x,
521 error_norms_met = _convergence_criterion->isSatisfied();
524 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
534 if (iteration >= _maxiter)
540 if (iteration > _maxiter)
542 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
551 return {error_norms_met, iteration};
563 if (type ==
"Picard")
567 return std::make_pair(
568 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
570 if (type ==
"Newton")
572 auto const recompute_jacobian =
581 "The damping factor for the Newon method must be positive, got "
585 auto const damping_reduction =
590 return std::make_pair(std::make_unique<ConcreteNLS>(
591 linear_solver, max_iter, recompute_jacobian,
592 damping, damping_reduction),
596 if (boost::iequals(type,
"PETScSNES"))
603 return std::make_pair(std::make_unique<ConcreteNLS>(
604 linear_solver, max_iter, std::move(prefix)),
609 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());