Reference documentation for deal.II version 9.5.0
\(\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\}}\)
Loading...
Searching...
No Matches
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
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
54namespace 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)
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_jacobian_system || solve_with_jacobian,
200 "solve_jacobian_system || solve_with_jacobian"));
201 }
202
203 // Create a new solver object:
204 int status = 0;
205 (void)status;
206
207 KINFree(&kinsol_mem);
208# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
209 status = SUNContext_Free(&kinsol_ctx);
210 AssertKINSOL(status);
211# endif
212
213# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
214 kinsol_mem = KINCreate();
215# else
216 // Same comment applies as in class constructor:
217 status =
218 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
219 &mpi_communicator,
220 &kinsol_ctx);
221 AssertKINSOL(status);
222
223 kinsol_mem = KINCreate(kinsol_ctx);
224# endif
225
226 status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
227 AssertKINSOL(status);
228
229 // helper function to create N_Vectors compatible with different versions
230 // of SUNDIALS
231 const auto make_compatible_nvector_view =
232# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
233 [](auto &v) { return internal::make_nvector_view(v); };
234# else
235 [this](auto &v) { return internal::make_nvector_view(v, kinsol_ctx); };
236# endif
237
238
239 VectorType ones;
240 // Prepare constant vector for scaling f or u only if we need it
241 if (!get_function_scaling || !get_solution_scaling)
242 {
243 reinit_vector(ones);
244 ones = 1.0;
246
247 auto u_scale = make_compatible_nvector_view(
248 get_solution_scaling ? get_solution_scaling() : ones);
249 auto f_scale = make_compatible_nvector_view(
250 get_function_scaling ? get_function_scaling() : ones);
251
252 auto solution = make_compatible_nvector_view(initial_guess_and_solution);
253
254 // This must be called before KINSetMAA
255 status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
256 AssertKINSOL(status);
257
258 // From the manual: this must be called BEFORE KINInit
259 status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
260 AssertKINSOL(status);
261
262 if (data.strategy == AdditionalData::fixed_point)
263 status = KINInit(
264 kinsol_mem,
265 /* wrap up the iteration_function() callback: */
266 [](N_Vector yy, N_Vector FF, void *user_data) -> int {
267 KINSOL<VectorType> &solver =
268 *static_cast<KINSOL<VectorType> *>(user_data);
269
270 auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
271 auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
272
274
276 solver.iteration_function,
277 solver.pending_exception,
278 *src_yy,
279 *dst_FF);
280
281 return err;
282 },
283 solution);
284 else
285 status = KINInit(
286 kinsol_mem,
287 /* wrap up the residual() callback: */
288 [](N_Vector yy, N_Vector FF, void *user_data) -> int {
289 KINSOL<VectorType> &solver =
290 *static_cast<KINSOL<VectorType> *>(user_data);
291
292 auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
293 auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
294
296
298 solver.residual, solver.pending_exception, *src_yy, *dst_FF);
299
300 return err;
301 },
302 solution);
303 AssertKINSOL(status);
304
305 status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
306 AssertKINSOL(status);
307
308 status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
309 AssertKINSOL(status);
310
311 status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
312 AssertKINSOL(status);
313
314 status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
315 AssertKINSOL(status);
316
317 status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
318 AssertKINSOL(status);
319
320 status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
321 AssertKINSOL(status);
322
323 status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
324 AssertKINSOL(status);
325
326 SUNMatrix J = nullptr;
327 SUNLinearSolver LS = nullptr;
328
329 if (solve_jacobian_system ||
330 solve_with_jacobian) // user assigned a function
331 // object to the solver slot
332 {
333 // Set the operations we care for in the sun_linear_solver object
334 // and attach it to the KINSOL object. The functions that will get
335 // called do not actually receive the KINSOL object, just the LS
336 // object, so we have to store a pointer to the current
337 // object in the LS object
338# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
339 LS = SUNLinSolNewEmpty();
340# else
341 LS = SUNLinSolNewEmpty(kinsol_ctx);
342# endif
343 LS->content = this;
344
345 LS->ops->gettype =
346 [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
347 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
348 };
349
350 LS->ops->free = [](SUNLinearSolver LS) -> int {
351 if (LS->content)
352 {
353 LS->content = nullptr;
354 }
355 if (LS->ops)
356 {
357 free(LS->ops);
358 LS->ops = nullptr;
359 }
360 free(LS);
361 LS = nullptr;
362 return 0;
363 };
364
365 LS->ops->solve = [](SUNLinearSolver LS,
366 SUNMatrix /*ignored*/,
367 N_Vector x,
368 N_Vector b,
369 realtype tol) -> int {
370 // Receive the object that describes the linear solver and
371 // unpack the pointer to the KINSOL object from which we can then
372 // get the 'reinit' and 'solve' functions.
373 const KINSOL<VectorType> &solver =
374 *static_cast<const KINSOL<VectorType> *>(LS->content);
375
376 // This is where we have to make a decision about which of the two
377 // signals to call. Let's first check the more modern one:
378 if (solver.solve_with_jacobian)
379 {
380 auto src_b = internal::unwrap_nvector_const<VectorType>(b);
381 auto dst_x = internal::unwrap_nvector<VectorType>(x);
382
384 solver.solve_with_jacobian,
385 solver.pending_exception,
386 *src_b,
387 *dst_x,
388 tol);
389
390 return err;
391 }
392 else
393 {
394 // User has not provided the modern callback, so the fact that
395 // we are here means that they must have given us something for
396 // the old signal. Check this.
397 Assert(solver.solve_jacobian_system, ExcInternalError());
398
399 // Allocate temporary (deal.II-type) dummy vectors
401 typename VectorMemory<VectorType>::Pointer src_ycur(mem);
402 typename VectorMemory<VectorType>::Pointer src_fcur(mem);
403
404 auto src_b = internal::unwrap_nvector_const<VectorType>(b);
405 auto dst_x = internal::unwrap_nvector<VectorType>(x);
406
407 // Call the user-provided setup function with these arguments.
408 // Note that Sundials 4.x and later no longer provide values for
409 // src_ycur and src_fcur, and so we simply pass dummy vector in.
410 // These vectors will have zero lengths because we don't reinit
411 // them above.
413 solver.solve_jacobian_system,
414 solver.pending_exception,
415 *src_ycur,
416 *src_fcur,
417 *src_b,
418 *dst_x);
419
420 return err;
421 }
422 };
423
424 // Even though we don't use it, KINSOL still wants us to set some
425 // kind of matrix object for the nonlinear solver. This is because
426 // if we don't set it, it won't call the functions that set up
427 // the matrix object (i.e., the argument to the 'KINSetJacFn'
428 // function below).
429# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
430 J = SUNMatNewEmpty();
431# else
432 J = SUNMatNewEmpty(kinsol_ctx);
433# endif
434 J->content = this;
435
436 J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
437 return SUNMATRIX_CUSTOM;
438 };
439
440 J->ops->destroy = [](SUNMatrix A) {
441 if (A->content)
442 {
443 A->content = nullptr;
444 }
445 if (A->ops)
446 {
447 free(A->ops);
448 A->ops = nullptr;
449 }
450 free(A);
451 A = nullptr;
452 };
453
454 // Now set the linear system and Jacobian objects in the solver:
455 status = KINSetLinearSolver(kinsol_mem, LS, J);
456 AssertKINSOL(status);
457
458 // Finally, if we were given a set-up function, tell KINSOL about
459 // it as well. The manual says that this must happen *after*
460 // calling KINSetLinearSolver
461 if (!setup_jacobian)
462 setup_jacobian = [](const VectorType &, const VectorType &) {
463 return 0;
464 };
465 status = KINSetJacFn(
466 kinsol_mem,
467 [](N_Vector u,
468 N_Vector f,
469 SUNMatrix /* ignored */,
470 void *user_data,
471 N_Vector /* tmp1 */,
472 N_Vector /* tmp2 */) {
473 // Receive the object that describes the linear solver and
474 // unpack the pointer to the KINSOL object from which we can then
475 // get the 'setup' function.
476 const KINSOL<VectorType> &solver =
477 *static_cast<const KINSOL<VectorType> *>(user_data);
478
479 auto ycur = internal::unwrap_nvector_const<VectorType>(u);
480 auto fcur = internal::unwrap_nvector<VectorType>(f);
481
482 // Call the user-provided setup function with these arguments:
484 solver.setup_jacobian, solver.pending_exception, *ycur, *fcur);
485 });
486 AssertKINSOL(status);
487 }
488
489 // Having set up all of the ancillary things, finally call the main KINSol
490 // function. One we return, check that what happened:
491 // - If we have a pending recoverable exception, ignore it if SUNDIAL's
492 // return code was zero -- in that case, SUNDIALS managed to indeed
493 // recover and we no longer need the exception
494 // - If we have any other exception, rethrow it
495 // - If no exception, test that SUNDIALS really did successfully return
496 //
497 // This all creates difficult exit paths from this function. We have to
498 // do some manual clean ups to get rid of the explicitly created
499 // temporary objects of this class. To avoid having to repeat the clean-up
500 // code on each exit path, we package it up and put the code into a
501 // ScopeExit object that is executed automatically on each such path
502 // out of this function.
503 Assert(pending_exception == nullptr, ExcInternalError());
504 status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
505
506 ScopeExit upon_exit([this, &J, &LS]() mutable {
507 if (J != nullptr)
508 SUNMatDestroy(J);
509 if (LS != nullptr)
510 SUNLinSolFree(LS);
511 KINFree(&kinsol_mem);
512 });
513
514 if (pending_exception)
515 {
516 try
517 {
518 std::rethrow_exception(pending_exception);
519 }
520 catch (const RecoverableUserCallbackError &exc)
521 {
522 pending_exception = nullptr;
523 if (status == 0)
524 /* just eat the exception */;
525 else
526 throw;
527 }
528 catch (...)
529 {
530 pending_exception = nullptr;
531 throw;
532 }
533 }
534 AssertKINSOL(status);
535
536 long nniters;
537 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
538 AssertKINSOL(status);
539
540 return static_cast<unsigned int>(nniters);
541 }
542
543
544
545 template <typename VectorType>
546 void
548 {
549 reinit_vector = [](VectorType &) {
550 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
551 };
552 }
553
554 template class KINSOL<Vector<double>>;
555 template class KINSOL<BlockVector<double>>;
556
559
560# ifdef DEAL_II_WITH_MPI
561
562# ifdef DEAL_II_WITH_TRILINOS
565# endif
566
567# ifdef DEAL_II_WITH_PETSC
568# ifndef PETSC_USE_COMPLEX
571# endif
572# endif
573
574# endif
575
576} // namespace SUNDIALS
577
579
580#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:547
MPI_Comm mpi_communicator
Definition kinsol.h:728
unsigned int solve(VectorType &initial_guess_and_solution)
Definition kinsol.cc:187
SUNContext kinsol_ctx
Definition kinsol.h:739
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:753
AdditionalData data
Definition kinsol.h:723
KINSOL(const AdditionalData &data=AdditionalData())
Definition kinsol.cc:134
#define DEAL_II_NAMESPACE_OPEN
Definition config.h:472
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition config.h:359
#define DEAL_II_NAMESPACE_CLOSE
Definition config.h:473
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & RecoverableUserCallbackError()
#define AssertThrow(cond, exc)
#define AssertKINSOL(code)
Definition kinsol.h:48
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