13 #include <boost/algorithm/string.hpp>
29 std::vector<GlobalVector*>
const& x,
30 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
32 if (!_compensate_non_equilibrium_initial_residuum)
39 _equation_system->assemble(x, x_prev, process_id);
40 _equation_system->getA(A);
41 _equation_system->getRhs(*x_prev[process_id], rhs);
50 std::vector<GlobalVector*>& x,
51 std::vector<GlobalVector*>
const& x_prev,
52 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
53 postIterationCallback,
57 auto& sys = *_equation_system;
62 std::vector<GlobalVector*> x_new{x};
67 bool error_norms_met =
false;
69 _convergence_criterion->preFirstIteration();
72 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
75 double time_dirichlet = 0.0;
78 time_iteration.
start();
80 timer_dirichlet.
start();
81 sys.computeKnownSolutions(*x_new[process_id], process_id);
82 sys.applyKnownSolutions(*x_new[process_id]);
83 time_dirichlet += timer_dirichlet.
elapsed();
85 sys.preIteration(iteration, *x_new[process_id]);
88 time_assembly.
start();
89 sys.assemble(x_new, x_prev, process_id);
91 sys.getRhs(*x_prev[process_id], rhs);
92 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
95 if (_r_neq !=
nullptr)
100 timer_dirichlet.
start();
101 sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
102 time_dirichlet += timer_dirichlet.
elapsed();
103 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
105 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
110 _convergence_criterion->checkResidual(res);
114 time_linear_solver.
start();
115 bool iteration_succeeded =
116 _linear_solver.solve(A, rhs, *x_new[process_id]);
117 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
119 if (!iteration_succeeded)
121 ERR(
"Picard: The linear solver failed.");
125 if (postIterationCallback)
127 postIterationCallback(iteration, x_new);
130 switch (sys.postIteration(*x_new[process_id]))
137 ERR(
"Picard: The postIteration() hook reported a "
138 "non-recoverable error.");
139 iteration_succeeded =
false;
147 "Picard: The postIteration() hook decided that this "
148 "iteration has to be repeated.");
156 if (!iteration_succeeded)
159 error_norms_met =
false;
165 error_norms_met =
true;
169 if (_convergence_criterion->hasDeltaXCheck())
174 _convergence_criterion->checkDeltaX(minus_delta_x,
178 error_norms_met = _convergence_criterion->isSatisfied();
184 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
194 if (iteration >= _maxiter)
200 if (iteration > _maxiter)
202 ERR(
"Picard: Could not solve the given nonlinear system within {:d} "
211 return {error_norms_met, iteration};
216 std::vector<GlobalVector*>
const& x,
217 std::vector<GlobalVector*>
const& x_prev,
int const process_id)
219 if (!_compensate_non_equilibrium_initial_residuum)
224 _equation_system->assemble(x, x_prev, process_id);
226 _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
230 std::vector<GlobalVector*>& x,
231 std::vector<GlobalVector*>
const& x_prev,
232 std::function<
void(
int, std::vector<GlobalVector*>
const&)>
const&
233 postIterationCallback,
234 int const process_id)
237 auto& sys = *_equation_system;
240 auto& minus_delta_x =
244 bool error_norms_met =
false;
250 _convergence_criterion->preFirstIteration();
253 for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
256 double time_dirichlet = 0.0;
259 time_iteration.
start();
261 timer_dirichlet.
start();
262 sys.computeKnownSolutions(*x[process_id], process_id);
263 sys.applyKnownSolutions(*x[process_id]);
264 time_dirichlet += timer_dirichlet.
elapsed();
266 sys.preIteration(iteration, *x[process_id]);
269 time_assembly.
start();
272 sys.assemble(x, x_prev, process_id);
276 ERR(
"Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
278 error_norms_met =
false;
279 iteration = _maxiter;
282 sys.getResidual(*x[process_id], *x_prev[process_id], res);
284 INFO(
"[time] Assembly took {:g} s.", time_assembly.
elapsed());
287 if (_r_neq !=
nullptr)
290 minus_delta_x.setZero();
292 timer_dirichlet.
start();
293 sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
294 time_dirichlet += timer_dirichlet.
elapsed();
295 INFO(
"[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
297 if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
299 _convergence_criterion->checkResidual(res);
303 time_linear_solver.
start();
304 bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
305 INFO(
"[time] Linear solver took {:g} s.", time_linear_solver.
elapsed());
307 if (!iteration_succeeded)
309 ERR(
"Newton: The linear solver failed.");
318 std::vector<GlobalVector*> x_new{x};
321 *x[process_id], _x_new_id);
322 LinAlg::axpy(*x_new[process_id], -_damping, minus_delta_x);
324 if (postIterationCallback)
326 postIterationCallback(iteration, x_new);
329 switch (sys.postIteration(*x_new[process_id]))
334 ERR(
"Newton: The postIteration() hook reported a "
335 "non-recoverable error.");
336 iteration_succeeded =
false;
340 "Newton: The postIteration() hook decided that this "
341 "iteration has to be repeated.");
354 if (!iteration_succeeded)
357 error_norms_met =
false;
363 error_norms_met =
true;
367 if (_convergence_criterion->hasDeltaXCheck())
370 _convergence_criterion->checkDeltaX(minus_delta_x,
374 error_norms_met = _convergence_criterion->isSatisfied();
377 INFO(
"[time] Iteration #{:d} took {:g} s.", iteration,
387 if (iteration >= _maxiter)
393 if (iteration > _maxiter)
395 ERR(
"Newton: Could not solve the given nonlinear system within {:d} "
404 return {error_norms_met, iteration};
416 if (type ==
"Picard")
420 return std::make_pair(
421 std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
423 if (type ==
"Newton")
430 "The damping factor for the Newon method must be positive, got "
436 return std::make_pair(
437 std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
441 if (boost::iequals(type,
"PETScSNES"))
448 return std::make_pair(std::make_unique<ConcreteNLS>(
449 linear_solver, max_iter, std::move(prefix)),
454 OGS_FATAL(
"Unsupported nonlinear solver type '{:s}'.", type.c_str());
459 if (_r_neq !=
nullptr)
467 if (_r_neq !=
nullptr)
void INFO(char const *fmt, Args const &... args)
void ERR(char const *fmt, Args const &... args)
Definition of the RunTime class.
T getConfigParameter(std::string const ¶m) const
double elapsed() const
Get the elapsed time in seconds.
void start()
Start the timer.
Global vector based on Eigen vector.
virtual GlobalMatrix & getMatrix(std::size_t &id)=0
Get an uninitialized matrix with the given id.
virtual void releaseMatrix(GlobalMatrix const &A)=0
virtual GlobalVector & getVector(std::size_t &id)=0
Get an uninitialized vector with the given id.
virtual void releaseVector(GlobalVector const &x)=0
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
std::pair< std::unique_ptr< NonlinearSolverBase >, NonlinearSolverTag > createNonlinearSolver(GlobalLinearSolver &linear_solver, BaseLib::ConfigTree const &config)
void axpy(PETScVector &y, double const a, PETScVector const &x)
void copy(PETScVector const &x, PETScVector &y)
void matMult(PETScMatrix const &A, PETScVector const &x, PETScVector &y)
void calculateNonEquilibriumInitialResiduum(std::vector< std::unique_ptr< ProcessData >> const &per_process_data, std::vector< GlobalVector * > process_solutions, std::vector< GlobalVector * > const &process_solutions_prev)
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Status of the non-linear solver.