Reference documentation for deal.II version GIT e22cb6a53e 2023-09-21 14:00:02+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
kinsol.cc
Go to the documentation of this file.
1 // ---------------------------------------------------------------------
2 //
3 // Copyright (C) 2017 - 2023 by the deal.II authors
4 //
5 // This file is part of the deal.II library.
6 //
7 // The deal.II library is free software; you can use it, redistribute
8 // it, and/or modify it under the terms of the GNU Lesser General
9 // Public License as published by the Free Software Foundation; either
10 // version 2.1 of the License, or (at your option) any later version.
11 // The full text of the license can be found in the file LICENSE.md at
12 // the top level directory of deal.II.
13 //
14 // ---------------------------------------------------------------------
15 
16 
17 
18 #include <deal.II/base/config.h>
19 
21 
22 #ifdef DEAL_II_WITH_SUNDIALS
23 
24 # include <deal.II/base/scope_exit.h>
25 # include <deal.II/base/utilities.h>
26 
30 # ifdef DEAL_II_WITH_TRILINOS
33 # endif
34 # ifdef DEAL_II_WITH_PETSC
37 # endif
38 
41 
42 // Make sure we #include the SUNDIALS config file...
43 # include <sundials/sundials_config.h>
44 // ...before the rest of the SUNDIALS files:
45 # include <kinsol/kinsol_direct.h>
46 # include <sunlinsol/sunlinsol_dense.h>
47 # include <sunmatrix/sunmatrix_dense.h>
48 
49 # include <iomanip>
50 # include <iostream>
51 
53 
54 namespace SUNDIALS
55 {
56  template <typename VectorType>
58  const SolutionStrategy &strategy,
59  const unsigned int maximum_non_linear_iterations,
60  const double function_tolerance,
61  const double step_tolerance,
62  const bool no_init_setup,
63  const unsigned int maximum_setup_calls,
64  const double maximum_newton_step,
65  const double dq_relative_error,
66  const unsigned int maximum_beta_failures,
67  const unsigned int anderson_subspace_size)
68  : strategy(strategy)
69  , maximum_non_linear_iterations(maximum_non_linear_iterations)
70  , function_tolerance(function_tolerance)
71  , step_tolerance(step_tolerance)
72  , no_init_setup(no_init_setup)
73  , maximum_setup_calls(maximum_setup_calls)
74  , maximum_newton_step(maximum_newton_step)
75  , dq_relative_error(dq_relative_error)
76  , maximum_beta_failures(maximum_beta_failures)
77  , anderson_subspace_size(anderson_subspace_size)
78  {}
79 
80 
81 
82  template <typename VectorType>
83  void
85  {
86  static std::string strategy_str("newton");
87  prm.add_parameter("Solution strategy",
88  strategy_str,
89  "Choose among newton|linesearch|fixed_point|picard",
91  "newton|linesearch|fixed_point|picard"));
92  prm.add_action("Solution strategy", [&](const std::string &value) {
93  if (value == "newton")
94  strategy = newton;
95  else if (value == "linesearch")
96  strategy = linesearch;
97  else if (value == "fixed_point")
98  strategy = fixed_point;
99  else if (value == "picard")
100  strategy = picard;
101  else
102  Assert(false, ExcInternalError());
103  });
104  prm.add_parameter("Maximum number of nonlinear iterations",
105  maximum_non_linear_iterations);
106  prm.add_parameter("Function norm stopping tolerance", function_tolerance);
107  prm.add_parameter("Scaled step stopping tolerance", step_tolerance);
108 
109  prm.enter_subsection("Newton parameters");
110  prm.add_parameter("No initial matrix setup", no_init_setup);
111  prm.add_parameter("Maximum iterations without matrix setup",
112  maximum_setup_calls);
113  prm.add_parameter("Maximum allowable scaled length of the Newton step",
114  maximum_newton_step);
115  prm.add_parameter("Relative error for different quotient computation",
116  dq_relative_error);
117  prm.leave_subsection();
118 
119  prm.enter_subsection("Linesearch parameters");
120  prm.add_parameter("Maximum number of beta-condition failures",
121  maximum_beta_failures);
122  prm.leave_subsection();
123 
124 
125  prm.enter_subsection("Fixed point and Picard parameters");
126  prm.add_parameter("Anderson acceleration subspace size",
127  anderson_subspace_size);
128  prm.leave_subsection();
129  }
130 
131 
132 
133  template <typename VectorType>
135  : KINSOL(data, MPI_COMM_SELF)
136  {}
137 
138 
139 
140  template <typename VectorType>
142  const MPI_Comm mpi_comm)
143  : data(data)
144  , mpi_communicator(mpi_comm)
145  , kinsol_mem(nullptr)
146 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
147  , kinsol_ctx(nullptr)
148 # endif
149  , pending_exception(nullptr)
150  {
152 
153 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
154  // SUNDIALS will always duplicate communicators if we provide them. This
155  // can cause problems if SUNDIALS is configured with MPI and we pass along
156  // MPI_COMM_SELF in a serial application as MPI won't be
157  // initialized. Hence, work around that by just not providing a
158  // communicator in that case.
159  const int status =
160  SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
162  &kinsol_ctx);
163  (void)status;
164  AssertKINSOL(status);
165 # endif
166  }
167 
168 
169 
170  template <typename VectorType>
172  {
173  KINFree(&kinsol_mem);
174 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
175  const int status = SUNContext_Free(&kinsol_ctx);
176  (void)status;
177  AssertKINSOL(status);
178 # endif
179 
180  Assert(pending_exception == nullptr, ExcInternalError());
181  }
182 
183 
184 
185  template <typename VectorType>
186  unsigned int
187  KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
188  {
189  // Make sure we have what we need
190  if (data.strategy == AdditionalData::fixed_point)
191  {
192  Assert(iteration_function,
193  ExcFunctionNotProvided("iteration_function"));
194  }
195  else
196  {
197  Assert(residual, ExcFunctionNotProvided("residual"));
198  Assert(solve_with_jacobian,
199  ExcFunctionNotProvided("solve_with_jacobian"));
200  }
201 
202  // Create a new solver object:
203  int status = 0;
204  (void)status;
205 
206  KINFree(&kinsol_mem);
207 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
208  status = SUNContext_Free(&kinsol_ctx);
209  AssertKINSOL(status);
210 # endif
211 
212 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
213  kinsol_mem = KINCreate();
214 # else
215  // Same comment applies as in class constructor:
216  status =
217  SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
218  &mpi_communicator,
219  &kinsol_ctx);
220  AssertKINSOL(status);
221 
222  kinsol_mem = KINCreate(kinsol_ctx);
223 # endif
224 
225  status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
226  AssertKINSOL(status);
227 
228  // helper function to create N_Vectors compatible with different versions
229  // of SUNDIALS
230  const auto make_compatible_nvector_view =
231 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
232  [](auto &v) { return internal::make_nvector_view(v); };
233 # else
234  [this](auto &v) { return internal::make_nvector_view(v, kinsol_ctx); };
235 # endif
236 
237 
238  VectorType ones;
239  // Prepare constant vector for scaling f or u only if we need it
240  if (!get_function_scaling || !get_solution_scaling)
241  {
242  reinit_vector(ones);
243  ones = 1.0;
244  }
245 
246  auto u_scale = make_compatible_nvector_view(
247  get_solution_scaling ? get_solution_scaling() : ones);
248  auto f_scale = make_compatible_nvector_view(
249  get_function_scaling ? get_function_scaling() : ones);
250 
251  auto solution = make_compatible_nvector_view(initial_guess_and_solution);
252 
253  // This must be called before KINSetMAA
254  status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
255  AssertKINSOL(status);
256 
257  // From the manual: this must be called BEFORE KINInit
258  status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
259  AssertKINSOL(status);
260 
261  if (data.strategy == AdditionalData::fixed_point)
262  status = KINInit(
263  kinsol_mem,
264  /* wrap up the iteration_function() callback: */
265  [](N_Vector yy, N_Vector FF, void *user_data) -> int {
266  KINSOL<VectorType> &solver =
267  *static_cast<KINSOL<VectorType> *>(user_data);
268 
269  auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
270  auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
271 
273 
275  solver.iteration_function,
276  solver.pending_exception,
277  *src_yy,
278  *dst_FF);
279 
280  return err;
281  },
282  solution);
283  else
284  status = KINInit(
285  kinsol_mem,
286  /* wrap up the residual() callback: */
287  [](N_Vector yy, N_Vector FF, void *user_data) -> int {
288  KINSOL<VectorType> &solver =
289  *static_cast<KINSOL<VectorType> *>(user_data);
290 
291  auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
292  auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
293 
294  Assert(solver.residual, ExcInternalError());
295 
297  solver.residual, solver.pending_exception, *src_yy, *dst_FF);
298 
299  return err;
300  },
301  solution);
302  AssertKINSOL(status);
303 
304  status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
305  AssertKINSOL(status);
306 
307  status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
308  AssertKINSOL(status);
309 
310  status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
311  AssertKINSOL(status);
312 
313  status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
314  AssertKINSOL(status);
315 
316  status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
317  AssertKINSOL(status);
318 
319  status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
320  AssertKINSOL(status);
321 
322  status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
323  AssertKINSOL(status);
324 
325  SUNMatrix J = nullptr;
326  SUNLinearSolver LS = nullptr;
327 
328  // user assigned a function object to the solver slot
329  if (solve_with_jacobian)
330  {
331  // Set the operations we care for in the sun_linear_solver object
332  // and attach it to the KINSOL object. The functions that will get
333  // called do not actually receive the KINSOL object, just the LS
334  // object, so we have to store a pointer to the current
335  // object in the LS object
336 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
337  LS = SUNLinSolNewEmpty();
338 # else
339  LS = SUNLinSolNewEmpty(kinsol_ctx);
340 # endif
341  LS->content = this;
342 
343  LS->ops->gettype =
344  [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
345  return SUNLINEARSOLVER_MATRIX_ITERATIVE;
346  };
347 
348  LS->ops->free = [](SUNLinearSolver LS) -> int {
349  if (LS->content)
350  {
351  LS->content = nullptr;
352  }
353  if (LS->ops)
354  {
355  free(LS->ops);
356  LS->ops = nullptr;
357  }
358  free(LS);
359  LS = nullptr;
360  return 0;
361  };
362 
363  LS->ops->solve = [](SUNLinearSolver LS,
364  SUNMatrix /*ignored*/,
365  N_Vector x,
366  N_Vector b,
367  realtype tol) -> int {
368  // Receive the object that describes the linear solver and
369  // unpack the pointer to the KINSOL object from which we can then
370  // get the 'reinit' and 'solve' functions.
371  const KINSOL<VectorType> &solver =
372  *static_cast<const KINSOL<VectorType> *>(LS->content);
373 
374  Assert(solver.solve_with_jacobian, ExcInternalError());
375 
376  auto src_b = internal::unwrap_nvector_const<VectorType>(b);
377  auto dst_x = internal::unwrap_nvector<VectorType>(x);
378 
380  solver.solve_with_jacobian,
381  solver.pending_exception,
382  *src_b,
383  *dst_x,
384  tol);
385 
386  return err;
387  };
388 
389  // Even though we don't use it, KINSOL still wants us to set some
390  // kind of matrix object for the nonlinear solver. This is because
391  // if we don't set it, it won't call the functions that set up
392  // the matrix object (i.e., the argument to the 'KINSetJacFn'
393  // function below).
394 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
395  J = SUNMatNewEmpty();
396 # else
397  J = SUNMatNewEmpty(kinsol_ctx);
398 # endif
399  J->content = this;
400 
401  J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
402  return SUNMATRIX_CUSTOM;
403  };
404 
405  J->ops->destroy = [](SUNMatrix A) {
406  if (A->content)
407  {
408  A->content = nullptr;
409  }
410  if (A->ops)
411  {
412  free(A->ops);
413  A->ops = nullptr;
414  }
415  free(A);
416  A = nullptr;
417  };
418 
419  // Now set the linear system and Jacobian objects in the solver:
420  status = KINSetLinearSolver(kinsol_mem, LS, J);
421  AssertKINSOL(status);
422 
423  // Finally, if we were given a set-up function, tell KINSOL about
424  // it as well. The manual says that this must happen *after*
425  // calling KINSetLinearSolver
426  if (!setup_jacobian)
427  setup_jacobian = [](const VectorType &, const VectorType &) {
428  return 0;
429  };
430  status = KINSetJacFn(
431  kinsol_mem,
432  [](N_Vector u,
433  N_Vector f,
434  SUNMatrix /* ignored */,
435  void *user_data,
436  N_Vector /* tmp1 */,
437  N_Vector /* tmp2 */) {
438  // Receive the object that describes the linear solver and
439  // unpack the pointer to the KINSOL object from which we can then
440  // get the 'setup' function.
441  const KINSOL<VectorType> &solver =
442  *static_cast<const KINSOL<VectorType> *>(user_data);
443 
444  auto ycur = internal::unwrap_nvector_const<VectorType>(u);
445  auto fcur = internal::unwrap_nvector<VectorType>(f);
446 
447  // Call the user-provided setup function with these arguments:
449  solver.setup_jacobian, solver.pending_exception, *ycur, *fcur);
450  });
451  AssertKINSOL(status);
452  }
453 
454  // Right before calling the main KINSol() call, allow expert users to
455  // perform additional setup operations on the KINSOL object.
456  if (custom_setup)
457  custom_setup(kinsol_mem);
458 
459 
460  // Having set up all of the ancillary things, finally call the main KINSol
461  // function. One we return, check that what happened:
462  // - If we have a pending recoverable exception, ignore it if SUNDIAL's
463  // return code was zero -- in that case, SUNDIALS managed to indeed
464  // recover and we no longer need the exception
465  // - If we have any other exception, rethrow it
466  // - If no exception, test that SUNDIALS really did successfully return
467  //
468  // This all creates difficult exit paths from this function. We have to
469  // do some manual clean ups to get rid of the explicitly created
470  // temporary objects of this class. To avoid having to repeat the clean-up
471  // code on each exit path, we package it up and put the code into a
472  // ScopeExit object that is executed automatically on each such path
473  // out of this function.
474  Assert(pending_exception == nullptr, ExcInternalError());
475  status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
476 
477  ScopeExit upon_exit([this, &J, &LS]() mutable {
478  if (J != nullptr)
479  SUNMatDestroy(J);
480  if (LS != nullptr)
481  SUNLinSolFree(LS);
482  KINFree(&kinsol_mem);
483  });
484 
485  if (pending_exception)
486  {
487  try
488  {
489  std::rethrow_exception(pending_exception);
490  }
491  catch (const RecoverableUserCallbackError &exc)
492  {
493  pending_exception = nullptr;
494  if (status == 0)
495  /* just eat the exception */;
496  else
497  throw;
498  }
499  catch (...)
500  {
501  pending_exception = nullptr;
502  throw;
503  }
504  }
505  AssertKINSOL(status);
506 
507  long nniters;
508  status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
509  AssertKINSOL(status);
510 
511  return static_cast<unsigned int>(nniters);
512  }
513 
514 
515 
516  template <typename VectorType>
517  void
519  {
520  reinit_vector = [](VectorType &) {
521  AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
522  };
523  }
524 
525  template class KINSOL<Vector<double>>;
526  template class KINSOL<BlockVector<double>>;
527 
530 
531 # ifdef DEAL_II_WITH_MPI
532 
533 # ifdef DEAL_II_WITH_TRILINOS
536 # endif
537 
538 # ifdef DEAL_II_WITH_PETSC
539 # ifndef PETSC_USE_COMPLEX
540  template class KINSOL<PETScWrappers::MPI::Vector>;
542 # endif
543 # endif
544 
545 # endif
546 
547 } // namespace SUNDIALS
548 
550 
551 #endif
void add_parameter(const std::string &entry, ParameterType &parameter, const std::string &documentation="", const Patterns::PatternBase &pattern= *Patterns::Tools::Convert< ParameterType >::to_pattern(), const bool has_to_be_set=false)
void add_action(const std::string &entry, const std::function< void(const std::string &value)> &action, const bool execute_action=true)
void enter_subsection(const std::string &subsection)
AdditionalData(const SolutionStrategy &strategy=linesearch, const unsigned int maximum_non_linear_iterations=200, const double function_tolerance=0.0, const double step_tolerance=0.0, const bool no_init_setup=false, const unsigned int maximum_setup_calls=0, const double maximum_newton_step=0.0, const double dq_relative_error=0.0, const unsigned int maximum_beta_failures=0, const unsigned int anderson_subspace_size=0)
Definition: kinsol.cc:57
void add_parameters(ParameterHandler &prm)
Definition: kinsol.cc:84
void set_functions_to_trigger_an_assert()
Definition: kinsol.cc:518
MPI_Comm mpi_communicator
Definition: kinsol.h:687
unsigned int solve(VectorType &initial_guess_and_solution)
Definition: kinsol.cc:187
SUNContext kinsol_ctx
Definition: kinsol.h:698
std::function< void(const VectorType &src, VectorType &dst)> residual
Definition: kinsol.h:441
std::function< void(const VectorType &src, VectorType &dst)> iteration_function
Definition: kinsol.h:457
std::exception_ptr pending_exception
Definition: kinsol.h:712
AdditionalData data
Definition: kinsol.h:682
KINSOL(const AdditionalData &data=AdditionalData())
Definition: kinsol.cc:134
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:477
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition: config.h:378
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:478
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & RecoverableUserCallbackError()
#define Assert(cond, exc)
Definition: exceptions.h:1616
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define AssertThrow(cond, exc)
Definition: exceptions.h:1705
#define AssertKINSOL(code)
Definition: kinsol.h:48
static const char A
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
int call_and_possibly_capture_exception(const F &f, std::exception_ptr &eptr, Args &&...args)
Definition: utilities.h:46
void free(T *&pointer)
Definition: cuda.h:97