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 time_dirichlet += timer_dirichlet.
elapsed();
350 sys.preIteration(iteration, *x[process_id]);
353 time_assembly.
start();
356 sys.assemble(x, x_prev, process_id);
360 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
362 error_norms_met =
false;
363 iteration = _maxiter;
366 sys.getResidual(*x[process_id], *x_prev[process_id], res);
368 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
371 if (_r_neq !=
nullptr)
374 minus_delta_x.setZero();
376 timer_dirichlet.
start();
377 sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
378 time_dirichlet += timer_dirichlet.
elapsed();
379 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
381 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
383 _convergence_criterion->checkResidual(res);
387 time_linear_solver.
start();
388 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
389 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
391 if (!iteration_succeeded)
393 ERR(
"Newton: The linear solver failed.");
402 std::vector<GlobalVector*> x_new{x};
405 *x[process_id], _x_new_id);
406 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
408 if (postIterationCallback)
410 postIterationCallback(iteration, x_new);
413 switch (sys.postIteration(*x_new[process_id]))
418 ERR(
"Newton: The postIteration() hook reported a "
419 "non-recoverable error.");
420 iteration_succeeded =
false;
424 "Newton: The postIteration() hook decided that this "
425 "iteration has to be repeated.");
438 if (!iteration_succeeded)
441 error_norms_met =
false;
447 error_norms_met =
true;
451 if (_convergence_criterion->hasDeltaXCheck())
454 _convergence_criterion->checkDeltaX(minus_delta_x,
458 error_norms_met = _convergence_criterion->isSatisfied();
461 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
471 if (iteration >= _maxiter)
477 if (iteration > _maxiter)
479 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
488 return {error_norms_met, iteration};
500 if (type ==
"Picard")
504 return std::make_pair(
505 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
507 if (type ==
"Newton")
514 "The damping factor for the Newon method must be positive, got "
520 return std::make_pair(
521 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
525 if (boost::iequals(type,
"PETScSNES"))
532 return std::make_pair(std::make_unique<ConcreteNLS>(
533 linear_solver, max_iter, std::move(prefix)),
538 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());