21#ifdef DEAL_II_WITH_SUNDIALS
29# ifdef DEAL_II_WITH_TRILINOS
33# ifdef DEAL_II_WITH_PETSC
41# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
42# include <arkode/arkode_impl.h>
43# include <sundials/sundials_config.h>
45# include <arkode/arkode_arkstep.h>
46# include <sunlinsol/sunlinsol_spgmr.h>
47# include <sunnonlinsol/sunnonlinsol_fixedpoint.h>
48# if DEAL_II_SUNDIALS_VERSION_LT(5, 0, 0)
55# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
68 template <
typename VectorType>
70 t_arkode_explicit_function(realtype tt,
76 ARKode<VectorType> &solver =
77 *
static_cast<ARKode<VectorType> *
>(user_data);
79 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
80 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
82 return solver.explicit_function(tt, *src_yy, *dst_yp);
87 template <
typename VectorType>
89 t_arkode_implicit_function(realtype tt,
95 ARKode<VectorType> &solver =
96 *
static_cast<ARKode<VectorType> *
>(user_data);
98 auto *src_yy = internal::unwrap_nvector_const<VectorType>(yy);
99 auto *dst_yp = internal::unwrap_nvector<VectorType>(yp);
101 return solver.implicit_function(tt, *src_yy, *dst_yp);
106# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
107 template <
typename VectorType>
109 t_arkode_setup_jacobian(ARKodeMem arkode_mem,
113 booleantype *jcurPtr,
119 ARKode<VectorType> &solver =
120 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
122 auto *src_ypred = internal::unwrap_nvector_const<VectorType>(ypred);
123 auto *src_fpred = internal::unwrap_nvector_const<VectorType>(fpred);
125 bool jcurPtr_tmp =
false;
126 int err = solver.setup_jacobian(convfail,
128 arkode_mem->ark_gamma,
132# if DEAL_II_SUNDIALS_VERSION_GTE(2, 0, 0)
133 *jcurPtr = jcurPtr_tmp ? SUNTRUE : SUNFALSE;
135 *jcurPtr = jcurPtr_tmp ? TRUE : FALSE;
143 template <
typename VectorType>
145 t_arkode_solve_jacobian(ARKodeMem arkode_mem,
154 ARKode<VectorType> &solver =
155 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
157 auto *dst = internal::unwrap_nvector<VectorType>(
b);
158 auto *src_ycur = internal::unwrap_nvector_const<VectorType>(ycur);
159 auto *src_fcur = internal::unwrap_nvector_const<VectorType>(fcur);
162 VectorType src = *dst;
164 int err = solver.solve_jacobian_system(arkode_mem->ark_tn,
165 arkode_mem->ark_gamma,
177 template <
typename VectorType>
179 t_arkode_setup_mass(ARKodeMem arkode_mem, N_Vector, N_Vector, N_Vector)
182 ARKode<VectorType> &solver =
183 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
184 int err = solver.setup_mass(arkode_mem->ark_tn);
190 template <
typename VectorType>
192 t_arkode_solve_mass(ARKodeMem arkode_mem,
202 ARKode<VectorType> &solver =
203 *
static_cast<ARKode<VectorType> *
>(arkode_mem->ark_user_data);
205 auto *dst = internal::unwrap_nvector<VectorType>(
b);
208 VectorType src = *dst;
210 int err = solver.solve_mass_system(src, *dst);
216 template <
typename VectorType>
218 t_arkode_jac_times_vec_function(N_Vector v,
227 ARKode<VectorType> &solver =
228 *
static_cast<ARKode<VectorType> *
>(user_data);
230 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
231 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
232 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
234 auto *dst_Jv = internal::unwrap_nvector<VectorType>(Jv);
236 return solver.jacobian_times_vector(*src_v, *dst_Jv, t, *src_y, *src_fy);
241 template <
typename VectorType>
243 t_arkode_jac_times_setup_function(realtype t,
249 ARKode<VectorType> &solver =
250 *
static_cast<ARKode<VectorType> *
>(user_data);
252 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
253 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
255 return solver.jacobian_times_setup(t, *src_y, *src_fy);
260 template <
typename VectorType>
262 t_arkode_prec_solve_function(realtype t,
273 ARKode<VectorType> &solver =
274 *
static_cast<ARKode<VectorType> *
>(user_data);
276 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
277 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
278 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
280 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
282 return solver.jacobian_preconditioner_solve(
283 t, *src_y, *src_fy, *src_r, *dst_z,
gamma, delta, lr);
288 template <
typename VectorType>
290 t_arkode_prec_setup_function(realtype t,
294 booleantype *jcurPtr,
299 ARKode<VectorType> &solver =
300 *
static_cast<ARKode<VectorType> *
>(user_data);
302 auto *src_y = internal::unwrap_nvector_const<VectorType>(y);
303 auto *src_fy = internal::unwrap_nvector_const<VectorType>(fy);
305 return solver.jacobian_preconditioner_setup(
306 t, *src_y, *src_fy, jok, *jcurPtr,
gamma);
311 template <
typename VectorType>
313 t_arkode_mass_times_vec_function(N_Vector v,
319 ARKode<VectorType> &solver =
320 *
static_cast<ARKode<VectorType> *
>(mtimes_data);
322 auto *src_v = internal::unwrap_nvector_const<VectorType>(v);
323 auto *dst_Mv = internal::unwrap_nvector<VectorType>(Mv);
325 return solver.mass_times_vector(t, *src_v, *dst_Mv);
330 template <
typename VectorType>
332 t_arkode_mass_times_setup_function(realtype t,
void *mtimes_data)
335 ARKode<VectorType> &solver =
336 *
static_cast<ARKode<VectorType> *
>(mtimes_data);
338 return solver.mass_times_setup(t);
343 template <
typename VectorType>
345 t_arkode_mass_prec_solve_function(realtype t,
353 ARKode<VectorType> &solver =
354 *
static_cast<ARKode<VectorType> *
>(user_data);
356 auto *src_r = internal::unwrap_nvector_const<VectorType>(r);
357 auto *dst_z = internal::unwrap_nvector<VectorType>(z);
359 return solver.mass_preconditioner_solve(t, *src_r, *dst_z, delta, lr);
364 template <
typename VectorType>
366 t_arkode_mass_prec_setup_function(realtype t,
void *user_data)
369 ARKode<VectorType> &solver =
370 *
static_cast<ARKode<VectorType> *
>(user_data);
372 return solver.mass_preconditioner_setup(t);
379 template <
typename VectorType>
383 , arkode_mem(nullptr)
387 , last_end_time(data.initial_time)
392 template <
typename VectorType>
396# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
397 ARKodeFree(&arkode_mem);
399 ARKStepFree(&arkode_mem);
402# ifdef DEAL_II_WITH_MPI
405 const int ierr = MPI_Comm_free(&communicator);
414 template <
typename VectorType>
420 data.initial_step_size);
422 return do_evolve_time(solution, time,
true);
427 template <
typename VectorType>
430 const double intermediate_time,
431 const bool reset_solver)
434 intermediate_time > last_end_time,
436 "The requested intermediate time is smaller than the last requested "
437 "intermediate time."));
439 const bool do_reset = reset_solver || arkode_mem ==
nullptr;
440 DiscreteTime time(last_end_time, intermediate_time, data.initial_step_size);
441 return do_evolve_time(solution, time, do_reset);
446 template <
typename VectorType>
461 auto solution_nvector = internal::make_nvector_view(solution);
466 double actual_next_time;
467# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
474 const auto status = ARKStepEvolve(arkode_mem,
502# if DEAL_II_SUNDIALS_VERSION_LT(4, 0, 0)
503 template <
typename VectorType>
506 const double current_time_step,
507 const VectorType &solution)
509 last_end_time = current_time;
511 ARKodeFree(&arkode_mem);
513 arkode_mem = ARKodeCreate();
518 Assert(explicit_function || implicit_function,
519 ExcFunctionNotProvided(
"explicit_function || implicit_function"));
523 auto initial_condition_nvector = internal::make_nvector_view(solution);
527 explicit_function ? &t_arkode_explicit_function<VectorType> :
nullptr,
528 implicit_function ? &t_arkode_implicit_function<VectorType> :
nullptr,
530 initial_condition_nvector);
533 if (get_local_tolerances)
535 const auto abs_tols = make_nvector_view(get_local_tolerances());
537 ARKodeSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
542 status = ARKodeSStolerances(arkode_mem,
543 data.relative_tolerance,
544 data.absolute_tolerance);
548 status = ARKodeSetInitStep(arkode_mem, current_time_step);
551 status = ARKodeSetUserData(arkode_mem,
this);
554 status = ARKodeSetStopTime(arkode_mem, data.final_time);
558 ARKodeSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
562 auto ARKode_mem =
static_cast<ARKodeMem
>(arkode_mem);
564 if (solve_jacobian_system)
566 status = ARKodeSetNewton(arkode_mem);
568 if (data.implicit_function_is_linear)
570 status = ARKodeSetLinear(
571 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
576 ARKode_mem->ark_lsolve = t_arkode_solve_jacobian<VectorType>;
579 ARKode_mem->ark_lsetup = t_arkode_setup_jacobian<VectorType>;
580# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
581 ARKode_mem->ark_setupNonNull =
true;
588 ARKodeSetFixedPoint(arkode_mem, data.maximum_non_linear_iterations);
593 if (solve_mass_system)
595 ARKode_mem->ark_msolve = t_arkode_solve_mass<VectorType>;
599 ARKode_mem->ark_msetup = t_arkode_setup_mass<VectorType>;
600# if DEAL_II_SUNDIALS_VERSION_LT(3, 0, 0)
601 ARKode_mem->ark_MassSetupNonNull =
true;
606 status = ARKodeSetOrder(arkode_mem, data.maximum_order);
610 custom_setup(arkode_mem);
615 template <
typename VectorType>
618 const double current_time_step,
619 const VectorType &solution)
621 last_end_time = current_time;
623 ARKStepFree(&arkode_mem);
630 auto initial_condition_nvector = internal::make_nvector_view(solution);
632 Assert(explicit_function || implicit_function,
633 ExcFunctionNotProvided(
"explicit_function || implicit_function"));
635 arkode_mem = ARKStepCreate(
636 explicit_function ? &t_arkode_explicit_function<VectorType> :
nullptr,
637 implicit_function ? &t_arkode_implicit_function<VectorType> :
nullptr,
639 initial_condition_nvector);
643 if (get_local_tolerances)
645 const auto abs_tols =
646 internal::make_nvector_view(get_local_tolerances());
648 ARKStepSVtolerances(arkode_mem, data.relative_tolerance, abs_tols);
653 status = ARKStepSStolerances(arkode_mem,
654 data.relative_tolerance,
655 data.absolute_tolerance);
660 status = ARKStepSetUserData(arkode_mem,
this);
663 setup_system_solver(solution);
665 setup_mass_solver(solution);
667 status = ARKStepSetInitStep(arkode_mem, current_time_step);
670 status = ARKStepSetStopTime(arkode_mem, data.final_time);
673 status = ARKStepSetOrder(arkode_mem, data.maximum_order);
677 custom_setup(arkode_mem);
682 template <
typename VectorType>
687 if (!implicit_function)
694 if (jacobian_times_vector)
697 if (solve_linearized_system)
700 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
701 solve_linearized_system);
702 sun_linear_solver = *linear_solver;
708 auto y_template = internal::make_nvector_view(solution);
710 SUNLinSol_SPGMR(y_template,
714 status = ARKStepSetLinearSolver(arkode_mem, sun_linear_solver,
nullptr);
717 ARKStepSetJacTimes(arkode_mem,
718 jacobian_times_setup ?
719 t_arkode_jac_times_setup_function<VectorType> :
721 t_arkode_jac_times_vec_function<VectorType>);
723 if (jacobian_preconditioner_solve)
725 status = ARKStepSetPreconditioner(
727 jacobian_preconditioner_setup ?
728 t_arkode_prec_setup_function<VectorType> :
730 t_arkode_prec_solve_function<VectorType>);
733 if (data.implicit_function_is_linear)
735 status = ARKStepSetLinear(
736 arkode_mem, data.implicit_function_is_time_independent ? 0 : 1);
742 auto y_template = internal::make_nvector_view(solution);
744 SUNNonlinearSolver fixed_point_solver =
745 SUNNonlinSol_FixedPoint(y_template,
746 data.anderson_acceleration_subspace);
748 status = ARKStepSetNonlinearSolver(arkode_mem, fixed_point_solver);
753 ARKStepSetMaxNonlinIters(arkode_mem, data.maximum_non_linear_iterations);
759 template <
typename VectorType>
766 if (mass_times_vector)
772 std::make_unique<internal::LinearSolverWrapper<VectorType>>(
774 sun_mass_linear_solver = *mass_solver;
778 auto y_template = internal::make_nvector_view(solution);
779 sun_mass_linear_solver =
780 SUNLinSol_SPGMR(y_template,
784 booleantype mass_time_dependent =
785 data.mass_is_time_independent ? SUNFALSE : SUNTRUE;
786 status = ARKStepSetMassLinearSolver(arkode_mem,
787 sun_mass_linear_solver,
789 mass_time_dependent);
792 ARKStepSetMassTimes(arkode_mem,
794 t_arkode_mass_times_setup_function<VectorType> :
796 t_arkode_mass_times_vec_function<VectorType>,
800 if (mass_preconditioner_solve)
802 status = ARKStepSetMassPreconditioner(
804 mass_preconditioner_setup ?
805 t_arkode_mass_prec_setup_function<VectorType> :
807 t_arkode_mass_prec_solve_function<VectorType>);
816 template <
typename VectorType>
820 reinit_vector = [](VectorType &) {
821 AssertThrow(
false, ExcFunctionNotProvided(
"reinit_vector"));
824 solver_should_restart = [](
const double, VectorType &) ->
bool {
831 template <
typename VectorType>
844# ifdef DEAL_II_WITH_MPI
846# ifdef DEAL_II_WITH_TRILINOS
852# ifdef DEAL_II_WITH_PETSC
853# ifndef PETSC_USE_COMPLEX
const auto & SundialsARKode
#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
void setup_mass_solver(const VectorType &solution)
void * get_arkode_memory() const
ARKode(const AdditionalData &data=AdditionalData(), const MPI_Comm &mpi_comm=MPI_COMM_WORLD)
int do_evolve_time(VectorType &solution, ::DiscreteTime &time, const bool do_reset)
void reset(const double t, const double h, const VectorType &y)
void set_functions_to_trigger_an_assert()
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)
#define DEAL_II_NAMESPACE_OPEN
#define DEAL_II_SUNDIALS_VERSION_LT(major, minor, patch)
#define DEAL_II_NAMESPACE_CLOSE
#define Assert(cond, exc)
#define AssertNothrow(cond, exc)
static ::ExceptionBase & ExcInternalError()
static ::ExceptionBase & ExcMessage(std::string arg1)
#define AssertThrow(cond, exc)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
MPI_Comm duplicate_communicator(const MPI_Comm &mpi_communicator)
long double gamma(const unsigned int n)