21#ifdef DEAL_II_WITH_SUNDIALS
29# ifdef DEAL_II_WITH_TRILINOS
33# ifdef DEAL_II_WITH_PETSC
42# include <arkode/arkode_arkstep.h>
43# include <sunlinsol/sunlinsol_spgmr.h>
44# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
52 template <
typename VectorType>
58 template <
typename VectorType>
77# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
84# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
96 template <
typename VectorType>
101# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
102 const int status = SUNContext_Free(&
arkode_ctx);
112 template <
typename VectorType>
118 data.initial_step_size);
125 template <
typename VectorType>
128 const double intermediate_time,
129 const bool reset_solver)
134 "The requested intermediate time is smaller than the last requested "
135 "intermediate time."));
137 const bool do_reset = reset_solver ||
arkode_mem ==
nullptr;
144 template <
typename VectorType>
188 double actual_next_time;
236 const auto status = ARKStepGetNumSteps(
arkode_mem, &n_steps);
245 template <
typename VectorType>
248 const double current_time_step,
249 const VectorType &solution)
255# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
265# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
298 void *user_data) ->
int {
318 void *user_data) ->
int {
339 initial_condition_nvector
356 ARKStepSVtolerances(
arkode_mem,
data.relative_tolerance, abs_tols);
362 data.relative_tolerance,
363 data.absolute_tolerance);
368 status = ARKStepSetUserData(
arkode_mem,
this);
375 status = ARKStepSetInitStep(
arkode_mem, current_time_step);
390 template <
typename VectorType>
408 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
411# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
428# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
430 SUNLinSol_SPGMR(y_template,
435 SUNLinSol_SPGMR(y_template,
441 status = ARKStepSetLinearSolver(
arkode_mem, sun_linear_solver,
nullptr);
444 auto jacobian_times_vector_callback = [](N_Vector v,
474 void *user_data) ->
int {
491 jacobian_times_vector_setup_callback :
492 ARKLsJacTimesSetupFn(
nullptr),
493 jacobian_times_vector_callback);
505 void *user_data) ->
int {
529 auto jacobian_solver_setup_callback =
536 void *user_data) ->
int {
557 jacobian_solver_setup_callback :
558 ARKLsPrecSetupFn(
nullptr),
559 solve_with_jacobian_callback);
562 if (
data.implicit_function_is_linear)
564 status = ARKStepSetLinear(
578# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
579 SUNNonlinearSolver fixed_point_solver =
580 SUNNonlinSol_FixedPoint(y_template,
581 data.anderson_acceleration_subspace);
583 SUNNonlinearSolver fixed_point_solver =
584 SUNNonlinSol_FixedPoint(y_template,
585 data.anderson_acceleration_subspace,
589 status = ARKStepSetNonlinearSolver(
arkode_mem, fixed_point_solver);
594 ARKStepSetMaxNonlinIters(
arkode_mem,
data.maximum_non_linear_iterations);
600 template <
typename VectorType>
613 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
617# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
632# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
633 sun_mass_linear_solver =
634 SUNLinSol_SPGMR(y_template,
638 sun_mass_linear_solver =
639 SUNLinSol_SPGMR(y_template,
647 data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
649 status = ARKStepSetMassLinearSolver(
arkode_mem,
650 sun_mass_linear_solver,
652 mass_time_dependent);
655 auto mass_matrix_times_vector_setup_callback =
665 auto mass_matrix_times_vector_callback = [](N_Vector v,
668 void *mtimes_data) ->
int {
686 mass_matrix_times_vector_setup_callback :
687 ARKLsMassTimesSetupFn(
nullptr),
688 mass_matrix_times_vector_callback,
694 auto mass_matrix_solver_setup_callback =
709 void *user_data) ->
int {
730 mass_matrix_solver_setup_callback :
731 ARKLsMassPrecSetupFn(
nullptr),
732 solve_with_mass_matrix_callback);
740 template <
typename VectorType>
751 template <
typename VectorType>
764# ifdef DEAL_II_WITH_MPI
766# ifdef DEAL_II_WITH_TRILINOS
772# ifdef DEAL_II_WITH_PETSC
773# ifndef PETSC_USE_COMPLEX
#define AssertARKode(code)
double get_next_step_size() const
void set_desired_next_step_size(const double time_step_size)
unsigned int get_step_number() const
double get_current_time() const
void set_next_step_size(const double time_step_size)
double get_next_time() const
double get_previous_step_size() const
double get_end_time() const
void setup_mass_solver(const VectorType &solution)
std::function< void(const double t)> mass_preconditioner_setup
std::function< void(const double t, const VectorType &y, VectorType &explicit_f)> explicit_function
std::exception_ptr pending_exception
std::function< void(const double t, const VectorType &y, const VectorType &fy, const VectorType &r, VectorType &z, const double gamma, const double tol, const int lr)> jacobian_preconditioner_solve
ARKode(const AdditionalData &data=AdditionalData())
void * get_arkode_memory() const
std::function< void(const double t, const VectorType &sol, const unsigned int step_number)> output_step
std::function< VectorType &()> get_local_tolerances
LinearSolveFunction< VectorType > solve_linearized_system
std::function< void(const double t, const VectorType &y, VectorType &res)> implicit_function
std::function< void(const double t, const VectorType &r, VectorType &z, const double tol, const int lr)> mass_preconditioner_solve
MPI_Comm mpi_communicator
void reset(const double t, const double h, const VectorType &y)
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > mass_solver
std::function< void(const double t, const VectorType &y, const VectorType &fy, const int jok, int &jcur, const double gamma)> jacobian_preconditioner_setup
std::function< bool(const double t, VectorType &sol)> solver_should_restart
unsigned int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
std::function< void(const double t)> mass_times_setup
void set_functions_to_trigger_an_assert()
std::unique_ptr< internal::LinearSolverWrapper< VectorType > > linear_solver
std::function< void(const double t, const VectorType &v, VectorType &Mv)> mass_times_vector
LinearSolveFunction< VectorType > solve_mass
std::function< void(const double t, const VectorType &y, const VectorType &fy)> jacobian_times_setup
unsigned int solve_ode(VectorType &solution)
void setup_system_solver(const VectorType &solution)
unsigned int solve_ode_incrementally(VectorType &solution, const double intermediate_time, const bool reset_solver=false)
std::function< void(void *arkode_mem)> custom_setup
std::function< void(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
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)
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)