Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.4.1
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
kinsol.cc
Go to the documentation of this file.
1// ---------------------------------------------------------------------
2//
3// Copyright (C) 2017 - 2022 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
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 : strategy(strategy)
68 , maximum_non_linear_iterations(maximum_non_linear_iterations)
69 , function_tolerance(function_tolerance)
70 , step_tolerance(step_tolerance)
71 , no_init_setup(no_init_setup)
72 , maximum_setup_calls(maximum_setup_calls)
73 , maximum_newton_step(maximum_newton_step)
74 , dq_relative_error(dq_relative_error)
75 , maximum_beta_failures(maximum_beta_failures)
76 , anderson_subspace_size(anderson_subspace_size)
77 {}
78
79
80
81 template <typename VectorType>
82 void
84 {
85 static std::string strategy_str("newton");
86 prm.add_parameter("Solution strategy",
87 strategy_str,
88 "Choose among newton|linesearch|fixed_point|picard",
90 "newton|linesearch|fixed_point|picard"));
91 prm.add_action("Solution strategy", [&](const std::string &value) {
92 if (value == "newton")
93 strategy = newton;
94 else if (value == "linesearch")
95 strategy = linesearch;
96 else if (value == "fixed_point")
97 strategy = fixed_point;
98 else if (value == "picard")
99 strategy = picard;
100 else
101 Assert(false, ExcInternalError());
102 });
103 prm.add_parameter("Maximum number of nonlinear iterations",
104 maximum_non_linear_iterations);
105 prm.add_parameter("Function norm stopping tolerance", function_tolerance);
106 prm.add_parameter("Scaled step stopping tolerance", step_tolerance);
107
108 prm.enter_subsection("Newton parameters");
109 prm.add_parameter("No initial matrix setup", no_init_setup);
110 prm.add_parameter("Maximum iterations without matrix setup",
111 maximum_setup_calls);
112 prm.add_parameter("Maximum allowable scaled length of the Newton step",
113 maximum_newton_step);
114 prm.add_parameter("Relative error for different quotient computation",
115 dq_relative_error);
116 prm.leave_subsection();
117
118 prm.enter_subsection("Linesearch parameters");
119 prm.add_parameter("Maximum number of beta-condition failures",
120 maximum_beta_failures);
121 prm.leave_subsection();
122
123
124 prm.enter_subsection("Fixed point and Picard parameters");
125 prm.add_parameter("Anderson acceleration subspace size",
126 anderson_subspace_size);
127 prm.leave_subsection();
128 }
129
130
131 namespace
132 {
133 template <typename VectorType>
134 int
135 residual_callback(N_Vector yy, N_Vector FF, void *user_data)
136 {
137 KINSOL<VectorType> &solver =
138 *static_cast<KINSOL<VectorType> *>(user_data);
140
141 typename VectorMemory<VectorType>::Pointer src_yy(mem);
142 solver.reinit_vector(*src_yy);
143
144 typename VectorMemory<VectorType>::Pointer dst_FF(mem);
145 solver.reinit_vector(*dst_FF);
146
147 internal::copy(*src_yy, yy);
148
149 int err = 0;
150 if (solver.residual)
151 err = solver.residual(*src_yy, *dst_FF);
152 else
153 Assert(false, ExcInternalError());
154
155 internal::copy(FF, *dst_FF);
156
157 return err;
158 }
159
160
161
162 template <typename VectorType>
163 int
164 iteration_callback(N_Vector yy, N_Vector FF, void *user_data)
165 {
166 KINSOL<VectorType> &solver =
167 *static_cast<KINSOL<VectorType> *>(user_data);
169
170 typename VectorMemory<VectorType>::Pointer src_yy(mem);
171 solver.reinit_vector(*src_yy);
172
173 typename VectorMemory<VectorType>::Pointer dst_FF(mem);
174 solver.reinit_vector(*dst_FF);
175
176 internal::copy(*src_yy, yy);
177
178 int err = 0;
179 if (solver.iteration_function)
180 err = solver.iteration_function(*src_yy, *dst_FF);
181 else
182 Assert(false, ExcInternalError());
183
184 internal::copy(FF, *dst_FF);
185
186 return err;
187 }
188
189
190
191 template <typename VectorType>
192 int
193 setup_jacobian_callback(N_Vector u,
194 N_Vector f,
195 SUNMatrix /* ignored */,
196 void *user_data,
197 N_Vector /* tmp1 */,
198 N_Vector /* tmp2 */)
199 {
200 // Receive the object that describes the linear solver and
201 // unpack the pointer to the KINSOL object from which we can then
202 // get the 'setup' function.
203 const KINSOL<VectorType> &solver =
204 *static_cast<const KINSOL<VectorType> *>(user_data);
205
206 // Allocate temporary (deal.II-type) vectors into which to copy the
207 // N_vectors
211 solver.reinit_vector(*ycur);
212 solver.reinit_vector(*fcur);
213
214 internal::copy(*ycur, u);
215 internal::copy(*fcur, f);
216
217 // Call the user-provided setup function with these arguments:
218 solver.setup_jacobian(*ycur, *fcur);
219
220 return 0;
221 }
222
223
224
225 template <typename VectorType>
226 int
227 solve_with_jacobian_callback(SUNLinearSolver LS,
228 SUNMatrix /*ignored*/,
229 N_Vector x,
230 N_Vector b,
231 realtype tol)
232 {
233 // Receive the object that describes the linear solver and
234 // unpack the pointer to the KINSOL object from which we can then
235 // get the 'reinit' and 'solve' functions.
236 const KINSOL<VectorType> &solver =
237 *static_cast<const KINSOL<VectorType> *>(LS->content);
238
239 // This is where we have to make a decision about which of the two
240 // signals to call. Let's first check the more modern one:
241 if (solver.solve_with_jacobian)
242 {
243 // Allocate temporary (deal.II-type) vectors into which to copy the
244 // N_vectors
248
249 solver.reinit_vector(*src_b);
250 solver.reinit_vector(*dst_x);
251
252 internal::copy(*src_b, b);
253
254 const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
255
256 internal::copy(x, *dst_x);
257
258 return err;
259 }
260 else
261 {
262 // User has not provided the modern callback, so the fact that we are
263 // here means that they must have given us something for the old
264 // signal. Check this.
265 Assert(solver.solve_jacobian_system, ExcInternalError());
266
267 // Allocate temporary (deal.II-type) vectors into which to copy the
268 // N_vectors
270 typename VectorMemory<VectorType>::Pointer src_ycur(mem);
271 typename VectorMemory<VectorType>::Pointer src_fcur(mem);
274
275 solver.reinit_vector(*src_b);
276 solver.reinit_vector(*dst_x);
277
278 internal::copy(*src_b, b);
279
280 // Call the user-provided setup function with these arguments. Note
281 // that Sundials 4.x and later no longer provide values for
282 // src_ycur and src_fcur, and so we simply pass dummy vector in.
283 // These vectors will have zero lengths because we don't reinit them
284 // above.
285 const int err =
286 solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
287
288 internal::copy(x, *dst_x);
289
290 return err;
291 }
292 }
293 } // namespace
294
295
296
297 template <typename VectorType>
299 : KINSOL(data, MPI_COMM_SELF)
300 {}
301
302
303
304 template <typename VectorType>
306 const MPI_Comm & mpi_comm)
307 : data(data)
308 , mpi_communicator(mpi_comm)
309 , kinsol_mem(nullptr)
311 , kinsol_ctx(nullptr)
312# endif
313 , solution(nullptr)
314 , u_scale(nullptr)
315 , f_scale(nullptr)
316 {
318
319# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
320 // SUNDIALS will always duplicate communicators if we provide them. This
321 // can cause problems if SUNDIALS is configured with MPI and we pass along
322 // MPI_COMM_SELF in a serial application as MPI won't be
323 // initialized. Hence, work around that by just not providing a
324 // communicator in that case.
325 const int status =
326 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
328 &kinsol_ctx);
329 (void)status;
330 AssertKINSOL(status);
331# endif
332 }
333
334
335
336 template <typename VectorType>
338 {
339 KINFree(&kinsol_mem);
340# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
341 const int status = SUNContext_Free(&kinsol_ctx);
342 (void)status;
343 AssertKINSOL(status);
344# endif
345 }
346
347
348
349 template <typename VectorType>
350 unsigned int
351 KINSOL<VectorType>::solve(VectorType &initial_guess_and_solution)
352 {
353 // Make sure we have what we need
354 if (data.strategy == AdditionalData::fixed_point)
355 {
356 Assert(iteration_function,
357 ExcFunctionNotProvided("iteration_function"));
358 }
359 else
360 {
361 Assert(residual, ExcFunctionNotProvided("residual"));
362 Assert(solve_jacobian_system || solve_with_jacobian,
363 ExcFunctionNotProvided(
364 "solve_jacobian_system || solve_with_jacobian"));
365 }
366
367 // Create a new solver object:
368 int status = 0;
369 (void)status;
370
371 KINFree(&kinsol_mem);
372# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
373 status = SUNContext_Free(&kinsol_ctx);
374 AssertKINSOL(status);
375# endif
376
377# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
378 kinsol_mem = KINCreate();
379# else
380 // Same comment applies as in class constructor:
381 status =
382 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? nullptr :
383 &mpi_communicator,
384 &kinsol_ctx);
385 AssertKINSOL(status);
386
387 kinsol_mem = KINCreate(kinsol_ctx);
388# endif
389
390 status = KINSetUserData(kinsol_mem, static_cast<void *>(this));
391 AssertKINSOL(status);
392
393 const auto system_size = initial_guess_and_solution.size();
394
395 // The solution is stored in solution. Here we take only a view of it.
396# ifdef DEAL_II_WITH_MPI
398 {
399 const IndexSet &is =
400 initial_guess_and_solution.locally_owned_elements();
401 const auto local_system_size = is.n_elements();
402
403 solution = N_VNew_Parallel(mpi_communicator,
404 local_system_size,
405 system_size
407 ,
408 kinsol_ctx
409# endif
410 );
411
412 u_scale = N_VNew_Parallel(mpi_communicator,
413 local_system_size,
414 system_size
416 ,
417 kinsol_ctx
418# endif
419 );
420 N_VConst_Parallel(1.e0, u_scale);
421
422 f_scale = N_VNew_Parallel(mpi_communicator,
423 local_system_size,
424 system_size
426 ,
427 kinsol_ctx
428# endif
429 );
430 N_VConst_Parallel(1.e0, f_scale);
431 }
432 else
433# endif
434 {
436 solution = N_VNew_Serial(system_size
438 ,
439 kinsol_ctx
440# endif
441 );
442 u_scale = N_VNew_Serial(system_size
444 ,
445 kinsol_ctx
446# endif
447 );
448 N_VConst_Serial(1.e0, u_scale);
449 f_scale = N_VNew_Serial(system_size
451 ,
452 kinsol_ctx
453# endif
454 );
455 N_VConst_Serial(1.e0, f_scale);
456 }
457
458 if (get_solution_scaling)
459 internal::copy(u_scale, get_solution_scaling());
460
461 if (get_function_scaling)
462 internal::copy(f_scale, get_function_scaling());
463
464 internal::copy(solution, initial_guess_and_solution);
465
466 // This must be called before KINSetMAA
467 status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
468 AssertKINSOL(status);
469
470 // From the manual: this must be called BEFORE KINInit
471 status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
472 AssertKINSOL(status);
473
474 if (data.strategy == AdditionalData::fixed_point)
475 status = KINInit(kinsol_mem, iteration_callback<VectorType>, solution);
476 else
477 status = KINInit(kinsol_mem, residual_callback<VectorType>, solution);
478 AssertKINSOL(status);
479
480 status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
481 AssertKINSOL(status);
482
483 status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
484 AssertKINSOL(status);
485
486 status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
487 AssertKINSOL(status);
488
489 status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
490 AssertKINSOL(status);
491
492 status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
493 AssertKINSOL(status);
494
495 status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
496 AssertKINSOL(status);
497
498 status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
499 AssertKINSOL(status);
500
501 SUNMatrix J = nullptr;
502 SUNLinearSolver LS = nullptr;
503
504 if (solve_jacobian_system ||
505 solve_with_jacobian) // user assigned a function object to the solver
506 // slot
507 {
508 // Set the operations we care for in the sun_linear_solver object
509 // and attach it to the KINSOL object. The functions that will get
510 // called do not actually receive the KINSOL object, just the LS
511 // object, so we have to store a pointer to the current
512 // object in the LS object
513# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
514 LS = SUNLinSolNewEmpty();
515# else
516 LS = SUNLinSolNewEmpty(kinsol_ctx);
517# endif
518 LS->content = this;
519
520 LS->ops->gettype =
521 [](SUNLinearSolver /*ignored*/) -> SUNLinearSolver_Type {
522 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
523 };
524
525 LS->ops->free = [](SUNLinearSolver LS) -> int {
526 if (LS->content)
527 {
528 LS->content = nullptr;
529 }
530 if (LS->ops)
531 {
532 free(LS->ops);
533 LS->ops = nullptr;
534 }
535 free(LS);
536 LS = nullptr;
537 return 0;
538 };
539
540 LS->ops->solve = solve_with_jacobian_callback<VectorType>;
541
542 // Even though we don't use it, KINSOL still wants us to set some
543 // kind of matrix object for the nonlinear solver. This is because
544 // if we don't set it, it won't call the functions that set up
545 // the matrix object (i.e., the argument to the 'KINSetJacFn'
546 // function below).
547# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
548 J = SUNMatNewEmpty();
549# else
550 J = SUNMatNewEmpty(kinsol_ctx);
551# endif
552 J->content = this;
553
554 J->ops->getid = [](SUNMatrix /*ignored*/) -> SUNMatrix_ID {
555 return SUNMATRIX_CUSTOM;
556 };
557
558 J->ops->destroy = [](SUNMatrix A) {
559 if (A->content)
560 {
561 A->content = nullptr;
562 }
563 if (A->ops)
564 {
565 free(A->ops);
566 A->ops = nullptr;
567 }
568 free(A);
569 A = nullptr;
570 };
571
572 // Now set the linear system and Jacobian objects in the solver:
573 status = KINSetLinearSolver(kinsol_mem, LS, J);
574 AssertKINSOL(status);
575
576 // Finally, if we were given a set-up function, tell KINSOL about
577 // it as well. The manual says that this must happen *after*
578 // calling KINSetLinearSolver
579 if (!setup_jacobian)
580 setup_jacobian = [](const VectorType &, const VectorType &) {
581 return 0;
582 };
583 status = KINSetJacFn(kinsol_mem, &setup_jacobian_callback<VectorType>);
584 AssertKINSOL(status);
585 }
586
587 // call to KINSol
588 status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
589 AssertKINSOL(status);
590
591 internal::copy(initial_guess_and_solution, solution);
592
593 // Free the vectors which are no longer used.
594# ifdef DEAL_II_WITH_MPI
596 {
597 N_VDestroy_Parallel(solution);
598 N_VDestroy_Parallel(u_scale);
599 N_VDestroy_Parallel(f_scale);
600 }
601 else
602# endif
603 {
604 N_VDestroy_Serial(solution);
605 N_VDestroy_Serial(u_scale);
606 N_VDestroy_Serial(f_scale);
607 }
608
609 long nniters;
610 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
611 AssertKINSOL(status);
612
613 if (J != nullptr)
614 SUNMatDestroy(J);
615 if (LS != nullptr)
616 SUNLinSolFree(LS);
617 KINFree(&kinsol_mem);
618
619 return static_cast<unsigned int>(nniters);
620 }
621
622 template <typename VectorType>
623 void
625 {
626 reinit_vector = [](VectorType &) {
627 AssertThrow(false, ExcFunctionNotProvided("reinit_vector"));
628 };
629 }
630
631 template class KINSOL<Vector<double>>;
632 template class KINSOL<BlockVector<double>>;
633
636
637# ifdef DEAL_II_WITH_MPI
638
639# ifdef DEAL_II_WITH_TRILINOS
642# endif
643
644# ifdef DEAL_II_WITH_PETSC
645# ifndef PETSC_USE_COMPLEX
648# endif
649# endif
650
651# endif
652
653} // namespace SUNDIALS
654
656
657#endif
size_type n_elements() const
Definition: index_set.h:1834
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)
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:56
void add_parameters(ParameterHandler &prm)
Definition: kinsol.cc:83
void set_functions_to_trigger_an_assert()
Definition: kinsol.cc:624
MPI_Comm mpi_communicator
Definition: kinsol.h:697
unsigned int solve(VectorType &initial_guess_and_solution)
Definition: kinsol.cc:351
SUNContext kinsol_ctx
Definition: kinsol.h:708
GrowingVectorMemory< VectorType > mem
Definition: kinsol.h:729
std::function< int(const VectorType &src, VectorType &dst)> residual
Definition: kinsol.h:431
std::function< void(VectorType &)> reinit_vector
Definition: kinsol.h:416
AdditionalData data
Definition: kinsol.h:692
KINSOL(const AdditionalData &data=AdditionalData())
Definition: kinsol.cc:298
std::function< int(const VectorType &src, VectorType &dst)> iteration_function
Definition: kinsol.h:448
#define DEAL_II_NAMESPACE_OPEN
Definition: config.h:442
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
Definition: config.h:347
#define DEAL_II_NAMESPACE_CLOSE
Definition: config.h:443
#define Assert(cond, exc)
Definition: exceptions.h:1473
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
Definition: exceptions.h:1583
#define AssertKINSOL(code)
Definition: kinsol.h:47
void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src)
Definition: copy.cc:136