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