22 #ifdef DEAL_II_WITH_SUNDIALS
30 # ifdef DEAL_II_WITH_TRILINOS
34 # ifdef DEAL_II_WITH_PETSC
43 # include <sundials/sundials_config.h>
45 # include <kinsol/kinsol_direct.h>
46 # include <sunlinsol/sunlinsol_dense.h>
47 # include <sunmatrix/sunmatrix_dense.h>
56 template <
typename VectorType>
59 const unsigned int maximum_non_linear_iterations,
60 const double function_tolerance,
61 const double step_tolerance,
62 const bool no_init_setup,
63 const unsigned int maximum_setup_calls,
64 const double maximum_newton_step,
65 const double dq_relative_error,
66 const unsigned int maximum_beta_failures,
67 const unsigned int anderson_subspace_size)
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)
82 template <
typename VectorType>
86 static std::string strategy_str(
"newton");
89 "Choose among newton|linesearch|fixed_point|picard",
91 "newton|linesearch|fixed_point|picard"));
92 prm.
add_action(
"Solution strategy", [&](
const std::string &value) {
93 if (value ==
"newton")
95 else if (value ==
"linesearch")
96 strategy = linesearch;
97 else if (value ==
"fixed_point")
98 strategy = fixed_point;
99 else if (value ==
"picard")
105 maximum_non_linear_iterations);
106 prm.
add_parameter(
"Function norm stopping tolerance", function_tolerance);
107 prm.
add_parameter(
"Scaled step stopping tolerance", step_tolerance);
112 maximum_setup_calls);
113 prm.
add_parameter(
"Maximum allowable scaled length of the Newton step",
114 maximum_newton_step);
115 prm.
add_parameter(
"Relative error for different quotient computation",
120 prm.
add_parameter(
"Maximum number of beta-condition failures",
121 maximum_beta_failures);
127 anderson_subspace_size);
133 template <
typename VectorType>
140 template <
typename VectorType>
142 const MPI_Comm mpi_comm)
144 , mpi_communicator(mpi_comm)
145 , kinsol_mem(nullptr)
147 , kinsol_ctx(nullptr)
149 , pending_exception(nullptr)
153 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
170 template <
typename VectorType>
173 KINFree(&kinsol_mem);
174 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
175 const int status = SUNContext_Free(&kinsol_ctx);
185 template <
typename VectorType>
190 if (data.strategy == AdditionalData::fixed_point)
192 Assert(iteration_function,
198 Assert(solve_with_jacobian,
206 KINFree(&kinsol_mem);
207 # if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
208 status = SUNContext_Free(&kinsol_ctx);
212 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
213 kinsol_mem = KINCreate();
217 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
222 kinsol_mem = KINCreate(kinsol_ctx);
225 status = KINSetUserData(kinsol_mem,
static_cast<void *
>(
this));
230 const auto make_compatible_nvector_view =
231 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
232 [](
auto &v) {
return internal::make_nvector_view(v); };
234 [
this](
auto &v) {
return internal::make_nvector_view(v, kinsol_ctx); };
240 if (!get_function_scaling || !get_solution_scaling)
246 auto u_scale = make_compatible_nvector_view(
247 get_solution_scaling ? get_solution_scaling() : ones);
248 auto f_scale = make_compatible_nvector_view(
249 get_function_scaling ? get_function_scaling() : ones);
251 auto solution = make_compatible_nvector_view(initial_guess_and_solution);
254 status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
258 status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
261 if (data.strategy == AdditionalData::fixed_point)
265 [](N_Vector yy, N_Vector FF,
void *user_data) ->
int {
269 auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
270 auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
287 [](N_Vector yy, N_Vector FF,
void *user_data) ->
int {
291 auto src_yy = internal::unwrap_nvector_const<VectorType>(yy);
292 auto dst_FF = internal::unwrap_nvector<VectorType>(FF);
304 status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
307 status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
310 status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
313 status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
316 status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
319 status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
322 status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
325 SUNMatrix J =
nullptr;
326 SUNLinearSolver LS =
nullptr;
329 if (solve_with_jacobian)
336 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
337 LS = SUNLinSolNewEmpty();
339 LS = SUNLinSolNewEmpty(kinsol_ctx);
344 [](SUNLinearSolver ) -> SUNLinearSolver_Type {
345 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
348 LS->ops->free = [](SUNLinearSolver LS) ->
int {
351 LS->content =
nullptr;
363 LS->ops->solve = [](SUNLinearSolver LS,
367 realtype tol) ->
int {
371 const KINSOL<VectorType> &solver =
372 *
static_cast<const KINSOL<VectorType> *
>(LS->content);
376 auto src_b = internal::unwrap_nvector_const<VectorType>(
b);
377 auto dst_x = internal::unwrap_nvector<VectorType>(x);
380 solver.solve_with_jacobian,
381 solver.pending_exception,
394 # if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
395 J = SUNMatNewEmpty();
397 J = SUNMatNewEmpty(kinsol_ctx);
401 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
402 return SUNMATRIX_CUSTOM;
405 J->ops->destroy = [](SUNMatrix
A) {
408 A->content =
nullptr;
420 status = KINSetLinearSolver(kinsol_mem, LS, J);
427 setup_jacobian = [](
const VectorType &,
const VectorType &) {
430 status = KINSetJacFn(
441 const KINSOL<VectorType> &solver =
442 *
static_cast<const KINSOL<VectorType> *
>(user_data);
444 auto ycur = internal::unwrap_nvector_const<VectorType>(u);
445 auto fcur = internal::unwrap_nvector<VectorType>(f);
449 solver.setup_jacobian, solver.pending_exception, *ycur, *fcur);
457 custom_setup(kinsol_mem);
475 status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
477 ScopeExit upon_exit([
this, &J, &LS]()
mutable {
482 KINFree(&kinsol_mem);
485 if (pending_exception)
489 std::rethrow_exception(pending_exception);
493 pending_exception =
nullptr;
501 pending_exception =
nullptr;
508 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
511 return static_cast<unsigned int>(nniters);
516 template <
typename VectorType>
520 reinit_vector = [](VectorType &) {
531 # ifdef DEAL_II_WITH_MPI
533 # ifdef DEAL_II_WITH_TRILINOS
538 # ifdef DEAL_II_WITH_PETSC
539 # 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 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)
void add_parameters(ParameterHandler &prm)
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
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 & ExcInternalError()
static ::ExceptionBase & RecoverableUserCallbackError()
#define Assert(cond, exc)
static ::ExceptionBase & ExcFunctionNotProvided(std::string arg1)
#define AssertThrow(cond, exc)
#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)