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 <kinsol/kinsol.h>
50# include <sunlinsol/sunlinsol_dense.h>
51# include <sunmatrix/sunmatrix_dense.h>
60 template <
typename VectorType>
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,
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)
88 template <
typename VectorType>
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")
101 else if (value ==
"linesearch")
102 strategy = linesearch;
103 else if (value ==
"fixed_point")
104 strategy = fixed_point;
105 else if (value ==
"picard")
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);
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",
126 prm.
add_parameter(
"Maximum number of beta-condition failures",
127 maximum_beta_failures);
133 anderson_subspace_size);
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;
163 template <
typename VectorType>
170 template <
typename VectorType>
174 , mpi_communicator(mpi_comm)
188# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
195# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
207 template <
typename VectorType>
211# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
222 template <
typename VectorType>
227 if (
data.strategy == AdditionalData::fixed_point)
229 Assert(iteration_function,
235 Assert(solve_with_jacobian,
244# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
250# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
259# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
278# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
279 [](
auto &v) {
return internal::make_nvector_view(v); };
281 [
this](
auto &v) {
return internal::make_nvector_view(v, kinsol_ctx); };
287 if (!get_function_scaling || !get_solution_scaling)
294 get_solution_scaling ? get_solution_scaling() :
ones);
296 get_function_scaling ? get_function_scaling() :
ones);
308# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
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."));
323 if (
data.strategy == AdditionalData::fixed_point)
331 auto src_yy = internal::unwrap_nvector_const<VectorType>(
yy);
332 auto dst_FF = internal::unwrap_nvector<VectorType>(
FF);
337 solver.iteration_function,
338 solver.pending_exception,
353 auto src_yy = internal::unwrap_nvector_const<VectorType>(
yy);
354 auto dst_FF = internal::unwrap_nvector<VectorType>(
FF);
359 solver.residual, solver.pending_exception, *
src_yy, *
dst_FF);
391 if (solve_with_jacobian)
398# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
413 LS->content =
nullptr;
438 auto src_b = internal::unwrap_nvector_const<VectorType>(b);
439 auto dst_x = internal::unwrap_nvector<VectorType>(
x);
442 solver.solve_with_jacobian,
443 solver.pending_exception,
456# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
470 A->content =
nullptr;
506 auto ycur = internal::unwrap_nvector_const<VectorType>(
u);
507 auto fcur = internal::unwrap_nvector<VectorType>(f);
511 solver.setup_jacobian, solver.pending_exception, *
ycur, *
fcur);
519 custom_setup(kinsol_mem);
547 if (pending_exception)
551 std::rethrow_exception(pending_exception);
555 pending_exception =
nullptr;
563 pending_exception =
nullptr;
573 return static_cast<unsigned int>(
nniters);
578 template <
typename VectorType>
582 reinit_vector = [](VectorType &) {
593# ifdef DEAL_II_WITH_MPI
595# ifdef DEAL_II_WITH_TRILINOS
600# ifdef DEAL_II_WITH_PETSC
601# 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)
OrthogonalizationStrategy
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)
void add_parameters(ParameterHandler &prm)
void set_functions_to_trigger_an_assert()
MPI_Comm mpi_communicator
unsigned int solve(VectorType &initial_guess_and_solution)
KINSOL(const AdditionalData &data=AdditionalData())
#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)
#define DEAL_II_ASSERT_UNREACHABLE()
#define AssertKINSOL(code)
std::vector< index_type > data
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)