22#ifdef DEAL_II_WITH_SUNDIALS
29# ifdef DEAL_II_WITH_TRILINOS
33# ifdef DEAL_II_WITH_PETSC
42# include <sundials/sundials_config.h>
44# include <kinsol/kinsol_direct.h>
45# include <sunlinsol/sunlinsol_dense.h>
46# include <sunmatrix/sunmatrix_dense.h>
55 template <
typename VectorType>
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)
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)
81 template <
typename VectorType>
85 static std::string strategy_str(
"newton");
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")
94 else if (value ==
"linesearch")
95 strategy = linesearch;
96 else if (value ==
"fixed_point")
97 strategy = fixed_point;
98 else if (value ==
"picard")
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);
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",
119 prm.
add_parameter(
"Maximum number of beta-condition failures",
120 maximum_beta_failures);
126 anderson_subspace_size);
133 template <
typename VectorType>
135 residual_callback(N_Vector yy, N_Vector FF,
void *user_data)
151 err = solver.
residual(*src_yy, *dst_FF);
162 template <
typename VectorType>
164 iteration_callback(N_Vector yy, N_Vector FF,
void *user_data)
191 template <
typename VectorType>
193 setup_jacobian_callback(N_Vector u,
203 const KINSOL<VectorType> &solver =
204 *
static_cast<const KINSOL<VectorType> *
>(user_data);
211 solver.reinit_vector(*ycur);
212 solver.reinit_vector(*fcur);
218 solver.setup_jacobian(*ycur, *fcur);
225 template <
typename VectorType>
236 const KINSOL<VectorType> &solver =
237 *
static_cast<const KINSOL<VectorType> *
>(LS->content);
241 if (solver.solve_with_jacobian)
249 solver.reinit_vector(*src_b);
250 solver.reinit_vector(*dst_x);
254 const int err = solver.solve_with_jacobian(*src_b, *dst_x, tol);
275 solver.reinit_vector(*src_b);
276 solver.reinit_vector(*dst_x);
286 solver.solve_jacobian_system(*src_ycur, *src_fcur, *src_b, *dst_x);
297 template <
typename VectorType>
304 template <
typename VectorType>
308 , mpi_communicator(mpi_comm)
309 , kinsol_mem(nullptr)
311 , kinsol_ctx(nullptr)
319# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
336 template <
typename VectorType>
339 KINFree(&kinsol_mem);
340# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
341 const int status = SUNContext_Free(&kinsol_ctx);
349 template <
typename VectorType>
354 if (data.strategy == AdditionalData::fixed_point)
356 Assert(iteration_function,
357 ExcFunctionNotProvided(
"iteration_function"));
361 Assert(residual, ExcFunctionNotProvided(
"residual"));
362 Assert(solve_jacobian_system || solve_with_jacobian,
363 ExcFunctionNotProvided(
364 "solve_jacobian_system || solve_with_jacobian"));
371 KINFree(&kinsol_mem);
372# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
373 status = SUNContext_Free(&kinsol_ctx);
377# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
378 kinsol_mem = KINCreate();
382 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
387 kinsol_mem = KINCreate(kinsol_ctx);
390 status = KINSetUserData(kinsol_mem,
static_cast<void *
>(
this));
393 const auto system_size = initial_guess_and_solution.size();
396# ifdef DEAL_II_WITH_MPI
400 initial_guess_and_solution.locally_owned_elements();
401 const auto local_system_size = is.
n_elements();
403 solution = N_VNew_Parallel(mpi_communicator,
412 u_scale = N_VNew_Parallel(mpi_communicator,
420 N_VConst_Parallel(1.e0, u_scale);
422 f_scale = N_VNew_Parallel(mpi_communicator,
430 N_VConst_Parallel(1.e0, f_scale);
436 solution = N_VNew_Serial(system_size
442 u_scale = N_VNew_Serial(system_size
448 N_VConst_Serial(1.e0, u_scale);
449 f_scale = N_VNew_Serial(system_size
455 N_VConst_Serial(1.e0, f_scale);
458 if (get_solution_scaling)
461 if (get_function_scaling)
467 status = KINSetNumMaxIters(kinsol_mem, data.maximum_non_linear_iterations);
471 status = KINSetMAA(kinsol_mem, data.anderson_subspace_size);
474 if (data.strategy == AdditionalData::fixed_point)
475 status = KINInit(kinsol_mem, iteration_callback<VectorType>, solution);
477 status = KINInit(kinsol_mem, residual_callback<VectorType>, solution);
480 status = KINSetFuncNormTol(kinsol_mem, data.function_tolerance);
483 status = KINSetScaledStepTol(kinsol_mem, data.step_tolerance);
486 status = KINSetMaxSetupCalls(kinsol_mem, data.maximum_setup_calls);
489 status = KINSetNoInitSetup(kinsol_mem, data.no_init_setup);
492 status = KINSetMaxNewtonStep(kinsol_mem, data.maximum_newton_step);
495 status = KINSetMaxBetaFails(kinsol_mem, data.maximum_beta_failures);
498 status = KINSetRelErrFunc(kinsol_mem, data.dq_relative_error);
501 SUNMatrix J =
nullptr;
504 if (solve_jacobian_system ||
513# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
514 LS = SUNLinSolNewEmpty();
516 LS = SUNLinSolNewEmpty(kinsol_ctx);
522 return SUNLINEARSOLVER_MATRIX_ITERATIVE;
528 LS->content =
nullptr;
540 LS->ops->solve = solve_with_jacobian_callback<VectorType>;
547# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
548 J = SUNMatNewEmpty();
550 J = SUNMatNewEmpty(kinsol_ctx);
554 J->ops->getid = [](SUNMatrix ) -> SUNMatrix_ID {
555 return SUNMATRIX_CUSTOM;
558 J->ops->destroy = [](SUNMatrix A) {
561 A->content =
nullptr;
573 status = KINSetLinearSolver(kinsol_mem, LS, J);
580 setup_jacobian = [](
const VectorType &,
const VectorType &) {
583 status = KINSetJacFn(kinsol_mem, &setup_jacobian_callback<VectorType>);
588 status = KINSol(kinsol_mem, solution, data.strategy, u_scale, f_scale);
594# ifdef DEAL_II_WITH_MPI
597 N_VDestroy_Parallel(solution);
598 N_VDestroy_Parallel(u_scale);
599 N_VDestroy_Parallel(f_scale);
604 N_VDestroy_Serial(solution);
605 N_VDestroy_Serial(u_scale);
606 N_VDestroy_Serial(f_scale);
610 status = KINGetNumNonlinSolvIters(kinsol_mem, &nniters);
617 KINFree(&kinsol_mem);
619 return static_cast<unsigned int>(nniters);
622 template <
typename VectorType>
626 reinit_vector = [](VectorType &) {
627 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
637# ifdef DEAL_II_WITH_MPI
639# ifdef DEAL_II_WITH_TRILINOS
644# ifdef DEAL_II_WITH_PETSC
645# ifndef PETSC_USE_COMPLEX
size_type n_elements() const
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)
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)
GrowingVectorMemory< VectorType > mem
std::function< int(const VectorType &src, VectorType &dst)> residual
std::function< void(VectorType &)> reinit_vector
KINSOL(const AdditionalData &data=AdditionalData())
std::function< int(const VectorType &src, VectorType &dst)> iteration_function
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_SUNDIALS_VERSION_GTE(major, minor, patch)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
static ::ExceptionBase & ExcInternalError()
#define AssertThrow(cond, exc)
#define AssertKINSOL(code)
void copy(TrilinosWrappers::MPI::Vector &dst, const N_Vector &src)