OGS 6.2.1-499-g3b941532c.dirty.20191012113459
NonlinearSolver.cpp
Go to the documentation of this file.
1 
11 #include "NonlinearSolver.h"
12 
13 #include <logog/include/logog.hpp>
14 
15 #include "BaseLib/ConfigTree.h"
16 #include "BaseLib/Error.h"
17 #include "BaseLib/RunTime.h"
18 #include "MathLib/LinAlg/LinAlg.h"
20 #include "ConvergenceCriterion.h"
21 
22 namespace NumLib
23 {
25  std::vector<GlobalVector*> const& x, int const process_id) const
26 {
27  _equation_system->assemble(x, process_id);
28 }
29 
31  std::vector<GlobalVector*>& x,
32  std::function<void(int, GlobalVector const&)> const& postIterationCallback,
33  int const process_id)
34 {
35  namespace LinAlg = MathLib::LinAlg;
36  auto& sys = *_equation_system;
37 
38  auto& A =
41  _rhs_id);
42 
43  std::vector<GlobalVector*> x_new{x};
44  x_new[process_id] =
46  LinAlg::copy(*x[process_id], *x_new[process_id]); // set initial guess
47 
48  bool error_norms_met = false;
49 
50  _convergence_criterion->preFirstIteration();
51 
52  int iteration = 1;
53  for (; iteration <= _maxiter;
54  ++iteration, _convergence_criterion->reset())
55  {
56  BaseLib::RunTime timer_dirichlet;
57  double time_dirichlet = 0.0;
58 
59  BaseLib::RunTime time_iteration;
60  time_iteration.start();
61 
62  timer_dirichlet.start();
63  sys.computeKnownSolutions(*x_new[process_id], process_id);
64  sys.applyKnownSolutions(*x_new[process_id]);
65  time_dirichlet += timer_dirichlet.elapsed();
66 
67  sys.preIteration(iteration, *x_new[process_id]);
68 
69  BaseLib::RunTime time_assembly;
70  time_assembly.start();
71  sys.assemble(x_new, process_id);
72  sys.getA(A);
73  sys.getRhs(rhs);
74  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
75 
76  timer_dirichlet.start();
77  sys.applyKnownSolutionsPicard(A, rhs, *x_new[process_id]);
78  time_dirichlet += timer_dirichlet.elapsed();
79  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
80 
81  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck()) {
82  GlobalVector res;
83  LinAlg::matMult(A, *x_new[process_id], res); // res = A * x_new
84  LinAlg::axpy(res, -1.0, rhs); // res -= rhs
85  _convergence_criterion->checkResidual(res);
86  }
87 
88  BaseLib::RunTime time_linear_solver;
89  time_linear_solver.start();
90  bool iteration_succeeded =
91  _linear_solver.solve(A, rhs, *x_new[process_id]);
92  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
93 
94  if (!iteration_succeeded)
95  {
96  ERR("Picard: The linear solver failed.");
97  }
98  else
99  {
100  if (postIterationCallback)
101  {
102  postIterationCallback(iteration, *x_new[process_id]);
103  }
104 
105  switch (sys.postIteration(*x_new[process_id]))
106  {
108  // Don't copy here. The old x might still be used further
109  // below. Although currently it is not.
110  break;
112  ERR("Picard: The postIteration() hook reported a "
113  "non-recoverable error.");
114  iteration_succeeded = false;
115  // Copy new solution to x.
116  // Thereby the failed solution can be used by the caller for
117  // debugging purposes.
118  LinAlg::copy(*x_new[process_id], *x[process_id]);
119  break;
121  INFO(
122  "Picard: The postIteration() hook decided that this "
123  "iteration has to be repeated.");
124  LinAlg::copy(
125  *x[process_id],
126  *x_new[process_id]); // throw the iteration result away
127  continue;
128  }
129  }
130 
131  if (!iteration_succeeded)
132  {
133  // Don't compute error norms, break here.
134  error_norms_met = false;
135  break;
136  }
137 
138  if (sys.isLinear()) {
139  error_norms_met = true;
140  } else {
141  if (_convergence_criterion->hasDeltaXCheck()) {
142  GlobalVector minus_delta_x(*x[process_id]);
143  LinAlg::axpy(minus_delta_x, -1.0,
144  *x_new[process_id]); // minus_delta_x = x - x_new
145  _convergence_criterion->checkDeltaX(minus_delta_x,
146  *x_new[process_id]);
147  }
148 
149  error_norms_met = _convergence_criterion->isSatisfied();
150  }
151 
152  // Update x s.t. in the next iteration we will compute the right delta x
153  LinAlg::copy(*x_new[process_id], *x[process_id]);
154 
155  INFO("[time] Iteration #%u took %g s.", iteration,
156  time_iteration.elapsed());
157 
158  if (error_norms_met)
159  {
160  break;
161  }
162 
163  // Avoid increment of the 'iteration' if the error norms are not met,
164  // but maximum number of iterations is reached.
165  if (iteration >= _maxiter)
166  {
167  break;
168  }
169  }
170 
171  if (iteration > _maxiter)
172  {
173  ERR("Picard: Could not solve the given nonlinear system within %u "
174  "iterations",
175  _maxiter);
176  }
177 
181 
182  return {error_norms_met, iteration};
183 }
184 
186  std::vector<GlobalVector*> const& x, int const process_id) const
187 {
188  _equation_system->assemble(x, process_id);
189  // TODO if the equation system would be reset to nullptr after each
190  // assemble() or solve() call, the user would be forced to set the
191  // equation every time and could not forget it.
192 }
193 
195  std::vector<GlobalVector*>& x,
196  std::function<void(int, GlobalVector const&)> const& postIterationCallback,
197  int const process_id)
198 {
199  namespace LinAlg = MathLib::LinAlg;
200  auto& sys = *_equation_system;
201 
203  _res_id);
204  auto& minus_delta_x =
206  _minus_delta_x_id);
207  auto& J =
209 
210  bool error_norms_met = false;
211 
212  // TODO be more efficient
213  // init minus_delta_x to the right size
214  LinAlg::copy(*x[process_id], minus_delta_x);
215 
216  _convergence_criterion->preFirstIteration();
217 
218  int iteration = 1;
219  for (; iteration <= _maxiter;
220  ++iteration, _convergence_criterion->reset())
221  {
222  BaseLib::RunTime timer_dirichlet;
223  double time_dirichlet = 0.0;
224 
225  BaseLib::RunTime time_iteration;
226  time_iteration.start();
227 
228  timer_dirichlet.start();
229  sys.computeKnownSolutions(*x[process_id], process_id);
230  sys.applyKnownSolutions(*x[process_id]);
231  time_dirichlet += timer_dirichlet.elapsed();
232 
233  sys.preIteration(iteration, *x[process_id]);
234 
235  BaseLib::RunTime time_assembly;
236  time_assembly.start();
237  sys.assemble(x, process_id);
238  sys.getResidual(*x[process_id], res);
239  sys.getJacobian(J);
240  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
241 
242  minus_delta_x.setZero();
243 
244  timer_dirichlet.start();
245  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
246  time_dirichlet += timer_dirichlet.elapsed();
247  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
248 
249  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
250  {
251  _convergence_criterion->checkResidual(res);
252  }
253 
254  BaseLib::RunTime time_linear_solver;
255  time_linear_solver.start();
256  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
257  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
258 
259  if (!iteration_succeeded)
260  {
261  ERR("Newton: The linear solver failed.");
262  }
263  else
264  {
265  // TODO could be solved in a better way
266  // cf.
267  // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
269  *x[process_id], _x_new_id);
270  LinAlg::axpy(x_new, -_damping, minus_delta_x);
271 
272  if (postIterationCallback)
273  {
274  postIterationCallback(iteration, x_new);
275  }
276 
277  switch(sys.postIteration(x_new))
278  {
280  break;
282  ERR("Newton: The postIteration() hook reported a "
283  "non-recoverable error.");
284  iteration_succeeded = false;
285  break;
287  INFO(
288  "Newton: The postIteration() hook decided that this "
289  "iteration"
290  " has to be repeated.");
291  // TODO introduce some onDestroy hook.
293  continue; // That throws the iteration result away.
294  }
295 
296  LinAlg::copy(x_new, *x[process_id]); // copy new solution to x
298  }
299 
300  if (!iteration_succeeded)
301  {
302  // Don't compute further error norms, but break here.
303  error_norms_met = false;
304  break;
305  }
306 
307  if (sys.isLinear()) {
308  error_norms_met = true;
309  } else {
310  if (_convergence_criterion->hasDeltaXCheck()) {
311  // Note: x contains the new solution!
312  _convergence_criterion->checkDeltaX(minus_delta_x,
313  *x[process_id]);
314  }
315 
316  error_norms_met = _convergence_criterion->isSatisfied();
317  }
318 
319  INFO("[time] Iteration #%u took %g s.", iteration,
320  time_iteration.elapsed());
321 
322  if (error_norms_met)
323  {
324  break;
325  }
326 
327  // Avoid increment of the 'iteration' if the error norms are not met,
328  // but maximum number of iterations is reached.
329  if (iteration >= _maxiter)
330  {
331  break;
332  }
333  }
334 
335  if (iteration > _maxiter)
336  {
337  ERR("Newton: Could not solve the given nonlinear system within %u "
338  "iterations",
339  _maxiter);
340  }
341 
345  minus_delta_x);
346 
347  return {error_norms_met, iteration};
348 }
349 
350 std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
351 createNonlinearSolver(GlobalLinearSolver& linear_solver,
352  BaseLib::ConfigTree const& config)
353 {
355  auto const type = config.getConfigParameter<std::string>("type");
357  auto const max_iter = config.getConfigParameter<int>("max_iter");
358 
359  if (type == "Picard") {
360  auto const tag = NonlinearSolverTag::Picard;
361  using ConcreteNLS = NonlinearSolver<tag>;
362  return std::make_pair(
363  std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
364  }
365  if (type == "Newton")
366  {
368  auto const damping = config.getConfigParameter<double>("damping", 1.0);
369  if (damping <= 0)
370  {
371  OGS_FATAL(
372  "The damping factor for the Newon method must be positive, got "
373  "%g.",
374  damping);
375  }
376  auto const tag = NonlinearSolverTag::Newton;
377  using ConcreteNLS = NonlinearSolver<tag>;
378  return std::make_pair(
379  std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
380  tag);
381  }
382  OGS_FATAL("Unsupported nonlinear solver type");
383 }
384 } // namespace NumLib
void matMult(Matrix const &A, Vector const &x, Vector &y)
Definition: LinAlg.h:114
virtual void releaseMatrix(GlobalMatrix const &A)=0
std::pair< std::unique_ptr< NonlinearSolverBase >, NonlinearSolverTag > createNonlinearSolver(GlobalLinearSolver &linear_solver, BaseLib::ConfigTree const &config)
virtual GlobalVector & getVector()=0
Get an uninitialized vector.
T getConfigParameter(std::string const &param) const
static NUMLIB_EXPORT MatrixProvider & provider
static NUMLIB_EXPORT VectorProvider & provider
Definition of the RunTime class.
NonlinearSolverTag
Tag used to specify which nonlinear solver will be used.
Definition: Types.h:20
Status of the non-linear solver.
virtual GlobalMatrix & getMatrix()=0
Get an uninitialized matrix.
Count the running time.
Definition: RunTime.h:31
void start()
Start the timer.
Definition: RunTime.h:35
#define OGS_FATAL(fmt,...)
Definition: Error.h:64
virtual void releaseVector(GlobalVector const &x)=0
void axpy(MatrixOrVector &y, double const a, MatrixOrVector const &x)
Computes .
Definition: LinAlg.h:58
void copy(MatrixOrVector const &x, MatrixOrVector &y)
Copies x to y.
Definition: LinAlg.h:37
double elapsed() const
Get the elapsed time after started.
Definition: RunTime.h:51