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);
128 std::vector<GlobalVector*>& x,
129 std::vector<GlobalVector*>
const& x_prev,
130 std::function<
void(
int,
bool, std::vector<GlobalVector*>
const&)>
const&
131 postIterationCallback,
132 int const process_id)
135 auto& sys = *_equation_system;
140 std::vector<GlobalVector*> x_new{x};
145 bool error_norms_met =
false;
147 _convergence_criterion->preFirstIteration();
150 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
153 double time_dirichlet = 0.0;
156 time_iteration.
start();
158 INFO(
"Iteration #{:d} started.", iteration);
159 timer_dirichlet.
start();
160 auto& x_new_process = *x_new[process_id];
162 sys.computeKnownSolutions(x_new_process, process_id);
163 sys.applyKnownSolutions(x_new_process);
164 time_dirichlet += timer_dirichlet.
elapsed();
166 sys.preIteration(iteration, x_new_process);
169 time_assembly.
start();
170 sys.assemble(x_new, x_prev, process_id);
172 sys.getRhs(*x_prev[process_id], rhs);
175 if (sys.requiresNormalization() &&
176 !_linear_solver.canSolveRectangular())
178 sys.getAandRhsNormalized(A, rhs);
180 "The equation system is rectangular, but the current linear "
181 "solver only supports square systems. "
182 "The system will be normalized, which lead to a squared "
183 "condition number and potential numerical issues. "
184 "It is recommended to use a solver that supports rectangular "
185 "equation systems for better numerical stability.");
188 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
191 if (_r_neq !=
nullptr)
196 timer_dirichlet.
start();
197 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
198 time_dirichlet += timer_dirichlet.
elapsed();
199 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
201 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
206 _convergence_criterion->checkResidual(res);
209 bool iteration_succeeded =
211 sys.linearSolverNeedsToCompute());
213 if (iteration_succeeded)
215 if (postIterationCallback)
217 postIterationCallback(iteration, error_norms_met, x_new);
220 switch (sys.postIteration(x_new_process))
227 ERR(
"Picard: The postIteration() hook reported a "
228 "non-recoverable error.");
229 iteration_succeeded =
false;
237 "Picard: The postIteration() hook decided that this "
238 "iteration has to be repeated.");
246 if (!iteration_succeeded)
249 error_norms_met =
false;
255 error_norms_met =
true;
259 if (_convergence_criterion->hasDeltaXCheck())
264 _convergence_criterion->checkDeltaX(minus_delta_x,
268 error_norms_met = _convergence_criterion->isSatisfied();
274 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
284 if (iteration >= _maxiter)
290 if (iteration > _maxiter)
292 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
301 return {error_norms_met, iteration};
306 std::vector<GlobalVector*>
const& x,
307 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
309 if (!_compensate_non_equilibrium_initial_residuum)
314 INFO(
"Calculate non-equilibrium initial residuum.");
316 _equation_system->assemble(x, x_prev, process_id);
318 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
323 const auto selected_global_indices =
324 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
325 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
327 _r_neq->set(selected_global_indices, zero_entries);
333 std::vector<GlobalVector*>& x,
334 std::vector<GlobalVector*>
const& x_prev,
335 std::function<
void(
int,
bool, std::vector<GlobalVector*>
const&)>
const&
336 postIterationCallback,
337 int const process_id)
340 auto& sys = *_equation_system;
343 auto& minus_delta_x =
347 bool error_norms_met =
false;
353 _convergence_criterion->preFirstIteration();
356#if !defined(USE_PETSC) && !defined(USE_LIS)
357 int next_iteration_inv_jacobian_recompute = 1;
359 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
362 double time_dirichlet = 0.0;
365 INFO(
"Iteration #{:d} started.", iteration);
366 time_iteration.
start();
368 timer_dirichlet.
start();
369 sys.computeKnownSolutions(*x[process_id], process_id);
370 time_dirichlet += timer_dirichlet.
elapsed();
372 sys.preIteration(iteration, *x[process_id]);
375 time_assembly.
start();
376 bool mpi_rank_assembly_ok =
true;
379 sys.assemble(x, x_prev, process_id);
383 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
385 error_norms_met =
false;
386 iteration = _maxiter;
387 mpi_rank_assembly_ok =
false;
393 sys.getResidual(*x[process_id], *x_prev[process_id], res);
395 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
398 if (_r_neq !=
nullptr)
401 minus_delta_x.setZero();
403 timer_dirichlet.
start();
404 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
405 time_dirichlet += timer_dirichlet.
elapsed();
406 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
408 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
410 _convergence_criterion->checkResidual(res);
414 time_linear_solver.
start();
415#if !defined(USE_PETSC) && !defined(USE_LIS)
417 if (iteration == next_iteration_inv_jacobian_recompute)
419 linear_solver_behaviour =
421 next_iteration_inv_jacobian_recompute =
422 next_iteration_inv_jacobian_recompute + _recompute_jacobian;
424 if (!_linear_solver.compute(J, linear_solver_behaviour))
426 ERR(
"Newton: The linear solver failed in the compute() step.");
428 bool iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
430 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
432 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
434 if (!iteration_succeeded)
436 ERR(
"Newton: The linear solver failed.");
445 std::vector<GlobalVector*> x_new{x};
448 *x[process_id], _x_new_id);
449 double damping = _damping;
450 if (_convergence_criterion->hasNonNegativeDamping())
452 damping = _convergence_criterion->getDampingFactor(
453 minus_delta_x, *x[process_id], _damping);
455 if (_damping_reduction)
460 std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
462 LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
464 if (postIterationCallback)
466 postIterationCallback(iteration, error_norms_met, x_new);
469 switch (sys.postIteration(*x_new[process_id]))
474 ERR(
"Newton: The postIteration() hook reported a "
475 "non-recoverable error.");
476 iteration_succeeded =
false;
480 "Newton: The postIteration() hook decided that this "
481 "iteration has to be repeated.");
494 if (!iteration_succeeded)
497 error_norms_met =
false;
503 error_norms_met =
true;
507 if (_convergence_criterion->hasDeltaXCheck())
510 _convergence_criterion->checkDeltaX(minus_delta_x,
514 error_norms_met = _convergence_criterion->isSatisfied();
517 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
527 if (iteration >= _maxiter)
533 if (iteration > _maxiter)
535 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
544 return {error_norms_met, iteration};
556 if (type ==
"Picard")
560 return std::make_pair(
561 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
563 if (type ==
"Newton")
565 auto const recompute_jacobian =
574 "The damping factor for the Newon method must be positive, got "
578 auto const damping_reduction =
583 return std::make_pair(std::make_unique<ConcreteNLS>(
584 linear_solver, max_iter, recompute_jacobian,
585 damping, damping_reduction),
589 if (boost::iequals(type,
"PETScSNES"))
596 return std::make_pair(std::make_unique<ConcreteNLS>(
597 linear_solver, max_iter, std::move(prefix)),
602 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());