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>
54 :
ARKode(data, MPI_COMM_SELF)
58 template <
typename VectorType>
66 , mpi_communicator(mpi_comm)
67 , last_end_time(data.initial_time)
68 , pending_exception(nullptr)
77# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
84# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
96 template <
typename VectorType>
99 ARKStepFree(&arkode_mem);
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);
120 return do_evolve_time(solution, time,
true);
125 template <
typename VectorType>
128 const double intermediate_time,
129 const bool reset_solver)
132 intermediate_time > last_end_time,
134 "The requested intermediate time is smaller than the last requested "
135 "intermediate time."));
137 const bool do_reset = reset_solver || arkode_mem ==
nullptr;
138 DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size);
139 return do_evolve_time(solution, time, do_reset);
144 template <
typename VectorType>
164 const int status = ARKStepSetStopTime(arkode_mem, time.
get_end_time());
169 auto solution_nvector = internal::make_nvector_view(solution
188 double actual_next_time;
189 const auto status = ARKStepEvolve(arkode_mem,
194 if (pending_exception)
198 std::rethrow_exception(pending_exception);
202 pending_exception =
nullptr;
210 pending_exception =
nullptr;
239 template <
typename VectorType>
242 const double current_time_step,
243 const VectorType &solution)
245 last_end_time = current_time;
249# if DEAL_II_SUNDIALS_VERSION_GTE(7, 0, 0)
250 status = SUNContext_Free(&arkode_ctx);
255 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ? SUN_COMM_NULL :
259# elif DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
260 status = SUNContext_Free(&arkode_ctx);
265 SUNContext_Create(mpi_communicator == MPI_COMM_SELF ?
nullptr :
273 ARKStepFree(&arkode_mem);
279 auto initial_condition_nvector = internal::make_nvector_view(solution
286 Assert(explicit_function || implicit_function,
292 void *user_data) ->
int {
297 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
298 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
312 void *user_data) ->
int {
317 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
318 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
328 arkode_mem = ARKStepCreate(explicit_function ? explicit_function_callback :
330 implicit_function ? implicit_function_callback :
333 initial_condition_nvector
334# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
341 if (get_local_tolerances)
343 const auto abs_tols = internal::make_nvector_view(get_local_tolerances()
350 ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
355 status = ARKStepSStolerances(arkode_mem,
356 data.relative_tolerance,
357 data.absolute_tolerance);
362 status = ARKStepSetUserData(arkode_mem,
this);
365 setup_system_solver(solution);
367 setup_mass_solver(solution);
369 status = ARKStepSetInitStep(arkode_mem, current_time_step);
372 status = ARKStepSetStopTime(arkode_mem, data.final_time);
375 status = ARKStepSetOrder(arkode_mem, data.maximum_order);
379 custom_setup(arkode_mem);
384 template <
typename VectorType>
389 if (!implicit_function)
396 if (jacobian_times_vector)
399 if (solve_linearized_system)
402 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
403 solve_linearized_system,
405# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
410 sun_linear_solver = *linear_solver;
416 auto y_template = internal::make_nvector_view(solution
422# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
424 SUNLinSol_SPGMR(y_template,
429 SUNLinSol_SPGMR(y_template,
435 status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver,
nullptr);
438 auto jacobian_times_vector_callback = [](N_Vector v,
449 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
450 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
451 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
453 auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
468 void *user_data) ->
int {
473 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
474 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
483 status = ARKStepSetJacTimes(arkode_mem,
484 jacobian_times_setup ?
485 jacobian_times_vector_setup_callback :
486 ARKLsJacTimesSetupFn(
nullptr),
487 jacobian_times_vector_callback);
489 if (jacobian_preconditioner_solve)
499 void *user_data) ->
int {
504 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
505 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
506 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
508 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
523 auto jacobian_solver_setup_callback =
530 void *user_data) ->
int {
535 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
536 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
549 status = ARKStepSetPreconditioner(arkode_mem,
550 jacobian_preconditioner_setup ?
551 jacobian_solver_setup_callback :
552 ARKLsPrecSetupFn(
nullptr),
553 solve_with_jacobian_callback);
556 if (data.implicit_function_is_linear)
558 status = ARKStepSetLinear(
559 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
565 auto y_template = internal::make_nvector_view(solution
572# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
573 SUNNonlinearSolver fixed_point_solver =
574 SUNNonlinSol_FixedPoint(y_template,
575 data.anderson_acceleration_subspace);
577 SUNNonlinearSolver fixed_point_solver =
578 SUNNonlinSol_FixedPoint(y_template,
579 data.anderson_acceleration_subspace,
583 status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
588 ARKStepSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
594 template <
typename VectorType>
601 if (mass_times_vector)
607 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
611# if DEAL_II_SUNDIALS_VERSION_GTE(6, 0, 0)
616 sun_mass_linear_solver = *mass_solver;
620 auto y_template = internal::make_nvector_view(solution
626# if DEAL_II_SUNDIALS_VERSION_LT(6, 0, 0)
627 sun_mass_linear_solver =
628 SUNLinSol_SPGMR(y_template,
632 sun_mass_linear_solver =
633 SUNLinSol_SPGMR(y_template,
641 data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
643 status = ARKStepSetMassLinearSolver(arkode_mem,
644 sun_mass_linear_solver,
646 mass_time_dependent);
649 auto mass_matrix_times_vector_setup_callback =
659 auto mass_matrix_times_vector_callback = [](N_Vector v,
662 void *mtimes_data) ->
int {
667 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
668 auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
678 status = ARKStepSetMassTimes(arkode_mem,
680 mass_matrix_times_vector_setup_callback :
681 ARKLsMassTimesSetupFn(
nullptr),
682 mass_matrix_times_vector_callback,
686 if (mass_preconditioner_solve)
688 auto mass_matrix_solver_setup_callback =
703 void *user_data) ->
int {
708 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
709 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
722 ARKStepSetMassPreconditioner(arkode_mem,
723 mass_preconditioner_setup ?
724 mass_matrix_solver_setup_callback :
725 ARKLsMassPrecSetupFn(
nullptr),
726 solve_with_mass_matrix_callback);
734 template <
typename VectorType>
738 solver_should_restart = [](
const double, VectorType &) ->
bool {
745 template <
typename VectorType>
758# ifdef DEAL_II_WITH_MPI
760# ifdef DEAL_II_WITH_TRILINOS
766# ifdef DEAL_II_WITH_PETSC
767# 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 &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
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
MPI_Comm mpi_communicator
void reset(const double t, const double h, const VectorType &y)
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< void(const double t)> mass_times_setup
void set_functions_to_trigger_an_assert()
std::function< void(const double t, const VectorType &v, VectorType &Mv)> mass_times_vector
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(const VectorType &v, VectorType &Jv, const double t, const VectorType &y, const VectorType &fy)> jacobian_times_vector
#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)