Reference documentation for deal.II version Git 191d06ed00 2021-05-11 21:15:49 -0400
\(\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 - 2021 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 #include <deal.II/base/config.h>
18 
20 
21 #ifdef DEAL_II_WITH_SUNDIALS
22 
23 # include <deal.II/base/utilities.h>
24 
28 # ifdef DEAL_II_WITH_TRILINOS
31 # endif
32 # ifdef DEAL_II_WITH_PETSC
35 # endif
36 
37 # include <deal.II/sundials/copy.h>
38 
39 // Make sure we #include the SUNDIALS config file...
40 # include <sundials/sundials_config.h>
41 // ...before the rest of the SUNDIALS files:
42 # include <kinsol/kinsol_direct.h>
43 # include <sunlinsol/sunlinsol_dense.h>
44 # include <sunmatrix/sunmatrix_dense.h>
45 
46 # include <iomanip>
47 # include <iostream>
48 
50 
51 namespace SUNDIALS
52 {
53  using namespace internal;
54 
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")
97  else if (value == "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",
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",
113  prm.add_parameter("Maximum allowable scaled length of the Newton step",
115  prm.add_parameter("Relative error for different quotient computation",
117  prm.leave_subsection();
118 
119  prm.enter_subsection("Linesearch parameters");
120  prm.add_parameter("Maximum number of beta-condition failures",
122  prm.leave_subsection();
123 
124 
125  prm.enter_subsection("Fixed point and Picard parameters");
126  prm.add_parameter("Anderson acceleration subspace size",
128  prm.leave_subsection();
129  }
130 
131 
132  namespace
133  {
134  template <typename VectorType>
135  int
136  residual_or_iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
137  {
138  KINSOL<VectorType> &solver =
139  *static_cast<KINSOL<VectorType> *>(user_data);
141 
142  typename VectorMemory<VectorType>::Pointer src_yy(mem);
143  solver.reinit_vector(*src_yy);
144 
145  typename VectorMemory<VectorType>::Pointer dst_FF(mem);
146  solver.reinit_vector(*dst_FF);
147 
148  copy(*src_yy, yy);
149 
150  int err = 0;
151  if (solver.residual)
152  err = solver.residual(*src_yy, *dst_FF);
153  else if (solver.iteration_function)
154  err = solver.iteration_function(*src_yy, *dst_FF);
155  else
156  Assert(false, ExcInternalError());
157 
158  copy(FF, *dst_FF);
159 
160  return err;
161  }
162 
163 
164 
165 # if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
166  template <typename VectorType>
167  int
168  setup_jacobian_callback(KINMem kinsol_mem)
169  {
170  KINSOL<VectorType> &solver =
171  *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
173 
174  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
175  solver.reinit_vector(*src_ycur);
176 
177  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
178  solver.reinit_vector(*src_fcur);
179 
180  copy(*src_ycur, kinsol_mem->kin_uu);
181  copy(*src_fcur, kinsol_mem->kin_fval);
182 
183  int err = solver.setup_jacobian(*src_ycur, *src_fcur);
184  return err;
185  }
186 
187 
188 
189  template <typename VectorType>
190  int
191  solve_with_jacobian_callback(KINMem kinsol_mem,
192  N_Vector x,
193  N_Vector b,
194  realtype *sJpnorm,
195  realtype *sFdotJp)
196  {
197  KINSOL<VectorType> &solver =
198  *static_cast<KINSOL<VectorType> *>(kinsol_mem->kin_user_data);
200 
201  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
202  solver.reinit_vector(*src_ycur);
203 
204  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
205  solver.reinit_vector(*src_fcur);
206 
207  copy(*src_ycur, kinsol_mem->kin_uu);
208  copy(*src_fcur, kinsol_mem->kin_fval);
209 
210  typename VectorMemory<VectorType>::Pointer src(mem);
211  solver.reinit_vector(*src);
212 
213  typename VectorMemory<VectorType>::Pointer dst(mem);
214  solver.reinit_vector(*dst);
215 
216  copy(*src, b);
217 
218  int err = solver.solve_jacobian_system(*src_ycur, *src_fcur, *src, *dst);
219  copy(x, *dst);
220 
221  *sJpnorm = N_VWL2Norm(b, kinsol_mem->kin_fscale);
222  N_VProd(b, kinsol_mem->kin_fscale, b);
223  N_VProd(b, kinsol_mem->kin_fscale, b);
224  *sFdotJp = N_VDotProd(kinsol_mem->kin_fval, b);
225 
226  return err;
227  }
228 
229 # else // SUNDIALS 5.0 or later
230 
231  template <typename VectorType>
232  int
233  setup_jacobian_callback(N_Vector u,
234  N_Vector f,
235  SUNMatrix /* ignored */,
236  void *user_data,
237  N_Vector /* tmp1 */,
238  N_Vector /* tmp2 */)
239  {
240  // Receive the object that describes the linear solver and
241  // unpack the pointer to the KINSOL object from which we can then
242  // get the 'setup' function.
243  const KINSOL<VectorType> &solver =
244  *static_cast<const KINSOL<VectorType> *>(user_data);
245 
246  // Allocate temporary (deal.II-type) vectors into which to copy the
247  // N_vectors
249  typename VectorMemory<VectorType>::Pointer ycur(mem);
250  typename VectorMemory<VectorType>::Pointer fcur(mem);
251  solver.reinit_vector(*ycur);
252  solver.reinit_vector(*fcur);
253 
254  copy(*ycur, u);
255  copy(*fcur, f);
256 
257  // Call the user-provided setup function with these arguments:
258  solver.setup_jacobian(*ycur, *fcur);
259 
260  return 0;
261  }
262 
263 
264 
265  template <typename VectorType>
266  int
267  solve_with_jacobian_callback(SUNLinearSolver LS,
268  SUNMatrix /*ignored*/,
269  N_Vector x,
270  N_Vector b,
271  realtype tol)
272  {
273  // Receive the object that describes the linear solver and
274  // unpack the pointer to the KINSOL object from which we can then
275  // get the 'reinit' and 'solve' functions.
276  const KINSOL<VectorType> &solver =
277  *static_cast<const KINSOL<VectorType> *>(LS->content);
278 
279  // This is where we have to make a decision about which of the two
280  // signals to call. Let's first check the more modern one:
281  if (solver.solve_with_jacobian)
282  {
283  // Allocate temporary (deal.II-type) vectors into which to copy the
284  // N_vectors
286  typename VectorMemory<VectorType>::Pointer src_b(mem);
287  typename VectorMemory<VectorType>::Pointer dst_x(mem);
288 
289  solver.reinit_vector(*src_b);
290  solver.reinit_vector(*dst_x);
291 
292  copy(*src_b, b);
293 
294  const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
295 
296  copy(x, *dst_x);
297 
298  return err;
299  }
300  else
301  {
302  // User has not provided the modern callback, so the fact that we are
303  // here means that they must have given us something for the old
304  // signal. Check this.
306 
307  // Allocate temporary (deal.II-type) vectors into which to copy the
308  // N_vectors
310  typename VectorMemory<VectorType>::Pointer src_ycur(mem);
311  typename VectorMemory<VectorType>::Pointer src_fcur(mem);
312  typename VectorMemory<VectorType>::Pointer src_b(mem);
313  typename VectorMemory<VectorType>::Pointer dst_x(mem);
314 
315  solver.reinit_vector(*src_b);
316  solver.reinit_vector(*dst_x);
317 
318  copy(*src_b, b);
319 
320  // Call the user-provided setup function with these arguments. Note
321  // that Sundials 4.x and later no longer provide values for
322  // src_ycur and src_fcur, and so we simply pass dummy vector in.
323  // These vectors will have zero lengths because we don't reinit them
324  // above.
325  const int err =
326  solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
327 
328  copy(x, *dst_x);
329 
330  return err;
331  }
332  }
333 
334 # endif
335  } // namespace
336 
337 
338 
339  template <typename VectorType>
341  const MPI_Comm & mpi_comm)
342  : data(data)
343  , kinsol_mem(nullptr)
344  , solution(nullptr)
345  , u_scale(nullptr)
346  , f_scale(nullptr)
348  MPI_COMM_SELF :
349  Utilities::MPI::duplicate_communicator(mpi_comm))
350  {
352  }
353 
354 
355 
356  template <typename VectorType>
358  {
359  if (kinsol_mem)
360  KINFree(&kinsol_mem);
361 
362 # ifdef DEAL_II_WITH_MPI
364  {
365  const int ierr = MPI_Comm_free(&communicator);
366  (void)ierr;
367  AssertNothrow(ierr == MPI_SUCCESS, ExcMPI(ierr));
368  }
369 # endif
370  }
371 
372 
373 
374  template <typename VectorType>
375  unsigned int
376  KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
377  {
378  unsigned int system_size = initial_guess_and_solution.size();
379 
380  // The solution is stored in
381  // solution. Here we take only a
382  // view of it.
383 # ifdef DEAL_II_WITH_MPI
385  {
386  const IndexSet is = initial_guess_and_solution.locally_owned_elements();
387  const unsigned int local_system_size = is.n_elements();
388 
389  solution =
390  N_VNew_Parallel(communicator, local_system_size, system_size);
391 
392  u_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
393  N_VConst_Parallel(1.e0, u_scale);
394 
395  f_scale = N_VNew_Parallel(communicator, local_system_size, system_size);
396  N_VConst_Parallel(1.e0, f_scale);
397  }
398  else
399 # endif
400  {
403  "Trying to use a serial code with a parallel vector."));
404  solution = N_VNew_Serial(system_size);
405  u_scale = N_VNew_Serial(system_size);
406  N_VConst_Serial(1.e0, u_scale);
407  f_scale = N_VNew_Serial(system_size);
408  N_VConst_Serial(1.e0, f_scale);
409  }
410 
413 
416 
417  copy(solution, initial_guess_and_solution);
418 
419  if (kinsol_mem)
420  KINFree(&kinsol_mem);
421 
422  kinsol_mem = KINCreate();
423 
424  int status =
425  KINInit(kinsol_mem, residual_or_iteration_callback<VectorType>, solution);
426  (void)status;
427  AssertKINSOL(status);
428 
429  status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
430  AssertKINSOL(status);
431 
432  status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
433  AssertKINSOL(status);
434 
435  status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
436  AssertKINSOL(status);
437 
438  status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
439  AssertKINSOL(status);
440 
441  status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
442  AssertKINSOL(status);
443 
444  status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
445  AssertKINSOL(status);
446 
447  status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
448  AssertKINSOL(status);
449 
450  status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
451  AssertKINSOL(status);
452 
453  status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
454  AssertKINSOL(status);
455 
456  status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
457  AssertKINSOL(status);
458 
459  SUNMatrix J = nullptr;
460  SUNLinearSolver LS = nullptr;
461 
462  if (solve_jacobian_system ||
463  solve_with_jacobian) // user assigned a function object to the solver
464  // slot
465  {
466 /* interface up to and including 4.0 */
467 # if DEAL_II_SUNDIALS_VERSION_LT(4, 1, 0)
468  auto KIN_mem = static_cast<KINMem>(kinsol_mem);
469  KIN_mem->kin_lsolve = solve_with_jacobian_callback<VectorType>;
470  if (setup_jacobian) // user assigned a function object to the Jacobian
471  // set-up slot
472  KIN_mem->kin_lsetup = setup_jacobian_callback<VectorType>;
473 
474 /* interface up to and including 4.1 */
475 # elif DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
476 
477  // deal.II does not currently have support for KINSOL in
478  // SUNDIALS 4.1. One could write this and update this section,
479  // but it does not seem worthwhile spending the time to
480  // interface with an old version of SUNDIAL given that the
481  // code below supports modern SUNDIAL versions just fine.
482  Assert(false, ExcNotImplemented());
483 
484 # else /* interface starting with SUNDIALS 5.0 */
485  // Set the operations we care for in the sun_linear_solver object
486  // and attach it to the KINSOL object. The functions that will get
487  // called do not actually receive the KINSOL object, just the LS
488  // object, so we have to store a pointer to the current
489  // object in the LS object
490  LS = SUNLinSolNewEmpty();
491  LS->content = this;
492 
493  LS->ops->gettype =
494  [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
495  return SUNLINEARSOLVER_MATRIX_ITERATIVE;
496  };
497 
498  LS->ops->free = [](SUNLinearSolver LS) -> int {
499  if (LS->content)
500  {
501  LS->content = nullptr;
502  }
503  if (LS->ops)
504  {
505  free(LS->ops);
506  LS->ops = nullptr;
507  }
508  free(LS);
509  LS = nullptr;
510  return 0;
511  };
512 
513  LS->ops->solve = solve_with_jacobian_callback<VectorType>;
514 
515  // Even though we don't use it, KINSOL still wants us to set some
516  // kind of matrix object for the nonlinear solver. This is because
517  // if we don't set it, it won't call the functions that set up
518  // the matrix object (i.e., the argument to the 'KINSetJacFn'
519  // function below).
520  J = SUNMatNewEmpty();
521  J->content = this;
522 
523  J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
524  return SUNMATRIX_CUSTOM;
525  };
526 
527  J->ops->destroy = [](SUNMatrix A) {
528  if (A->content)
529  {
530  A->content = nullptr;
531  }
532  if (A->ops)
533  {
534  free(A->ops);
535  A->ops = nullptr;
536  }
537  free(A);
538  A = nullptr;
539  };
540 
541  // Now set the linear system and Jacobian objects in the solver:
542  status = KINSetLinearSolver(kinsol_mem, LS, J);
543  AssertKINSOL(status);
544 
545  // Finally, if we were given a set-up function, tell KINSOL about
546  // it as well. The manual says that this must happen *after*
547  // calling KINSetLinearSolver
548  if (setup_jacobian)
549  {
550  status =
551  KINSetJacFn(kinsol_mem, &setup_jacobian_callback<VectorType>);
552  AssertKINSOL(status);
553  }
554 # endif
555  }
556  else
557  {
558  J = SUNDenseMatrix(system_size, system_size);
559  LS = SUNDenseLinearSolver(u_scale, J);
560  status = KINDlsSetLinearSolver(kinsol_mem, LS, J);
561  AssertKINSOL(status);
562  }
563 
566  Assert(residual, ExcFunctionNotProvided("residual"));
567 
570  Assert(iteration_function, ExcFunctionNotProvided("iteration_function"));
571 
572  // call to KINSol
573  status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
574  AssertKINSOL(status);
575 
576  copy(initial_guess_and_solution, solution);
577 
578  // Free the vectors which are no longer used.
579 # ifdef DEAL_II_WITH_MPI
581  {
582  N_VDestroy_Parallel(solution);
583  N_VDestroy_Parallel(u_scale);
584  N_VDestroy_Parallel(f_scale);
585  }
586  else
587 # endif
588  {
589  N_VDestroy_Serial(solution);
590  N_VDestroy_Serial(u_scale);
591  N_VDestroy_Serial(f_scale);
592  }
593 
594  long nniters;
595  status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
596  AssertKINSOL(status);
597 
598  SUNMatDestroy(J);
599  SUNLinSolFree(LS);
600  KINFree(&kinsol_mem);
601 
602  return static_cast<unsigned int>(nniters);
603  }
604 
605  template <typename VectorType>
606  void
608  {
609  reinit_vector = [](VectorType &) {
610  AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
611  };
612  }
613 
614  template class KINSOL<Vector<double>>;
615  template class KINSOL<BlockVector<double>>;
616 
619 
620 # ifdef DEAL_II_WITH_MPI
621 
622 # ifdef DEAL_II_WITH_TRILINOS
625 # endif
626 
627 # ifdef DEAL_II_WITH_PETSC
628 # ifndef PETSC_USE_COMPLEX
629  template class KINSOL<PETScWrappers::MPI::Vector>;
631 # endif
632 # endif
633 
634 # endif
635 
636 } // namespace SUNDIALS
637 
639 
640 #endif
#define AssertNothrow(cond, exc)
Definition: exceptions.h:1528
N_Vector solution
Definition: kinsol.h:717
unsigned int maximum_non_linear_iterations
Definition: kinsol.h:330
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
std::function< int(const VectorType &src, VectorType &dst)> residual
Definition: kinsol.h:446
unsigned int maximum_setup_calls
Definition: kinsol.h:366
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)
std::function< int(const VectorType &src, VectorType &dst)> iteration_function
Definition: kinsol.h:463
KINSOL(const AdditionalData &data=AdditionalData(), const MPI_Comm &mpi_comm=MPI_COMM_WORLD)
Definition: kinsol.cc:340
void set_functions_to_trigger_an_assert()
Definition: kinsol.cc:607
#define AssertThrow(cond, exc)
Definition: exceptions.h:1575
AdditionalData data
Definition: kinsol.h:707
SolutionStrategy strategy
Definition: kinsol.h:325
N_Vector u_scale
Definition: kinsol.h:722
std::function< VectorType &()> get_function_scaling
Definition: kinsol.h:674
MPI_Comm communicator
Definition: kinsol.h:732
unsigned int anderson_subspace_size
Definition: kinsol.h:397
N_Vector f_scale
Definition: kinsol.h:727
std::function< int(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
Definition: kinsol.h:618
void add_parameters(ParameterHandler &prm)
Definition: kinsol.cc:84
void enter_subsection(const std::string &subsection)
#define Assert(cond, exc)
Definition: exceptions.h:1465
GrowingVectorMemory< VectorType > mem
Definition: kinsol.h:737
std::function< VectorType &()> get_solution_scaling
Definition: kinsol.h:658
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:396
unsigned int maximum_beta_failures
Definition: kinsol.h:389
std::function< int(const VectorType &ycur, const VectorType &fcur, const VectorType &rhs, VectorType &dst)> solve_jacobian_system
Definition: kinsol.h:575
void * kinsol_mem
Definition: kinsol.h:712
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
std::function< void(VectorType &)> reinit_vector
Definition: kinsol.h:431
#define AssertKINSOL(code)
Definition: kinsol.h:49
static const char A
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:160
void add_action(const std::string &entry, const std::function< void(const std::string &value)> &action)
SUNLinearSolver SUNLinSolNewEmpty()
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:395
static ::ExceptionBase & ExcNotImplemented()
std::function< int(const VectorType &current_u, const VectorType &current_f)> setup_jacobian
Definition: kinsol.h:511
void free(T *&pointer)
Definition: cuda.h:97
unsigned int solve(VectorType &initial_guess_and_solution)
Definition: kinsol.cc:376
size_type n_elements() const
Definition: index_set.h:1832
void copy(const T *begin, const T *end, U *dest)
static ::ExceptionBase & ExcInternalError()