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);
173 if (sys.requiresNormalization() &&
174 !_linear_solver.canSolveRectangular())
176 sys.getAandRhsNormalized(A, rhs);
178 "The equation system is rectangular, but the current linear "
179 "solver only supports square systems. "
180 "The system will be normalized, which lead to a squared "
181 "condition number and potential numerical issues. "
182 "It is recommended to use a solver that supports rectangular "
183 "equation systems for better numerical stability.");
186 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
189 if (_r_neq !=
nullptr)
194 timer_dirichlet.
start();
195 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
196 time_dirichlet += timer_dirichlet.
elapsed();
197 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
199 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
204 _convergence_criterion->checkResidual(res);
207 bool iteration_succeeded =
209 sys.linearSolverNeedsToCompute());
211 if (iteration_succeeded)
213 if (postIterationCallback)
215 postIterationCallback(iteration, x_new);
218 switch (sys.postIteration(x_new_process))
225 ERR(
"Picard: The postIteration() hook reported a "
226 "non-recoverable error.");
227 iteration_succeeded =
false;
235 "Picard: The postIteration() hook decided that this "
236 "iteration has to be repeated.");
244 if (!iteration_succeeded)
247 error_norms_met =
false;
253 error_norms_met =
true;
257 if (_convergence_criterion->hasDeltaXCheck())
262 _convergence_criterion->checkDeltaX(minus_delta_x,
266 error_norms_met = _convergence_criterion->isSatisfied();
272 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
282 if (iteration >= _maxiter)
288 if (iteration > _maxiter)
290 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
299 return {error_norms_met, iteration};
304 std::vector<GlobalVector*>
const& x,
305 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
307 if (!_compensate_non_equilibrium_initial_residuum)
312 INFO(
"Calculate non-equilibrium initial residuum.");
314 _equation_system->assemble(x, x_prev, process_id);
316 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
321 const auto selected_global_indices =
322 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
323 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
325 _r_neq->set(selected_global_indices, zero_entries);
331 std::vector<GlobalVector*>& x,
332 std::vector<GlobalVector*>
const& x_prev,
333 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
334 postIterationCallback,
335 int const process_id)
338 auto& sys = *_equation_system;
341 auto& minus_delta_x =
345 bool error_norms_met =
false;
351 _convergence_criterion->preFirstIteration();
354 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
357 double time_dirichlet = 0.0;
360 time_iteration.
start();
362 timer_dirichlet.
start();
363 sys.computeKnownSolutions(*x[process_id], process_id);
364 time_dirichlet += timer_dirichlet.
elapsed();
366 sys.preIteration(iteration, *x[process_id]);
369 time_assembly.
start();
370 int mpi_rank_assembly_ok = 1;
371 int mpi_all_assembly_ok = 0;
374 sys.assemble(x, x_prev, process_id);
378 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
380 error_norms_met =
false;
381 iteration = _maxiter;
382 mpi_rank_assembly_ok = 0;
385 if (mpi_all_assembly_ok == 0)
389 sys.getResidual(*x[process_id], *x_prev[process_id], res);
391 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
394 if (_r_neq !=
nullptr)
397 minus_delta_x.setZero();
399 timer_dirichlet.
start();
400 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
401 time_dirichlet += timer_dirichlet.
elapsed();
402 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
404 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
406 _convergence_criterion->checkResidual(res);
410 time_linear_solver.
start();
411 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
412 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
414 if (!iteration_succeeded)
416 ERR(
"Newton: The linear solver failed.");
425 std::vector<GlobalVector*> x_new{x};
428 *x[process_id], _x_new_id);
429 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
431 if (postIterationCallback)
433 postIterationCallback(iteration, x_new);
436 switch (sys.postIteration(*x_new[process_id]))
441 ERR(
"Newton: The postIteration() hook reported a "
442 "non-recoverable error.");
443 iteration_succeeded =
false;
447 "Newton: The postIteration() hook decided that this "
448 "iteration has to be repeated.");
461 if (!iteration_succeeded)
464 error_norms_met =
false;
470 error_norms_met =
true;
474 if (_convergence_criterion->hasDeltaXCheck())
477 _convergence_criterion->checkDeltaX(minus_delta_x,
481 error_norms_met = _convergence_criterion->isSatisfied();
484 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
494 if (iteration >= _maxiter)
500 if (iteration > _maxiter)
502 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
511 return {error_norms_met, iteration};
523 if (type ==
"Picard")
527 return std::make_pair(
528 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
530 if (type ==
"Newton")
537 "The damping factor for the Newon method must be positive, got "
543 return std::make_pair(
544 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
548 if (boost::iequals(type,
"PETScSNES"))
555 return std::make_pair(std::make_unique<ConcreteNLS>(
556 linear_solver, max_iter, std::move(prefix)),
561 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());