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