OGS 6.1.0-1721-g6382411ad
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  postIterationCallback(iteration, x_new);
98 
99  switch (sys.postIteration(x_new))
100  {
102  // Don't copy here. The old x might still be used further
103  // below. Although currently it is not.
104  break;
106  ERR("Picard: The postIteration() hook reported a "
107  "non-recoverable error.");
108  iteration_succeeded = false;
109  // Copy new solution to x.
110  // Thereby the failed solution can be used by the caller for
111  // debugging purposes.
112  LinAlg::copy(x_new, x);
113  break;
115  INFO(
116  "Picard: The postIteration() hook decided that this "
117  "iteration has to be repeated.");
118  LinAlg::copy(x, x_new); // throw the iteration result away
119  continue;
120  }
121  }
122 
123  if (!iteration_succeeded)
124  {
125  // Don't compute error norms, break here.
126  error_norms_met = false;
127  break;
128  }
129 
130  if (sys.isLinear()) {
131  error_norms_met = true;
132  } else {
133  if (_convergence_criterion->hasDeltaXCheck()) {
134  GlobalVector minus_delta_x(x);
135  LinAlg::axpy(minus_delta_x, -1.0,
136  x_new); // minus_delta_x = x - x_new
137  _convergence_criterion->checkDeltaX(minus_delta_x, x_new);
138  }
139 
140  error_norms_met = _convergence_criterion->isSatisfied();
141  }
142 
143  // Update x s.t. in the next iteration we will compute the right delta x
144  LinAlg::copy(x_new, x);
145 
146  INFO("[time] Iteration #%u took %g s.", iteration,
147  time_iteration.elapsed());
148 
149  if (error_norms_met)
150  break;
151 
152  // Avoid increment of the 'iteration' if the error norms are not met,
153  // but maximum number of iterations is reached.
154  if (iteration >= _maxiter)
155  break;
156  }
157 
158  if (iteration > _maxiter)
159  {
160  ERR("Picard: Could not solve the given nonlinear system within %u "
161  "iterations",
162  _maxiter);
163  }
164 
168 
169  return {error_norms_met, iteration};
170 }
171 
173  GlobalVector const& x) const
174 {
175  _equation_system->assemble(x);
176  // TODO if the equation system would be reset to nullptr after each
177  // assemble() or solve() call, the user would be forced to set the
178  // equation every time and could not forget it.
179 }
180 
182  GlobalVector& x,
183  std::function<void(int, GlobalVector const&)> const&
184  postIterationCallback)
185 {
186  namespace LinAlg = MathLib::LinAlg;
187  auto& sys = *_equation_system;
188 
190  _res_id);
191  auto& minus_delta_x =
193  _minus_delta_x_id);
194  auto& J =
196 
197  bool error_norms_met = false;
198 
199  // TODO be more efficient
200  // init minus_delta_x to the right size
201  LinAlg::copy(x, minus_delta_x);
202 
203  _convergence_criterion->preFirstIteration();
204 
205  int iteration = 1;
206  for (; iteration <= _maxiter;
207  ++iteration, _convergence_criterion->reset())
208  {
209  BaseLib::RunTime timer_dirichlet;
210  double time_dirichlet = 0.0;
211 
212  BaseLib::RunTime time_iteration;
213  time_iteration.start();
214 
215  timer_dirichlet.start();
216  sys.computeKnownSolutions(x);
217  sys.applyKnownSolutions(x);
218  time_dirichlet += timer_dirichlet.elapsed();
219 
220  sys.preIteration(iteration, x);
221 
222  BaseLib::RunTime time_assembly;
223  time_assembly.start();
224  sys.assemble(x);
225  sys.getResidual(x, res);
226  sys.getJacobian(J);
227  INFO("[time] Assembly took %g s.", time_assembly.elapsed());
228 
229  minus_delta_x.setZero();
230 
231  timer_dirichlet.start();
232  sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
233  time_dirichlet += timer_dirichlet.elapsed();
234  INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet);
235 
236  if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
237  _convergence_criterion->checkResidual(res);
238 
239  BaseLib::RunTime time_linear_solver;
240  time_linear_solver.start();
241  bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
242  INFO("[time] Linear solver took %g s.", time_linear_solver.elapsed());
243 
244  if (!iteration_succeeded)
245  {
246  ERR("Newton: The linear solver failed.");
247  }
248  else
249  {
250  // TODO could be solved in a better way
251  // cf.
252  // http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecWAXPY.html
253  auto& x_new =
255  LinAlg::axpy(x_new, -_damping, minus_delta_x);
256 
257  if (postIterationCallback)
258  postIterationCallback(iteration, x_new);
259 
260  switch(sys.postIteration(x_new))
261  {
263  break;
265  ERR("Newton: The postIteration() hook reported a "
266  "non-recoverable error.");
267  iteration_succeeded = false;
268  break;
270  INFO(
271  "Newton: The postIteration() hook decided that this "
272  "iteration"
273  " has to be repeated.");
274  // TODO introduce some onDestroy hook.
276  continue; // That throws the iteration result away.
277  }
278 
279  LinAlg::copy(x_new, x); // copy new solution to x
281  }
282 
283  if (!iteration_succeeded)
284  {
285  // Don't compute further error norms, but break here.
286  error_norms_met = false;
287  break;
288  }
289 
290  if (sys.isLinear()) {
291  error_norms_met = true;
292  } else {
293  if (_convergence_criterion->hasDeltaXCheck()) {
294  // Note: x contains the new solution!
295  _convergence_criterion->checkDeltaX(minus_delta_x, x);
296  }
297 
298  error_norms_met = _convergence_criterion->isSatisfied();
299  }
300 
301  INFO("[time] Iteration #%u took %g s.", iteration,
302  time_iteration.elapsed());
303 
304  if (error_norms_met)
305  break;
306 
307  // Avoid increment of the 'iteration' if the error norms are not met,
308  // but maximum number of iterations is reached.
309  if (iteration >= _maxiter)
310  break;
311  }
312 
313  if (iteration > _maxiter)
314  {
315  ERR("Newton: Could not solve the given nonlinear system within %u "
316  "iterations",
317  _maxiter);
318  }
319 
323  minus_delta_x);
324 
325  return {error_norms_met, iteration};
326 }
327 
328 std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
329 createNonlinearSolver(GlobalLinearSolver& linear_solver,
330  BaseLib::ConfigTree const& config)
331 {
333  auto const type = config.getConfigParameter<std::string>("type");
335  auto const max_iter = config.getConfigParameter<int>("max_iter");
336 
337  if (type == "Picard") {
338  auto const tag = NonlinearSolverTag::Picard;
339  using ConcreteNLS = NonlinearSolver<tag>;
340  return std::make_pair(
341  std::make_unique<ConcreteNLS>(linear_solver, max_iter), tag);
342  }
343  if (type == "Newton")
344  {
346  auto const damping = config.getConfigParameter<double>("damping", 1.0);
347  if (damping <= 0)
348  {
349  OGS_FATAL(
350  "The damping factor for the Newon method must be positive, got "
351  "%g.",
352  damping);
353  }
354  auto const tag = NonlinearSolverTag::Newton;
355  using ConcreteNLS = NonlinearSolver<tag>;
356  return std::make_pair(
357  std::make_unique<ConcreteNLS>(linear_solver, max_iter, damping),
358  tag);
359  }
360  OGS_FATAL("Unsupported nonlinear solver type");
361 }
362 }
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:71
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