21#ifdef DEAL_II_WITH_SUNDIALS
29# ifdef DEAL_II_WITH_TRILINOS
33# ifdef DEAL_II_WITH_PETSC
42# include <sundials/sundials_config.h>
44# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
45# include <kinsol/kinsol_ls.h>
47# include <kinsol/kinsol_direct.h>
49# include <sunlinsol/sunlinsol_dense.h>
50# include <sunmatrix/sunmatrix_dense.h>
59 template <
typename VectorType>
87 template <
typename VectorType>
91 static std::string strategy_str(
"newton");
94 "Choose among newton|linesearch|fixed_point|picard",
96 "newton|linesearch|fixed_point|picard"));
98 if (
value ==
"newton")
100 else if (
value ==
"linesearch")
102 else if (
value ==
"fixed_point")
104 else if (
value ==
"picard")
118 prm.
add_parameter(
"Maximum allowable scaled length of the Newton step",
120 prm.
add_parameter(
"Relative error for different quotient computation",
125 prm.
add_parameter(
"Maximum number of beta-condition failures",
134 static std::string orthogonalization_str(
"modified_gram_schmidt");
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")
147 else if (
value ==
"inverse_compact")
149 else if (
value ==
"classical_gram_schmidt")
151 else if (
value ==
"delayed_classical_gram_schmidt")
162 template <
typename VectorType>
169 template <
typename VectorType>
187# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
194# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
206 template <
typename VectorType>
210# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
211 const int status = SUNContext_Free(&
kinsol_ctx);
221 template <
typename VectorType>
243# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
249# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
258# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
271 status = KINSetUserData(
kinsol_mem,
static_cast<void *
>(
this));
276 const auto make_compatible_nvector_view =
277# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
292 auto u_scale = make_compatible_nvector_view(
294 auto f_scale = make_compatible_nvector_view(
297 auto solution = make_compatible_nvector_view(initial_guess_and_solution);
300 status = KINSetNumMaxIters(
kinsol_mem,
data.maximum_non_linear_iterations);
307# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
309 status = KINSetOrthAA(
kinsol_mem,
data.anderson_qr_orthogonalization);
313 data.anderson_qr_orthogonalization ==
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."));
326 [](N_Vector yy, N_Vector FF,
void *user_data) ->
int {
348 [](N_Vector yy, N_Vector FF,
void *user_data) ->
int {
371 status = KINSetMaxSetupCalls(
kinsol_mem,
data.maximum_setup_calls);
377 status = KINSetMaxNewtonStep(
kinsol_mem,
data.maximum_newton_step);
380 status = KINSetMaxBetaFails(
kinsol_mem,
data.maximum_beta_failures);
386 SUNMatrix J =
nullptr;
397# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
398 LS = SUNLinSolNewEmpty();
406 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
412 LS->content =
nullptr;
424 LS->ops->solve = [](SUNLinearSolver LS,
441 solver.solve_with_jacobian,
442 solver.pending_exception,
455# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
456 J = SUNMatNewEmpty();
462 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
463 return SUNMATRIX_CUSTOM;
466 J->ops->destroy = [](SUNMatrix A) {
469 A->content =
nullptr;
481 status = KINSetLinearSolver(kinsol_mem, LS, J);
488 setup_jacobian = [](
const VectorType &,
const VectorType &) {
491 status = KINSetJacFn(
502 const KINSOL<VectorType> &solver =
503 *
static_cast<const KINSOL<VectorType> *
>(user_data);
510 solver.setup_jacobian, solver.pending_exception, *ycur, *fcur);
518 custom_setup(kinsol_mem);
536 status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
538 ScopeExit upon_exit([
this, &J, &LS]()
mutable {
543 KINFree(&kinsol_mem);
546 if (pending_exception)
550 std::rethrow_exception(pending_exception);
554 pending_exception =
nullptr;
562 pending_exception =
nullptr;
569 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
572 return static_cast<unsigned int>(nniters);
577 template <
typename VectorType>
592# ifdef DEAL_II_WITH_MPI
594# ifdef DEAL_II_WITH_TRILINOS
599# ifdef DEAL_II_WITH_PETSC
600# ifndef PETSC_USE_COMPLEX
void add_parameter(const std::string &entry, ParameterType ¶meter, 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)
unsigned int maximum_setup_calls
double maximum_newton_step
SolutionStrategy strategy
unsigned int maximum_beta_failures
OrthogonalizationStrategy
@ delayed_classical_gram_schmidt
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)
OrthogonalizationStrategy anderson_qr_orthogonalization
void add_parameters(ParameterHandler &prm)
unsigned int maximum_non_linear_iterations
unsigned int anderson_subspace_size
double function_tolerance
std::function< VectorType &()> get_solution_scaling
void set_functions_to_trigger_an_assert()
MPI_Comm mpi_communicator
unsigned int solve(VectorType &initial_guess_and_solution)
std::function< void(const VectorType &src, VectorType &dst)> residual
std::function< void(const VectorType &src, VectorType &dst)> iteration_function
std::exception_ptr pending_exception
std::function< void(const VectorType &rhs, VectorType &dst, const double tolerance)> solve_with_jacobian
std::function< VectorType &()> get_function_scaling
std::function< void(VectorType &)> reinit_vector
KINSOL(const AdditionalData &data=AdditionalData())
NVectorView< VectorType > make_nvector_view(VectorType &vector, SUNContext nvector_context)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
#define DEAL_II_NAMESPACE_CLOSE
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)
const bool IsBlockVector< VectorType >::value
#define DEAL_II_ASSERT_UNREACHABLE()
#define AssertKINSOL(code)
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)
VectorType * unwrap_nvector(N_Vector v)
const VectorType * unwrap_nvector_const(N_Vector v)