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