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