89 std::vector<GlobalVector*>
const& x,
90 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
92 if (!_compensate_non_equilibrium_initial_residuum)
97 INFO(
"Calculate non-equilibrium initial residuum.");
101 _equation_system->assemble(x, x_prev, process_id);
102 _equation_system->getA(A);
103 _equation_system->getRhs(*x_prev[process_id], rhs);
113 const auto selected_global_indices =
114 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
116 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
117 _r_neq->set(selected_global_indices, zero_entries);
126 std::vector<GlobalVector*>& x,
127 std::vector<GlobalVector*>
const& x_prev,
128 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
129 postIterationCallback,
130 int const process_id)
133 auto& sys = *_equation_system;
138 std::vector<GlobalVector*> x_new{x};
143 bool error_norms_met =
false;
145 _convergence_criterion->preFirstIteration();
148 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
151 double time_dirichlet = 0.0;
154 time_iteration.
start();
156 timer_dirichlet.
start();
157 auto& x_new_process = *x_new[process_id];
159 sys.computeKnownSolutions(x_new_process, process_id);
160 sys.applyKnownSolutions(x_new_process);
161 time_dirichlet += timer_dirichlet.
elapsed();
163 sys.preIteration(iteration, x_new_process);
166 time_assembly.
start();
167 sys.assemble(x_new, x_prev, process_id);
169 sys.getRhs(*x_prev[process_id], rhs);
170 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
173 if (_r_neq !=
nullptr)
178 timer_dirichlet.
start();
179 sys.applyKnownSolutionsPicard(A, rhs, x_new_process);
180 time_dirichlet += timer_dirichlet.
elapsed();
181 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
183 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
188 _convergence_criterion->checkResidual(res);
191 bool iteration_succeeded =
193 sys.linearSolverNeedsToCompute());
195 if (iteration_succeeded)
197 if (postIterationCallback)
199 postIterationCallback(iteration, x_new);
202 switch (sys.postIteration(x_new_process))
209 ERR(
"Picard: The postIteration() hook reported a "
210 "non-recoverable error.");
211 iteration_succeeded =
false;
219 "Picard: The postIteration() hook decided that this "
220 "iteration has to be repeated.");
228 if (!iteration_succeeded)
231 error_norms_met =
false;
237 error_norms_met =
true;
241 if (_convergence_criterion->hasDeltaXCheck())
246 _convergence_criterion->checkDeltaX(minus_delta_x,
250 error_norms_met = _convergence_criterion->isSatisfied();
256 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
266 if (iteration >= _maxiter)
272 if (iteration > _maxiter)
274 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
283 return {error_norms_met, iteration};
288 std::vector<GlobalVector*>
const& x,
289 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
291 if (!_compensate_non_equilibrium_initial_residuum)
296 INFO(
"Calculate non-equilibrium initial residuum.");
298 _equation_system->assemble(x, x_prev, process_id);
300 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
305 const auto selected_global_indices =
306 _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
307 std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
309 _r_neq->set(selected_global_indices, zero_entries);
315 std::vector<GlobalVector*>& x,
316 std::vector<GlobalVector*>
const& x_prev,
317 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
318 postIterationCallback,
319 int const process_id)
322 auto& sys = *_equation_system;
325 auto& minus_delta_x =
329 bool error_norms_met =
false;
335 _convergence_criterion->preFirstIteration();
338 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
341 double time_dirichlet = 0.0;
344 time_iteration.
start();
346 timer_dirichlet.
start();
347 sys.computeKnownSolutions(*x[process_id], process_id);
348 sys.applyKnownSolutions(*x[process_id]);
349 time_dirichlet += timer_dirichlet.
elapsed();
351 sys.preIteration(iteration, *x[process_id]);
354 time_assembly.
start();
357 sys.assemble(x, x_prev, process_id);
361 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
363 error_norms_met =
false;
364 iteration = _maxiter;
367 sys.getResidual(*x[process_id], *x_prev[process_id], res);
369 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
372 if (_r_neq !=
nullptr)
375 minus_delta_x.setZero();
377 timer_dirichlet.
start();
378 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
379 time_dirichlet += timer_dirichlet.
elapsed();
380 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
382 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
384 _convergence_criterion->checkResidual(res);
388 time_linear_solver.
start();
389 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
390 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
392 if (!iteration_succeeded)
394 ERR(
"Newton: The linear solver failed.");
403 std::vector<GlobalVector*> x_new{x};
406 *x[process_id], _x_new_id);
407 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
409 if (postIterationCallback)
411 postIterationCallback(iteration, x_new);
414 switch (sys.postIteration(*x_new[process_id]))
419 ERR(
"Newton: The postIteration() hook reported a "
420 "non-recoverable error.");
421 iteration_succeeded =
false;
425 "Newton: The postIteration() hook decided that this "
426 "iteration has to be repeated.");
439 if (!iteration_succeeded)
442 error_norms_met =
false;
448 error_norms_met =
true;
452 if (_convergence_criterion->hasDeltaXCheck())
455 _convergence_criterion->checkDeltaX(minus_delta_x,
459 error_norms_met = _convergence_criterion->isSatisfied();
462 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
472 if (iteration >= _maxiter)
478 if (iteration > _maxiter)
480 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
489 return {error_norms_met, iteration};
501 if (type ==
"Picard")
505 return std::make_pair(
506 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
508 if (type ==
"Newton")
515 "The damping factor for the Newon method must be positive, got "
521 return std::make_pair(
522 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
526 if (boost::iequals(type,
"PETScSNES"))
533 return std::make_pair(std::make_unique<ConcreteNLS>(
534 linear_solver, max_iter, std::move(prefix)),
539 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());